examples/scf_water_631gd.viv

# Restricted Hartree-Fock for WATER (H2O) in 6-31G(d) [6-31G*], with a general
# McMurchie-Davidson integral engine over s, p, AND d Cartesian Gaussians.
#
# A real polyatomic calculation with a polarized split-valence basis (19 basis
# functions, variable contraction lengths, 6 Cartesian d functions on oxygen).
# At the (stretched) Crawford geometry, E(SCF) = -76.0192 hartree; the basis
# series STO-3G (-74.942) > 6-31G (-75.936) > 6-31G* (-76.019) descends
# monotonically toward the HF limit (~-76.067), exactly as variational theory
# requires. The s/p engine is anchored by the exact STO-3G match in
# scf_water.viv (-74.9420787).

pi = 3.141592653589793
# Max primitives per contracted function; shorter contractions pad with
# coefficient 0 (skipped in the contraction loops below).
n_prim = 6

# =====================================================================
# McMurchie-Davidson integral engine (molecule-independent).
# =====================================================================

# Boys function F_n(T) = integral_0^1 t^(2n) exp(-T t^2) dt.
def boys_n(nn, T) {
    if T < 1.0e-12 {
        return 1.0 / (2.0 * nn + 1.0)
    }
    return simpson(fn(t) -> pow(t, 2.0 * nn) * exp(-T * t * t), 0.0, 1.0, 100)
}

# Hermite expansion coefficient E_t^{ij} (1D), exponents a,b, Qx = Ax - Bx.
def hermite_E(i, j, t, Qx, a, b) {
    if i < 0.0 { return 0.0 }
    if j < 0.0 { return 0.0 }
    if t < 0.0 { return 0.0 }
    if t > i + j { return 0.0 }
    p = a + b
    if i == 0.0 {
        if j == 0.0 {
            mu = a * b / p
            return exp(-mu * Qx * Qx)
        }
        xpb = a / p * Qx
        return (0.5 / p * hermite_E(i, j - 1.0, t - 1.0, Qx, a, b)
             + xpb * hermite_E(i, j - 1.0, t, Qx, a, b)
             + (t + 1.0) * hermite_E(i, j - 1.0, t + 1.0, Qx, a, b))
    }
    xpa = -b / p * Qx
    return (0.5 / p * hermite_E(i - 1.0, j, t - 1.0, Qx, a, b)
         + xpa * hermite_E(i - 1.0, j, t, Qx, a, b)
         + (t + 1.0) * hermite_E(i - 1.0, j, t + 1.0, Qx, a, b))
}

# Hermite Coulomb integral R_{tuv}^n.
def hermite_R(t, u, v, nn, p, PCx, PCy, PCz, RPC2) {
    if t == 0.0 {
        if u == 0.0 {
            if v == 0.0 {
                return pow(-2.0 * p, nn) * boys_n(nn, p * RPC2)
            }
            res = PCz * hermite_R(t, u, v - 1.0, nn + 1.0, p, PCx, PCy, PCz, RPC2)
            if v > 1.0 {