Notebooks · lib/gto_md

Angular momentum, by McMurchie and Davidson

lib/gto does s functions, which is hydrogen. Oxygen has p functions, and handling them is the difference between a toy and something that can put water on the screen. This is that difference.

The idea is to stop integrating Cartesian Gaussians directly and expand their product in Hermite Gaussians instead, because a Hermite Gaussian integrates against 1/r in closed form. Two things carry the whole scheme:

The expansion coefficients E(i, j, t). A product of two one-dimensional Cartesian Gaussians of powers i and j is a sum of i + j + 1 Hermite Gaussians, and E says how much of each. E(0, 0, 0) is the Gaussian prefactor, and two recursions reach the rest: one lowers i, the other lowers j. Overlap is E(i, j, 0) times sqrt(pi/p), per dimension, and kinetic energy is three E values combined.

The Hermite Coulomb integrals R(t, u, v, n). These are what 1/r becomes. R(0, 0, 0, n) is (-2p)^n F_n(pR^2), where F_n is the Boys function, and raising any index costs one more recursion step. Nuclear attraction is a sum of Es against Rs; a two-electron integral is a sum of Es against Es against Rs, with a sign that alternates on the second pair.

So the Boys function is needed at every order up to the total angular momentum, rather than only at zero. It comes from its series for small arguments and its asymptotic form for large ones, and where to change over is a measured question. The largest relative gap between the two, over orders 0 to 5:

x       20      3.6e-5
x       30      9.3e-9
x       40      1.5e-12
x       50      1.2e-15

So the changeover is at 50, where they agree to the last bit. At 30, which is where this started, an integral with an argument near the join would have been wrong in the ninth digit and the energy would have followed.

Everything is held against PySCF through lib/qcref, one integral at a time and then the energy: the Hartree-Fock energies of water and of hydrogen fluoride land within 10⁻¹⁰ hartree of the ones PySCF printed there (WATER_ERHF and HF_ERHF).

Used with import gto_md.

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 47–750

from gto import boys0, GTO_DIM, gto_norm_s, GTO_PI, gto_rhf, gto_rhf_full, gto_worst
from qcref import (
    water631g_bf_pow, water631g_bf_shell, water631g_charges, water631g_coords, WATER631G_ENUC,
    WATER631G_ERHF, WATER631G_NATM, WATER631G_NBF, WATER631G_NOCC, WATER631G_NPRIM, water631g_S,
    water631g_shell_atom, water631g_shell_coef, water631g_shell_expo, water631g_shell_nprim,
    water631g_T, water631g_V, water_bf_pow, water_bf_shell, water_charges, water_coords,
    WATER_ENUC, WATER_ERHF, water_ERI, WATER_NATM, WATER_NBF, WATER_NOCC, WATER_NPRIM, water_S,
    water_shell_atom, water_shell_coef, water_shell_expo, water_shell_nprim, water_T, water_V,
    hf_bf_pow, hf_bf_shell, hf_charges, hf_coords, HF_ENUC, HF_ERHF, hf_ERI, HF_NATM, HF_NBF,
    HF_NOCC, HF_NPRIM, hf_S, hf_shell_atom, hf_shell_coef, hf_shell_expo, hf_shell_nprim, hf_T,
    hf_V,
)

MD_BOYS_SWITCH = 50.0

#| The most terms the Boys series may take below the switch; it stops as soon
#| as a term no longer changes the sum.
MD_BOYS_MAX_TERMS = 300

#| How many SCF iterations the water energies run.
MD_WATER_SCF_ITERS = 60

#| The Boys function `F_n(x) = int_0^1 t^2n exp(-x t^2) dt`, at any order.
#|
#| Small `x`: the series `exp(-x) sum_k (2x)^k / (2n+1)(2n+3)...(2n+2k+1)`,
#| summed until a term stops mattering -- the terms grow before they shrink, and
#| the peak is near `k = x`, so a fixed count is wrong at the top of the range.
#| Large `x`: the integrand
#| has decayed long before `t = 1`, so the limit `(2n-1)!!/2^(n+1) sqrt(pi/x^(2n+1))`
#| takes over. The test below checks they agree across the join.
def boys_n(n, x) {
    if x < 1.0e-12 {
        return 1.0 / (2.0 * n + 1.0)
    }
    if x > MD_BOYS_SWITCH {
        dfac = 1.0
        for k to n {
            dfac = dfac * (2.0 * k + 1.0)
        }
        return dfac / pow(2.0, n + 1.0) * sqrt(GTO_PI / pow(x, 2.0 * n + 1.0))
    }
    term = 1.0 / (2.0 * n + 1.0)
    total = term
    k = 0
    while k < MD_BOYS_MAX_TERMS {
        term = term * 2.0 * x / (2.0 * n + 2.0 * k + 3.0)
        total = total + term
        if term < 1.0e-17 * total {
            break
        }
        k = k + 1
    }
    return exp(-x) * total
}

