Electron repulsion integrals by Rys quadrature

This is the third ERI algorithm in lib/hf, next to McMurchie-Davidson and Head-Gordon-Pople. It shares no code with either: it never builds Hermite expansions or Obara-Saika intermediates, and it does not call the Boys function. Its tests use Boys as the oracle for the quadrature rule.

The idea. A primitive ERI over Cartesian Gaussians reduces to

(ab|cd) = K * int_0^1 Ix(t) Iy(t) Iz(t) exp(-x t^2) dt

where x = rho |P - Q|^2 (P and Q are the product centers of the bra and the ket, rho = pq/(p+q)), K collects the Gaussian prefactors, and each I is a polynomial in u = t^2 of degree at most L/2, with L = la + lb + lc + ld. The Boys function is the special case where all three I are 1. A polynomial of degree L/2 in u is integrated exactly by an n-point Gauss rule for the weight exp(-x u) / (2 sqrt(u)) on [0, 1] once 2n - 1 >= L/2, so

n = floor(L/2) + 1

roots suffice: 1 for (ss|ss), 3 for (pp|pp), 5 for (dd|dd), 7 for (ff|ff). The rule for weight x is the Rys quadrature; its moments are the Boys values,

sum_i w_i u_i^k = int_0^1 t^(2k) exp(-x t^2) dt = F_k(x),   k < 2n,

which is what the first test checks.

Roots and weights. The rule changes with x, and every primitive quartet has its own x, so the roots are computed on the spot:

1. Discretize the measure. An 80-point Gauss-Legendre rule on [-1, 1], folded to its 40 positive nodes, integrates t^(2k) exp(-x t^2) on [0, 1]. Above x = 80 the integrand is negligible past t = sqrt(80/x), so the nodes shrink onto [0, sqrt(80/x)] and the integrand the rule sees is never sharper than exp(-80 v^2). The dropped tail is at most 7.7 × 10⁻²¹ of F_13, the highest moment 7 roots use (the regularized upper incomplete gamma Q(13.5, 80), from mpmath).
2. Stieltjes. On the discrete measure, build the monic orthogonal polynomials by the three-term recurrence pi_(k+1)(u) = (u - alpha_k) pi_k(u) - beta_k pi_(k-1)(u), reading alpha_k and beta_k off discrete inner products.
3. Golub-Welsch. The roots are the eigenvalues of the Jacobi matrix with alpha on the diagonal and sqrt(beta) beside it, and each weight is beta_0 (the total mass, F_0) times the squared first component of its unit eigenvector. The eigenproblem is tridiagonal and at most 7 by 7, so an implicit QL iteration that tracks only the first row of the eigenvector matrix does it.

Measured here, over the 16 seeded x of the first test (10⁻⁴ to 200): the worst relative error in F_k, k < 2n, against lib/hf/boys is 8.9 × 10⁻¹⁶ at n = 1 and 4.7 × 10⁻¹⁵ at n = 7. A numpy prototype of the same three steps chose the constants: with 72 nodes and a cutoff of 100 its worst error (n <= 5, x up to 10⁵) was 3.9 × 10⁻¹², and with 80 nodes and a cutoff of 80 it was 6.2 × 10⁻¹⁴ (n <= 7, x up to 10⁶, numpy's eigensolver rather than this QL). The node count and the cutoff move together.

Against PySCF, the worst absolute gap is 9.2 × 10⁻¹⁵ over the 512 prims picks (every class through (ff|ff)) and 2.2 × 10⁻¹⁴ over the 4186 water 6-31G values. Shifting C00 by one part in 10⁷ makes the prims test fail at 4.0e-9.

The 2-D integrals. At a root u, each Cartesian direction separates. The integral G(i, k) of (x - Ax)^i (x - Cx)^k against the root's 2-D Gaussian obeys

G(i+1, k) = C00 G(i, k) + i B10 G(i-1, k) + k B00 G(i, k-1)
G(i, k+1) = D00 G(i, k) + k B01 G(i, k-1) + i B00 G(i-1, k)

with G(0, 0) = 1 and

B00 = u / (2(p+q))
B10 = (1 - q u/(p+q)) / (2p),   B01 = (1 - p u/(p+q)) / (2q)
C00 = (P - A) - q u/(p+q) (P - Q)
D00 = (Q - C) + p u/(p+q) (P - Q)

G puts all the bra's angular momentum on A and the ket's on C. Writing x - B = (x - A) + (A - B) and expanding the power moves it back:

I(a, b; c, d) = sum_i sum_j C(b,i) (A-B)^(b-i) C(d,j) (C-D)^(d-j) G(a+i, c+j)

Then (ab|cd) = K sum_r w_r Ix Iy Iz, with K = 2 pi^(5/2) / (p q sqrt(p+q)) exp(-ab/p |AB|^2) exp(-cd/q |CD|^2).

Normalization. Everything is in the library's unit-Cartesian convention: each primitive carries cn from bas_contract, and each component its bas_nc factor.

Validated range. rys_applicable(L) is 1 for L <= RYS_LTOT_MAX = 12, (ff|ff), the highest class lib/hfref/prims holds. Every prims case agrees with PySCF there, and so does every water 6-31G value (tests below). The coverage counters use it to decide which classes Rys may vote on.

Used with import hf.rys_eri.