Refactor n_electrons out of psp correction eval

Author: mfherbst

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/terms/psp_correction.jl

@@ -24,16 +24,11 @@ Compute the correction term for properly modelling the interaction of the pseudo
 core with the compensating background charge induced by the `Ewald` term.
 """
 function energy_psp_correction(lattice::AbstractMatrix{T}, atoms, atom_groups) where {T}
-    psp_groups = [group for group in atom_groups if atoms[first(group)] isa ElementPsp]
-    isempty(psp_groups) && return zero(T)
-
+    correction_per_cell_and_electron = sum(atom_groups) do group
+        length(group) * eval_psp_energy_correction(T, atoms[first(group)])::T
+    end
     n_electrons::Int = n_electrons_from_atoms(atoms)
-    correction_per_cell = sum(
-        length(group) * eval_psp_energy_correction(T, atoms[first(group)].psp, n_electrons)
-        for group in psp_groups
-    )
-
-    correction_per_cell / compute_unit_cell_volume(lattice)
+    correction_per_cell_and_electron * n_electrons / compute_unit_cell_volume(lattice)
 end
 function energy_psp_correction(model::Model)
     energy_psp_correction(model.lattice, model.atoms, model.atom_groups)

test/PspHgh.jl

@@ -195,20 +195,19 @@ end
 
     reg_param = 1e-6  # divergent integral, needs regularization
     p_small = 1e-6    # We are interested in p→0 term
-    function integrand(psp, n_electrons, r)
+    function integrand(psp, r)
         # Difference of potential of point-like atom (what is assumed in Ewald)
         # versus actual structure of the pseudo potential
         coulomb = -psp.Zion / r
-        diff = n_electrons * (eval_psp_local_real(psp, r) - coulomb)
+        diff = (eval_psp_local_real(psp, r) - coulomb)
         4π * diff * exp(-reg_param * r) * sin(p_small*r) / p_small * r
     end
 
-    n_electrons = 20
     family = PseudoFamily("cp2k.nc.sr.lda.v0_1.semicore.gth")
     for element in (:Au, :Ba)
         psp = load_psp(family, element)
-        reference = quadgk(r -> integrand(psp, n_electrons, r), 0, Inf)[1]
-        @test reference ≈ eval_psp_energy_correction(psp, n_electrons) atol=1e-2
+        reference = quadgk(r -> integrand(psp, r), 0, Inf)[1]
+        @test reference ≈ eval_psp_energy_correction(psp) atol=1e-3
     end
 end
 
@@ -224,7 +223,7 @@ end
         psp = load_psp(family, element)
         coulomb = -4π * psp.Zion / p_small^2
         reference = eval_psp_local_fourier(psp, p_small) - coulomb
-        @test reference ≈ eval_psp_energy_correction(psp, 1) atol=1e-3
+        @test reference ≈ eval_psp_energy_correction(psp) atol=1e-3
     end
 end
 

test/PspUpf.jl

@@ -104,9 +104,8 @@ end
     for psp_pair in mPspUpf.gth_pseudos
         upf = psp_pair.upf
         gth = psp_pair.gth
-        n_electrons = 3
-        reference_gth = eval_psp_energy_correction(gth, n_electrons)
-        @test reference_gth ≈ eval_psp_energy_correction(upf, n_electrons) atol=1e-3 rtol=1e-3
+        reference_gth = eval_psp_energy_correction(gth)
+        @test reference_gth ≈ eval_psp_energy_correction(upf) atol=1e-3 rtol=1e-3
     end
 end
 
@@ -199,7 +198,7 @@ end
     for psp in values(mPspUpf.upf_pseudos)
         coulomb = -4π * (psp.Zion) / p_small^2
         reference = eval_psp_local_fourier(psp, p_small) - coulomb
-        @test reference ≈ eval_psp_energy_correction(psp, 1) atol=1e-2
+        @test reference ≈ eval_psp_energy_correction(psp) atol=1e-2
     end
 end