Add Hubbard corrections and DFT+U

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"

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",

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)
 

examples/hubbard.jl

@@ -0,0 +1,58 @@
+# # 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=32, kgrid=[2, 2, 2])
+scfres = self_consistent_field(basis; tol=1e-10, ρ=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
+U = 10u"eV"
+manifold = OrbitalManifold(; species=:Ni, label="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=32, kgrid=[2, 2, 2] )
+scfres = self_consistent_field(basis; tol=1e-10, ρ=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)

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]
@@ -158,12 +159,12 @@ function plot_pdos(basis::PlaneWaveBasis{T}, eigenvalues, ψ; iatom, label=nothi
     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)])
+                         [OrbitalManifold(;iatom, 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[σ])

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

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 = 2*l+2 # We need at least 2*l+1 equations, we add one for numerical stability
+    B = Matrix{T}(undef, 2*l+1, neq)
+    A = Matrix{T}(undef, 2*l+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] = DFTK.ylm_real(l, m, r0)
+            A[m+l+1, n] = DFTK.ylm_real(l, m, r)
+        end
+    end
+    κ = cond(A)
+    @assert κ < 100.0 "The Wigner matrix computation is badly conditioned. κ(A)=$(κ)"
+    B / A
+end

src/postprocess/band_structure.jl

@@ -15,6 +15,7 @@ All kwargs not specified below are passed to [`diagonalize_all_kblocks`](@ref):
                                kgrid::Union{AbstractKgrid,AbstractKgridGenerator};
                                n_bands=default_n_bands_bandstructure(basis.model),
                                n_extra=3, ρ=nothing, τ=nothing, εF=nothing,
+                               occupation=nothing, nhubbard=nothing,
                                eigensolver=lobpcg_hyper, tol=1e-3, seed=nothing,
                                kwargs...)
     # kcoords are the kpoint coordinates in fractional coordinates
@@ -36,7 +37,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; ρ, τ, nhubbard, occupation)
     eigres = diagonalize_all_kblocks(eigensolver, ham, n_bands + n_extra;
                                      n_conv_check=n_bands, tol, kwargs...)
     if !eigres.converged
@@ -68,7 +69,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...)
+    nhubbard = haskey(scfres, :nhubbard) ? scfres.nhubbard : nothing
+    compute_bands(scfres.basis, kgrid; 
+                  scfres.ρ, τ, nhubbard, scfres.occupation, 
+                  scfres.εF, n_bands, kwargs...)
 end
 
 """

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])   
 
@@ -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))

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,
+    nhubbard=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, τ, nhubbard) = 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, τ, nhubbard, eigenvalues, εF)
 
         # Diagonalize `ham` to get the new state
         nextstate = next_density(ham, nbandsalg, fermialg; eigensolver, ψ, eigenvalues,
@@ -180,16 +181,22 @@ Overview of parameters:
         if any(needs_τ, basis.terms)
             τ = compute_kinetic_energy_density(basis, ψ, occupation)
         end
+        for term in basis.terms
+            if isa(term, DFTK.TermHubbard)
+                nhubbard = compute_nhubbard(term.manifold, basis, ψ, occupation;
+                                            projectors=term.P, labels=term.labels)
+            end
+        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, τ, nhubbard, α=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, τ, nhubbard, eigenvalues, εF)
         end
         history_Etot = vcat(info.history_Etot, energies.total)
         history_Δρ = vcat(info.history_Δρ, norm(Δρ) * sqrt(basis.dvol))
@@ -211,7 +218,7 @@ Overview of parameters:
         ρnext, info_next
     end
 
-    info_init = (; ρin=ρ, τ, ψ, occupation=nothing, eigenvalues=nothing, εF=nothing,
+    info_init = (; ρin=ρ, τ, nhubbard, ψ, occupation=nothing, eigenvalues=nothing, εF=nothing,
                    n_iter=0, n_matvec=0, timedout=false, converged=false,
                    history_Etot=T[], history_Δρ=T[])
 
@@ -221,13 +228,13 @@ 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, τ, nhubbard, ψ, occupation, eigenvalues, εF, converged) = info
+    energies, ham = energy_hamiltonian(basis, ψ, occupation; ρ=ρout, τ, nhubbard, 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.n_iter, info.n_matvec, ψ, nhubbard, info.diagonalization, stage=:finalize,
                 info.history_Δρ, info.history_Etot, info.timedout, mixing, seed,
                 runtime_ns=time_ns() - start_ns, algorithm="SCF")
     callback(scfres)

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_nhubbard(nhubbard::Array{Matrix{Complex{T}}}, 
+                             model, symmetries, positions) 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!
+
+    nspins = size(nhubbard, 1)
+    natoms = size(nhubbard, 2)
+    nsym = length(symmetries)
+    # We extract the l value from the manifold size per atom (2l+1)
+    l = Int((size(nhubbard[1, 1, 1], 1)-1)/2)
+    ldim = 2*l+1
+
+    # Initialize the nhubbard 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(l, Wcart)
+        for σ in 1:nspins, iatom in 1:natoms
+            sym_atom = find_symmetry_preimage(positions, positions[iatom], symmetry)
+            ns[σ, iatom, iatom] .+= WigD' * nhubbard[σ, sym_atom, sym_atom] * WigD
+        end
+    end
+    ns ./= nsym
+    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.nhubbard, eigenvalues, scfres.εF)
     @assert energies.total ≈ scfres.energies.total
     new_scfres = (; basis=basis_unfolded, ψ, ham, eigenvalues, occupation)
     merge(scfres, new_scfres)