Multi-centre one-electron integrals

Quantum Chemistry

Up to this point every integral has lived at a single centre. The helium calculation had one nucleus, the basis sat on top of it, and the nuclear attraction was just $- Z / r$ evaluated against same-origin Gaussians. H₂ breaks all of that. The basis functions live on two different atoms, the nuclear-attraction operator runs over both nuclei in turn, and even the overlap between two basis functions involves Gaussians with different centres. This page builds the multi-centre one-electron integrals that the molecule needs.

The Gaussian product theorem

Two normalized $1 s$ Gaussians at different centres multiply into a single Gaussian at a new centre:

e^{- α ∣ r - A ∣^{2}} e^{- β ∣ r - B ∣^{2}} = K_{A B} e^{- (α + β) ∣ r - P ∣^{2}},

where the new centre is the weighted midpoint

P = \frac{α A + β B}{α + β},

and the Gaussian-product factor — the suppression from the separation of the two centres — is

K_{A B} = exp (- \frac{α β}{α + β} ∣ A - B ∣^{2}) .

This is the trick that makes Gaussians practical. Slater orbitals don't have an analogous identity — the product of two Slater functions at different centres is some ugly thing without a closed form. Boys saw in 1950 that with Gaussians, every multi-centre integral collapses to an integral over a single Gaussian at a single point, with an analytic prefactor that captures the geometry. We're going to lean on this for every formula on the page.

Overlap and kinetic

Apply the Gaussian-product theorem to the overlap integral and the radial integral evaluates to $(π / (α + β))^{3/2}$ by the standard three-dimensional Gaussian normalization. Including the primitive normalization factors, the result for normalized 1s primitives is

S (α, A; β, B) = (\frac{4 α β}{( α + β ) ^{2}})^{3/4} exp (- \frac{α β}{α + β} ∣ A - B ∣^{2}) .

The same-centre limit $A = B$ recovers the formula we used in step 3. The kinetic-energy integral is obtained by acting $- \frac{1}{2} \nabla^{2}$ on the second Gaussian and computing the result in closed form:

T (α, A; β, B) = \frac{α β}{α + β} [3 - \frac{2 α β}{α + β} ∣ A - B ∣^{2}] S .

Same-centre limit $∣ A - B ∣^{2} = 0$ gives $T = 3 α β / (α + β) S$ , again the formula from step 3.

The Boys function

Nuclear attraction needs us to integrate $1/∣ r - C ∣$ against a Gaussian. There is no elementary closed form for this, but a one-dimensional integral representation does exist: the Boys function

F_{0} (x) = \int_{0}^{1} e^{- x t^{2}} d t = \frac{1}{2} \frac{π}{x} erf (x) (x > 0),

with limits $F_{0} (0) = 1$ and $F_{0} (x) \sim π / (4 x)$ as $x \to \infty$ . The error-function form is fine for moderate $x$ ; for $x$ very close to zero it has a removable $0/0$ that we sidestep with the Taylor expansion $F_{0} (x) = 1 - x /3 + x^{2} /10 - \dots$ . We'll see this function come up again on the next page in the multi-centre two-electron integral; it's the single function that encodes Coulomb's law inside a Gaussian basis. Higher-angular- momentum bases involve $F_{n} (x)$ for $n > 0$ ; we only need $n = 0$ .

Nuclear attraction

With $F_{0}$ in hand, the nuclear-attraction integral over two normalized 1s Gaussians and a nucleus at $C$ with charge $Z$ becomes

V_{C} (α, A; β, B) = - Z \cdot 2 \frac{α + β}{π} S \cdot F_{0} ((α + β) ∣ P - C ∣^{2}) .

The Boys argument is the squared distance from the Gaussian-product centre to the attracting nucleus, scaled by the combined exponent. Same-centre cases ( $P = C$ ) drop the Boys factor to $1$ and recover the step-3 expression. For the molecule we sum $V_{C}$ over every nucleus, and that sum is the total nuclear-attraction matrix element for the pair of Gaussians.

From primitives to AOs

H₂ in STO-3G has two atomic orbitals, one $1 s$ per hydrogen, each built from three primitive Gaussians sharing that atom's centre. Six primitives total, two AOs total. Every primitive matrix is $6 \times 6$ ; every AO matrix is $2 \times 2$ . The contraction is block-diagonal in atom: the AO on atom $a$ is

