Many-electron wavefunctions

Quantum Chemistry

Up to here the program has been a one-electron calculator. To get to H₂ we need it to handle two electrons, which means facing two facts at once: the wavefunction has to be antisymmetric under exchange, and the two electrons see each other through $1/ r_{12}$ . This page introduces the antisymmetric structure (the Slater determinant) and writes the two-electron integral that the structure forces into the energy. We test the new machinery on helium — same single-centre geometry as hydrogen, one extra electron — before turning to H₂ with its two centres.

Why a Hartree product is wrong

The naive guess for a two-electron wavefunction is the product of two single-particle states:

ψ_{Hartree} (r_{1}, r_{2}) = ϕ (r_{1}) ϕ (r_{2}) .

This is wrong, and the way it's wrong is structural rather than quantitative. Electrons are fermions, and any two-electron wavefunction must change sign when we swap the two electrons: $ψ (r_{1}, r_{2}) = - ψ (r_{2}, r_{1})$ . The Hartree product is symmetric under that swap, not antisymmetric. It cannot describe fermions, full stop. Pauli exclusion (no two electrons in the same state) is one consequence of the antisymmetry requirement; the deeper requirement is the antisymmetry itself.

The Slater determinant

Antisymmetry is what determinants do — swap two columns and a determinant changes sign. So we put the single-particle states into the columns of a matrix and take the determinant:

ψ (r_{1}, r_{2}) = \frac{1}{2} det (χ_{1} (x_{1}) χ_{1} (x_{2}) χ_{2} (x_{1}) χ_{2} (x_{2}))

The $χ_{i}$ are spin orbitals — products of a spatial orbital and a spin function ( $α$ or $β$ ), and $x_{i}$ stands for the combined spatial-and-spin coordinates of electron $i$ . Expand the determinant:

ψ = \frac{1}{2} [χ_{1} (x_{1}) χ_{2} (x_{2}) - χ_{2} (x_{1}) χ_{1} (x_{2})] .

Swap $x_{1} \leftrightarrow x_{2}$ and the expression flips sign. The Pauli requirement is built in. Generalizing to $N$ electrons gives an $N \times N$ determinant of spin orbitals — the construction is the same.

Closed-shell with one spatial orbital

The simplest non-trivial case is two electrons in the same spatial orbital with opposite spins — the closed-shell singlet. Both electrons occupy the same spatial function $ϕ$ , so the spin orbitals are $χ_{1} = ϕ α$ and $χ_{2} = ϕβ$ . The Slater determinant becomes

ψ = ϕ (r_{1}) ϕ (r_{2}) \cdot \frac{1}{2} [α (1) β (2) - β (1) α (2)],

a symmetric spatial part times an antisymmetric spin part — the classic singlet structure. Apply the Slater-Condon rules (or expand $⟨ ψ ∣ \hat{H} ∣ ψ ⟩$ by hand and use spin orthogonality $⟨ α ∣ β ⟩ = 0$ ) and the energy collapses to two terms:

E = 2 ⟨ ϕ ∣ \hat{h} ∣ ϕ ⟩ + (ϕϕ ∣ ϕϕ)

The first term is twice the one-electron energy of the spatial orbital — both electrons sit in $ϕ$ , so each contributes the same $h$ -integral. The second term is the only two-electron integral that survives. Notation: $\hat{h} = - \frac{1}{2} \nabla^{2} - Z / r$ is the one-electron Hamiltonian (kinetic plus nuclear attraction); the chemist-notation Coulomb integral is

(ϕϕ ∣ ϕϕ) = \iint ∣ ϕ (r_{1}) ∣^{2} \frac{1}{r _{12}} ∣ ϕ (r_{2}) ∣^{2} d r_{1} d r_{2} .

The exchange integral $(ϕ α ϕβ ∣ ϕβ ϕ α)$ that one might expect from antisymmetry vanishes here: it carries a factor of $⟨ α ∣ β ⟩$ , and that's zero. Exchange will reappear with a vengeance once we have two spatial orbitals with the same spin (page 7), but for closed-shell same-orbital it's silent.

The two-electron integral over Gaussians

For our STO-3G basis the contracted Coulomb integral is a sum over primitive integrals:

(ϕϕ ∣ ϕϕ) = ijk l \sum d_{i} d_{j} d_{k} d_{l} (g_{i} g_{j} ∣ g_{k} g_{l}) .

For four 1s primitive Gaussians at the same centre, the primitive integral has a closed form:

(g_{α} g_{β} ∣ g_{γ} g_{δ}) = \frac{2 ( 16 α β γ δ ) ^{3/4}}{π ( α + β + γ + δ ) ( α + β ) ( γ + δ )} .

This is the same-centre, all-s-type specialization of the general Boys-function expression that handles arbitrary centres and angular momenta. We'll need the more general version for H₂; for the He calculation today, this is enough.

The code (new this page)

One new integral function (coulomb_g), the helium-specific STO-3G constants, and a function that assembles the closed-shell energy $E = 2 h + J$ from the matrix machinery we already have.

