Add Hubbard corrections and DFT+U

examples/collinear_magnetism.jl

@@ -109,7 +109,7 @@ plot_dos(bands_666)
 # is sufficient.
 # We can clearly see that the origin of this spin-polarization traces back to the 
 # 3D orbital contribution if we look at the corresponding projected density of states (PDOS).
-plot_pdos(bands_666; iatoms=[1], label="3D")
+plot_pdos(bands_666; iatom=1, label="3D")
 
 # Similarly the band structure shows clear differences between both spin components.
 using Unitful

examples/dos.jl

@@ -28,5 +28,5 @@ plot_dos(scfres)
 plot_ldos(scfres; n_points=100, ldos_xyz=[:, 10, 10])
 
 # Plot the projected DOS
-p = plot_pdos(scfres; iatoms=[1], label="3S", εrange=(-0.3, 0.5))
-plot_pdos(scfres; p, colors=[:red], iatoms=[1], label="3P", εrange=(-0.3, 0.5))
+p = plot_pdos(scfres; iatom=1, label="3S", εrange=(-0.3, 0.5))
+plot_pdos(scfres; p, colors=[:red], iatom=1, label="3P", εrange=(-0.3, 0.5))

examples/hubbard.jl

@@ -21,7 +21,8 @@ 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)
+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))
@@ -33,11 +34,17 @@ band_gap = bands.eigenvalues[1][lowest_unocc_band] - bands.eigenvalues[1][lowest
 width = 5.0u"eV"
 εrange = (εF - austrip(width), εF + austrip(width))
 p = plot_dos(bands; εrange, colors=[:red, :red])
-plot_pdos(bands; p, iatoms=[1], label="3D", colors=[:yellow, :orange], εrange)
+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
+# 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 a manifold in the mathematical sense.
 U = 10u"eV"
-manifold = OrbitalManifold([1,3], Ni.psp, 2, 1)
+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)],
@@ -55,4 +62,4 @@ band_gap = bands_hub.eigenvalues[1][lowest_unocc_band] - bands_hub.eigenvalues[1
 ε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, iatoms=[1], label="3D", colors=[:green, :purple], εrange)
+plot_pdos(bands_hub; p, iatom=1, label="3D", colors=[:green, :purple], εrange)

ext/DFTKPlotsExt.jl

@@ -142,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, ψ; iatoms=nothing, 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,
@@ -155,17 +155,17 @@ function plot_pdos(basis::PlaneWaveBasis{T}, eigenvalues, ψ; iatoms=nothing, la
     n_spin = basis.model.n_spin_components
     to_unit = ustrip(auconvert(unit, 1.0))
 
-    species = isnothing(iatoms) ? "all atoms" : "atoms $(iatoms) ($(basis.model.atoms[iatoms].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())
     p = Plots.plot(p; kwargs...)
     spinlabels = spin_components(basis.model)
-    labels = basis.terms[findfirst(term -> isa(term, DFTK.TermHubbard), basis.terms)].labels
     pdos = DFTK.sum_pdos(compute_pdos(εs, basis, ψ, eigenvalues;
-                                      positions, temperature, smearing), 
-                         [OrbitalManifold(basis.model.atoms, labels; iatoms, 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[σ])

src/postprocess/band_structure.jl

@@ -71,7 +71,7 @@ function compute_bands(scfres::NamedTuple,
     τ = haskey(scfres, :τ) ? scfres.τ : nothing
     hubbard_n = haskey(scfres, :hubbard_n) ? scfres.hubbard_n : nothing
     compute_bands(scfres.basis, kgrid; 
-                  scfres.ρ, τ, hubbard_n, scfres.occupation, 
+                  scfres.ρ, τ, hubbard_n, scfres.occupation,
                   scfres.εF, n_bands, kwargs...)
 end
 

src/postprocess/dos.jl

@@ -136,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
@@ -212,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(is_on_manifold(orb, manifold) for manifold in manifolds)
+            if any(filt(orb) for filt in projector_filters)
                 pdos[:, σ] += pdos_res.pdos[:, j, σ]
             end
         end

src/pseudo/PspUpf.jl

@@ -175,14 +175,6 @@ count_n_pswfc_radial(psp::PspUpf, l) = length(psp.r2_pswfcs[l+1])
 
 pswfc_label(psp::PspUpf, i, l) = psp.pswfc_labels[l+1][i]
 
-function pswfc_indices(psp::PspUpf, label)
-    for l in 0:psp.lmax, n in 1:DFTK.count_n_pswfc_radial(psp, l)
-        if (DFTK.pswfc_label(psp, n, l) == label)
-            return (; l, n)
-        end
-    end
-end
-
 function eval_psp_pswfc_real(psp::PspUpf, i, l, r::T)::T where {T<:Real}
     psp.r2_pswfcs_interp[l+1][i](r) / r^2  # TODO if r is below a threshold, return zero
 end

src/scf/self_consistent_field.jl

@@ -181,11 +181,11 @@ 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)
-                hubbard_n = compute_hubbard_n(term.manifold, basis, ψ, occupation;
-                                              projectors=term.P, labels=term.labels)
-            end
+        ihubbard = findfirst(t -> t isa TermHubbard, basis.terms)
+        if !isnothing(ihubbard)
+            term = basis.terms[ihubbard]
+            hubbard_n = compute_hubbard_n(term.manifold, basis, ψ, occupation;
+                                          projectors=term.P, labels=term.labels)
         end
 
         # Update info with results gathered so far
@@ -233,7 +233,7 @@ Overview of parameters:
 
     # Callback is run one last time with final state to allow callback to clean up
     scfres = (; ham, basis, energies, converged, nbandsalg.occupation_threshold,
-                ρ=ρout, τ, hubbard_n, α=damping, eigenvalues, occupation, εF, 
+                ρ=ρ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")

src/symmetry.jl

@@ -425,23 +425,23 @@ end
 """
 Symmetrize the Hubbard occupation matrix according to the l quantum number of the manifold.
 """
-function symmetrize_hubbard_n(hubbard_n::Array{Matrix{Complex{T}}}, 
-                              model, symmetries, positions) where {T}
+function symmetrize_hubbard_n(manifold::OrbitalManifold,
+                              hubbard_n::Array{Matrix{Complex{T}}},
+                              model, symmetries) 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)
     nsym = length(symmetries)
-    # We extract the l value from the manifold size per atom (2l+1)
-    l = div(size(hubbard_n[1, 1, 1], 1)-1,2)
-    ldim = 2*l+1
+    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(l, Wcart)
+        WigD = wigner_d_matrix(manifold.l, Wcart)
         for σ in 1:nspins, iatom in 1:natoms
             sym_atom = find_symmetry_preimage(positions, positions[iatom], symmetry)
             ns[σ, iatom, iatom] .+= WigD' * hubbard_n[σ, sym_atom, sym_atom] * WigD