Notebooks · lib/hf/md_eri

Electron repulsion integrals by McMurchie-Davidson, one shell quartet at a time

The two-electron integral over four Cartesian Gaussians is, after the Hermite expansion of both products (see lib/hf/hermite),

(ab|cd) = 2 pi^(5/2) / (p q sqrt(p + q))
          * sum_h  E_ab(h)  sum_h'  (-1)^(t'+u'+v') E_cd(h')  R(h + h')

where h = (t, u, v) runs over the Hermite functions of the bra pair, h' over those of the ket pair, and R is taken at P - Q with exponent alpha = pq / (p + q). The sign comes from moving the derivatives of the ket's Hermite Gaussians over to the bra's: d/dQ = -d/dP on a function of P - Q.

Written as two matrix products. For one shell pair, stack the E values of every primitive pair side by side: E_ab is [components, primitive pairs x Hermite functions], with the contraction coefficients, the component factors and the Gaussian product prefactor folded in. For one shell quartet, the signed R values for a batch of bra primitive pairs against every ket primitive pair form a matrix BR, with a row per ket (primitive pair, Hermite function) and a column per bra (primitive pair, Hermite function). Then

Z_b = E_cd @ BR_b        [ket components, bra primitive-Hermite]
(ab|cd) = E_ab @ Z^T     [bra components, ket components]

where Z^T stacks the Z_b^T of every batch, and both sums over primitives happen inside the products. The batch is as many bra primitive pairs as keep BR near MDERI_BATCH = 4096 entries: all 36 at once for a contracted (ss|ss) quartet of oxygen 1s shells, one at a time for a (dd|dd) quartet with 10 primitives per shell, where a single BR over both primitive lists would be 3500 x 3500. Each primitive quartet costs one Boys evaluation and one R table, shared by every component, and the component work goes to @, which skips the zeros of E (an entry is zero whenever t passes i + j on its axis) and runs the same sum in the same order in both engines.

The shell-pair list is built once per basis: for every pair A >= B, the primitive-pair exponents p, centers P and the stacked E_ab. All of it lives in one flat array (mderi_store), because Viv has no records.

Schwarz screening. (ab|cd)^2 <= (ab|ab)(cd|cd), so with Q_AB = sqrt(max (ab|ab)) over the components of a shell pair, a quartet with Q_AB Q_CD below the threshold contributes less than the threshold to every integral it holds. mderi_screened skips those quartets (threshold MDERI_SCREEN_TOL = 10⁻¹⁴), and mderi_full computes every one; the two agree to the threshold.

One flat quartet list. The build loops over a list of (bra pair, ket pair) rows with parallel for. Each row writes a disjoint set of packed-8 slots, so the rows are independent, and the default (non-OpenMP) build runs them in order with the same result.

One integral at a time. mderi_prim1 computes a single (ab|cd) over four primitive components: E only for those powers, R only inside the box they need (herm_rfill_box), and a plain Hermite sum. It is for spot checks where a full block is wasted work, and it shares no matrix code with the block path, so the tests hold the two against each other.

Measured cost (native, cc -O2, Apple M2 Pro). A full (dd|dd) block of single primitives, 1296 integrals, takes 12.7 us; a (pp|pp) block, 81 integrals, 1.1 us: about 10 and 14 ns per integral. One (xx,xx|xx,xx) value through mderi_prim1 takes 0.67 us and one (x,x|x,x) 0.34 us, against 9.67 us and 1.11 us for the legacy lib/gto_md. The water 6-31G build (4186 values, 41,050 primitive quartets) takes 2.7 ms against the legacy 53.8 ms.

PySCF is not always the better reference. For the prims cases, PySCF's value and a 50-digit mpmath evaluation differ by more than 5 × 10⁻¹³ relative in 7 of 1536 values, by up to 3.5 × 10⁻¹¹ for an (ff|ff) pair far apart, where this module is within 5.3 × 10⁻¹⁴ of mpmath. The tests use the mpmath value for those 7 (mderi_prims_fix).

Output is packed-8 (lib/hf/eri_store) over unit-normalized Cartesian AOs.

Used with import hf.md_eri.

Look inside this program

Atlas
The program’s map: every function and which calls which.
There’s no recorded run of this program yet.

Shift+Enter in a cell runs it and moves to the next. Running a cell runs every cell above it too, as one program.

[1]lines 80–872

from hf.basis import (
    AT_X, AT_Y, AT_Z, bas_contract, bas_nao_cart, bas_nao_sph, bas_nc, bas_ncart, BAS_NDIM,
    bas_nr, bas_powers, PR_EXP, SH_ATOM, SH_CART0, SH_L, SH_NCART, SH_NPRIM, SH_NSPH, SH_P0,
    SH_SPH0,
)
import hf.boys
from hf.eri_store import es_checksums, es_len8, es_sample_gap
from hf.hermite import (
    herm_boys_fill, herm_efill, herm_esize, herm_grid, HERM_GRID_XMAX, herm_nh, herm_rfill,
    herm_rfill_box, herm_tuv,
)

