Electron repulsion integrals by the Head-Gordon-Pople recurrences

The two-electron integral over four contracted Cartesian Gaussians,

(ab|cd) = double integral a(1) b(1) (1/r12) c(2) d(2),

is the costly part of every Hartree-Fock build: water in 6-31G* already has 18,145 of them that differ by symmetry. This module computes them by the scheme of Head-Gordon and Pople (J. Chem. Phys. 89, 5777, 1988), which splits the work in two:

1. A vertical recurrence (Obara and Saika's, with the auxiliary index m) builds [e0|f0]^(m) for one primitive quartet: all the angular momentum of the bra sits on A (e runs up to la + lb) and all of the ket on C. These are summed over the primitives.
2. A horizontal recurrence then moves power from A to B and from C to D on the contracted sums:

(a, b+1_i | = (a+1_i, b | + AB_i (a, b |

It needs no exponents, so it runs once per shell quartet instead of once per primitive quartet. That is the point of the scheme: the expensive step works on a short list of targets, and the cheap step does the fan-out.

The vertical step. With zeta = a + b, eta = c + d, the product centers P and Q, W = (zeta P + eta Q)/(zeta + eta) and rho = zeta eta / (zeta + eta):

[00|00]^(m) = 2 pi^(5/2) / (zeta eta sqrt(zeta + eta))
              exp(-ab/zeta AB^2 - cd/eta CD^2) F_m(rho PQ^2)

[e+1_i,0|f0]^(m) = PA_i [e0|f0]^(m) + WP_i [e0|f0]^(m+1)
    + e_i/(2 zeta) ([e-1_i,0|f0]^(m) - rho/zeta [e-1_i,0|f0]^(m+1))
    + f_i/(2(zeta+eta)) [e0|f-1_i,0]^(m+1)

and the same with the roles of bra and ket swapped for raising f: QC_i, WQ_i, f_i/(2 eta), rho/eta. The code first raises e with f = 0, where the last term vanishes, and then raises f for every e it still needs. F_m comes from boys_values in lib/hf/boys.viv; nothing else in the module touches a special function.

Components. Every Cartesian component with total power up to L has one index: l(l+1)(l+2)/6 for the components below l, then PySCF's order in the shell (lx counting down, then ly). hgp_ctab(L) tabulates the powers and the neighbors one step up and down along each axis. lib/hf/os_onebody has the same table under its own prefix; this module keeps its own copy so that the ERI family can be loaded, tested and replaced without the one-body one.

Normalization and storage. Contracted coefficients come from bas_contract, and each component gets bas_nc(lx, ly, lz) at the end, so the values are PySCF's int2e_cart divided by F for each index. hgp_eri fills the packed-8 vector of lib/hf/eri_store.viv, one shell quartet per symmetry class, with no screening.

Used with import hf.hgp_eri.