examples/scf_sto3g.viv

# Molecule-driven restricted Hartree-Fock in STO-3G.
#
# The basis is assembled from an element library (H, He, C, N, O, F) given only
# the atoms; the McMurchie-Davidson engine handles the s/p integrals. This is the
# reusable driver behind the molecule database (tests/molecule_db.rs), which
# substitutes the MOL_BEGIN..MOL_END block per molecule and checks the energy
# against known values. Default molecule = water (E = -74.9420787, Crawford).

pi = 3.141592653589793
n_prim = 3

# =====================================================================
# 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 {
                res += (v - 1.0) * hermite_R(t, u, v - 2.0, nn + 1.0, p, PCx, PCy, PCz, RPC2)
            }
            return res
        }
        res = PCy * hermite_R(t, u - 1.0, v, nn + 1.0, p, PCx, PCy, PCz, RPC2)