This is a requirement for integration with legacy codes

The angular grids are precomputed and stored statically.

I wonder whether it would make sense to also precompute the primitive quadratures, since this could be done in arbitrary precision with e.g. sympy.

Storing the weights and nodes for $n$ points costs $2 n$ floating point numbers. This means that storing all rules from $n=1$ to $N$ costs $(N+1)^2 - N -1$.

It might be a good idea to to have this functionality anyway, since one should test the accuracy of the computed nodes and weights by comparison to accurate reference values.

In Int. J. Quantum Chem. 34, 55 (1988), Savin employs a partitioning where the atomic weight is given by

$$ w_A({\bf r}) = z_A \zeta_A^3 \exp(-\zeta_A r) $$

where $\zeta_A = Z_A$ "is a reasonable choice".

This partitioning is already available in CRYSTAL.

The other quadratures in quadrature/ have consistent file names, but the Gauss-Chebyshev rules have shorter names. They should be renamed, e.g. gausscheby1.hpp -> gausschebyshev1.hpp.

Add the following weight partition schemes:

1. Becke (w/ and w/o atomic size adjustments https://doi.org/10.1063/1.454033
2. Strattman-Scuseria-Frisch https://doi.org/10.1016/0009-2614(96)00600-8
3. Laqua-Kussmann-Ochsenfeld https://doi.org/10.1063/1.5049435
4. Delley https://doi.org/10.1063/1.458452 #60
5. Franchini-Philipsen-Visscher https://doi.org/10.1002/jcc.23323 #69

All quadratures seem to be on the same footing, separate into (1) priimtive quadratures (GC, GL, etc), (2) radial quadratures and (3) angular quadratures.

The tests are currently using Catch2. However, many Linux distributions have already switched to Catch3.

The migration guide is available at
https://github.com/catchorg/Catch2/blob/devel/docs/migrate-v2-to-v3.md

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.

According to the Catch2 documentation

We strongly recommend against using Approx when writing new code. You should be using floating point matchers instead.

The tests should accordingly use

• WithinAbs(double target, double margin),
• WithinRel(FloatingPoint target, FloatingPoint eps) and
• WithinULP(FloatingPoint target, uint64_t maxUlpDiff)

Baker, Andzelm, Scheiner, and Delley describe a radial grid in eq (12) of J. Chem. Phys. 101, 8894 (1994), which is not currently available.

Womersley has published t designs ranging up to 16 k points.

https://web.maths.unsw.edu.au/~rsw/Sphere/EffSphDes/sf.html

(Sorry, hit enter too soon.) Anyway, to be brief, I think you're mentioning IntegratorXX::integratorxx but exporting IntegratorXX::IntegratorXX.

https://github.com/wavefunction91/IntegratorXX/blob/master/cmake/IntegratorXXConfig.cmake.in#L48 is the lowercase. Otherwise, easy enough to add at psi4/psi4#3041

Shizgal, Ho, and Yang propose a Gauss-Maxwell quadrature in J. Math. Chem. 55, 413 (2017).

It is not immediately obvious what this scheme actually is.

addendum: it looks like a much better reference for the actual nodes and weights is J. Comput. Phys. 41, 309 (1981)

Lindh, Malmqvist and Gagliardi describe a radial quadrature rule for Gaussian basis sets in Theor. Chem. Acc. 106, 178 (2001). This would be nice to have as well.

The source files don't mention any license at all.

The license should be included at the start of every source file.

Per discussion in #27 (comment), the "canonical" implementation of the MK quadrature uses a Uniform Trapezoid base quadrature. However Lehtola, Marques J. Chem. Phys. 157, 174114 (2022) points out that accuracy of this scheme can be improved via using Gauss-Chebyshev (second kind) as the base quadrature.

Slater's rules offer a simple way to compute somewhat realistic atomic electron densities. They are also used in the formation of some quadrature grids. The rules used in Molpro have been summarized in J. Chem. Phys. 149, 044103 (2018); in addition, one needs to assign formal occupations to the orbitals.

#35 pointed out that the bounds transformations from the "natural" bounds of a quadrature should be decoupled from the quadrature implementation to reduce code replication. This will require two steps:

• Remove existing bounds internal transformations in the GC quadratures
  • fixed in #27 SHA 44e6c97
• Add accessors to the "natural" bounds of a particular quadrature to dispatch subsequent bounds transformers.

For some weight functions, (e.g. $\ln^2$ of Gill and Chien), there are not simple methods for the computation of quadrature nodes and weights. The Golub-Welsch algorithm provides a systematic way of numerically calculating these quantities for arbitrary weight functions (given the ability to compute moments of specified order).

It requires the solution of a Hankel-type tridiagonal eigenvalue problem, which tends to be numerically unstable for large orders, but having it as an option would lower the barrier to entry for testing out new quadratures and an important validation tool moving forward.

Requires:

• Integration of LAPACK++ to solve the EVP
• Tabulate coefficients for the classical polynomials

Franchini et al proposed an atomic partitioning in Journal of Computational Chemistry 2013, 34, 1819–1827

$$ p_i(\mathbf{r}) = \frac{\mathcal{P}_i(\mathbf{r})}{\sum_j\mathcal{P}_j(\mathbf{r})} $$

where the atomic partition function has a very simple form:

$$ \mathcal{P}_i(\mathbf{r}) = \eta_i \frac{\mathrm{exp}(-2r_i)}{r_i^3} $$

The same partitioning was also used in their paper on Coulomb fitting J. Chem. Theory Comput. 2014, 10, 1994.

Hirshfeld partitioning could be easily implemented in IntegratorXX. Here, the weight function is defined by the ratio of electron densities

$$ w_A({\bf r}) = \frac {n_A({\bf r})} {\sum_B n_B({\bf r})} $$

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.

The weight used in gausscheby2 are incorrect.

There are two problems here. Compared to the equations given on Wikipedia

• the sine squared term is missing, and
• the transformation to unit weight - which is present in gausscheby1 - is also missing.

Delley published high order integration schemes on the unit sphere in J. Comput. Chem. 17, 1152 (1996).

These schemes are similar to Lebedev's schemes, but are determined to higher numerical precision. I have also heard that Delley's grids do not suffer from the problem of negative weights, which affects some Lebedev grids.

I have asked Delley by email whether he would be willing to contribute his grids.

We should provide generators for the "Standard Grids" commonly employed in quantum chemistry:

• SG-0
• SG-1
• SG-2
• SG-3

Mieke Peels' PhD thesis states that they have obtained new angular integration grids

The new integration grids have been generated to integrate polynomials, and therefore spherical harmonics, up through order 𝐿 = 105. Many of the new grids contain fewer points than the Lebedev grids. None of the grids contain negative weights, and the weight spreads are generally smaller than for the Lebedev grids.

The thesis also states that

In the future, grid generators containing optimized seed points will be released open source as part of the angular integration grid generator (aigg) program, which would allow for the new grids to be implemented in existing quantum chemistry software packages.

I have asked Gerald Knizia about these grids in October 2022, and in March 2023, but have not received a reply. I did contact Peels as well, who said I should be in touch with Knizia about the code.

My assumption is that Gerald is simply too busy, so if someone runs into him at a conference, please ask about this.

#22 fixed the GaussChebyshev2 implementation, but it still needs a unit test.

Messages such as the following are less-than-meaningful:

/home/dbwy/Development/IntegratorXX/test/gausslegendre.cxx:42824: FAILED:
  REQUIRE_THAT( wgt[i], Catch::Matchers::WithinAbs(ref_wgt[i],w_tolerance) )
with expansion:
  0.0000203866 is within 0.0 of 0.1000203866

The following information should be printed in a sane manner

1. The point index where the failure occurred
2. Which rule generated the error should be clearly noted
3. The observed difference between calculated and reference values (in scientific notation)
4. Comparison epsilon should be printed in scientific notation

#23 added formatting files, the corresponding GH action should be enabled

https://github.com/wavefunction91/MACIS/blob/master/.github/workflows/clang_format_lint.yaml

Delley's paper on numerical atomic orbital calculations, J. Chem. Phys. 92, 508 (1990), states employing a partition function (eq. 3)

$$ p_A = \frac { g_A({\bf r} - {\bf R}_A) } { \sum_B g_B ({\bf r} - {\bf R}_B) } $$

where "preferred choices are the functions"

• $g_A = \rho_{A}^2(r)$,
• $g_A = \rho_{A}^2(r) / r^2$, and
• $g_A = \rho_{A}(r) [\exp(r_0/r) - 1 - r_0/r] $

where $\rho_A$ is the atom's electron density.

Delley writes

These partition functions automatically depend on atom sizes via $p_A$ and have given consistently good results for all sorts of compounds including heteronuclear compounds containing heavy transition metals and light elements including hydrogen. The last function leads to continuous vanishing of all partition functions and all their derivatives on all nuclei where they are not centered.

These schemes are clearly similar to Becke's, and should be quite easy to implement. The only thing necessary are just the atomic electron densities. The adaptive Molpro grid wavefunction91/GauXC#51 implementation of J. Chem. Phys. 157, 234106 (2022) employed Slater atomic densities:

In Ref. 79, these are defined for elements up to $n = 6$ with effective principal numbers of $n^\star = 3.7$, 4.0, 4.2 for $n = 4, 5, 6$. The implementation in Molpro uses $n^\star = 4.4$ for seventh row elements based on a simple linear extrapolation of the $n^\star$ values given by Leach.79 The advantage of using Slater’s densities over tabulated numerical densities is that they are given analytically and can be easily computed at any given grid point.

NWChem implements an error function based partitioning described in J. Chem. Phys. 152, 184102 (2020) on page 8

$$ w_A({\bf r}) = \prod_{B \neq A} \frac 1 2 [ 1 - \text{erf } \mu_{AB}'] $$

where

$$ \mu_{AB}' = \frac 1 \alpha \frac {\mu_{AB}} { [1 - \mu_{AB}^2]^{n}} $$

and

$$ \mu_{AB} = \frac {\mathbf{r}_A - \mathbf{r}_B} {|\mathbf{r}_A - \mathbf{r}_B|} $$

In particular, each of the primitive quadratures are guaranteed to produce exact integrals for polynomials of a particular order subject to a particular weight function. This can be tested directly by randomly generating polnomials for the unit tests.

For the radial quadratures, each ansatz was designed to integrate function of a specific form (e.g. radial integrands). These tests should be updated to test these classes of functions.

The Gauss-Legendre routine only computes half of the nodes and weights, since they are symmetric.

However, I think the Gauss-Chebyshev rules 1 and 2 (and its modification) also generate symmetric nodes and weights, allowing this symmetry to be taken advantage of.

At

      // TODO: Tranform on bound                                                                                                                                                              
      points[i] = pt;
      weights[i] = wgt;

What do you mean by this?

#46 includes tests for low orders of nodes.

I think it would be good to also include tests for higher numbers of nodes, since quadrature grids often use even several hundreds of points.

However, not every number of points needs to be tested. I think that going from 100, it might be reasonable to test rules every 20 points, (100, 120, 140, 160, 180, 200), then maybe every 25 points (225, 250, 275, 300, 325, 350, 375, 400), then every 50 points (450, 500, 550, 600).

Reducing the spacing of the tests is a big reduction storage space, and the cost to generate the tests; even the 100-point rules take some seconds to generate on my laptop with SymPy...

The Gauss-Chebyshev rules 1 and 2 are found in SciPy; one just has to modify the weights, which can be done in the generator.

Weber, Daul and Baltensperger discuss adaptive radial grids using sinc quadrature in Computer Physics Communications 163, 133 (2004). They also present equations for STOs for using Lindh's scheme from Theor Chem Acc 106, 178 (2001).

I don't know if this will be useful, but I am opening the ticket as a note.

Ahrens and Beylkin have described rotationally invariant quadratures on the sphere in Proc. R. Soc. A (2009) 465, 3103–3125.

The code should check that the point and weight arrays have the same size. This check should be configurable and possible to turn off for production level.

An improvement over #87 is to use minimal basis non-integer Slater-type orbital wave functions. Such wave functions have been reported by Koga et al in Int. J. Quantum Chem. 62, 1 (1997).

Unfortunately, this paper does not specify what the employed orbital occupations are. However, the occupations are likely the same as in other works reporting (differences to) numerical Hartree-Fock energies; occupations have been given in e.g. Atomic Data and Nuclear Data Tables 95 (2009) 836

Self explanatory

The Gauss-Legendre quadrature routine hardcodes an absolute tolerance eps = 3.0e-11. The origin of this value is not documented, and it also seems a bit too large to me. Since $x\in[-1,1]$, eps should be of the order of machine epsilon of the employed data type.

Another thing that needs to be fixed is to employ a for loop instead of the current while loop which risks hanging. Instead, one should take e.g. a maximum of 100 Newton iterations, and then throw an error if convergence was not achieved.

I am not sure if these are already implemented, but...

1. Murray et al describe a scheme which is exact for all spherical harmonics up to angular momentum $L$: Gauss-Legendre quadrature is used in $\theta$ with $N_\theta =(L + 1)/2$ points, and "the simple Gauss scheme" for $\phi$ quadrature with $N_\phi = L + 1$ points

2. Treutler and Ahlrichs describe a more compact scheme with Lobatto polynomials with $N_\theta = (L+3)/2$ and $N_\phi = L$ except for $\cos \theta = \pm 1$ for which $N_\phi=1$ and for $\cos \theta = 0$ for which $N_\phi = L+1$. They also discuss reducing the number of points in phi near the pole by determining the number of points from $(\sin \theta_i)^{N_\phi} < \epsilon$.

Add Becke's original radial transform

$$ r = R \frac {1-x} {1+x} $$

as a quadrature method

Gauss-Lobatto quadrature should be added as well

The Gauss-Jacobi quadrature integrates functions of the form

$$ \int_{-1}^1 f(x) (1-x)^\alpha (1+x)^\beta\mathrm{d}x \approx \sum_{i=1}^n w_i f(x_i) $$

where ${x_i}$ are the nodes of $n$-th order Jacobi polynomial, $P^{(\alpha,\beta)}_n$. The Chebyshev and Legendre polynomials are special cases of the polynomials which already have quadrature rules generated.

Hayle and Trefethen discuss quadrature formulas from conformal maps in SIAM J. Numer. Anal. 46, 930 (2008), which may prove useful.