Modularize radial rules

[susilehtola]

At present, the implementations of the radial rules hardcode the primitive quadrature rule, and the radial size adjustment. This leads to code duplication, and limits the reusability of the code.

The radial rules should be constructed from an array of nodes and weights. The size adjustment, which consists of scaling the radial nodes by R and the weights by R^3 should likewise happen in a dedicated routine.

[wavefunction91]

I think we're working with a different set of vernacular here - what exactly do you mean by "primitive quadrature"?

For example, for the MHL quadrature, I would view the trapezoidal rule to be the "primitive quadrature" (i.e. the "base" quadrature), but it's not immediately obvious if that's what you mean.

Do you just mean Jacobian inclusion?

Can you perhaps provide a concrete example that exemplifies this issue?

[susilehtola]

Yes, that is exactly what I mean. The radial quadratures are composed of two things:

1. the primitive quadrature rule, e.g. trapezoidal quadrature, Gauss-Legendre quadrature, Gauss-Chebyshev quadrature, and
2. the radial transform, which is the key difference between the various rules.

Our recent study J. Chem. Phys. 157, 174114 (2022) examined some combinations of various radial transforms and primitive quadratures. The Ahlrichs M4 grid with Chebyshev quadrature afforded the best convergence.

[susilehtola]

Psi4 has a pretty nice implementation, see https://github.com/psi4/psi4/blob/afd50aa42403e7980b5e2f1fdf5d8be0f239275a/psi4/src/psi4/libfock/cubature.cc#L2396

[wavefunction91]

Ah OK, so essentially a desirable scheme would be to wrap up the radial quadrature generation into a unified API. That should be doable.

I'm not a fan of the runtime nature of the Psi4 implementation (function pointers, string enumeration, etc), since it opens up the possibility of ABI incompatibility and silent errors. I'd prefer to keep things compile time for the discrete "standard" use cases. There should be a template metaprogramming solution here, I'll give it some thought.

If there's a desire to enable runtime flexibility as an option we should be able wrap it in a special type (e.g. that stores std::function's to allow hooking up to free C/F*/Python functions). Would this be a desirable feature for your purposes?

[susilehtola]

No, we can definitely restrict to compile time only. A std::function interface would be awesome, though! What I have in my mind is a RadialTransform class that defines r and the Jacobian dr, since these are all what you need to compute the nodes and weights.