Fix Fermi level variation at zero temperature

src/occupation.jl

@@ -174,20 +174,21 @@ end
 
 
 """
-Estimate a Fermi level by assuming (a) integer occupation and (b) an equal number of bands
-is filled for each k-point. This is largely the setting for occupation without temperature.
-Note, that while this function can be used in cases with spin or with temperature, there
-is no guarantee that the Fermi-level estimated by this function *is* the actual Fermi level.
-Therefore this function should generally only be used as a starting point for other routines.
+Find the HOMO and LUMO levels and their local (k-point, band) indices, assuming
+(a) integer occupation and (b) an equal number of bands is filled for each k-point.
+Returns `(; HOMO, LUMO, ik_HOMO, n_HOMO, ik_LUMO, n_LUMO)`.
+`LUMO`, `ik_LUMO`, `n_LUMO` are `nothing` if no unoccupied bands are available.
+The `ik_*` indices refer to this MPI rank's local k-points.
 """
-function guess_fermi_level_intocc_(basis::PlaneWaveBasis, eigenvalues)
+function find_HOMO_LUMO(basis::PlaneWaveBasis, eigenvalues)
     filled_occ = filled_occupation(basis.model)
     n_spin = basis.model.n_spin_components
     n_fill = div(basis.model.n_electrons, n_spin * filled_occ, RoundUp)
 
     # Highest occupied energy level
     HOMO = maximum([εk[n_fill] for εk in eigenvalues])
     HOMO = mpi_max(HOMO, basis.comm_kpts)
+    ik_HOMO = argmax(ik -> eigenvalues[ik][n_fill], eachindex(eigenvalues))
 
     # Lowest unoccupied energy level: not all k-points might have at least n_fill+1
     # energy levels so we have to take care of that by specifying init to minimum
@@ -197,6 +198,27 @@ function guess_fermi_level_intocc_(basis::PlaneWaveBasis, eigenvalues)
     LUMO = mpi_min(LUMO, basis.comm_kpts)
 
     if LUMO == typemax(HOMO)
+        (; HOMO, LUMO=nothing, ik_HOMO, n_HOMO=n_fill,
+         ik_LUMO=nothing, n_LUMO=nothing)
+    else
+        ik_LUMO = argmin(eachindex(eigenvalues)) do ik
+            minimum(eigenvalues[ik][n_fill+1:end]; init=typemax(HOMO))
+        end
+        n_LUMO = n_fill + argmin(eigenvalues[ik_LUMO][n_fill+1:end])
+        (; HOMO, LUMO, ik_HOMO, n_HOMO=n_fill, ik_LUMO, n_LUMO)
+    end
+end
+
+"""
+Estimate a Fermi level by assuming (a) integer occupation and (b) an equal number of bands
+is filled for each k-point. This is largely the setting for occupation without temperature.
+Note, that while this function can be used in cases with spin or with temperature, there
+is no guarantee that the Fermi-level estimated by this function *is* the actual Fermi level.
+Therefore this function should generally only be used as a starting point for other routines.
+"""
+function guess_fermi_level_intocc_(basis::PlaneWaveBasis, eigenvalues)
+    (; HOMO, LUMO) = find_HOMO_LUMO(basis, eigenvalues)
+    if isnothing(LUMO)
         HOMO + 1  # Just to make sure the εF is a sane number and above HOMO
     else
         (HOMO + LUMO) / 2

src/response/chi0.jl

