Implement double-exponential radial quadrature

[susilehtola]

The double-exponential radial quadrature of Takahasi and Mori, also known as tanh-sinc quadrature, was suggested by Mitani in Theor. Chem. Acc. 130, 645 (2011) and Theor. Chem. Acc. 131, 1169 (2012) for density functional quadrature.

This quadrature is used in the SG-2 and SG-3 standard grids, for instance.

[wavefunction91]

Great, if you have the closed-form for these quadratures, it should be an easy add.

This also brings up a point I've been meaning to address - we should have stock generators for SG-X spherical quadratures. I'll open an issue to track.

[susilehtola]

This also brings up a point I've been meaning to address - we should have stock generators for SG-X spherical quadratures. I'll open an issue to track.

Yeah, although I think those grids also go hand-in-hand with specific radial quadratures and atomic partitionings. The latter can probably be changed, the former shouldn't. But the whole idea about a standard quadrature library is that these quadratures can be made available to any code.

[wavefunction91]

@susilehtola I wanted to hash some of this out before spending time on implementing it. Let's take DE1 as an example

$$ r(x) = \exp( \alpha \sinh(x) ) $$

$$ r'(x) = \alpha \cosh(x)\ r(x) $$

The integral is then approximated (I'm neglecting the $r^2$ spherical Jacobian due to the code conventions)

$$ \int_0^\infty f(r)\mathrm{d}r = \sum_{i=-\infty}^{+\infty} w_i r'(x_i) f(r(x_i)) $$

where $x_i$ is taken to be the uniform trapezoid with grid spacing $h$ (i.e. $w_i = h$). This expression is an infinite sum that must be truncated in practice.

In the paper, they do this by selecting a $R_\mathrm{min}$ and $R_\mathrm{max}$ and truncating the sum to include $x_i\in [r^{-1}(R_\mathrm{min}),\ r^{-1}(R_\mathrm{max})]$. For each of the DE grids, the radial transformation is analytically invertible for $r>0$, e.g.

$$ r^{-1}(R) = \mathrm{asinh}\left(\frac{\ln R}{\alpha}\right) $$

where asinh has an STL implementation.

My understanding of how this should work is the following:

1. Select a sensible $\alpha$, the paper has some defaults, we can make it configurable defaulting to these values
2. Determine the $x$-range via selecting a sensible $R_\mathrm{min}$ and $R_\mathrm{max}$ s.t. $x_\mathrm{min/max} = r^{-1}(R_\mathrm{min/max})$
3. Generate a UniformTrapezoid base quadrature on $x_i \in [x_\mathrm{min}, x_\mathrm{max}]$ with a specified number of points
4. Apply the DE-X transformation to $x_i$ ala the modular design introduced in Modularize Radial Quadrature Generation #27

Make sense?

[wavefunction91]

For some undocumented $\alpha$ (but likely $\alpha = \frac{\pi}{2}$ per the paper: edit: confirmed numerically), a select set of nodes and weights are provided in boost::math