Add Hubbard corrections and DFT+U

docs/make.jl

@@ -38,6 +38,7 @@ PAGES = [
         "examples/convergence_study.jl",
         "examples/pseudopotentials.jl",
         "examples/supercells.jl",
+        "examples/hubbard.jl",
         "examples/gaas_surface.jl",
         "examples/graphene.jl",
         "examples/geometry_optimization.jl",

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 Printf
+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, O = ElementPsp(:Ni, pseudopotentials), 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 NSCF without the Hubbard term
+model = model_DFT(lattice, atoms, positions; temperature=5e-3,
+                  functionals = PBE(), magnetic_moments=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(2, 2, 2); ρ=scfres.ρ);
+band_gap = bands.eigenvalues[1][25] - bands.eigenvalues[1][24]
+
+# 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, temperature=2e-3, 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 = DFTK.OrbitalManifold(;species=:Ni, label="3D")
+
+# Run SCF
+model = model_DFT(lattice, atoms, positions; extra_terms=[DFTK.Hubbard(manifold, U)],
+                  functionals = PBE(), temperature=5e-3, magnetic_moments=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 NSCF
+bands_hub = compute_bands(scfres, MonkhorstPack(2, 2, 2); ρ=scfres.ρ)
+band_gap = bands_hub.eigenvalues[1][26] - bands_hub.eigenvalues[1][25]
+
+# 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; 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,7 +159,7 @@ 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;

src/DFTK.jl

@@ -111,6 +111,7 @@ export AtomicNonlocal
 export Ewald
 export PspCorrection
 export Entropy
+export Hubbard
 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 hcat(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, kwargs...)
     # kcoords are the kpoint coordinates in fractional coordinates
     if isnothing(ρ)
@@ -34,7 +35,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
@@ -66,7 +67,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
+        count_n_pswfc(psp) # 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),
@@ -157,12 +158,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,
@@ -178,16 +179,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).nhubbard
+            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,
                        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))
@@ -209,7 +216,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[])
 
@@ -219,13 +226,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,
                 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]