Notebooks · lib/hf/hermite

Hermite tables for McMurchie-Davidson integrals, built bottom-up

Every integral in Hartree-Fock over Gaussians goes through the same move. The product of two Cartesian Gaussians on centers A and B, exponents a and b, is a Gaussian on the point P = (aA + bB)/p with exponent p = a + b, times a polynomial. McMurchie and Davidson write that polynomial as a sum of Hermite Gaussians, which are derivatives of the plain Gaussian with respect to P:

x_A^i x_B^j exp(-a x_A^2 - b x_B^2) = sum_t E(i, j, t) d^t/dP_x^t exp(-p x_P^2)

one axis at a time. E(i, j, t) is nonzero only for t <= i + j.

Why bother: a derivative with respect to P passes straight through an integral over r. So the overlap keeps only the t = 0 term, E(i, j, 0) sqrt(pi/p), and the Coulomb integrals become derivatives of one function of P. For 1/|r - C| that function is (2 pi / p) F_0(p |P - C|^2), and its derivatives are

R(t, u, v) = d^t/dX^t d^u/dY^u d^v/dZ^v  F_0(alpha (X^2 + Y^2 + Z^2))

evaluated at (X, Y, Z) = P - C, with alpha = p for a nucleus and alpha = pq/(p + q) for a pair of pairs.

E, bottom-up. Start from E(0, 0, 0) = exp(-mu X_AB^2), where mu = ab/p and X_AB = A_x - B_x; that is the Gaussian product prefactor for one axis, so the product over x, y and z is K_AB. Then raise i and j one step at a time:

E(i+1, j, t) = E(i, j, t-1) / (2p) + X_PA E(i, j, t) + (t+1) E(i, j, t+1)
E(i, j+1, t) = E(i, j, t-1) / (2p) + X_PB E(i, j, t) + (t+1) E(i, j, t+1)

with X_PA = -b X_AB / p and X_PB = a X_AB / p. Every entry is computed once. The legacy code in lib/gto_md recursed downward for each entry it wanted, which recomputed the same values about 3^(i+j) times; at (dd|dd) that was most of its time.

R, bottom-up. The derivatives come from an auxiliary index n:

R(n; 0, 0, 0)   = (-2 alpha)^n F_n(alpha R^2)
R(n; t+1, u, v) = t R(n+1; t-1, u, v) + X R(n+1; t, u, v)

and the same with u, Y and v, Z. The answer is R(0; t, u, v). Level n needs only level n + 1, so two cubes of side L + 1 hold everything, alternating by the parity of n, and one Boys evaluation at the top order L serves every entry.

Used with import hf.hermite.

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 52–594

from hf.basis import BAS_NDIM, bas_odfact
from hf.boys import (
    BOYS_GRID_HDR, boys_grid_new, boys_series_one, boys_switch, BOYS_TAYLOR_K, boys_values,
)

HERM_PI = 3.141592653589793
# the top of the Boys grid an integral build makes: past every order's switch
# (34.4 at most, lib/hf/boys), so the grid spans all of the series side
HERM_GRID_XMAX = 40.0

#| Hermite functions with `t + u + v <= L`.
def herm_nh(L) {
    return (L + 1) * (L + 2) * (L + 3) / 6
}

#| The Hermite index list for total order `L`, `[nh, 3]`: `t` counts up, then
#| `u`, then `v`. Every table in this family uses this order.
def herm_tuv(L) {
    H = zeros(herm_nh(L), BAS_NDIM)
    h = 0
    for t to L + 1 {
        for u to L - t + 1 {
            for v to L - t - u + 1 {
                H[h, 0] = t
                H[h, 1] = u
                H[h, 2] = v
                h = h + 1
            }
        }
    }
    return H
}

#| Entries of a flat E table for one axis with powers up to `la` and `lb`.
def herm_esize(la, lb) {
    return (la + 1) * (lb + 1) * (la + lb + 1)
}

#| Fill the flat 1-D table `E(i, j, t)`, `i <= la`, `j <= lb`, `t <= i + j`,
#| into `E` from entry `e0` on, at `e0 + (i (lb + 1) + j) (la + lb + 1) + t`.
#| `X` is `A - B` on this axis. The table is overwritten, so a caller can
#| reuse one buffer for every primitive pair.
def herm_efill(E: Arr, e0, la, lb, a, b, X) {
    #: the table starts from E(0, 0, 0) = exp(-ab/p X^2)
    nt = la + lb + 1
    ld = lb + 1
    for k to herm_esize(la, lb) {
        E[e0 + k] = 0.0
    }
    p = a + b
    o2p = 0.5 / p
    xpa = -b * X / p
    xpb = a * X / p
    E[e0] = exp(-a * b / p * X * X)
    # raise i along j = 0
    for i to la {
        src = e0 + (i * ld) * nt
        dst = e0 + ((i + 1) * ld) * nt
        for t to i + 2 {
            v = 0.0

from hf.basis import BAS_NDIM, bas_odfact
from hf.boys import (
    BOYS_GRID_HDR, boys_grid_new, boys_series_one, boys_switch, BOYS_TAYLOR_K, boys_values,
)

HERM_PI = 3.141592653589793
# the top of the Boys grid an integral build makes: past every order's switch
# (34.4 at most, lib/hf/boys), so the grid spans all of the series side
HERM_GRID_XMAX = 40.0

#| Hermite functions with `t + u + v <= L`.
def herm_nh(L) {
    return (L + 1) * (L + 2) * (L + 3) / 6
}

#| The Hermite index list for total order `L`, `[nh, 3]`: `t` counts up, then
#| `u`, then `v`. Every table in this family uses this order.
def herm_tuv(L) {
    H = zeros(herm_nh(L), BAS_NDIM)
    h = 0
    for t to L + 1 {
        for u to L - t + 1 {
            for v to L - t - u + 1 {
                H[h, 0] = t
                H[h, 1] = u
                H[h, 2] = v
                h = h + 1
            }
        }
    }
    return H
}

#| Entries of a flat E table for one axis with powers up to `la` and `lb`.
def herm_esize(la, lb) {
    return (la + 1) * (lb + 1) * (la + lb + 1)
}

#| Fill the flat 1-D table `E(i, j, t)`, `i <= la`, `j <= lb`, `t <= i + j`,
#| into `E` from entry `e0` on, at `e0 + (i (lb + 1) + j) (la + lb + 1) + t`.
#| `X` is `A - B` on this axis. The table is overwritten, so a caller can
#| reuse one buffer for every primitive pair.
def herm_efill(E: Arr, e0, la, lb, a, b, X) {
    #: the table starts from E(0, 0, 0) = exp(-ab/p X^2)
    nt = la + lb + 1
    ld = lb + 1
    for k to herm_esize(la, lb) {
        E[e0 + k] = 0.0
    }
    p = a + b
    o2p = 0.5 / p
    xpa = -b * X / p
    xpb = a * X / p
    E[e0] = exp(-a * b / p * X * X)
    # raise i along j = 0
    for i to la {
        src = e0 + (i * ld) * nt
        dst = e0 + ((i + 1) * ld) * nt
        for t to i + 2 {
            v = 0.0