Multi-centre two-electron integrals

Quantum Chemistry

The last piece of the engine. Hartree-Fock needs the four-index electron-repulsion tensor over basis functions, and our basis now lives on two atoms, which means each of the four indices in $(g_{i} g_{j} ∣ g_{k} g_{l})$ can sit at either centre. We have to evaluate every multi-centre case. Same trick as the one-electron integrals: the Gaussian product theorem turns the integrand into a product of two single-centre Gaussians, and Coulomb's law against those Gaussians collapses to a single Boys-function call.

The four-centre formula

Apply the Gaussian-product theorem twice — once to the pair on electron 1, once to the pair on electron 2:

g_{a} (r_{1} - A) g_{b} (r_{1} - B) = K_{A B} g_{p} (r_{1} - P), p = α + β,

g_{c} (r_{2} - C) g_{d} (r_{2} - D) = K_{C D} g_{q} (r_{2} - Q), q = γ + δ .

The four-centre integral becomes a Coulomb integral between two single Gaussians at $P$ and $Q$ , with the two prefactors out front. Carrying out the radial integrals and bookkeeping the normalization gives

(ab ∣ c d) = \frac{2 π ^{5/2}}{p q p + q} K_{A B} K_{C D} F_{0} (\frac{p q}{p + q} ∣ P - Q ∣^{2})

times the four primitive normalization factors $(2 α / π)^{3/4}$ . The same-centre limit ( $A = B = C = D$ ) drops the $K$ 's to $1$ and the Boys argument to $0$ , recovering the same-centre formula from step 4. The general formula handles every geometry with the same single function call.

There are six floating-point exponentials and one Boys evaluation per integral. With six primitives we have $6^{4} = 1296$ calls to make, and the contraction down to AO indices is another $2^{4} \times block size^{4} = 16 \times 81 = 1296$ multiply-and-accumulate operations. For STO-3G H₂ the whole tensor takes a few milliseconds. For larger systems this is the dominant cost — every production quantum-chemistry code has decades of effort poured into making this loop fast (Schwarz screening, density fitting, symmetry, integral direct methods).

Symmetry: 8-fold equivalence

A real ERI computation exploits that $(ab ∣ c d)$ is invariant under several index permutations:

(ab ∣ c d) = (ba ∣ c d) = (ab ∣ d c) = (ba ∣ d c) = (c d ∣ ab) = (d c ∣ ab) = (c d ∣ ba) = (d c ∣ ba) .

Real-orbital ERIs have an 8-fold symmetry that reduces the distinct-integral count by a factor of 8. Our code doesn't exploit it — at this scale the loop is fast enough — but every production code does, and the saving compounds with basis size.

Putting H₂ together

Once the AO-basis Hamiltonian and ERI tensor are in hand, every single thing the SCF cycle needs is available, and the SCF code from page 7 is generic — it doesn't care whether the matrices came from one atom or two. We just feed it the H₂ matrices.

One more thing: the molecular total energy. SCF gives us the electronic energy at the given nuclear positions. The total Born-Oppenheimer energy adds the classical nuclear-nuclear repulsion:

E_{HF} (R) = E_{elec} (R) + A < B \sum \frac{Z _{A} Z _{B}}{∣ R _{A} - R _{B} ∣} .

For H₂ with both $Z = 1$ , that's just $1/ R$ . Without it, the energy diverges to $- \infty$ as $R \to 0$ — the electrons happily collapse onto a single nucleus and lower the energy forever. The nuclear repulsion is what creates a finite equilibrium bond length.

The code (new this page)

The four-centre primitive integral, the routine that builds the full primitive ERI tensor, the AO contraction, and a wrapper that runs SCF on H₂ at a given bond length and adds the nuclear repulsion.

def coulomb_4c(a, A, b, B, c, C, d_, D):
    """(g_a g_b | g_c g_d) — four-centre electron-repulsion integral
    over normalized 1s primitive Gaussians.

        (ab|cd) = int int g_a(r1-A) g_b(r1-B) (1/r12) g_c(r2-C) g_d(r2-D) dr1 dr2

    Closed form via two applications of the Gaussian-product theorem
    plus the Boys function (Szabo & Ostlund A.41):
    """
    p  = a + b
    q  = c + d_
    P  = (a * A + b * B) / p
    Q  = (c * C + d_ * D) / q
    R2_AB = np.sum((A - B) ** 2)
    R2_CD = np.sum((C - D) ** 2)
    R2_PQ = np.sum((P - Q) ** 2)

    K_ab = np.exp(-a * b * R2_AB / p)
    K_cd = np.exp(-c * d_ * R2_CD / q)
    arg  = p * q / (p + q) * R2_PQ

    norm = (2.0 * a / np.pi) ** 0.75 * (2.0 * b / np.pi) ** 0.75 \
         * (2.0 * c / np.pi) ** 0.75 * (2.0 * d_ / np.pi) ** 0.75

    return (norm
            * 2.0 * np.pi ** 2.5 / (p * q * np.sqrt(p + q))
            * K_ab * K_cd
            * boys_F0(arg))


def build_eri_tensor_2c(alphas, centres):
    """Four-index ERI tensor over a multi-centre primitive basis."""
    n = len(alphas)
    ERI = np.zeros((n, n, n, n))
    for i in range(n):
        for j in range(n):
            for k in range(n):
                for l in range(n):
                    ERI[i, j, k, l] = coulomb_4c(
                        alphas[i], centres[i],
                        alphas[j], centres[j],
                        alphas[k], centres[k],
                        alphas[l], centres[l],
                    )
    return ERI


