Modularize Radial Quadrature Generation

[wavefunction91]

Unify implementation of radial quadrature generation. Closes #12.

[wavefunction91]

@susilehtola This should add everything that we need to unify the implementations. Dox incoming.

I've written the traits for TA, but I can't consolidate the Chebyshev roots with the closed form expression given in https://doi.org/10.1002/jcc.10211 (which also agrees with the Psi4 implementation as far as I can tell). For order n Chebyshev quadrature, the points are sinusoidally generated with spacing pi/n, where as these implementations generate on pi/(n+1). Unclear how to resolve this. Is there a reason that the TA radial qudature of order n uses the Chebyshev nodes of order n+1? We just need to be able form the logic to relate the dimension of the base quadrature with that of the transformed quadrature.

[susilehtola]

I've written the traits for TA, but I can't consolidate the Chebyshev roots with the closed form expression given in https://doi.org/10.1002/jcc.10211 (which also agrees with the Psi4 implementation as far as I can tell). For order n Chebyshev quadrature, the points are sinusoidally generated with spacing pi/n, where as these implementations generate on pi/(n+1). Unclear how to resolve this. Is there a reason that the TA radial qudature of order n uses the Chebyshev nodes of order n+1? We just need to be able form the logic to relate the dimension of the base quadrature with that of the transformed quadrature.

I can't say I fully understand what you are saying, but yes, it looks like there may be an issue with the Psi4 code. I've opened up a ticket psi4/psi4#3004

[susilehtola]

Treutler and Ahlrichs clearly state that

These mappings have been combined with Chebyshev quadratures of the first (T1) and second (T2) kind, where $-1 \le x \le 1$, and a third formula based on the weight function $\sqrt{x/(1-x)}, 0 \le x \le 1$, here denoted T3, which is a variant of T1 (for formulas of T1, T2, and T3, see p. 889 of Ref. 1).

Their Ref. 1 is Abramowicz and Stegun, and page 889 gives the following rules:

1. Chebyshev polynomials of the first kind

$$ \int_{-1}^{1} \frac {f(x)} {\sqrt{1-x^2}} {\rm d}x = \sum_{i=1}^{n} w_i f(x_i) + R_{n} $$

with abscissas and weights

$$x_i = \cos \frac {(2i-1) \pi} {2n} $$

$$w_i = \frac {\pi} {n}$$

2. Chebyshev polynomials of the second kind

$$ \int_{-1}^{1} f(x) \sqrt{1-x^2} {\rm d}x = \sum_{i=1}^{n} w_i f(x_i) + R_{n} $$

with abscissas and weights

$$x_i = \cos \frac {i \pi} {n+1} $$

$$w_i = \frac {\pi} {n+1} \sin^2 \frac {i \pi} {n+1}$$

3. The third rule is

$$ \int_{0}^{1} f(x) \sqrt{\frac {x} {1-x}} {\rm d}x = \sum_{i=1}^{n} w_i f(x_i) + R_{n} $$

with abscissas and weights

$$x_i = \cos^2 \frac {(2i-1) \pi} {2(2n+1)} $$

$$w_i = \frac {2 \pi} {2n+1} x_i$$

The rules T1 and T2 are in perfect agreement with the current equations in https://github.com/wavefunction91/IntegratorXX/blob/master/include/integratorxx/quadratures/gausscheby1.hpp and https://github.com/wavefunction91/IntegratorXX/blob/master/include/integratorxx/quadratures/gausscheby2.hpp.

As we discuss in J. Chem. Phys. 157, 174114 (2022), Gill and Chien's equations are wrong for Mura-Knowles quadrature (they write it with trapezoidal quadrature instead of Chebychev quadrature), so I would not rely on their equations being correct for Treutler-Ahlrichs, either...

[wavefunction91]

Thanks for looking in to it - glad to know it's not just something I'm missing in my understanding.

Gill and Chien's equations are wrong for Mura-Knowles quadrature (they write it with trapezoidal quadrature instead of Chebychev quadrature)

That's what we do too... Interesting, @hjjvandam what does NWChem do? I was able to get (within epsilon) identical grids to NWChem at one point. If there's just an honest discrepancy in the lit / community software, we can expose both via the new modular API.

[loriab]

• should https://github.com/wavefunction91/IntegratorXX/blob/master/include/integratorxx/quadratures/gausscheby2.hpp#L26 be scaling by 1/sqrt( 1-x_i^2) ?
• https://github.com/wavefunction91/IntegratorXX/blob/master/include/integratorxx/quadratures/gausscheby1.hpp#L71 denominator quite right: formula is pi/2n and code is 1/(2n pi)

I looked at Susi's equations from Abramowicz a little closer, and the below is what I would do for psi4. The Type1 call is new and the Type2 call isn't actually changing code. I think the scaling for unit form actually is present, but it gets folded into the weight. Of course that doesn't explain why IntegratorXX and Psi4 are producing different values. Could it be something about Integrator looping from i=0 and filling arrays weights[i] with i+1 and Psi looping from i=1 and filling arrays weights[i-1] with i?

Basically, I don't have an answer, so putting forward status from another set of eyes.

void RadialGridMgr::getChebychevRootsType1(int n, double r[], double w[]) {
    double piOver2N = M_PI / (2 * n);
    for (int i = 1; i <= n; i++) {
        double x = cos((2 * i - 1.0) * piOver2N);
        r[i - 1] = x;
        w[i - 1] =
            2 * piOver2N * sqrt(1.0 - x * x);  
            // raw weight = pi / n = 2 * piOver2N
            // but scaling for unit form by *= sqrt(1 - x^2), final weight = 2 * piOver2N * sqrt(1 - x^2) = 2 * piOver2N * sin((2 * i - 1.0) * piOver2N)
    }
}

