diff --git a/src/occupation.jl b/src/occupation.jl index bc3fda24dd..341f2c4cfe 100644 --- a/src/occupation.jl +++ b/src/occupation.jl @@ -174,13 +174,13 @@ 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) @@ -188,6 +188,7 @@ function guess_fermi_level_intocc_(basis::PlaneWaveBasis, eigenvalues) # 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 diff --git a/src/response/chi0.jl b/src/response/chi0.jl index 5e60b7dbae..c512b23d26 100644 --- a/src/response/chi0.jl +++ b/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 diff --git a/test/chi0.jl b/test/chi0.jl index 616f7015c5..5be8f88315 100644 --- a/test/chi0.jl +++ b/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