# ---- two-electron integral over s-Gaussians at the same centre ----

def coulomb_g(a, b, c, d):
    """(g_a g_b | g_c g_d) for normalized same-centre 1s primitives.

    Closed form (Szabo & Ostlund A.41 specialized to a single centre and
    s-type angular momentum, with all four exponents non-negative):

        (ab|cd) = 2 (16 a b c d)^(3/4) / [sqrt(pi (a+b+c+d)) * (a+b)(c+d)]
    """
    num = 2 * (16 * a * b * c * d) ** 0.75
    den = np.sqrt(np.pi * (a + b + c + d)) * (a + b) * (c + d)
    return num / den


# ---- STO-3G basis for helium ----
# Same contraction coefficients as hydrogen; exponents scaled to the
# effective zeta of helium (~1.69, by Slater's rules).
STO_3G_HE_ALPHAS = np.array([6.36242139, 1.15892300, 0.31364979])
STO_3G_HE_COEFFS = np.array([0.15432897, 0.53532814, 0.44463454])


def helium_energy_sto3g():
    """HF energy of the He atom: closed-shell singlet, one spatial orbital.

    With both electrons in the same spatial orbital and opposite spins,
    the Slater-determinant energy reduces to
        E = 2 <phi|h|phi> + (phi phi | phi phi)
    where the exchange contribution vanishes because <alpha|beta> = 0.
    """
    a = STO_3G_HE_ALPHAS
    d = STO_3G_HE_COEFFS
    n = len(a)
    Z = 2  # helium nucleus

    # one-electron primitive matrices
    S = np.array([[overlap_g(a[i], a[j])               for j in range(n)] for i in range(n)])
    T = np.array([[kinetic_g(a[i], a[j])               for j in range(n)] for i in range(n)])
    V = np.array([[nuclear_attraction_g(a[i], a[j], Z) for j in range(n)] for i in range(n)])

    # contracted basis function: phi = sum_i d_i g_i.
    # h matrix element <phi|h|phi> = d.T (T+V) d.
    # the contracted phi may not be exactly unit-normalized; divide through.
    norm_sq = d @ S @ d
    h_phi   = (d @ (T + V) @ d) / norm_sq

    # two-electron integral (phi phi | phi phi) = sum_{ijkl} d_i d_j d_k d_l (ij|kl)
    J_phi = 0.0
    for i in range(n):
        for j in range(n):
            for k in range(n):
                for l in range(n):
                    J_phi += d[i] * d[j] * d[k] * d[l] * coulomb_g(a[i], a[j], a[k], a[l])
    J_phi /= norm_sq ** 2

    return 2 * h_phi + J_phi

The four-deep loop in helium_energy_sto3g is the structure that becomes $O (N^{4})$ when the basis grows beyond a single contracted function. For three primitives this is 81 integrals — fine. For a real basis (hundreds of functions), the same loop becomes the dominant cost of any HF calculation, and is why two-electron-integral evaluation is what every quantum chemistry code spends most of its time on.

Run it

Hydrogen atom:
  Slater zeta=1   : E = -0.500000  (exact)
  STO-3G basis    : E = -0.466582

Helium atom (closed-shell HF, single spatial orbital):
  STO-3G basis    : E = -2.807784
  exact           : E = -2.903724
  HF correlation gap: +0.095940

Two new lines worth examining. The He STO-3G energy comes out at $- 2.808$ Hartree; the experimental ground-state energy of helium is $- 2.904$ . The 96 mHa gap between them is the correlation energy — the bit of the true wavefunction that a single Slater determinant cannot capture, no matter how good the basis. HF describes each electron as moving in the average field of the other; correlation is the part of the motion where the electrons actually dodge each other instantaneously. We won't fix that in this series. Recovering correlation is the entire subject of post-HF methods.

Note that the variational principle still holds: every method we build will give an energy strictly above the exact value. That's not a bug; it's the floor that tells us when we've stopped improving. The next page derives the Hartree-Fock equations themselves — we've used the closed-shell formula $E = 2 h + J$ for a single fixed orbital here, but for H₂ the orbital itself becomes something we have to find, and the equations that determine it are the Hartree-Fock equations.

The full file

"""qc.py — Computational QC from scratch.
Step 4: helium. Slater determinants, antisymmetry, and the first
two-electron integral.
"""

import numpy as np
from scipy.optimize import minimize_scalar


# ---- analytic Slater 1s integrals (kept from steps 1-2 as a reference) ----

def kinetic_1s(zeta):
    return 0.5 * zeta ** 2


def nuclear_attraction_1s(zeta):
    return -zeta


def hydrogen_energy(zeta):
    return kinetic_1s(zeta) + nuclear_attraction_1s(zeta)


def optimise_hydrogen():
    result = minimize_scalar(hydrogen_energy, bracket=(0.5, 1.5))
    return result.x, result.fun


# ---- primitive 1s Gaussian integrals (same-centre, nucleus at origin) ----

