Consider using SphericalBesselTransforms.jl for log-grid pseudos

[azadoks]

I just submitted SphericalBesselTransforms.jl the registrator.

It's a fairly direct translation of pySBT which implements the algorithm described in Talman, Comp. Phys. Comms. 180, 332 (2009) for the SBT/Hankel transform of functions on a logarithmic radial mesh in O(Nlog(N)) time using FFTs.

The implementation is fairly minimal, and the interface is not super ergonomic at the moment. I spent some time with the profiler to make sure it's pretty fast, so the typing might be a bit too strict.

An implementation in DFTK would probably also require some interpolation because the transform returns g(k) on a logarithmic grid. For AD support, we'd also need to write a rule.

I wrote up a quick example here which does the transforms for HGH pseudos and compares with the analytic reference. It seems to do quite well, but I should probably also compare to simpson integration / QuadGK.jl.

[mfherbst]

Related to #986

[mfherbst]

Not sure how much this is still a performance bottle neck right now (@abussy may know), but could be worth considering, indeed.

[azadoks]

I just ran a little bit of timing on the Silicon / PseudoDojo example, and indeed the time in hankel is negligible. Things might get more expensive when ultrasoft/PAW are implemented because there are so many augmentation charges, etc., but I can't imagine that it would get significantly worse.

One other benefit to the FFT-SBT is that it can also be more stable / accurate than direct integration (see the second plot in pySBT's README).
I'll eventually run some tests and see if this ever affects ONCVPSPs.

[abussy]

Performance wise, hankel is negligible in most cases. There is however one specific case where it is not: when Dual numbers are involved. In stress and response calculations, this becomes a major bottleneck.

For example, in the Xc term instantiation, the core charge density is calculated:

DFTK.jl/src/pseudo/PspUpf.jl

Lines 256 to 260 in b6880b8

	function eval_psp_density_core_fourier(psp::PspUpf, p::T) where {T<:Real}
	rgrid = @view psp.rgrid[1:psp.ircut]
	r2_ρcore = @view psp.r2_ρcore[1:psp.ircut]
	return hankel(rgrid, r2_ρcore, 0, p)
	end

Each individual integration is fast, but there are 1000s taking place. I will try and parallelize it for the GPU in the spirit of #1163 , but CPU stress/response calculations of large systems will remain very expensive. There might be a benefit there.

[antoine-levitt]

The nice thing about having a better algorithm is hopefully we don't need to optimize it as much :)

[Technici4n]

I remember the AtomicLocal instantiation being quite slow, not sure if it was with duals or not. Generally I feel like basis construction ought to be faster.