Add Hubbard corrections and DFT+U

src/scf/self_consistent_field.jl

@@ -183,9 +183,7 @@ Overview of parameters:
         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)
+            hubbard_n = compute_hubbard_n(basis.terms[ihubbard], basis, ψ, occupation)
         end
 
         # Update info with results gathered so far

src/terms/hubbard.jl

@@ -68,42 +68,31 @@ function extract_manifold(manifold::OrbitalManifold, projectors, labels)
 end
 
 """
-    compute_hubbard_n(manifold, basis, ψ, occupation; [projectors, labels, positions])
+    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.
+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`.
-
-Overview of inputs:
-- `manifold`: OrbitalManifold with the atomic orbital type to define the Hubbard manifold.
-- `occupation`: Occupation matrix for the bands.
-- `projectors` (kwarg): Vector of projection matrices. For each matrix, each column corresponds
-                        to a different atomic orbital projector, as specified in labels.
-- `labels` (kwarg): Vector of NamedTuples. Each projectors[ik][:,iproj] column has all relevant 
-                    chemical information stored in the corresponding labels[iproj] NamedTuple.
-
-Overview of outputs:
-- `hubbard_n`: 3-tensor of matrices. See above for details.
+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(manifold::OrbitalManifold,
+function compute_hubbard_n(term::TermHubbard,
                            basis::PlaneWaveBasis{T},
-                           ψ, occupation;
-                           projectors, labels) where {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(projectors, labels, manifold)
+    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

test/hamiltonian_consistency.jl

@@ -49,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
 
@@ -60,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ε)
@@ -91,26 +102,36 @@ end
     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())
     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
-    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
         a = 6

test/hubbard.jl

@@ -91,10 +91,8 @@ end
              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(manifold, scfres_nosym.basis, 
-                                                 scfres_nosym.ψ, scfres_nosym.occupation;
-                                                 projectors=term_hub.P,
-                                                 labels=term_hub.labels)
+             nhub_nosym = DFTK.compute_hubbard_n(term_hub, scfres_nosym.basis,
+                                                 scfres_nosym.ψ, scfres_nosym.occupation)
              @test n_hub ≈ nhub_nosym
         end
    end