#| `E(i, j, t)`: how much of the `t`-th Hermite Gaussian is in the product of
#| two one-dimensional Cartesian Gaussians of powers `i` and `j`, separated by
#| `Q` and with exponents `a` and `b`.
def herm_e(i, j, t, Q, a, b) {
    p = a + b

from gto import boys0, GTO_DIM, gto_norm_s, GTO_PI, gto_rhf, gto_rhf_full, gto_worst
from qcref import (
    water631g_bf_pow, water631g_bf_shell, water631g_charges, water631g_coords, WATER631G_ENUC,
    WATER631G_ERHF, WATER631G_NATM, WATER631G_NBF, WATER631G_NOCC, WATER631G_NPRIM, water631g_S,
    water631g_shell_atom, water631g_shell_coef, water631g_shell_expo, water631g_shell_nprim,
    water631g_T, water631g_V, water_bf_pow, water_bf_shell, water_charges, water_coords,
    WATER_ENUC, WATER_ERHF, water_ERI, WATER_NATM, WATER_NBF, WATER_NOCC, WATER_NPRIM, water_S,
    water_shell_atom, water_shell_coef, water_shell_expo, water_shell_nprim, water_T, water_V,
    hf_bf_pow, hf_bf_shell, hf_charges, hf_coords, HF_ENUC, HF_ERHF, hf_ERI, HF_NATM, HF_NBF,
    HF_NOCC, HF_NPRIM, hf_S, hf_shell_atom, hf_shell_coef, hf_shell_expo, hf_shell_nprim, hf_T,
    hf_V,
)

MD_BOYS_SWITCH = 50.0

#| The most terms the Boys series may take below the switch; it stops as soon
#| as a term no longer changes the sum.
MD_BOYS_MAX_TERMS = 300

#| How many SCF iterations the water energies run.
MD_WATER_SCF_ITERS = 60

#| The Boys function `F_n(x) = int_0^1 t^2n exp(-x t^2) dt`, at any order.
#|
#| Small `x`: the series `exp(-x) sum_k (2x)^k / (2n+1)(2n+3)...(2n+2k+1)`,
#| summed until a term stops mattering -- the terms grow before they shrink, and
#| the peak is near `k = x`, so a fixed count is wrong at the top of the range.
#| Large `x`: the integrand
#| has decayed long before `t = 1`, so the limit `(2n-1)!!/2^(n+1) sqrt(pi/x^(2n+1))`
#| takes over. The test below checks they agree across the join.
def boys_n(n, x) {
    if x < 1.0e-12 {
        return 1.0 / (2.0 * n + 1.0)
    }
    if x > MD_BOYS_SWITCH {
        dfac = 1.0
        for k to n {
            dfac = dfac * (2.0 * k + 1.0)
        }
        return dfac / pow(2.0, n + 1.0) * sqrt(GTO_PI / pow(x, 2.0 * n + 1.0))
    }
    term = 1.0 / (2.0 * n + 1.0)
    total = term
    k = 0
    while k < MD_BOYS_MAX_TERMS {
        term = term * 2.0 * x / (2.0 * n + 2.0 * k + 3.0)
        total = total + term
        if term < 1.0e-17 * total {
            break
        }
        k = k + 1
    }
    return exp(-x) * total
}

#| `E(i, j, t)`: how much of the `t`-th Hermite Gaussian is in the product of
#| two one-dimensional Cartesian Gaussians of powers `i` and `j`, separated by
#| `Q` and with exponents `a` and `b`.
def herm_e(i, j, t, Q, a, b) {
    p = a + b