@@ -333,7 +333,7 @@ function compute_δocc!(δocc, basis::PlaneWaveBasis{T}, ψ, εF, ε, δHψ, δt
         end
         D = mpi_sum(D, basis.comm_kpts)
 
-        if isnothing(model.εF)  # εF === nothing means that Fermi level is fixed by model
+        if isnothing(model.εF)  # εF === nothing means that Fermi level is not fixed by model
             # Compute δεF from δ ∑ occ = 0…
             δocc_tot = mpi_sum(sum(basis.kweights .* sum.(δocc)), basis.comm_kpts)
             δεF = -δocc_tot / D
@@ -344,6 +344,20 @@ function compute_δocc!(δocc, basis::PlaneWaveBasis{T}, ψ, εF, ε, δHψ, δt
                 δocc[ik][n] -= fpnk * δεF / temperature
             end
         end
+    elseif isnothing(model.εF)
+        # Effective insulator (zero temperature or large gap): occupations don't change,
+        # but εF = (HOMO + LUMO) / 2 shifts with the eigenvalues.
+        # We reuse find_HOMO_LUMO to identify the relevant (k, band) indices,
+        # consistent with guess_fermi_level_intocc_.
+        (; ik_HOMO, n_HOMO, ik_LUMO, n_LUMO) = find_HOMO_LUMO(basis, ε)
+        if !isnothing(ik_LUMO)
+            δε_HOMO = real(dot(ψ[ik_HOMO][:, n_HOMO], δHψ[ik_HOMO][:, n_HOMO]))
+            δε_LUMO = real(dot(ψ[ik_LUMO][:, n_LUMO], δHψ[ik_LUMO][:, n_LUMO]))
+            # In MPI, each rank computes δε at its local argmax/argmin;
+            # the global extremum picks the correct derivative.
+            δεF = (mpi_max(δε_HOMO, basis.comm_kpts) +
+                   mpi_min(δε_LUMO, basis.comm_kpts)) / 2
+        end
     end
 
     (; δocc, δεF)
@@ -435,7 +449,7 @@ end
 
 """
 Compute the orbital and occupation changes as a result of applying the ``χ_0`` superoperator
-to the Hamiltonian change `δH` represented by the matrix-vector products `δHψ`. 
+to the Hamiltonian change `δH` represented by the matrix-vector products `δHψ`.
 """
 @views @timing function apply_χ0_4P(ham, ψ, occupation, εF, eigenvalues, δHψ;
                                     δtemperature=zero(eltype(ham.basis)),
@@ -469,16 +483,17 @@ to the Hamiltonian change `δH` represented by the matrix-vector products `δHψ
     bandtol_minus_q_occ = [bandtol_occ[k_to_k_minus_q[ik]] for ik in 1:length(basis.kpoints)]
     @assert bandtolalg.occupation_threshold == occupation_threshold
 
-    # First we compute δoccupation. We only need to do this for the actually occupied
-    # orbitals. So we make a fresh array padded with zeros, but only alter the elements
-    # corresponding to the occupied orbitals. (Note both compute_δocc! and compute_δψ!
-    # assume that the first array argument has already been initialised to zero).
+    # First we compute δoccupation and δεF. We pass the full eigenvalues, ψ and δHψ
+    # (not just the occupied subset) so that compute_δocc! can compute δεF for
+    # effective insulators, where εF = (HOMO + LUMO) / 2 requires LUMO data.
+    # For the metallic/finite-temperature path the extra empty bands contribute ≈ 0
+    # (exponentially small occupation derivatives). Both compute_δocc! and compute_δψ!
+    # assume that the first array argument has already been initialised to zero.
     # For phonon calculations when q ≠ 0, we do not use δoccupation, and compute directly
     # the full perturbation δψ.
     δoccupation = zero.(occupation)
     if iszero(q)
-        δocc_occ = [δoccupation[ik][maskk] for (ik, maskk) in enumerate(mask_occ)]
-        (; δεF) = compute_δocc!(δocc_occ, basis, ψ_occ, εF, ε_occ, δHψ_minus_q_occ, δtemperature)
+        (; δεF) = compute_δocc!(δoccupation, basis, ψ, εF, eigenvalues, δHψ, δtemperature)
     else
         # When δH is not periodic, δH ψnk is a Bloch wave at k+q and ψnk at k,
         # so that δεnk = <ψnk|δH|ψnk> = 0 and there is no occupation shift

test/chi0.jl

@@ -1,5 +1,6 @@
 @testmodule Chi0 begin
 using Test
+using ComponentArrays
 using DFTK
 using DFTK: mpi_mean!
 using LinearAlgebra
@@ -33,7 +34,9 @@ function test_chi0(testcase; symmetries=false, temperature=0, spin_polarization=
         ham0  = energy_hamiltonian(basis, nothing, nothing; ρ=ρ0).ham
         nbandsalg = is_εF_fixed ? FixedBands(; n_bands_converge=6) : AdaptiveBands(model)
         res = DFTK.next_density(ham0, nbandsalg; tol, eigensolver)
-        scfres = (; ham=ham0, res...)
+        scfres = (; basis, ham=ham0, res...)
+
+        bandtolalg = DFTK.BandtolBalanced(scfres)
 
         # create external small perturbation εδV
         n_spin = model.n_spin_components
@@ -46,24 +49,34 @@ function test_chi0(testcase; symmetries=false, temperature=0, spin_polarization=
             @test δV_sym ≈ δV
         end
 
-        function compute_ρ_FD(ε)
+        function compute_finite_perturbation(ε)
             term_builder = basis -> DFTK.TermExternal(ε * δV)
             model = model_DFT(testcase.lattice, testcase.atoms, testcase.positions;
                               functionals, model_kwargs..., extra_terms=[term_builder])
             basis = PlaneWaveBasis(model; basis_kwargs...)
             ham = energy_hamiltonian(basis, nothing, nothing; ρ=ρ0).ham
-            res = DFTK.next_density(ham, nbandsalg; tol, eigensolver)
-            res.ρout
+            (; ρout, occupation, εF) = DFTK.next_density(ham, nbandsalg; tol, eigensolver)
+            ComponentArray(; ρ=ρout, occupation, εF)
         end
 
         # middle point finite difference for more precision
-        ρ1 = compute_ρ_FD(-ε)
-        ρ2 = compute_ρ_FD(ε)
-        diff_findiff = (ρ2 - ρ1) / (2ε)
+        res1 = compute_finite_perturbation(-ε)
+        res2 = compute_finite_perturbation(+ε)
+        diff_findiff = (res2 - res1) / (2ε)
 
         # Test apply_χ0 and compare against finite differences
-        diff_applied_χ0 = apply_χ0(scfres, δV).δρ
-        @test norm(diff_findiff - diff_applied_χ0) < atol
+        diff_applied_χ0 = apply_χ0(scfres, δV)
+        @test norm(diff_findiff.ρ - diff_applied_χ0.δρ) < atol
+
+        # Test occupation and Fermi-level sensitivities
+        (; ψ, occupation, eigenvalues, εF, occupation_threshold) = scfres
+        q = zero(Vec3{eltype(ham0.basis)})
+        δHψ = DFTK.multiply_ψ_by_blochwave(basis, ψ, δV, q)
+        diff_applied_χ0_4P = DFTK.apply_χ0_4P(ham0, ψ, occupation, εF, eigenvalues, δHψ;
+                                              occupation_threshold, bandtolalg, q)
+        maximumabs(x) = maximum(abs, x)
+        @test maximum(maximumabs, diff_applied_χ0_4P.δoccupation - diff_findiff.occupation) < atol
+        @test abs(diff_applied_χ0_4P.δεF - diff_findiff.εF) < atol
 
         # Test apply_χ0 without extra bands
         ψ_occ, occ_occ = DFTK.select_occupied_orbitals(basis, scfres.ψ, scfres.occupation;
@@ -72,12 +85,13 @@ function test_chi0(testcase; symmetries=false, temperature=0, spin_polarization=
 
         diff_applied_χ0_noextra = apply_χ0(scfres.ham, ψ_occ, occ_occ, scfres.εF, ε_occ, δV;
                                            scfres.occupation_threshold).δρ
-        @test norm(diff_applied_χ0_noextra - diff_applied_χ0) < atol
+        @test norm(diff_applied_χ0_noextra - diff_applied_χ0.δρ) < atol
 
         # just to cover it here
         if temperature > 0
             D = compute_dos(res.εF, basis, res.eigenvalues)
             LDOS = compute_ldos(res.εF, basis, res.eigenvalues, res.ψ)
+            @test abs(sum(LDOS) * basis.dvol - sum(D)) < 1e-12
         end
 
         if !symmetries
@@ -86,7 +100,7 @@ function test_chi0(testcase; symmetries=false, temperature=0, spin_polarization=
             if compute_full_χ0
                 χ0 = compute_χ0(ham0)
                 diff_computed_χ0 = reshape(χ0 * vec(δV), basis.fft_size..., n_spin)
-                @test norm(diff_findiff - diff_computed_χ0) < atol
+                @test norm(diff_findiff.ρ - diff_computed_χ0) < atol
             end
 
             # Test that apply_χ0 is self-adjoint