def contract_eri_to_ao(ERI_prim, atom_of, n_ao=2, d=STO_3G_COEFFS):
    """Contract the four-index primitive ERI tensor down to AO indices."""
    ERI_ao = np.zeros((n_ao, n_ao, 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 c in range(n_ao):
                ic = np.where(atom_of == c)[0]
                for D in range(n_ao):
                    iD = np.where(atom_of == D)[0]
                    for ii, i in enumerate(ia):
                        for jj, j in enumerate(ib):
                            for kk, k in enumerate(ic):
                                for ll, l in enumerate(iD):
                                    ERI_ao[a, b, c, D] += (
                                        d[ii] * d[jj] * d[kk] * d[ll] * ERI_prim[i, j, k, l]
                                    )
    return ERI_ao


def h2_eri_tensor(R):
    """AO-basis ERI tensor for H2 at bond length R, with AO normalization."""
    alphas, centres, atom_of = h2_basis(R)
    ERI_prim = build_eri_tensor_2c(alphas, centres)
    ERI_ao = contract_eri_to_ao(ERI_prim, atom_of)

    # AO normalization factor (same N_a we used on page 8)
    S_prim = np.array([[overlap_2c(alphas[i], centres[i], alphas[j], centres[j])
                        for j in range(len(alphas))] for i in range(len(alphas))])
    S_ao = contract_to_ao(S_prim, atom_of)
    norms = np.sqrt(np.diag(S_ao))
    # divide each ERI index by its AO norm
    return ERI_ao / np.einsum('a,b,c,d->abcd', norms, norms, norms, norms)


def h2_scf(R):
    """Run SCF for H2 at bond length R. Returns (E_total, eps, C)."""
    S, T, V = h2_one_electron_matrices(R)
    h   = T + V
    ERI = h2_eri_tensor(R)
    E_elec, eps, C, _ = scf(h, S, ERI, n_occ=1, verbose=False)
    E_nuc = 1.0 * 1.0 / R                       # nuclear-nuclear repulsion
    return E_elec + E_nuc, eps, C

The five-deep nested loop in contract_eri_to_ao is the brute-force version of the AO contraction. A real code uses einsum with a sparse contraction matrix and runs in a fraction of the time, but for two atoms the brute-force version is more readable.

Run it

H2 at R = 1.4 bohr — symmetry-distinct AO ERIs (chemist notation):
  (11|11) = 0.774606    (one-centre on atom A)
  (11|22) = 0.569676    (Coulomb between AOs on different atoms)
  (12|12) = 0.297029    (exchange-type)
  (11|12) = 0.444108    (three-centre)

H2 SCF at R = 1.4 bohr:
  electronic energy : E_elec = -1.831000
  nuclear repulsion : 1/R    = +0.714286
  total HF energy   : E_HF   = -1.116714
  occupied orbital  : eps    = -0.578203

The four printed ERIs are the symmetry-distinct AO integrals for homonuclear H₂. $(11 ∣ 11)$ is the on-atom Coulomb integral — same shape as helium's. The $(11 ∣ 22)$ integral is the Coulomb repulsion between two charge clouds at different atoms; at $R = 1.4$ it's about $0.57$ , smaller than the on-atom value because the charge clouds are further apart on average. The exchange-type $(12 ∣ 12)$ involves products of orbitals straddling the bond, and is smaller again. The three-centre $(11 ∣ 12)$ appears in the Fock matrix off-diagonals.

SCF converges in a handful of iterations. The electronic energy at $R = 1.4$ is $- 1.831$ Hartree; nuclear repulsion is $1/1.4 = 0.714$ ; total HF energy is $- 1.117$ Hartree. This matches Szabo & Ostlund's canonical worked example for H₂ in STO-3G, which is nice corroboration that all six pages of integrals and iteration came out right. The experimental binding energy of H₂ relative to two free H atoms is about $0.174$ Hartree (4.75 eV). We have one more page to go, in which we scan $R$ across a range and watch the binding curve take shape.

The $- 0.5782$ Hartree occupied-orbital energy is what Koopmans' theorem identifies with the H₂ ionisation potential. Experiment puts H₂'s IP at $0.585$ Hartree (15.43 eV); HF gets it within one percent, despite this being one of the smallest possible basis sets. The reason it works so well: H₂ at equilibrium is a closed-shell, weakly correlated, single-bond system — basically the cleanest test case there is for a method built on a single Slater determinant.

How does this scale to larger molecules?

The number of two-electron integrals scales as $N^{4}$ with basis size $N$ , which is what dominates the cost of a real Hartree-Fock calculation. The fact that many of those integrals are tiny (because the basis functions are far apart, so the Gaussian product factor $K_{A B}$ is exponentially small) gets exploited by Schwarz screening: skip integrals whose Schwarz upper bound is below a tolerance. That reduces the effective scaling toward $N^{2}$ for large systems. Density fitting, integral direct algorithms, and the resolution of the identity all push further. None of this changes the underlying formula we just implemented.

Why does the orbital energy approximate the ionisation potential?

Koopmans' theorem: in a frozen-orbital approximation (the remaining N-1 electrons keep the same orbitals after removal), the energy difference between the N-electron and (N-1)-electron states is exactly the orbital energy $- ε_{H OMO}$ . The approximation works well when the cation is well-described by the same orbitals — true for closed-shell stable molecules, less true for systems with significant orbital relaxation.

What's missing from this calculation?

Three things, in roughly increasing severity. (1) Basis-set incompleteness: STO-3G is small. cc-pVTZ would lower the energy by tens of millihartree. (2) Electron correlation: HF misses about 40 mHa of correlation energy in H₂; post-HF methods (CISD, CCSD, CCSD(T)) recover most of it. (3) Static correlation at long bond length: HF dissociates H₂ incorrectly because at large $R$ the ground state is a superposition of two open-shell determinants, and a single-determinant method can't represent that. We'll see this defect concretely in the binding curve on the next page.