void RadialGridMgr::getChebychevRootsType2(int n, double r[], double w[]) {
    double piOverNPlusOne = M_PI / (n + 1);
    for (int i = 1; i <= n; i++) {
        double x = cos(i * piOverNPlusOne);
        r[i - 1] = x;
        w[i - 1] =
            piOverNPlusOne * sqrt(1.0 - x * x);  
            // raw weight = piOverNPlusOne * sin^2 (i * piOverNPlusOne) = piOverNPlusOne * (1 - cos^2 (i * piOverNPlusOne)) = piOverNPlusOne * (1 - x^2)
            // but scaling for unit form by /= sqrt(1 - x^2), final weight = piOverNplusOne * sqrt(1 - x^2) = piOverNPlusOne * sin(i * piOverNPlusOne)
    }                        

[wavefunction91]

• should https://github.com/wavefunction91/IntegratorXX/blob/master/include/integratorxx/quadratures/gausscheby2.hpp#L26 be scaling by 1/sqrt( 1-x_i^2) ?
• https://github.com/wavefunction91/IntegratorXX/blob/master/include/integratorxx/quadratures/gausscheby1.hpp#L71 denominator quite right: formula is pi/2n and code is 1/(2n pi)

There were bugs in the GC quadratures, fixed in #32 (+ added unit tests). These were place holders before we started the refactor, now they've been validated.

so putting forward status from another set of eyes.

Those seem to agree with what's currently getting PR-ed (modulo trig optimizations). I think this makes things consistent between the two codes for the base quadratures, would be easy enough to check.

The remaining question is how to validated that the generated Treutler-Ahlrichs (TA) quadratures are consistent. As @susilehtola pointed out, the Gill paper might have the wrong final expressions, so I may have been chasing my tail a bit trying to get things to agree. Do we have a third party source to validate the final TA quadratures?

[loriab]

Ok, so Psi4 is believed correct at the moment? Whew, but let me know if that changes. I don't have any particular sources on the TA quadratures.

[wavefunction91]

If the quadrature generated on this line is meant to be Chebyshev of the second kind, then yes, this is correct.

[loriab]

If the quadrature generated on this line is meant to be Chebyshev of the second kind, then yes, this is correct.

Yes, that line is second kind -- I'll fix up the naming and documentation.

[susilehtola]

The remaining question is how to validated that the generated Treutler-Ahlrichs (TA) quadratures are consistent. As @susilehtola pointed out, the Gill paper might have the wrong final expressions, so I may have been chasing my tail a bit trying to get things to agree. Do we have a third party source to validate the final TA quadratures?

I can get in touch with the Turbomole developers, since the original implementation of Treutler and Ahlrichs is in that code... but it does now look like the expressions are correct.

[wavefunction91]

@susilehtola I was able to confirm that generating the quadratures via GC-2 -> TA-M4 gave the same thing as the hand coded variant that I got from the Gill paper. Might be good to get some external validation, but it's likely not a priority.

[wavefunction91]

I've released that we'll need to expose some runtime configurability in this modular design. The grids that are exposed right now are relatively simple (the only real configurable parameter is a linear scaling factor), but the DE grids require a bit more configurability (deciding where to cut off the infinite sum, etc). The specification of alpha for TA could fall under this type of configurability.

This kind of configurability will also be useful for folks trying to optimize grids programmatically.

[hjjvandam]

Hi David, NWChem does the following: fmn = rm/(1d0+nr) do 10 i = 1,nr qi = dble(i)/(nr+1d0) ri = -alpha*log(1.D0 - qi**rm) wi = fmn*alpha*(ri*ri)/(1.D0 - qi**rm)*qi**(rm-1.d0) r(i) = ri w(i) = wi 10 continue where rm=3.0 and alpha=3.3*range where range is the maximum radius covered by the radial grid, and nr is the number of radial grid points. This code lives in nwchem/src/nwdft/xc/xc_murakn.F. Best wishes, Huub

…

----------------------------------------------------------------------------------------------------- Hubertus van Dam, 631-344-6020, ***@***.******@***.***> Brookhaven National Laboratory From: David Williams-Young ***@***.***> Date: Friday, July 7, 2023 at 1:59 PM To: wavefunction91/IntegratorXX ***@***.***> Cc: Van Dam, Hubertus ***@***.***>, Mention ***@***.***> Subject: Re: [wavefunction91/IntegratorXX] Modularize Radial Quadrature Generation (PR #27) Thanks for looking in to it - glad to know it's not just something I'm missing in my understanding. Gill and Chien's equations are wrong for Mura-Knowles quadrature (they write it with trapezoidal quadrature instead of Chebychev quadrature) That's what we do too... Interesting, @hjjvandam<https://github.com/hjjvandam> what does NWChem do? I was able to get (within epsilon) identical grids to NWChem at one point. If there's just an honest discrepancy in the lit / community software, we can expose both via the new modular API. — Reply to this email directly, view it on GitHub<#27 (comment)>, or unsubscribe<https://github.com/notifications/unsubscribe-auth/ABDS7HX2A6UNPY2J23ERTYTXPBFAJANCNFSM6AAAAAA2BFFRJU>. You are receiving this because you were mentioned.Message ID: ***@***.***>

[wavefunction91]

Thanks @hjjvandam, that's definitely the uniform trapezoid variant.

[wavefunction91]

@susilehtola Anything in particular that would be beneficial to roll up into this PR? I think this is more or less g2g.

[susilehtola]

@susilehtola Anything in particular that would be beneficial to roll up into this PR? I think this is more or less g2g.

I think this is g2g, especially since I want the changes to the Chebyshev routines in master... The modularization can continue in other pulls if necessary