ϕ_{a} (r) = N_{a} i \in a \sum d_{i} g_{i} (r - R_{a}),

where $N_{a}$ renormalises the contraction so that $⟨ ϕ_{a} ∣ ϕ_{a} ⟩ = 1$ exactly (the published STO-3G coefficients give a contraction that is only approximately normalized). For the two-by-two AO matrix we just contract the primitive matrix and divide by $N_{a} N_{b}$ on the matrix element $(a, b)$ .

The code (new this page)

The Boys function, the three multi-centre primitive integrals, and a routine that builds H₂'s $(S, T, V)$ in the AO basis at a given bond length.

from scipy.special import erf


# ---- Boys function F_0 (the Coulomb-from-Gaussian engine) ----

def boys_F0(x):
    """F_0(x) = int_0^1 exp(-x t^2) dt.

    Closed form: 0.5 * sqrt(pi/x) * erf(sqrt(x)) for x > 0.
    For x -> 0 use the Taylor series F_0(x) = 1 - x/3 + x^2/10 - ...
    so that the small-x branch is numerically stable.
    """
    if x < 1e-8:
        return 1.0 - x / 3.0 + x * x / 10.0
    return 0.5 * np.sqrt(np.pi / x) * erf(np.sqrt(x))


# ---- two-centre primitive 1s integrals ----

def overlap_2c(a, A, b, B):
    """<g_a(r-A) | g_b(r-B)> for normalized 1s primitive Gaussians."""
    R2 = np.sum((A - B) ** 2)
    p  = a + b
    return (4 * a * b / p ** 2) ** 0.75 * np.exp(-a * b * R2 / p)


def kinetic_2c(a, A, b, B):
    """<g_a | -1/2 nabla^2 | g_b>  (normalized 1s primitives)."""
    R2 = np.sum((A - B) ** 2)
    p  = a + b
    return a * b / p * (3.0 - 2.0 * a * b * R2 / p) * overlap_2c(a, A, b, B)


def nuclear_attraction_2c(a, A, b, B, C, Z):
    """<g_a(r-A) | -Z/|r-C| | g_b(r-B)>  (normalized 1s primitives)."""
    p  = a + b
    P  = (a * A + b * B) / p           # Gaussian product centre
    R2 = np.sum((P - C) ** 2)
    return -Z * 2.0 * np.sqrt(p / np.pi) * overlap_2c(a, A, b, B) * boys_F0(p * R2)


# ---- H2 in the STO-3G basis: AO matrices at a given bond length ----

def h2_basis(R):
    """Return (alphas, centres, atom_of) for H2 in STO-3G primitive basis."""
    Apos = np.array([0.0, 0.0, 0.0])
    Bpos = np.array([R,   0.0, 0.0])
    alphas  = np.concatenate([STO_3G_H_ALPHAS, STO_3G_H_ALPHAS])
    centres = np.array([Apos, Apos, Apos, Bpos, Bpos, Bpos])
    atom_of = np.array([0, 0, 0, 1, 1, 1])
    return alphas, centres, atom_of


def contract_to_ao(M_prim, atom_of, n_ao=2, d=STO_3G_COEFFS):
    """Block-contract a primitive matrix into the AO basis.

    Each AO is sum_i d_i g_i over the primitives belonging to one atom.
    """
    M_ao = np.zeros((n_ao, n_ao))
    for a in range(n_ao):
        ia = np.where(atom_of == a)[0]
        for b in range(n_ao):
            ib = np.where(atom_of == b)[0]
            for ii, i in enumerate(ia):
                for jj, j in enumerate(ib):
                    M_ao[a, b] += d[ii] * d[jj] * M_prim[i, j]
    return M_ao


