diff --git a/Project.toml b/Project.toml index 87b2c213ed..63633872ea 100644 --- a/Project.toml +++ b/Project.toml @@ -77,7 +77,7 @@ ASEconvert = "0.1" AbstractFFTs = "1" Aqua = "0.8.5" Artifacts = "1" -AtomsBase = "0.5" +AtomsBase = "0.5.2" AtomsBuilder = "0.2" AtomsCalculators = "0.2.3" AtomsIO = "0.3" diff --git a/docs/make.jl b/docs/make.jl index 4280db725e..55aa0a8402 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -42,6 +42,7 @@ PAGES = [ "examples/graphene.jl", "examples/geometry_optimization.jl", "examples/energy_cutoff_smearing.jl", + "examples/hubbard.jl", ], "Response and properties" => [ "examples/elastic_constants.jl", diff --git a/docs/src/features.md b/docs/src/features.md index fd2820cee2..b543594325 100644 --- a/docs/src/features.md +++ b/docs/src/features.md @@ -8,6 +8,7 @@ and runs out of the box on Linux, windows and macOS, see [Installation](@ref). ## Standard methods and models - DFT models (LDA, GGA, meta-GGA): Any functional from the [libxc](https://libxc.gitlab.io/) library +- DFT+U and Hubbard corrections - **Norm-conserving pseudopotentials**: Goedecker-type (GTH) or numerical (in UPF pseudopotential format), see [Pseudopotentials](@ref). @@ -25,7 +26,7 @@ and runs out of the box on Linux, windows and macOS, see [Installation](@ref). ## Ground-state properties and post-processing - Total energy, forces, stresses -- Density of states (DOS), local density of states (LDOS) +- Density of states (DOS), local density of states (LDOS), projected density of states (PDOS) - Band structures - Easy access to all intermediate quantities (e.g. density, Bloch waves) diff --git a/examples/hubbard.jl b/examples/hubbard.jl new file mode 100644 index 0000000000..720711ca32 --- /dev/null +++ b/examples/hubbard.jl @@ -0,0 +1,65 @@ +# # Hubbard correction (DFT+U) +# In this example, we'll plot the DOS and projected DOS of Nickel Oxide +# with and without the Hubbard term correction. + +using DFTK +using PseudoPotentialData +using Unitful +using UnitfulAtomic +using Plots + +# Define the geometry and pseudopotential +a = 7.9 # Nickel Oxide lattice constant in Bohr +lattice = a * [[ 1.0 0.5 0.5]; + [ 0.5 1.0 0.5]; + [ 0.5 0.5 1.0]] +pseudopotentials = PseudoFamily("dojo.nc.sr.pbe.v0_4_1.standard.upf") +Ni = ElementPsp(:Ni, pseudopotentials) +O = ElementPsp(:O, pseudopotentials) +atoms = [Ni, O, Ni, O] +positions = [zeros(3), ones(3) / 4, ones(3) / 2, ones(3) * 3 / 4] +magnetic_moments = [2, 0, -1, 0] + +# First, we run an SCF and band computation without the Hubbard term +model = model_DFT(lattice, atoms, positions; temperature=5e-3, + functionals=PBE(), magnetic_moments) +basis = PlaneWaveBasis(model; Ecut=20, kgrid=[2, 2, 2]) +scfres = self_consistent_field(basis; tol=1e-6, ρ=guess_density(basis, magnetic_moments)) +bands = compute_bands(scfres, MonkhorstPack(4, 4, 4)) +lowest_unocc_band = findfirst(ε -> ε-bands.εF > 0, bands.eigenvalues[1]) +band_gap = bands.eigenvalues[1][lowest_unocc_band] - bands.eigenvalues[1][lowest_unocc_band-1] + +# Then we plot the DOS and the PDOS for the relevant 3D (pseudo)atomic projector +εF = bands.εF +width = 5.0u"eV" +εrange = (εF - austrip(width), εF + austrip(width)) +p = plot_dos(bands; εrange, colors=[:red, :red]) +plot_pdos(bands; p, iatom=1, label="3D", colors=[:yellow, :orange], εrange) + +# To perform and Hubbard computation, we have to define the Hubbard manifold and associated constant. +# +# In DFTK there are a few ways to construct the `OrbitalManifold`. +# Here, we will apply the Hubbard correction on the 3D orbital of all nickel atoms. +# +# Note that "manifold" is the standard term used in the literature for the set of atomic orbitals +# used to compute the Hubbard correction, but it is not meant in the mathematical sense. +U = 10u"eV" +manifold = OrbitalManifold(atoms, Ni, "3D") + +# Run SCF with a DFT+U setup, notice the `extra_terms` keyword argument, setting up the Hubbard +U term. +model = model_DFT(lattice, atoms, positions; extra_terms=[Hubbard(manifold, U)], + functionals=PBE(), temperature=5e-3, magnetic_moments) +basis = PlaneWaveBasis(model; Ecut=20, kgrid=[2, 2, 2]) +scfres = self_consistent_field(basis; tol=1e-6, ρ=guess_density(basis, magnetic_moments)) + +# Run band computation +bands_hub = compute_bands(scfres, MonkhorstPack(4, 4, 4)) +lowest_unocc_band = findfirst(ε -> ε-bands_hub.εF > 0, bands_hub.eigenvalues[1]) +band_gap = bands_hub.eigenvalues[1][lowest_unocc_band] - bands_hub.eigenvalues[1][lowest_unocc_band-1] + +# With the electron localization introduced by the Hubbard term, the band gap has now opened, +# reflecting the experimental insulating behaviour of Nickel Oxide. +εF = bands_hub.εF +εrange = (εF - austrip(width), εF + austrip(width)) +p = plot_dos(bands_hub; p, colors=[:blue, :blue], εrange) +plot_pdos(bands_hub; p, iatom=1, label="3D", colors=[:green, :purple], εrange) diff --git a/ext/DFTKPlotsExt.jl b/ext/DFTKPlotsExt.jl index fad3ba9b5c..d88d94f472 100644 --- a/ext/DFTKPlotsExt.jl +++ b/ext/DFTKPlotsExt.jl @@ -76,6 +76,7 @@ end function plot_dos(basis, eigenvalues; εF=nothing, unit=u"hartree", temperature=basis.model.temperature, smearing=basis.model.smearing, + colors=[:blue, :red], p=nothing, εrange=default_band_εrange(eigenvalues; εF), n_points=1000, kwargs...) # TODO Should also split this up into one stage doing the DOS computation # and one stage doing the DOS plotting (like now for the bands.) @@ -87,9 +88,9 @@ function plot_dos(basis, eigenvalues; εF=nothing, unit=u"hartree", # Constant to convert from AU to the desired unit to_unit = ustrip(auconvert(unit, 1.0)) - p = Plots.plot(; kwargs...) + isnothing(p) && (p = Plots.plot()) + p = Plots.plot(p; kwargs...) spinlabels = spin_components(basis.model) - colors = [:blue, :red] Dεs = compute_dos.(εs, Ref(basis), Ref(eigenvalues); smearing, temperature) for σ = 1:n_spin D = [Dσ[σ] for Dσ in Dεs] @@ -141,7 +142,7 @@ function plot_ldos(basis, eigenvalues, ψ; εF=nothing, unit=u"hartree", end plot_ldos(scfres; kwargs...) = plot_ldos(scfres.basis, scfres.eigenvalues, scfres.ψ; scfres.εF, kwargs...) -function plot_pdos(basis::PlaneWaveBasis{T}, eigenvalues, ψ; iatom, label=nothing, +function plot_pdos(basis::PlaneWaveBasis{T}, eigenvalues, ψ; iatom=nothing, label=nothing, positions=basis.model.positions, εF=nothing, unit=u"hartree", temperature=basis.model.temperature, @@ -154,16 +155,17 @@ function plot_pdos(basis::PlaneWaveBasis{T}, eigenvalues, ψ; iatom, label=nothi n_spin = basis.model.n_spin_components to_unit = ustrip(auconvert(unit, 1.0)) - species = isnothing(iatom) ? "all atoms" : "atom $(iatom) ($(basis.model.atoms[iatom].species))" + species = isnothing(iatom) ? "all atoms" : "atom $(iatom) ($(basis.model.atoms[iatom].species))" orb_name = isnothing(label) ? "all orbitals" : label # Plot pdos - isnothing(p) && (p = Plots.plot(; kwargs...)) + isnothing(p) && (p = Plots.plot()) p = Plots.plot(p; kwargs...) spinlabels = spin_components(basis.model) pdos = DFTK.sum_pdos(compute_pdos(εs, basis, ψ, eigenvalues; - positions, temperature, smearing), - [DFTK.OrbitalManifold(;iatom, label)]) + positions, temperature, smearing), + [l -> ((isnothing(iatom) || l.iatom == iatom) + && (isnothing(label) || l.label == label))]) for σ = 1:n_spin plot_label = n_spin > 1 ? "$(species) $(orb_name) $(spinlabels[σ]) spin" : "$(species) $(orb_name)" Plots.plot!(p, (εs .- eshift) .* to_unit, pdos[:, σ]; label=plot_label, color=colors[σ]) diff --git a/src/DFTK.jl b/src/DFTK.jl index a812986bbb..85cfd850e0 100644 --- a/src/DFTK.jl +++ b/src/DFTK.jl @@ -112,6 +112,8 @@ export AtomicNonlocal export Ewald export PspCorrection export Entropy +export Hubbard +export OrbitalManifold export Magnetic export PairwisePotential export Anyonic diff --git a/src/common/spherical_harmonics.jl b/src/common/spherical_harmonics.jl index ec46b871eb..287a8521bb 100644 --- a/src/common/spherical_harmonics.jl +++ b/src/common/spherical_harmonics.jl @@ -47,3 +47,33 @@ function ylm_real(l::Integer, m::Integer, rvec::AbstractVector{T}) where {T} error("The case l = $l and m = $m is not implemented") end + +""" + This function returns the Wigner D matrix for real spherical harmonics, + for a given l and orthogonal matrix, solving a randomized linear system. + Such D matrix gives the decomposition of a spherical harmonic after application + of an orthogonal matrix back into the basis of spherical harmonics. + + Yₗₘ₁(Wr) = Σₘ₂ D(l,R̂)ₘ₁ₘ₂ * Yₗₘ₂(r) +""" +function wigner_d_matrix(l::Integer, Wcart::AbstractMatrix{T}) where {T} + if l == 0 # In this case no computation is needed + return [one(T);;] + end + rng = Xoshiro(1234) + neq = 2l+2 # We need at least 2l+1 equations, we add one for numerical stability + B = Matrix{T}(undef, 2l+1, neq) + A = Matrix{T}(undef, 2l+1, neq) + for n in 1:neq + r = rand(rng, T, 3) + r = r / norm(r) + r0 = Wcart * r + for m in -l:l + B[m+l+1, n] = ylm_real(l, m, r0) + A[m+l+1, n] = ylm_real(l, m, r) + end + end + κ = cond(A) + @assert κ < 100.0 "The Wigner matrix computation is badly conditioned. κ(A)=$(κ)" + B / A +end \ No newline at end of file diff --git a/src/postprocess/band_structure.jl b/src/postprocess/band_structure.jl index 8e264da302..d991a3405f 100644 --- a/src/postprocess/band_structure.jl +++ b/src/postprocess/band_structure.jl @@ -14,9 +14,9 @@ All kwargs not specified below are passed to [`diagonalize_all_kblocks`](@ref): @timing function compute_bands(basis::PlaneWaveBasis, kgrid::Union{AbstractKgrid,AbstractKgridGenerator}; n_bands=default_n_bands_bandstructure(basis.model), - n_extra=3, ρ=nothing, τ=nothing, εF=nothing, - eigensolver=lobpcg_hyper, tol=1e-3, seed=nothing, - kwargs...) + n_extra=3, ρ=nothing, τ=nothing, hubbard_n=nothing, + εF=nothing, eigensolver=lobpcg_hyper, tol=1e-3, + seed=nothing, kwargs...) # kcoords are the kpoint coordinates in fractional coordinates if isnothing(ρ) if any(t isa TermNonlinear for t in basis.terms) @@ -36,7 +36,7 @@ All kwargs not specified below are passed to [`diagonalize_all_kblocks`](@ref): # Create new basis with new kpoints bs_basis = PlaneWaveBasis(basis, kgrid) - ham = Hamiltonian(bs_basis; ρ, τ) + ham = Hamiltonian(bs_basis; ρ, τ, hubbard_n) eigres = diagonalize_all_kblocks(eigensolver, ham, n_bands + n_extra; n_conv_check=n_bands, tol, kwargs...) if !eigres.converged @@ -68,7 +68,10 @@ function compute_bands(scfres::NamedTuple, kgrid::Union{AbstractKgrid,AbstractKgridGenerator}; n_bands=default_n_bands_bandstructure(scfres), kwargs...) τ = haskey(scfres, :τ) ? scfres.τ : nothing - compute_bands(scfres.basis, kgrid; scfres.ρ, τ, scfres.εF, n_bands, kwargs...) + hubbard_n = haskey(scfres, :hubbard_n) ? scfres.hubbard_n : nothing + compute_bands(scfres.basis, kgrid; + scfres.ρ, τ, hubbard_n, + scfres.εF, n_bands, kwargs...) end """ diff --git a/src/postprocess/dos.jl b/src/postprocess/dos.jl index 3edf643ec4..cf5493272a 100644 --- a/src/postprocess/dos.jl +++ b/src/postprocess/dos.jl @@ -119,31 +119,6 @@ function compute_pdos(εs, bands; kwargs...) compute_pdos(εs, bands.basis, bands.ψ, bands.eigenvalues; kwargs...) end -""" -Structure for manifold choice and projectors extraction. - -Overview of fields: -- `iatom`: Atom position in the atoms array. -- `species`: Chemical Element as in ElementPsp. -- `label`: Orbital name, e.g.: "3S". - -All fields are optional, only the given ones will be used for selection. -Can be called with an orbital NamedTuple and returns a boolean - stating whether the orbital belongs to the manifold. -""" -@kwdef struct OrbitalManifold - iatom ::Union{Int64, Nothing} = nothing - species ::Union{Symbol, AtomsBase.ChemicalSpecies, Nothing} = nothing - label ::Union{String, Nothing} = nothing -end -function (s::OrbitalManifold)(orb) - iatom_match = isnothing(s.iatom) || (s.iatom == orb.iatom) - # See JuliaMolSim/AtomsBase.jl#139 why both species equalities are needed - species_match = isnothing(s.species) || (s.species == orb.species) || (orb.species == s.species) - label_match = isnothing(s.label) || (s.label == orb.label) - iatom_match && species_match && label_match -end - """ atomic_orbital_projectors(basis; [isonmanifold, positions]) @@ -161,8 +136,8 @@ The projectors are computed by decomposition into a form factor multiplied by a Overview of inputs: - `positions` : Positions of the atoms in the unit cell. Default is model.positions. -- `isonmanifold` (opt) : (see notes below) OrbitalManifold struct to select only a subset of orbitals - for the computation. +- `isonmanifold` (opt) : (see notes below) A function, typically a lambda, + to select projectors to include in the pdos. Overview of outputs: - `projectors`: Vector of matrices of projectors @@ -191,6 +166,7 @@ function atomic_orbital_projectors(basis::PlaneWaveBasis{T}; labels = [] for (iatom, atom) in enumerate(basis.model.atoms) psp = atom.psp + @assert count_n_pswfc(psp) > 0 # We need this to check if we have any atomic orbital projector for l in 0:psp.lmax, n in 1:DFTK.count_n_pswfc_radial(psp, l) label = DFTK.pswfc_label(psp, n, l) if !isonmanifold((; iatom, atom.species, label)) @@ -236,23 +212,24 @@ function atomic_orbital_projections(basis::PlaneWaveBasis{T}, ψ; end """ - sum_pdos(pdos_res, manifolds) + sum_pdos(pdos_res, projector_filters) This function extracts and sums up all the PDOSes, directly from the output of the `compute_pdos` function, - that match any of the manifolds. + that match any of the filters. Overview of inputs: - `pdos_res`: Whole output from compute_pdos. -- `manifolds`: Vector of OrbitalManifolds to select the desired projectors pdos. +- `projector_filters`: Vector of functions, typically lambdas, + to select projectors to include in the pdos. Overview of outputs: - `pdos`: Vector containing the pdos(ε). """ -function sum_pdos(pdos_res, manifolds::AbstractVector) +function sum_pdos(pdos_res, projector_filters::AbstractVector) pdos = zeros(Float64, length(pdos_res.εs), size(pdos_res.pdos, 3)) for σ in 1:size(pdos_res.pdos, 3) for (j, orb) in enumerate(pdos_res.projector_labels) - if any(manifold(orb) for manifold in manifolds) + if any(filt(orb) for filt in projector_filters) pdos[:, σ] += pdos_res.pdos[:, j, σ] end end diff --git a/src/pseudo/NormConservingPsp.jl b/src/pseudo/NormConservingPsp.jl index ccff5bc874..b3088195a2 100644 --- a/src/pseudo/NormConservingPsp.jl +++ b/src/pseudo/NormConservingPsp.jl @@ -221,3 +221,13 @@ count_n_pswfc(psp::NormConservingPsp, l) = count_n_pswfc_radial(psp, l) * (2l + function count_n_pswfc(psp::NormConservingPsp) sum(l -> count_n_pswfc(psp, l), 0:psp.lmax; init=0)::Int end + +function find_pswfc(psp::NormConservingPsp, label::String) + for l = 0:psp.lmax, i = 1:count_n_pswfc_radial(psp, l) + if pswfc_label(psp, i, l) == label + return (; l, i) + end + end + error("Could not find pseudo atomic orbital with label $label " + * "in pseudopotential $(psp.identifier).") +end diff --git a/src/scf/self_consistent_field.jl b/src/scf/self_consistent_field.jl index 699935aaf7..5a11b2eb27 100644 --- a/src/scf/self_consistent_field.jl +++ b/src/scf/self_consistent_field.jl @@ -131,6 +131,7 @@ Overview of parameters: basis::PlaneWaveBasis{T}; ρ=guess_density(basis), τ=any(needs_τ, basis.terms) ? zero(ρ) : nothing, + hubbard_n=nothing, ψ=nothing, tol=1e-6, is_converged=ScfConvergenceDensity(tol), @@ -159,12 +160,12 @@ Overview of parameters: # We do density mixing in the real representation # TODO support other mixing types function fixpoint_map(ρin, info) - (; ψ, occupation, eigenvalues, εF, n_iter, converged, timedout, τ) = info + (; ψ, occupation, eigenvalues, εF, n_iter, converged, timedout, τ, hubbard_n) = info n_iter += 1 # Note that ρin is not the density of ψ, and the eigenvalues # are not the self-consistent ones, which makes this energy non-variational - energies, ham = energy_hamiltonian(basis, ψ, occupation; ρ=ρin, τ, eigenvalues, εF) + energies, ham = energy_hamiltonian(basis, ψ, occupation; ρ=ρin, τ, hubbard_n, eigenvalues, εF) # Diagonalize `ham` to get the new state nextstate = next_density(ham, nbandsalg, fermialg; eigensolver, ψ, eigenvalues, @@ -180,16 +181,20 @@ Overview of parameters: if any(needs_τ, basis.terms) τ = compute_kinetic_energy_density(basis, ψ, occupation) end + ihubbard = findfirst(t -> t isa TermHubbard, basis.terms) + if !isnothing(ihubbard) + hubbard_n = compute_hubbard_n(basis.terms[ihubbard], basis, ψ, occupation) + end # Update info with results gathered so far info_next = (; ham, basis, converged, stage=:iterate, algorithm="SCF", - ρin, τ, α=damping, n_iter, nbandsalg.occupation_threshold, + ρin, τ, hubbard_n, α=damping, n_iter, nbandsalg.occupation_threshold, seed, runtime_ns=time_ns() - start_ns, nextstate..., diagonalization=[nextstate.diagonalization]) # Compute the energy of the new state if compute_consistent_energies - (; energies) = energy(basis, ψ, occupation; ρ=ρout, τ, eigenvalues, εF) + (; energies) = energy(basis, ψ, occupation; ρ=ρout, τ, hubbard_n, eigenvalues, εF) end history_Etot = vcat(info.history_Etot, energies.total) history_Δρ = vcat(info.history_Δρ, norm(Δρ) * sqrt(basis.dvol)) @@ -211,7 +216,7 @@ Overview of parameters: ρnext, info_next end - info_init = (; ρin=ρ, τ, ψ, occupation=nothing, eigenvalues=nothing, εF=nothing, + info_init = (; ρin=ρ, τ, hubbard_n, ψ, occupation=nothing, eigenvalues=nothing, εF=nothing, n_iter=0, n_matvec=0, timedout=false, converged=false, history_Etot=T[], history_Δρ=T[]) @@ -221,15 +226,15 @@ Overview of parameters: # We do not use the return value of solver but rather the one that got updated by fixpoint_map # ψ is consistent with ρout, so we return that. We also perform a last energy computation # to return a correct variational energy - (; ρin, ρout, τ, ψ, occupation, eigenvalues, εF, converged) = info - energies, ham = energy_hamiltonian(basis, ψ, occupation; ρ=ρout, τ, eigenvalues, εF) + (; ρin, ρout, τ, hubbard_n, ψ, occupation, eigenvalues, εF, converged) = info + energies, ham = energy_hamiltonian(basis, ψ, occupation; ρ=ρout, τ, hubbard_n, eigenvalues, εF) # Callback is run one last time with final state to allow callback to clean up scfres = (; ham, basis, energies, converged, nbandsalg.occupation_threshold, - ρ=ρout, τ, α=damping, eigenvalues, occupation, εF, info.n_bands_converge, - info.n_iter, info.n_matvec, ψ, info.diagonalization, stage=:finalize, - info.history_Δρ, info.history_Etot, info.timedout, mixing, seed, - runtime_ns=time_ns() - start_ns, algorithm="SCF") + ρ=ρout, τ, hubbard_n, α=damping, eigenvalues, occupation, εF, + info.n_bands_converge, info.n_iter, info.n_matvec, ψ, info.diagonalization, + stage=:finalize, info.history_Δρ, info.history_Etot, info.timedout, mixing, + seed, runtime_ns=time_ns() - start_ns, algorithm="SCF") callback(scfres) scfres end diff --git a/src/symmetry.jl b/src/symmetry.jl index 24ce4ed332..af22a2bdbf 100644 --- a/src/symmetry.jl +++ b/src/symmetry.jl @@ -305,7 +305,7 @@ function accumulate_over_symmetries!(ρaccu, ρin, basis::PlaneWaveBasis{T}, sym # Looping over symmetries inside of map! on G vectors allow for a single GPU kernel launch map!(ρaccu, Gs) do G acc = zero(complex(T)) - # Explicit loop over indicies because AMDGPU does not support zip() in map! + # Explicit loop over indices because AMDGPU does not support zip() in map! for i_symm in 1:n_symm invS = symm_invS[i_symm] τ = symm_τ[i_symm] @@ -422,6 +422,35 @@ function symmetrize_forces(basis::PlaneWaveBasis, forces) symmetrize_forces(basis.model, forces; basis.symmetries) end +""" +Symmetrize the Hubbard occupation matrix according to the l quantum number of the manifold. +""" +function symmetrize_hubbard_n(model, manifold::OrbitalManifold, + hubbard_n::Array{Matrix{Complex{T}}}; + symmetries, tol_symmetry=SYMMETRY_TOLERANCE) where {T} + # For now we apply symmetries only on nII terms, not on cross-atom terms (nIJ) + # WARNING: To implement +V this will need to be changed! + + positions = model.positions[manifold.iatoms] + nspins = size(hubbard_n, 1) + natoms = size(hubbard_n, 2) + ldim = 2*manifold.l+1 + + # Initialize the hubbard_n matrix + ns = [zeros(Complex{T}, ldim, ldim) for _ in 1:nspins, _ in 1:natoms, _ in 1:natoms] + for symmetry in symmetries + Wcart = model.lattice * symmetry.W * model.inv_lattice + WigD = wigner_d_matrix(manifold.l, Wcart) + for σ in 1:nspins, iatom in 1:natoms + sym_atom = find_symmetry_preimage(positions, positions[iatom], symmetry; + tol_symmetry) + ns[σ, iatom, iatom] .+= WigD' * hubbard_n[σ, sym_atom, sym_atom] * WigD + end + end + ns ./= length(symmetries) + ns +end + """" Convert a `basis` into one that doesn't use BZ symmetry. This is mainly useful for debug purposes (e.g. in cases we don't want to @@ -494,7 +523,7 @@ function unfold_bz(scfres) eigenvalues = unfold_array(scfres.basis, basis_unfolded, scfres.eigenvalues, false) occupation = unfold_array(scfres.basis, basis_unfolded, scfres.occupation, false) energies, ham = energy_hamiltonian(basis_unfolded, ψ, occupation; - scfres.ρ, eigenvalues, scfres.εF) + scfres.ρ, scfres.hubbard_n, eigenvalues, scfres.εF) @assert energies.total ≈ scfres.energies.total new_scfres = (; basis=basis_unfolded, ψ, ham, eigenvalues, occupation) merge(scfres, new_scfres) diff --git a/src/terms/hubbard.jl b/src/terms/hubbard.jl new file mode 100644 index 0000000000..d14303790e --- /dev/null +++ b/src/terms/hubbard.jl @@ -0,0 +1,199 @@ +using LinearAlgebra +using Random + +""" +Structure for Hubbard manifold choice and projectors extraction. + +"Manifold" is the standard name used in the literature to refer +to the set of atomic orbitals used to compute the Hubbard correction. +It is to be noted that this is not meant in the mathematical sense of "manifold". + +Overview of fields: +- `psp`: Pseudopotential containing the atomic orbital projectors +- `iatoms`: Atom indices that are part of the manifold. +- `l`: Angular momentum quantum number of the manifold. +- `i`: Index of the atomic orbital within the given l. + +See also the convenience constructors, to construct a manifold more easily. +""" +struct OrbitalManifold + psp::NormConservingPsp + iatoms::Vector{Int64} + l::Int64 + i::Int64 +end +function OrbitalManifold(atoms::Vector{<:Element}, atom::ElementPsp, label::AbstractString) + OrbitalManifold(atom, findall(at -> at === atom, atoms), label) +end +function OrbitalManifold(atom::ElementPsp, iatoms::Vector{Int64}, label::AbstractString) + OrbitalManifold(atom.psp, iatoms, label) +end +function OrbitalManifold(psp::NormConservingPsp, iatoms::Vector{Int64}, label::AbstractString) + (; l, i) = find_pswfc(psp, label) + OrbitalManifold(psp, iatoms, l, i) +end + +function check_hubbard_manifold(manifold::OrbitalManifold, model::Model) + for atom in model.atoms[manifold.iatoms] + atom isa ElementPsp || error("Orbital manifold elements must have a psp.") + atom.psp === manifold.psp || error("Orbital manifold psp $(manifold.psp.identifier) " * + "does not match the psp of atom $atom") + end + isempty(manifold.iatoms) && error("Orbital manifold has no atoms.") + # Tricky: make sure that the iatoms are consistent with the symmetries of the model, + # i.e. that the manifold would be a valid atom group + atom_positions = model.positions[manifold.iatoms] + for symop in model.symmetries + W, w = symop.W, symop.w + for coord in atom_positions + # If all elements of a difference in diffs is integer, then + # W * coord + w and pos are equivalent lattice positions + if !any(c -> is_approx_integer(W * coord + w - c; atol=SYMMETRY_TOLERANCE), + atom_positions) + error("Inconsistency between orbital manifold and model symmetries: " * + "Cannot map the atom at position $coord to another atom of the manifold " * + "under the symmetry operation (W, w):\n($W, $w)") + end + end + end + + nothing +end + +function extract_manifold(manifold::OrbitalManifold, projectors, labels) + # We extract the labels only for orbitals belonging to the manifold + proj_indices = findall(orb -> (orb.iatom ∈ manifold.iatoms + && orb.l == manifold.l + && orb.n == manifold.i), labels) + isempty(proj_indices) && @warn "Projector for $(manifold) not found." + manifold_labels = labels[proj_indices] + manifold_projectors = map(enumerate(projectors)) do (ik, projk) + projk[:, proj_indices] + end + (; manifold_labels, manifold_projectors) +end + +@doc raw""" +Hubbard energy, following the Dudarev et al. (1998) rotationally invariant formalism: +```math +1/2 Σ_{σI} U * Tr[hubbard_n[σ,I,I] * (1 - hubbard_n[σ,I,I])] +``` +""" +struct Hubbard{T} + manifold::OrbitalManifold + U::T + function Hubbard(manifold::OrbitalManifold, U) + U = austrip(U) + new{typeof(U)}(manifold, U) + end +end +function (hubbard::Hubbard)(basis::AbstractBasis) + check_hubbard_manifold(hubbard.manifold, basis.model) + + projs, labels = atomic_orbital_projectors(basis) + labels, projectors_matrix = extract_manifold(hubbard.manifold, projs, labels) + TermHubbard(hubbard.manifold, hubbard.U, projectors_matrix, labels) +end + +struct TermHubbard{T, PT, L} <: Term + manifold::OrbitalManifold + U::T + P::PT + labels::L +end + +@timing "ene_ops: hubbard" function ene_ops(term::TermHubbard, + basis::PlaneWaveBasis{T}, + ψ, occupation; hubbard_n=nothing, + kwargs...) where {T} + if isnothing(hubbard_n) + return (; E=zero(T), ops=[NoopOperator(basis, kpt) for kpt in basis.kpoints]) + end + proj = term.P + + filled_occ = filled_occupation(basis.model) + natoms = length(term.manifold.iatoms) + n_spin = basis.model.n_spin_components + nproj_atom = size(hubbard_n[1,1,1], 1) # This is the number of projectors per atom, namely 2l+1 + # For the ops we have to reshape hubbard_n to match the NonlocalOperator structure, using a block diagonal form + nhub = [zeros(Complex{T}, nproj_atom*natoms, nproj_atom*natoms) for _ in 1:n_spin] + E = zero(T) + for σ in 1:n_spin, iatom in 1:natoms + proj_range = (1+nproj_atom*(iatom-1)):(nproj_atom*iatom) + nhub[σ][proj_range, proj_range] = hubbard_n[σ, iatom, iatom] + E += filled_occ * 1/T(2) * term.U * + real(tr(hubbard_n[σ, iatom, iatom] * (I - hubbard_n[σ, iatom, iatom]))) + end + ops = [NonlocalOperator(basis, kpt, proj[ik], 1/T(2) * term.U * (I - 2*nhub[kpt.spin])) + for (ik,kpt) in enumerate(basis.kpoints)] + return (; E, ops) +end + +""" + compute_hubbard_n(term::TermHubbard, basis, ψ, occupation) + +Computes a matrix hubbard_n of size (n_spin, natoms, natoms), where each entry hubbard_n[σ, iatom, jatom] +contains the submatrix of the occupation matrix corresponding to the projectors +of atom iatom and atom jatom, with dimensions determined by the number of projectors for each atom. +The atoms and orbitals are defined by the manifold tuple. + + hubbard_n[σ, iatom, jatom][m1, m2] = Σₖ₍ₛₚᵢₙ₎Σₙ weights[ik, ibnd] * ψₙₖ' * Pᵢₘ₁ * Pᵢₘ₂' * ψₙₖ + +where n or ibnd is the band index, ``weights[ik, ibnd] = kweights[ik] * occupation[ik, ibnd]`` +and ``Pᵢₘ₁`` is the pseudoatomic orbital projector for atom i and orbital m₁ +(just the magnetic quantum number, since l is fixed, as is usual in the literature). +For details on the projectors see `atomic_orbital_projectors`. +""" +function compute_hubbard_n(term::TermHubbard, + basis::PlaneWaveBasis{T}, + ψ, occupation) where {T} + filled_occ = filled_occupation(basis.model) + n_spin = basis.model.n_spin_components + + manifold = term.manifold + manifold_atoms = manifold.iatoms + natoms = length(manifold_atoms) + l = manifold.l + projectors = reshape_hubbard_proj(term.P, term.labels, manifold) + hubbard_n = Array{Matrix{Complex{T}}}(undef, n_spin, natoms, natoms) + for σ in 1:n_spin + for idx in 1:natoms, jdx in 1:natoms + hubbard_n[σ, idx, jdx] = zeros(Complex{T}, 2l+1, 2l+1) + for ik = krange_spin(basis, σ) + j_projection = ψ[ik]' * projectors[ik][jdx] # <ψ|ϕJ> + i_projection = projectors[ik][idx]' * ψ[ik] # <ϕI|ψ> + # Sums over the bands, dividing by filled_occ to deal + # with the physical two spin channels separately + hubbard_n[σ, idx, jdx] .+= (basis.kweights[ik] * i_projection * + diagm(occupation[ik]/filled_occ) * j_projection) + end + hubbard_n[σ, idx, jdx] = mpi_sum(hubbard_n[σ, idx, jdx], basis.comm_kpts) + end + end + hubbard_n = symmetrize_hubbard_n(basis.model, manifold, hubbard_n; basis.symmetries) +end + +""" +This function reshapes for each kpoint the projectors matrix to a vector of matrices, +taking only the columns corresponding to orbitals in the manifold and splitting them +into different matrices, one for each atom. Columns in the same matrix differ only in +the value of the magnetic quantum number m of the corresponding orbitals. +""" +function reshape_hubbard_proj(projectors::Vector{Matrix{Complex{T}}}, + labels, manifold) where {T} + natoms = length(manifold.iatoms) + l = manifold.l + @assert all(label -> label.l == manifold.l, labels) "$(labels)" + @assert length(labels) == natoms * (2l+1) + p_I = [Vector{Matrix{Complex{T}}}(undef, natoms) for _ = 1:length(projectors)] + for (idx, iatom) in enumerate(manifold.iatoms) + for i = 1:2l+1:length(labels) + iatom != labels[i].iatom && continue + for (ik, projk) in enumerate(projectors) + p_I[ik][idx] = projk[:, i:i+2l] + end + end + end + + p_I +end diff --git a/src/terms/terms.jl b/src/terms/terms.jl index 11b9b77970..292ad59454 100644 --- a/src/terms/terms.jl +++ b/src/terms/terms.jl @@ -64,6 +64,7 @@ include("ewald.jl") include("psp_correction.jl") include("entropy.jl") include("pairwise.jl") +include("hubbard.jl") include("magnetic.jl") breaks_symmetries(::Magnetic) = true diff --git a/src/terms/xc.jl b/src/terms/xc.jl index b9ead3cdd3..fae8429dd5 100644 --- a/src/terms/xc.jl +++ b/src/terms/xc.jl @@ -84,7 +84,7 @@ function xc_potential_real(term::TermXc, basis::PlaneWaveBasis{T}, ψ, occupatio max_ρ_derivs = maximum(max_required_derivative, term.functionals) density = LibxcDensities(basis, max_ρ_derivs, ρ, τ) - if !isnothing(term.ρcore) && !isnothing(τ) + if !isnothing(term.ρcore) && needs_τ(term) negative_α = @views any(1:n_spin) do iσ # α = (τ - τ_W) / τ_unif should be positive with τ_W = |∇ρ|² / 8ρ # equivalently, check 8ρτ - |∇ρ|² ≥ 0 diff --git a/src/workarounds/forwarddiff_rules.jl b/src/workarounds/forwarddiff_rules.jl index c9295085f9..a11f98cee1 100644 --- a/src/workarounds/forwarddiff_rules.jl +++ b/src/workarounds/forwarddiff_rules.jl @@ -275,7 +275,8 @@ end # This has to be changed whenever the scfres structure changes (; ham, basis=basis_dual, energies, ρ, eigenvalues, occupation, εF, ψ, - scfres.τ, # TODO make τ also differentiable for meta-GGA DFPT + # TODO make τ and hubbard_n also differentiable for meta-GGA/DFT+U DFPT + scfres.τ, scfres.hubbard_n, # non-differentiable metadata: response=getfield.(δresults, :info_gmres), scfres.converged, scfres.occupation_threshold, scfres.α, scfres.n_iter, diff --git a/test/hamiltonian_consistency.jl b/test/hamiltonian_consistency.jl index d633ae4b4f..edbf2f248e 100644 --- a/test/hamiltonian_consistency.jl +++ b/test/hamiltonian_consistency.jl @@ -22,14 +22,16 @@ function test_matrix_repr_operator(hamk, ψk; atol=1e-8) end function test_consistency_term(term; rtol=1e-4, atol=1e-8, ε=1e-6, kgrid=[1, 2, 3], - kshift=[0, 1, 0]/2, lattice=testcase.lattice, Ecut=10, - spin_polarization=:none) + kshift=[0, 1, 0]/2, lattice=testcase.lattice, + atom=nothing, Ecut=10, spin_polarization=:none) sspol = spin_polarization != :none ? " ($spin_polarization)" : "" xc = term isa Xc ? "($(first(term.functionals)))" : "" @testset "$(typeof(term))$xc $sspol" begin n_dim = 3 - count(iszero, eachcol(lattice)) - Si = n_dim == 3 ? ElementPsp(14, load_psp(testcase.psp_gth)) : ElementCoulomb(:Si) - atoms = [Si, Si] + if isnothing(atom) + atom = n_dim == 3 ? ElementPsp(14, load_psp(testcase.psp_gth)) : ElementCoulomb(:Si) + end + atoms = [atom, atom] model = Model(lattice, atoms, testcase.positions; terms=[term], spin_polarization, symmetries=true) basis = PlaneWaveBasis(model; Ecut, kgrid=MonkhorstPack(kgrid; kshift)) @@ -47,7 +49,11 @@ function test_consistency_term(term; rtol=1e-4, atol=1e-8, ε=1e-6, kgrid=[1, 2, compute_density(basis, ψ, occupation) end τ = compute_kinetic_energy_density(basis, ψ, occupation) - E0, ham = energy_hamiltonian(basis, ψ, occupation; ρ, τ) + hubbard_n = nothing + if term isa Hubbard + hubbard_n = DFTK.compute_hubbard_n(only(basis.terms), basis, ψ, occupation) + end + E0, ham = energy_hamiltonian(basis, ψ, occupation; ρ, τ, hubbard_n) @assert length(basis.terms) == 1 @@ -58,7 +64,14 @@ function test_consistency_term(term; rtol=1e-4, atol=1e-8, ε=1e-6, kgrid=[1, 2, compute_density(basis, ψ_trial, occupation) end τ_trial = compute_kinetic_energy_density(basis, ψ_trial, occupation) - (; energies) = energy_hamiltonian(basis, ψ_trial, occupation; ρ=ρ_trial, τ=τ_trial) + hubbard_n_trial = nothing + if term isa Hubbard + hubbard_n_trial = DFTK.compute_hubbard_n(only(basis.terms), basis, + ψ_trial, occupation) + end + (; energies) = energy_hamiltonian(basis, ψ_trial, occupation; + ρ=ρ_trial, τ=τ_trial, + hubbard_n=hubbard_n_trial) energies.total end diff = (compute_E(ε) - compute_E(-ε)) / (2ε) @@ -73,6 +86,9 @@ function test_consistency_term(term; rtol=1e-4, atol=1e-8, ε=1e-6, kgrid=[1, 2, end diff_predicted = mpi_sum(diff_predicted, basis.comm_kpts) + # Make sure that we don't accidentally test 0 == 0 + @test abs(diff) > atol + err = abs(diff - diff_predicted) @test err < rtol * abs(E0.total) || err < atol end @@ -83,25 +99,39 @@ end @testitem "Hamiltonian consistency" setup=[TestCases, HamConsistency] tags=[:dont_test_mpi] begin using DFTK using LinearAlgebra - using .HamConsistency: test_consistency_term + using .HamConsistency: test_consistency_term, testcase test_consistency_term(Kinetic()) - test_consistency_term(AtomicLocal()) - test_consistency_term(AtomicNonlocal()) test_consistency_term(ExternalFromReal(X -> cos(X[1]))) test_consistency_term(ExternalFromFourier(X -> abs(norm(X)))) test_consistency_term(LocalNonlinearity(ρ -> ρ^2)) test_consistency_term(Hartree()) - test_consistency_term(Ewald()) - test_consistency_term(PspCorrection()) - test_consistency_term(Xc([:lda_xc_teter93])) - test_consistency_term(Xc([:lda_xc_teter93]), spin_polarization=:collinear) - test_consistency_term(Xc([:gga_x_pbe]), spin_polarization=:collinear) - test_consistency_term(Xc([:mgga_x_tpss])) - test_consistency_term(Xc([:mgga_x_scan])) - test_consistency_term(Xc([:mgga_c_scan]), spin_polarization=:collinear) - test_consistency_term(Xc([:mgga_x_b00])) - test_consistency_term(Xc([:mgga_c_b94]), spin_polarization=:collinear) + let + Si = ElementPsp(14, load_psp(testcase.psp_upf)) + test_consistency_term(Hubbard(OrbitalManifold(Si, [1, 2], "3P"), 0.01), atom=Si) + test_consistency_term(Hubbard(OrbitalManifold(Si, [1, 2], "3P"), 0.01), atom=Si, + spin_polarization=:collinear) + end + # Disabled since the energy is constant, and the test guards against 0 differences + # test_consistency_term(Ewald()) + # test_consistency_term(PspCorrection()) + for psp in [testcase.psp_gth, testcase.psp_upf] + Si = ElementPsp(14, load_psp(psp)) + test_consistency_term(AtomicLocal(), atom=Si) + test_consistency_term(AtomicNonlocal(), atom=Si) + test_consistency_term(Xc([:lda_xc_teter93]), atom=Si) + test_consistency_term(Xc([:lda_xc_teter93]), atom=Si, spin_polarization=:collinear) + test_consistency_term(Xc([:gga_x_pbe]), atom=Si, spin_polarization=:collinear) + # TODO: for use_nlcc=true need to fix consistency for meta-GGA with NLCC + # (see JuliaMolSim/DFTK.jl#1180) + test_consistency_term(Xc([:mgga_x_tpss]; use_nlcc=false), atom=Si) + test_consistency_term(Xc([:mgga_x_scan]; use_nlcc=false), atom=Si) + test_consistency_term(Xc([:mgga_c_scan]; use_nlcc=false), atom=Si, + spin_polarization=:collinear) + test_consistency_term(Xc([:mgga_x_b00]; use_nlcc=false), atom=Si) + test_consistency_term(Xc([:mgga_c_b94]; use_nlcc=false), atom=Si, + spin_polarization=:collinear) + end let a = 6 diff --git a/test/hubbard.jl b/test/hubbard.jl new file mode 100644 index 0000000000..ab246788b6 --- /dev/null +++ b/test/hubbard.jl @@ -0,0 +1,99 @@ +@testitem "Test Wigner matrices" setup=[TestCases] begin + using DFTK + using PseudoPotentialData + using LinearAlgebra + + @testset "Wigner Identity" begin + # Identity + Id = [1.0 0 0; 0 1 0; 0 0 1] + D = DFTK.wigner_d_matrix(1, Id) + @test D ≈ I + D = DFTK.wigner_d_matrix(2, Id) + @test D ≈ I + end + @testset "Wigner Inversion" begin + # This reverts all p orbitals, sends all d orbitals in themselves + Inv = -[1.0 0 0; 0 1 0; 0 0 1] + D = DFTK.wigner_d_matrix(1, Inv) + @test D ≈ -I + D = DFTK.wigner_d_matrix(2, Inv) + @test D ≈ I + end + @testset "Wigner invert x and y" begin + # This keeps pz, dz2, dx2-y2 and dxy unchanged, changes sign to all others + A3 = [1.0 0 0; 0 -1 0; 0 0 -1] + D3p = [-1.0 0 0; 0 -1 0; 0 0 1] + D3d = [-1.0 0 0 0 0; 0 1 0 0 0; 0 0 1 0 0; 0 0 0 -1 0; 0 0 0 0 1] + D = DFTK.wigner_d_matrix(1, A3) + @test D ≈ D3p + D = DFTK.wigner_d_matrix(2, A3) + @test D ≈ D3d + end + @testset "Wigner swap x and y" begin + # This sends: px <-> py, dxz <-> dyz, dx2-y2 -> -(dx2-y2) and keeps the other fixed + A3 = [0.0 1 0; 1 0 0; 0 0 1] + D3p = [0.0 0 1; 0 1 0; 1 0 0] + D3d = [1.0 0 0 0 0; 0 0 0 1 0; 0 0 1 0 0; 0 1 0 0 0; 0 0 0 0 -1] + D = DFTK.wigner_d_matrix(1, A3) + @test D ≈ D3p + D = DFTK.wigner_d_matrix(2, A3) + @test D ≈ D3d + end +end + +@testitem "Test Hubbard U term in Nickel Oxide" setup=[TestCases] begin + using DFTK + using PseudoPotentialData + using Unitful + using UnitfulAtomic + using LinearAlgebra + + a = 7.9 # Bohr + lattice = a * [[ 1.0 0.5 0.5]; + [ 0.5 1.0 0.5]; + [ 0.5 0.5 1.0]] + pseudopotentials = PseudoFamily("dojo.nc.sr.pbe.v0_4_1.standard.upf") + Ni = ElementPsp(:Ni, pseudopotentials) + O = ElementPsp(:O, pseudopotentials) + atoms = [Ni, O, Ni, O] + positions = [[0.0, 0.0, 0.0], + [0.25, 0.25, 0.25], + [0.5, 0.5, 0.5], + [0.75, 0.75, 0.75]] + magnetic_moments = [2, 0, -1, 0] + + # Hubbard parameters + U = 10u"eV" + manifold = OrbitalManifold(atoms, Ni, "3D") + + model = model_DFT(lattice, atoms, positions; + extra_terms=[Hubbard(manifold, U)], + temperature=0.01, functionals=PBE(), + smearing=DFTK.Smearing.Gaussian(), magnetic_moments=magnetic_moments) + basis = PlaneWaveBasis(model; Ecut = 15, kgrid = [2, 2, 2]) + ρ0 = guess_density(basis, magnetic_moments) + scfres = self_consistent_field(basis; tol=1e-10, ρ=ρ0) + + @testset "Test Energy results" begin + # The reference values are obtained with first released version + # of the Hubbard code and are in good agreement with Quantum Espresso + ref = -354.907446880021 + e_total = scfres.energies.total + @test e_total ≈ ref + ref_hub = 0.17629078433258719 + @test scfres.energies.Hubbard ≈ ref_hub + end + # The unfolding of the kpoints is not supported with MPI + if mpi_nprocs() == 1 + @testset "Test symmetry consistency" begin + n_hub = scfres.hubbard_n + scfres_nosym = DFTK.unfold_bz(scfres) + term_idx = findfirst(term -> isa(term, DFTK.TermHubbard), + scfres_nosym.basis.terms) + term_hub = scfres_nosym.basis.terms[term_idx] + nhub_nosym = DFTK.compute_hubbard_n(term_hub, scfres_nosym.basis, + scfres_nosym.ψ, scfres_nosym.occupation) + @test n_hub ≈ nhub_nosym + end + end +end