Notebooks · lib/gto

Gaussian-basis quantum chemistry, held against PySCF

Molecular electronic structure from the Gaussians up: contracted s-type basis functions, the four integrals a Hartree-Fock calculation needs, and the self-consistent field that turns them into an energy. Every number is checked against PySCF's own, which tools/pyscf_reference.py writes into lib/qcref with enough digits to round-trip. What separates PySCF's answer from this one is the arithmetic in between.

The four integrals, for primitives exp(-a|r-A|^2):

Overlap. Two Gaussians multiply to a third, at P = (aA + bB)/(a + b), scaled by K = exp(-ab/(a+b) |A-B|^2). Integrating that gives K (pi/p)^(3/2) and nothing harder.
Kinetic. The Laplacian of a Gaussian is a Gaussian, so this is the overlap times ab/p (3 - 2ab/p |A-B|^2).
Nuclear attraction. 1/|r-C| is not Gaussian, and this is where the Boys function F0(x) = sqrt(pi/4x) erf(sqrt(x)) arrives. It is the one special function this module calls, and it is built from erf.
Two-electron. Two Gaussian products, one Boys function between them. For n basis functions there are n^4 of these, which is why everything past a small molecule is about not computing them.

Then the self-consistent field: the Fock matrix depends on the density, which comes from the Fock matrix. Start with a zero density, solve F C = S C e in an orthonormal basis found by Cholesky, rebuild the density, repeat.

Only s functions, which means hydrogen. Oxygen has p, and the recursion that handles angular momentum is the next rung.

Used with import gto.

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 36–471

from linalg import inverse
from qcref import (
    h2_bf_atom, h2_charges, h2_coords, H2_ENUC, H2_ERHF, h2_ERI, h2_long_bf_atom,
    h2_long_charges, h2_long_coords, H2_LONG_ENUC, H2_LONG_ERHF, H2_LONG_NATM, H2_LONG_NBF,
    H2_LONG_NOCC, H2_NATM, H2_NBF, H2_NOCC, h2_S, h2_T, h2_V, h4_bf_atom, h4_charges, h4_coords,
    H4_ENUC, H4_ERHF, H4_NATM, H4_NBF, H4_NOCC, h4_S, h4_T, h4_V, NPRIM, sto3g_h_coef,
    sto3g_h_expo,
)

GTO_PI = 3.141592653589793

#| Space has three dimensions: a position, and a Cartesian power, has this many
#| components.
GTO_DIM = 3

#| How many SCF iterations `hydrogens_energy` runs; a chain of STO-3G
#| hydrogens has settled long before.
GTO_SCF_ITERS = 60

#| A primitive s Gaussian's normalization, `(2a/pi)^(3/4)`.
def gto_norm_s(a) {
    return pow(2.0 * a / GTO_PI, 0.75)
}

#| The Boys function of order zero. Both the nuclear attraction and the
#| two-electron integrals get their distance dependence from it, and it is the
#| only thing here that is not elementary. Near zero it is taken from its series,
#| because `erf(sqrt(x))/sqrt(x)` is 0/0 there.
def boys0(x) {
    if x < 1.0e-12 {
        return 1.0 - x / 3.0
    }
    return 0.5 * sqrt(GTO_PI / x) * erf(sqrt(x))
}

#| The squared distance between rows `i` and `j` of a coordinate table.
def gto_dist2(R: Arr, i, j) {
    total = 0.0
    for k to GTO_DIM {
        d = R[i, k] - R[j, k]
        total = total + d * d
    }
    return total
}

#| The contracted coefficients: the tabulated ones times each primitive's own
#| normalization, and then the whole contraction scaled so that it overlaps
#| itself exactly once.
#|
#| That last step is not decoration. The tabulated STO-3G coefficients are
#| rounded to eight figures, so the contraction they define is normalized only
#| to about 1e-8. PySCF renormalizes and so must anything held against it: the
#| overlap matrix differs in the ninth digit otherwise. The Hartree-Fock energy
#| does not, because a common scale on every basis function cancels between the
#| integrals and the density, which is exactly the kind of error that hides.
def gto_contraction(expo: Arr, coef: Arr, nprim) {
    out = zeros(nprim)
    for i to nprim {
        out[i] = coef[i] * gto_norm_s(expo[i])
    }

from linalg import inverse
from qcref import (
    h2_bf_atom, h2_charges, h2_coords, H2_ENUC, H2_ERHF, h2_ERI, h2_long_bf_atom,
    h2_long_charges, h2_long_coords, H2_LONG_ENUC, H2_LONG_ERHF, H2_LONG_NATM, H2_LONG_NBF,
    H2_LONG_NOCC, H2_NATM, H2_NBF, H2_NOCC, h2_S, h2_T, h2_V, h4_bf_atom, h4_charges, h4_coords,
    H4_ENUC, H4_ERHF, H4_NATM, H4_NBF, H4_NOCC, h4_S, h4_T, h4_V, NPRIM, sto3g_h_coef,
    sto3g_h_expo,
)

GTO_PI = 3.141592653589793

#| Space has three dimensions: a position, and a Cartesian power, has this many
#| components.
GTO_DIM = 3

#| How many SCF iterations `hydrogens_energy` runs; a chain of STO-3G
#| hydrogens has settled long before.
GTO_SCF_ITERS = 60

#| A primitive s Gaussian's normalization, `(2a/pi)^(3/4)`.
def gto_norm_s(a) {
    return pow(2.0 * a / GTO_PI, 0.75)
}

#| The Boys function of order zero. Both the nuclear attraction and the
#| two-electron integrals get their distance dependence from it, and it is the
#| only thing here that is not elementary. Near zero it is taken from its series,
#| because `erf(sqrt(x))/sqrt(x)` is 0/0 there.
def boys0(x) {
    if x < 1.0e-12 {
        return 1.0 - x / 3.0
    }
    return 0.5 * sqrt(GTO_PI / x) * erf(sqrt(x))
}

#| The squared distance between rows `i` and `j` of a coordinate table.
def gto_dist2(R: Arr, i, j) {
    total = 0.0
    for k to GTO_DIM {
        d = R[i, k] - R[j, k]
        total = total + d * d
    }
    return total
}

#| The contracted coefficients: the tabulated ones times each primitive's own
#| normalization, and then the whole contraction scaled so that it overlaps
#| itself exactly once.
#|
#| That last step is not decoration. The tabulated STO-3G coefficients are
#| rounded to eight figures, so the contraction they define is normalized only
#| to about 1e-8. PySCF renormalizes and so must anything held against it: the
#| overlap matrix differs in the ninth digit otherwise. The Hartree-Fock energy
#| does not, because a common scale on every basis function cancels between the
#| integrals and the density, which is exactly the kind of error that hides.
def gto_contraction(expo: Arr, coef: Arr, nprim) {
    out = zeros(nprim)
    for i to nprim {
        out[i] = coef[i] * gto_norm_s(expo[i])
    }