def h2_one_electron_matrices(R):
    """Return AO-basis (S, T, V) for H2 at bond length R (bohr)."""
    alphas, centres, atom_of = h2_basis(R)
    nuclei = [(np.array([0.0, 0.0, 0.0]), 1),
              (np.array([R,   0.0, 0.0]), 1)]
    n = len(alphas)

    S_prim = np.zeros((n, n))
    T_prim = np.zeros((n, n))
    V_prim = np.zeros((n, n))
    for i in range(n):
        for j in range(n):
            S_prim[i, j] = overlap_2c(alphas[i], centres[i], alphas[j], centres[j])
            T_prim[i, j] = kinetic_2c(alphas[i], centres[i], alphas[j], centres[j])
            for C, Z in nuclei:
                V_prim[i, j] += nuclear_attraction_2c(alphas[i], centres[i],
                                                      alphas[j], centres[j], C, Z)

    S = contract_to_ao(S_prim, atom_of)
    T = contract_to_ao(T_prim, atom_of)
    V = contract_to_ao(V_prim, atom_of)

    # AO normalization: each contracted AO should satisfy <phi|phi> = 1.
    norms = np.sqrt(np.diag(S))
    Norm  = np.outer(norms, norms)
    return S / Norm, T / Norm, V / Norm

A note on the structure: all geometry-dependent code now uses explicit centre vectors. The atom-only integrals from earlier pages still work — they're just the special case $A = B$ — and we keep them around for the helium calculation. The molecular code uses only the new _2c versions.

Run it

H2 at R = 1.4 bohr (STO-3G AO basis):

  S =
    [[ 1.000000  0.659318]
     [ 0.659318  1.000000]]

  T =
    [[ 0.760032  0.236455]
     [ 0.236455  0.760032]]

  V =
    [[-1.880441 -1.194835]
     [-1.194835 -1.880441]]

  h = T + V =
    [[-1.120409 -0.958380]
     [-0.958380 -1.120409]]

Three things to read off these matrices. (1) The off-diagonal overlap $S_{12} = 0.6593$ is the AO overlap of two hydrogen 1s functions separated by 1.4 bohr — large, because Gaussians at this distance still have substantial mass between them. The diagonals are $1$ by construction (AO normalization). (2) The kinetic matrix is symmetric and positive on the diagonal, with a smaller positive off-diagonal — the "kinetic coupling" between the two AOs.

(3) The nuclear-attraction matrix is where the molecule's geometry shows up most clearly. The diagonal $V_{11}$ is the energy of an electron in $ϕ_{A}$ attracted to both nuclei: most of the $- 1.880$ Hartree comes from the on-atom $A$ nucleus, with a smaller contribution from the off-atom $B$ . The off-diagonal $V_{12}$ is the cross-attraction between the two AOs summed over both nuclei, which is what builds the bonding interaction in H₂: when $V_{12}$ is more negative than $T_{12}$ is positive, the bonding orbital drops below the antibonding orbital and the molecule is stable. On page 9 we'll add the two-electron repulsion that fights against this bonding tendency.

Why does the Gaussian product theorem hold but a Slater product theorem doesn't?

Because Gaussians factorize into independent Cartesian components: $e^{- α r^{2}} = e^{- α x^{2}} e^{- α y^{2}} e^{- α z^{2}}$ . The product of Gaussians at different centres also factorizes, and one-dimensional Gaussian products are again one-dimensional Gaussians. Slater orbitals $e^{- ζ r}$ involve the Euclidean radius $r$ , which doesn't factorize into anything tractable when the centre changes. The price we pay for using Gaussians instead of Slater orbitals is that individual Gaussians have the wrong shape (no cusp, wrong tail); the price we pay for using Slater orbitals would be that you can't analytically integrate any non-trivial molecule. Quantum chemistry chose the first option in 1950 and never looked back.

What's the physical meaning of the Gaussian product centre P?

It's the centre of mass of the joint charge distribution $g_{α} (r - A) g_{β} (r - B)$ , with the "mass" of each Gaussian equal to its exponent. Steeper Gaussians (larger $α$ ) pull the product centre toward themselves because their decay length is shorter. In practice $P$ is just a convenient bookkeeping centre; nothing physical lives at it.

Why does the Boys function appear?

The Coulomb potential $1/ r$ has the Laplace representation $1/ r = (2/ π) \int_{0}^{\infty} e^{- t^{2} r^{2}} d t$ . Substituting that into the nuclear-attraction integral converts the nasty $1/ r$ into another Gaussian, which multiplies cleanly with the existing Gaussians (Gaussian product theorem again), integrates analytically over space, and leaves a one-dimensional integral over $t$ . After a change of variables, that one-dimensional integral is exactly $F_{0}$ . Higher-angular-momentum integrals produce $F_{n}$ for $n > 0$ , which can be computed from $F_{0}$ by recursion. The Boys family is the entire engine of Gaussian-basis Coulomb integrals.