Possible flag-based interface to reduce duplication

[azadoks]

Proposal

While working to integrate PseudoPotentialIO into DFTK, there were a few points where having a dispatch flag-based interface in PPIO would help to reduce code duplication, namely

1. When calculating atomic form-factors for constructing a density superposition, the code is identical except for the function used to calculate the Fourier transform of the atomic density quantity. This affects density guesses and non-linear core-correction.
2. When building atomic orbital-like projection vectors for the non-local potential term, wavefunction guesses, DFT+U, and projected quantities, everything is identical except for the fourier evaluator for the projector (beta / chi / future PAW orbital-like quantities)

Here, I define a set of AbstractPsPQuantity structs:

abstract type AbstractPsPQuantity end

struct PsPProjector <: AbstractPsPQuantity end
struct BetaProjector <: PsPProjector end
struct ChiProjector <: PsPProjector end

struct AtomicChargeDensity <: AbstractPsPQuantity end
struct CoreChargeDensity <: AtomicChargeDensity end
struct ValenceChargeDensity <: AtomicChargeDensity end

struct PsPPotential <: AbstractPsPQuantity end
struct LocalPotential <: PsPPotential end

struct PsPCoupling <: AbstractPsPQuantity end
struct BetaCoupling <: PsPCoupling end
struct AugmentationCoupling <: PsPCoupling end

struct AugmentationFunction <: AbstractPsPQuantity end

abstract type EvaluationSpace end
struct FourierSpace <: EvaluationSpace end
struct RealSpace <: EvaluationSpace end

# TODO
abstract type LocalPotentialCorrection end
struct CoulombCorrection <: LocalPotentialCorrection end
struct ErfCoulombCorrection <: LocalPotentialCorrection end

which are used to select dispatches for a variety of functions:

# Whether the given pseudopotential provides the quantity
has_quantity(::AbstractPsPQuantity, ::AbstractPsP)
# Get the internal representation of a quantity, if available
get_quantity(::AbstractPsPQuantity, ::AbstractPsP)
# Get the cutoff radius (minimum or maximum) of a given quantity
cutoff_radius(::AbstractPsPQuantity, ::AbstractPsP; f=minimum, tol=nothing)
# Number of radial functions
n_radials(::PsPProjector, ::AbstractPsP)
# Number of angular functions = (2l + 1) * n_radials
n_angulars(::PsPProjector, ::AbstractPsP)
# Return a function, interpolator, or other callable which evaluates the quantity
# in the evaluation space
psp_quantity_evaluator(::EvaluationSpace, ::AbstractPsPQuantity, ::AbstractPsP)

This reduces code duplication to some extent in PPIO and should allow for additional quantities to be added without much effort.

Additional features and changes

• Add qe_simpson, cp90_simpson, abinit_corrected_trapezoid integrators
• Rename dotprod integrator to rectangle
• Truncate UPF quantities at the same index as their corresponding PSP8 quantities in UPF-PSP8 agreement tests; significantly tighten tolerances
• Polish and reorganize load/list/show functions; add unit tests
• Add identifier field to PsPFiles (defaults to filename when loaded from a file)
• Make checksum and identifier optional fields in AbstractPsPs; this information may not be available / correct if the psp is constructed by hand or modified, e.g. in the DFTK HGH sensitivity test / example
• Remove Base.== and Base.hash dispatches from AbstractPsP now that checksum is optional
• Reorganize PsPFile constructors to have two common dispatches: io::IO and path::AbstractString (which calls down to io::IO and passes along the filename)
• More minor changes

TODOs

• Add LocalPotentialCorrection flags
• Implement QE-style erf(r) Z / r local potential correction
• Add r as an argument to all the integrators (cp90_simpson needs it to avoid r=0, future integrators might use it as well)
• Derive and explain the integrators (we know the true derivative of the radial mesh dr and use it rather than finite differences r[i+1] - r[i] or similar)
• Write long-form documentation for the new reduced interface surface which explains the different dispatches
• Ignore the DC component of Fourier quantities in tests
• Look at SIRIUS implementation of the Fourier transform of the augmentation charge -- it's different than the (shot in the dark, probably incorrect) implementation that we have

[codecov]

Codecov Report

Attention: 75 lines in your changes are missing coverage. Please review.

Comparison is base (8dd0f06) 82.45% compared to head (4d6b119) 88.13%.

❗ Current head 4d6b119 differs from pull request most recent head 0d70da7. Consider uploading reports for the commit 0d70da7 to get more accurate results

Files	Patch %	Lines
src/psp/psp.jl	10.71%	25 Missing ⚠️
src/common/quadrature.jl	77.92%	17 Missing ⚠️
src/file/file.jl	25.00%	9 Missing ⚠️
src/io/show.jl	86.79%	7 Missing ⚠️
src/psp/ultrasoft.jl	30.00%	7 Missing ⚠️
src/io/load.jl	89.28%	3 Missing ⚠️
src/psp/hgh.jl	94.73%	3 Missing ⚠️
src/file/psp8.jl	94.73%	2 Missing ⚠️
src/common/mesh.jl	90.00%	1 Missing ⚠️
src/psp/numeric.jl	98.76%	1 Missing ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##             main      #11      +/-   ##
==========================================
+ Coverage   82.45%   88.13%   +5.68%     
==========================================
  Files          22       25       +3     
  Lines        1533     1610      +77     
==========================================
+ Hits         1264     1419     +155     
+ Misses        269      191      -78     

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

[azadoks]

Some comments on integration:

• The trapezoid and rectangular integrators here are currently incorrect
• Non-uniform grids require that Δr be calculated from psp.r; using the grid derivative psp.dr leads to gaps / overlaps in integration
• Agreement between rectangle/trapezoid/simpson quadrature w.r.t Gauss quadrature of cubic spline-interpolated properties improves if the property on a non-uniform mesh is resampled on a linear mesh (sometimes by orders of magnitude!)

[mfherbst]

Quick comment on the dispatches:

• I usually put the container (i.e. AbstractPsP) first and then the representation of the quantity or the input data.
• When it comes to algorithmic selection, it's usually better to do this via kwargs rather then positional args (imagine a hierarchy of calls. Then you can just pass the kwargs through from top to bottom without needing to know their order or what they are. This means once this is pulled through all the way to DFTK, you could add new algorithmic variations in e.g. PSPIO and directly use them in DFTK without the need for additional modifications.

Otherwise, there is always the risk of overdoing it. I trust you on this, but I just wanted to point out alternatives to introducing too many types. Sometimes properties can also be elegantly distinguished by plain symbols. This has disadvantages (e.g. multiple dispatch does not work), but has advantages, too (e.g. seamless integration with getproperty or similar'

I'll do a more careful review during the review (if I get raound).

[mfherbst]

Quick comments.

[mfherbst]

This feels a bit like a too overwhelming method zoo. People will generally not remember all this. Let's brainstorm if we can simplify by show functions on clever types or so.

[mfherbst]

Nice docs, super valuable.

[mfherbst]

Simpson is a function. Breaks a bit with the usual pattern where algorithm selection is done by types. We should at least capitalise it (to be able to switch to a struct). Lets discuss.

[azadoks]

Yeah, I was also thinking about this and more generally about the integration interface yesterday when I was doing some testing on the integrators. You're right that these should be type flags as well. Could take some inspiration from NumericalIntegration.jl, though I'm not sure I'd want to have it as a dependency.

[mfherbst]

Regarding the kwarg thing: Inside pspio we probable want to deconstruct (via a dispatch) to positional args, but let's discuss if this mskes sense.