Various small fixes and interface improvements (#1179)

Author: mfherbst

docs/src/features.md

@@ -32,7 +32,8 @@ and runs out of the box on Linux, windows and macOS, see [Installation](@ref).
 ## Response and response properties
 - Density-functional perturbation theory (DFPT)
 - Integration of DFPT with **algorithmic differentiation**,
-  e.g. [Polarizability using automatic differentiation](@ref)
+  e.g. [Elastic constants](@ref),
+  [Polarizability using automatic differentiation](@ref)
 - [Phonon computations](@ref) *(preliminary implementation)*
 
 ## Unique features

docs/src/guide/introductory_resources.md

@@ -67,21 +67,22 @@ see [Publications](@ref).
   Covers topics such as DFT, pseudos, SCF, response, ...
 
 ## Recordings
+- [Algorithmic differentiation (AD) for plane-wave DFT](https://www.youtube.com/watch?v=g6j1beYSWV4) by M. F. Herbst:
+  45-min talk at the Institute for Pure and Applied Mathematics (IPAM)
+  at UCLA discussing the algorithmic differentiation techniques in DFTK.
+
 - [DFTK.jl: 5 years of a multidisciplinary electronic-structure code](https://www.youtube.com/watch?v=ox_j2zKOuIk) by M. F. Herbst:
   30-min talk at JuliaCon 2024 providing the state of DFTK 5 years after the project was started.
   [Slides](https://michael-herbst.com/talks/2024.07.12_5years_DFTK.pdf),
   [Pluto notebook](https://michael-herbst.com/talks/2024.07.12_5years_DFTK.html)
 
-- [Julia for Materials Modelling](https://www.youtube.com/watch?v=dujepKxxxkg) by M. F. Herbst:
+- [Julia for Materials Modelling](https://www.youtube.com/watch?v=dujepKxxxkg)
+  by M. F. Herbst (from 2023):
   One-hour talk providing an overview of materials modelling tools for Julia.
   Key features of DFTK are highlighted as part of the talk.
+  Starts to become a little outdated
   [Pluto notebooks](https://mfherbst.github.io/julia-for-materials/)
 
-- [DFTK: A Julian approach for simulating electrons in solids](https://www.youtube.com/watch?v=-RomkxjlIcQ) by M. F. Herbst:
-  Pre-recorded talk for JuliaCon 2020.
-  Assumes no knowledge about DFT and gives the broad picture of DFTK. Starts to become a little outdated.
-  [Slides](https://michael-herbst.com/talks/2020.07.29_juliacon_dftk.pdf).
-
 - [Juliacon 2021 DFT workshop](https://www.youtube.com/watch?v=HvpPMWVm8aw) by M. F. Herbst:
   Three-hour workshop session at the 2021 Juliacon providing a mathematical look on
   DFT, SCF solution algorithms as well as the integration of DFTK into the Julia

src/elements.jl

@@ -35,6 +35,9 @@ pseudofamily(::Element) = nothing
 has_core_density(::Element) = false
 # The preceding functions are fallback implementations that should be altered as needed.
 
+eval_psp_energy_correction(T, ::Element) = zero(T)
+eval_psp_energy_correction(psp::Element) = eval_psp_energy_correction(Float64, psp)
+
 # Fall back to the Gaussian table for Elements without pseudopotentials
 function valence_charge_density_fourier(el::Element, p::T)::T where {T <: Real}
     gaussian_valence_charge_density_fourier(el, p)
@@ -152,6 +155,7 @@ AtomsBase.mass(el::ElementPsp)    = el.mass
 AtomsBase.species(el::ElementPsp) = el.species
 charge_ionic(el::ElementPsp)      = charge_ionic(el.psp)
 has_core_density(el::ElementPsp)  = has_core_density(el.psp)
+eval_psp_energy_correction(T, el::ElementPsp) = eval_psp_energy_correction(T, el.psp)
 
 function local_potential_fourier(el::ElementPsp, p::T) where {T <: Real}
     p == 0 && return zero(T)  # Compensating charge background
@@ -237,6 +241,7 @@ function local_potential_fourier(el::ElementCohenBergstresser, p::T) where {T <:
     psq_pi = Int(round(p^2 / (2π / el.lattice_constant)^2, digits=2))
     T(get(el.V_sym, psq_pi, 0.0))
 end
+# TODO Strictly speaking needs a eval_psp_energy_correction
 
 
 #
@@ -269,3 +274,4 @@ end
 function local_potential_fourier(el::ElementGaussian, p::Real)
     -el.α * exp(- (p * el.L)^2 / 2)  # = ∫_ℝ³ V(x) exp(-ix⋅p) dx
 end
+# TODO Strictly speaking needs a eval_psp_energy_correction

src/pseudo/NormConservingPsp.jl

@@ -18,7 +18,7 @@ abstract type NormConservingPsp end
 # eval_psp_projector_fourier(psp, i, l, p::Real)
 # eval_psp_local_real(psp, r::Real)
 # eval_psp_local_fourier(psp, p::Real)
-# eval_psp_energy_correction(T::Type, psp, n_electrons::Integer)
+# eval_psp_energy_correction(T::Type, psp)
 
 #### Optional methods:
 # eval_psp_density_valence_real(psp, r::Real)
@@ -86,15 +86,15 @@ eval_psp_local_fourier(psp::NormConservingPsp, p::AbstractVector) =
     eval_psp_local_fourier(psp, norm(p))
 
 @doc raw"""
-    eval_psp_energy_correction([T=Float64,] psp, n_electrons)
+    eval_psp_energy_correction([T=Float64,] psp::NormConservingPsp)
+    eval_psp_energy_correction([T=Float64,] element::Element)
 
-Evaluate the energy correction to the Ewald electrostatic interaction energy of one unit
-cell, which is required compared the Ewald expression for point-like nuclei. `n_electrons`
-is the number of electrons per unit cell. This defines the uniform compensating background
-charge, which is assumed here.
-
-Notice: The returned result is the *energy per unit cell* and not the energy per volume.
-To obtain the latter, the caller needs to divide by the unit cell volume.
+Evaluate the energy correction to the Ewald electrostatic interaction energy per unit
+of uniform negative charge. This is the correction required to account for the fact that
+the Ewald expression assumes point-like nuclei and not nuclei of the shape induced by
+the pseudopotential. The compensating background charge assumed for this expression is
+scaled to ``1``. Therefore multiplying by the number of electrons and dividing by the unit
+cell volume yields the energy correction per volume for the DFT simulation.
 
 The energy correction is defined as the limit of the Fourier-transform of the local
 potential as ``p \to 0``, using the same correction as in the Fourier-transform of the local
@@ -103,12 +103,14 @@ potential:
 \lim_{p \to 0} 4π N_{\rm elec} ∫_{ℝ_+} (V(r) - C(r)) \frac{\sin(p·r)}{p·r} r^2 dr + F[C(r)]
  = 4π N_{\rm elec} ∫_{ℝ_+} (V(r) + Z/r) r^2 dr.
 ```
+where as discussed above the implementation is expected to return the result
+for ``N_{\rm elec} = 1``.
 """
 function eval_psp_energy_correction end
 # by default, no correction, see PspHgh implementation and tests
 # for details on what to implement
-eval_psp_energy_correction(psp::NormConservingPsp, n_electrons) =
-    eval_psp_energy_correction(Float64, psp, n_electrons)
+
+eval_psp_energy_correction(psp::NormConservingPsp) = eval_psp_energy_correction(Float64, psp)
 
 """
     eval_psp_density_valence_real(psp, r)

src/pseudo/PspHgh.jl

@@ -218,7 +218,7 @@ function eval_psp_projector_real(psp::PspHgh, i, l, r::T) where {T <: Real}
     sqrt(T(2)) * r^(l + 2(i - 1)) * exp(-r^2 / 2rp^2) / rp^(l + ired) / sqrt(gamma(l + ired))
 end
 
-function eval_psp_energy_correction(T, psp::PspHgh, n_electrons)
+function eval_psp_energy_correction(T, psp::PspHgh)
     # By construction we need to compute the DC component of the difference
     # of the Coulomb potential (-Z/G^2 in Fourier space) and the pseudopotential
     # i.e. -4πZ/(ΔG)^2 -  eval_psp_local_fourier(psp, ΔG) for ΔG → 0. This is:
@@ -228,5 +228,5 @@ function eval_psp_energy_correction(T, psp::PspHgh, n_electrons)
 
     # Multiply by number of electrons and 4π (spherical Hankel prefactor)
     # to get energy per unit cell
-    4T(π) * n_electrons * difference_DC
+    4T(π) * difference_DC
 end

src/pseudo/PspUpf.jl

@@ -232,10 +232,10 @@ function eval_psp_density_core_fourier(psp::PspUpf, p::T) where {T<:Real}
     return hankel(rgrid, r2_ρcore, 0, p)
 end
 
-function eval_psp_energy_correction(T, psp::PspUpf, n_electrons)
+function eval_psp_energy_correction(T, psp::PspUpf)
     rgrid = @view psp.rgrid[1:psp.ircut]
     vloc = @view psp.vloc[1:psp.ircut]
-    4T(π) * n_electrons * simpson(rgrid) do i, r
+    4T(π) * simpson(rgrid) do i, r
         r * (r * vloc[i] - -psp.Zion)
     end
 end

src/response/chi0.jl

@@ -575,6 +575,7 @@ function construct_bandtol(Bandtol::Type, basis::PlaneWaveBasis, ψ, occupation:
     # Distributing the error equally across all k-points leads to (with w = sqrt(Ω / Ng))
     #   ‖z_{nk}‖ ≤ sqrt(Ω / Ng) / (‖K v_i‖ sqrt(Nocck) ‖Re(F⁻¹ Φk)‖_{2,∞} * 2f_{nk} Nk wk)
     # If we bound ‖Re(F⁻¹ Φk)‖_{2,∞} from below this is sqrt(Nocc / Ω).
+    # See also section SM6 and Table SM4 in 2505.02319.
     #
     # Note that the kernel term ||K v_i|| of 2505.02319 is dropped here as it purely arises
     # from the rescaling of the RHS performed in apply_χ0 above. Consequently the function
@@ -598,10 +599,13 @@ function construct_bandtol(Bandtol::Type, scfres::NamedTuple; kwargs...)
 end
 
 function adaptive_bandtol_orbital_term_(::Type{BandtolGuaranteed}, basis, kpt, ψk, mask_k)
-    # Orbital term ‖Re(F⁻¹ Φk)‖_{2,∞} for thik k-point
+    # Orbital term ‖F⁻¹ Φk‖_{2,∞} for thik k-point
+    # Note that compared to the paper we deliberately do not take the real part,
+    # since taking the real part represents an additional approximation
+    # (thus making the strategy less guaranteed)
     row_sums_squared = sum(mask_k) do n
         ψnk_real = @views ifft(basis, kpt, ψk[:, n])
-        abs2.(real.(ψnk_real))
+        abs2.(ψnk_real)
     end
     sqrt(maximum(row_sums_squared))
 end

src/scf/anderson.jl

@@ -103,6 +103,8 @@ and returns ``xₙ₊₁``.
     M = stack(Pfxs) .- vec(Pfxₙ)  # Mᵢⱼ = (Pfxⱼ)ᵢ - (Pfxₙ)ᵢ
     # We need to solve 0 = M' Pfxₙ + M'M βs <=> βs = - (M'M)⁻¹ M' Pfxₙ
 
+    # TODO If memory pressure is high, automatically drop old iterates here
+    #      (Note: This quickly happens in GPU scenarios)
     # Ensure the condition number of M stays below maxcond, else prune the history
     Mfac = qr(M)
     while size(M, 2) > 1 && cond(Mfac.R) > anderson.maxcond

src/scf/mixing.jl

@@ -33,18 +33,20 @@ mix_density(::SimpleMixing, ::PlaneWaveBasis, δF; kwargs...) = δF
 @doc raw"""
 Kerker mixing: ``J^{-1} ≈ \frac{|G|^2}{k_{TF}^2 + |G|^2}``
 where ``k_{TF}`` is the Thomas-Fermi wave vector. For spin-polarized calculations
-by default the spin density is not preconditioned. Unless a non-default value
-for ``ΔDOS_Ω`` is specified. This value should roughly be the expected difference in density
-of states (per unit volume) between spin-up and spin-down.
+by default the spin density is not preconditioned unless a non-default value
+for `ΔDOS_Ω` is specified. This value should roughly be the expected difference in density
+of states (per unit volume) between spin-up and spin-down. Notably setting
+`ΔDOS_Ω = kTF^2 / 4π` disables acting on the ``β`` spin channel completely (as if the
+DOS on ``β`` spin was zero).
 
 Notes:
 - Abinit calls ``1/k_{TF}`` the dielectric screening length (parameter *dielng*)
 """
 @kwdef struct KerkerMixing <: Mixing
     # Default kTF parameter suggested by Kresse, Furthmüller 1996 (kTF=1.5Å⁻¹)
     # DOI 10.1103/PhysRevB.54.11169
-    kTF::Real    = 0.8  # == sqrt(4π (DOS_α + DOS_β)) / Ω
-    ΔDOS_Ω::Real = 0.0  # == (DOS_α - DOS_β) / Ω
+    kTF::Real    = 0.8  # == sqrt(4π (DOS_α + DOS_β) / Ω)
+    ΔDOS_Ω::Real = 0.0  # == (DOS_α - DOS_β) / Ω; set == kTF^2/4π to disable acting on β density
 end
 
 @timing "KerkerMixing" function mix_density(mixing::KerkerMixing, basis::PlaneWaveBasis,

src/terms/operators.jl

@@ -153,7 +153,7 @@ function apply!(Hψ, op::DivAgradOperator, ψ;
         Hψ.fourier .-= im .* G_plus_k[α] .* A∇ψ ./ 2
     end
 end
-# TODO Implement  Matrix(op::DivAgrad)
+# TODO Implement  Matrix(op::DivAgradOperator)
 
 
 # Optimize RFOs by combining terms that can be combined