Add Hubbard corrections and DFT+U

src/DFTK.jl

@@ -111,6 +111,7 @@ export AtomicNonlocal
 export Ewald
 export PspCorrection
 export Entropy
+export Hubbard
 export Magnetic
 export PairwisePotential
 export Anyonic

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,9 @@ 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
+    occupation = scfres.occupation
+    compute_bands(scfres.basis, kgrid; scfres.ρ, τ, nhubbard, 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 !iszero(size(psp.r2_pswfcs[1], 1)) "FATAL ERROR: No Atomic projector found within the provided PseudoPotential."
         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,21 @@ Overview of parameters:
         if any(needs_τ, basis.terms)
             τ = compute_kinetic_energy_density(basis, ψ, occupation)
         end
+        for (iterm, term) in enumerate(basis.terms)
+            if typeof(term)==DFTK.TermHubbard{Vector{Matrix{ComplexF64}}, Vector{Any}}
+                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 +215,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 +225,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]