MDERI_PI = 3.141592653589793
MDERI_SCREEN_TOL = 1.0e-14
MDERI_HDR = 5
MDERI_TOP = 2
# entries per primitive pair in a pack, (p, Px, Py, Pz), and per pair of
# shell groups in the store's directory, (offset, first shell of each group,
# shells in each group, the two group numbers, 0)
MDERI_PP_LEN = 4
MDERI_DIR_LEN = 8
MDERI_D_SA = 1
MDERI_D_SB = 2
MDERI_D_NA = 3
MDERI_D_NB = 4
MDERI_D_GA = 5
MDERI_D_GB = 6
MDERI_BATCH = 4096
# the four functions of one integral (ab|cd), and the columns of the row
# mderi_prim1 reads for each: (exponent, x, y, z, lx, ly, lz)
MDERI_QUARTET = 4
MDERI_PQ_LEN = 7
# a row of mderi_prims_fix: (case, pick, kind, value)
MDERI_FIX_ROW = 4
# the highest l in lib/hfref/prims: f
MDERI_T_LMAX = 3

#| Entries of one pack for groups of `ca` and `cb` shells of angular
#| momentum `la` and `lb` with `npp` primitive pairs: a 5-number header
#| `[la, lb, npp, nab, nh]`, `(p, Px, Py, Pz)` per primitive pair, then
#| `E_ab` as `[nab, npp * nh]`, `nab = ca cb ncart(la) ncart(lb)` rows.
def mderi_packsize(la, lb, npp, ca, cb) {
    return MDERI_HDR + MDERI_PP_LEN * npp + ca * cb * bas_ncart(la) * bas_ncart(lb) * npp * herm_nh(la + lb)
}

#| Fill the pack for the group pair of `nga` shells from `sa` and `ngb`
#| from `sb` at `ST[off ..]`. `cn` is from `bas_contract`. The shells of a
#| group share their primitives' exponents, so the primitive pairs, their
#| centers and the Hermite expansions are the first shells'; each row of
#| `E_ab` carries one shell pair's contraction coefficients, rows ordered
#| `((a nb_grp + b) ncart_a + i) ncart_b + j` for shell `a` of the first
#| group, `b` of the second, components `i` and `j`.
def mderi_fillpack(ST: Arr, off, atoms: Arr, shells: Arr, prims: Arr, cn: Arr, sa, nga, sb, ngb) {
    la = shells[sa, SH_L]
    lb = shells[sb, SH_L]
    na = bas_ncart(la)
    nb = bas_ncart(lb)
    npa = round(shells[sa, SH_NPRIM])
    npb = round(shells[sb, SH_NPRIM])
    npp = npa * npb

from hf.basis import (
    AT_X, AT_Y, AT_Z, bas_contract, bas_nao_cart, bas_nao_sph, bas_nc, bas_ncart, BAS_NDIM,
    bas_nr, bas_powers, PR_EXP, SH_ATOM, SH_CART0, SH_L, SH_NCART, SH_NPRIM, SH_NSPH, SH_P0,
    SH_SPH0,
)
import hf.boys
from hf.eri_store import es_checksums, es_len8, es_sample_gap
from hf.hermite import (
    herm_boys_fill, herm_efill, herm_esize, herm_grid, HERM_GRID_XMAX, herm_nh, herm_rfill,
    herm_rfill_box, herm_tuv,
)

MDERI_PI = 3.141592653589793
MDERI_SCREEN_TOL = 1.0e-14
MDERI_HDR = 5
MDERI_TOP = 2
# entries per primitive pair in a pack, (p, Px, Py, Pz), and per pair of
# shell groups in the store's directory, (offset, first shell of each group,
# shells in each group, the two group numbers, 0)
MDERI_PP_LEN = 4
MDERI_DIR_LEN = 8
MDERI_D_SA = 1
MDERI_D_SB = 2
MDERI_D_NA = 3
MDERI_D_NB = 4
MDERI_D_GA = 5
MDERI_D_GB = 6
MDERI_BATCH = 4096
# the four functions of one integral (ab|cd), and the columns of the row
# mderi_prim1 reads for each: (exponent, x, y, z, lx, ly, lz)
MDERI_QUARTET = 4
MDERI_PQ_LEN = 7
# a row of mderi_prims_fix: (case, pick, kind, value)
MDERI_FIX_ROW = 4
# the highest l in lib/hfref/prims: f
MDERI_T_LMAX = 3

#| Entries of one pack for groups of `ca` and `cb` shells of angular
#| momentum `la` and `lb` with `npp` primitive pairs: a 5-number header
#| `[la, lb, npp, nab, nh]`, `(p, Px, Py, Pz)` per primitive pair, then
#| `E_ab` as `[nab, npp * nh]`, `nab = ca cb ncart(la) ncart(lb)` rows.
def mderi_packsize(la, lb, npp, ca, cb) {
    return MDERI_HDR + MDERI_PP_LEN * npp + ca * cb * bas_ncart(la) * bas_ncart(lb) * npp * herm_nh(la + lb)
}

#| Fill the pack for the group pair of `nga` shells from `sa` and `ngb`
#| from `sb` at `ST[off ..]`. `cn` is from `bas_contract`. The shells of a
#| group share their primitives' exponents, so the primitive pairs, their
#| centers and the Hermite expansions are the first shells'; each row of
#| `E_ab` carries one shell pair's contraction coefficients, rows ordered
#| `((a nb_grp + b) ncart_a + i) ncart_b + j` for shell `a` of the first
#| group, `b` of the second, components `i` and `j`.
def mderi_fillpack(ST: Arr, off, atoms: Arr, shells: Arr, prims: Arr, cn: Arr, sa, nga, sb, ngb) {
    la = shells[sa, SH_L]
    lb = shells[sb, SH_L]
    na = bas_ncart(la)
    nb = bas_ncart(lb)
    npa = round(shells[sa, SH_NPRIM])
    npb = round(shells[sb, SH_NPRIM])
    npp = npa * npb