def overlap_g(a, b):
    return (4 * a * b / (a + b) ** 2) ** 0.75


def kinetic_g(a, b):
    return 3 * a * b / (a + b) * overlap_g(a, b)


def nuclear_attraction_g(a, b, Z=1):
    return -Z * 2 * np.sqrt((a + b) / np.pi) * overlap_g(a, b)


def coulomb_g(a, b, c, d):
    """(g_a g_b | g_c g_d) for normalized same-centre 1s primitives."""
    num = 2 * (16 * a * b * c * d) ** 0.75
    den = np.sqrt(np.pi * (a + b + c + d)) * (a + b) * (c + d)
    return num / den


# ---- STO-3G ----

STO_3G_H_ALPHAS  = np.array([3.42525091, 0.62391373, 0.16885540])
STO_3G_HE_ALPHAS = np.array([6.36242139, 1.15892300, 0.31364979])
STO_3G_COEFFS    = np.array([0.15432897, 0.53532814, 0.44463454])


def hydrogen_energy_sto3g():
    a = STO_3G_H_ALPHAS
    d = STO_3G_COEFFS
    n = len(a)
    S = np.array([[overlap_g(a[i], a[j])            for j in range(n)] for i in range(n)])
    T = np.array([[kinetic_g(a[i], a[j])            for j in range(n)] for i in range(n)])
    V = np.array([[nuclear_attraction_g(a[i], a[j]) for j in range(n)] for i in range(n)])
    H = T + V
    return (d @ H @ d) / (d @ S @ d)


def helium_energy_sto3g():
    """HF energy of the He atom: closed-shell singlet, one spatial orbital."""
    a = STO_3G_HE_ALPHAS
    d = STO_3G_COEFFS
    n = len(a)
    Z = 2

    S = np.array([[overlap_g(a[i], a[j])               for j in range(n)] for i in range(n)])
    T = np.array([[kinetic_g(a[i], a[j])               for j in range(n)] for i in range(n)])
    V = np.array([[nuclear_attraction_g(a[i], a[j], Z) for j in range(n)] for i in range(n)])

    norm_sq = d @ S @ d
    h_phi   = (d @ (T + V) @ d) / norm_sq

    J_phi = 0.0
    for i in range(n):
        for j in range(n):
            for k in range(n):
                for l in range(n):
                    J_phi += d[i] * d[j] * d[k] * d[l] * coulomb_g(a[i], a[j], a[k], a[l])
    J_phi /= norm_sq ** 2

    return 2 * h_phi + J_phi


if __name__ == "__main__":
    E_h_slater = optimise_hydrogen()[1]
    E_h_sto3g  = hydrogen_energy_sto3g()
    E_he_sto3g = helium_energy_sto3g()

    print("Hydrogen atom:")
    print(f"  Slater zeta=1   : E = {E_h_slater:+.6f}  (exact)")
    print(f"  STO-3G basis    : E = {E_h_sto3g:+.6f}")
    print()
    print("Helium atom (closed-shell HF, single spatial orbital):")
    print(f"  STO-3G basis    : E = {E_he_sto3g:+.6f}")
    print(f"  exact           : E = -2.903724")
    print(f"  HF correlation gap: {E_he_sto3g - (-2.903724):+.6f}")

Why does antisymmetry matter if we only have one spatial orbital?

Because the Pauli principle says you can't put two electrons in the same spin orbital. With one spatial orbital you can put one electron in $ϕ α$ and one in $ϕβ$ — different spin orbitals. The antisymmetric Slater determinant is what enforces this; without it, nothing would distinguish the two electrons or prevent them from sharing a single spin orbital. The fact that the energy formula simplifies to $E = 2 h + J$ is downstream of antisymmetry; it's only because we wrote down a Slater determinant that the two-electron part collapses to the single Coulomb integral.

What is the "correlation energy" exactly?

It is the difference between the true ground-state energy and the best possible Hartree-Fock energy in the limit of a complete basis. HF approximates the wavefunction as a single Slater determinant; the true wavefunction is in general a sum of many Slater determinants. The energy difference between the single-determinant and many-determinant pictures is the correlation energy. For He it's about 42 mHa in the basis-complete limit; the 96 mHa we see here also includes basis incompleteness (the STO-3G basis can't represent the exact 1s shape). The two sources of error are conventionally separated as "basis-set error" and "correlation error."

Why does the exchange integral vanish for closed-shell same-orbital?

Mechanically: the exchange integral $(ϕ α ϕβ ∣ ϕβ ϕ α)$ involves $⟨ α ∣ β ⟩$ on each electron coordinate when you separate spatial and spin parts, and $⟨ α ∣ β ⟩ = 0$ . Physically: exchange between electrons of opposite spin doesn't happen because they are distinguishable by spin and the antisymmetrization of the determinant isn't doing any work for opposite-spin pairs. Same-spin pairs are the ones where exchange has bite. With only one spatial orbital and a closed shell, every pair of electrons is opposite-spin, so exchange is silent.