The Hartree-Fock equations

On the previous page the helium energy formula $E=2h_{\varphi \varphi }+(\varphi \varphi \mid \varphi \varphi )$ was derived for a fixed orbital $\varphi$, which was in turn the STO-3G contraction we wrote down by hand. That's not yet a calculation — it's an evaluation. To do quantum chemistry we have to let the orbital itself become the unknown and ask: of all possible $\varphi$, which one minimizes the energy? The Hartree-Fock equations are the answer to that question.

The variational problem on the determinant

Take the energy of an $N$-electron Slater determinant as a functional of the spin orbitals $\{{\chi }_i\}$:

The first sum is the one-electron part (kinetic + nuclear attraction); the second is the two-electron Coulomb minus exchange. We want to minimize $E$ with respect to the ${\chi }_i$, with the constraint that they remain orthonormal: $⟨{\chi }_i\mid {\chi }_j⟩={\delta }_{ij}$. Standard Lagrangian move:

Vary $\mathcal{L}$ with respect to each ${\chi }_i^{\ast }$ and set the variation to zero. The matrix ${\epsilon }_{ij}$ is Hermitian, so we can rotate the spin orbitals among themselves to diagonalize it. The variational equations then collapse to a one-electron eigenvalue problem:

The Fock operator $\hat{F}$ is the one-electron Hamiltonian plus an averaged interaction with all the other occupied spin orbitals:

where the Coulomb operator ${\hat{J}}_j$ acts on a state $\chi$ as $({\hat{J}}_j\chi )({\mathbf{x}}_1)=[\int |{\chi }_j({\mathbf{x}}_2){|}^2/r_{12}d{\mathbf{x}}_2]\chi ({\mathbf{x}}_1)$ — the classical Coulomb potential from the charge density of orbital ${\chi }_j$ — and the exchange operator ${\hat{K}}_j$ is the same thing with the orbitals swapped, a purely quantum object with no classical analogue.

The thing to notice: $\hat{F}$ depends on the ${\chi }_j$ we are trying to find. We can't just set up the operator and diagonalize once; the operator is a functional of its own eigenvectors. This is the self-consistency problem, the "self-consistent" in self-consistent field. Pages 6 and 7 are about how to solve it.

Closed-shell simplification

For a closed-shell ground state — every occupied spatial orbital has both spin partners present — the spin sums collapse. Each spatial orbital ${\varphi }_a$ contributes twice as a Coulomb source (once via ${\varphi }_a\alpha$, once via ${\varphi }_a\beta$) but only once as an exchange source against any given target spin (because the opposite-spin partner has zero spin overlap, so its $\hat{K}$ contribution vanishes). The result is the closed-shell Fock operator:

For the closed-shell helium of page 4 there is exactly one occupied spatial orbital $\varphi$, so $\hat{F}=\hat{h}+2{\hat{J}}_{\varphi }-{\hat{K}}_{\varphi }$. The matrix element $⟨\varphi \mid \hat{F}\mid \varphi ⟩$ works out to $h_{\varphi \varphi }+(\varphi \varphi \mid \varphi \varphi )$ (the same-orbital exchange integral cancels half of the same-orbital Coulomb), which is the orbital energy $\epsilon$. We can read off the relation between total and orbital energies:

The total energy is not twice the orbital energy. If you just doubled $\epsilon$ you would count the electron-electron repulsion twice — once when each electron sits in the averaged field of the other, double-counting the single repulsion event. Subtracting $(\varphi \varphi \mid \varphi \varphi )$ fixes the over-count. This double-counting structure is going to appear in every HF energy expression we write from now on.

The Fock matrix in a basis

Pick a finite basis of single-particle functions $\{g_i\}$ (for us, the three primitive Gaussians) and write the orbital as $\varphi ={\sum }_i{c}_i{g}_i$. The Fock matrix in that basis is

where $(ij\mid kl)=(g_i{g}_j\mid g_k{g}_l)$ is a primitive electron-repulsion integral and the density matrix in the closed-shell case is

The factor of 2 absorbs the spin pair. The total energy in this matrix form is the symmetric average $E=\frac{1}{2}{\sum }_{ij}{P}_{ij}(h_{ij}+F_{ij})$. Once the four-index tensor $(ij\mid kl)$ is in hand, the rest is just contractions over indices.

The code (new this page)

Three new helpers — the four-index electron-repulsion tensor, the closed-shell Fock-matrix builder, and the energy-from-density expression — plus a sanity check that pushes the helium calculation of page 4 through the new machinery and recovers the same number.

# ---- four-index electron-repulsion tensor over primitive Gaussians ----

def build_eri_tensor(alphas):
    """ERI[i,j,k,l] = (g_i g_j | g_k g_l) for same-centre 1s primitives."""
    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_g(alphas[i], alphas[j], alphas[k], alphas[l])
    return ERI


# ---- closed-shell Fock matrix ----

def build_fock(h, ERI, P):
    """F_ij = h_ij + sum_{kl} P_kl [(ij|kl) - 1/2 (ik|jl)].

    Closed-shell: density matrix P is built from the doubly-occupied
    orbitals as P_kl = 2 sum_a c_ka c_la (the factor of 2 is the spin pair).
    """
    J = np.einsum('ijkl,kl->ij', ERI, P)        # Coulomb part
    K = np.einsum('ikjl,kl->ij', ERI, P)        # exchange part
    return h + J - 0.5 * K


def hf_energy_from_density(P, h, F):
    """Closed-shell HF energy:  E = 1/2 sum_ij P_ij (h_ij + F_ij)."""
    return 0.5 * np.sum(P * (h + F))


# ---- helium re-done through the Fock-matrix machinery ----

def helium_via_fock():
    """Verify: building F from the STO-3G density gives the same energy."""
    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)])
    h   = T + V
    ERI = build_eri_tensor(a)

    # normalize the STO-3G contraction so that c.T S c = 1
    c = d / np.sqrt(d @ S @ d)
    P = 2.0 * np.outer(c, c)        # closed-shell density: one occupied orbital

    F = build_fock(h, ERI, P)
    E   = hf_energy_from_density(P, h, F)
    eps = c @ F @ c                  # orbital energy of the (fixed) HOMO

    return E, eps

A note on the einsum calls. 'ijkl,kl->ij' contracts the last two indices of the four-index tensor against the density matrix to make the Coulomb-type matrix $J_{ij}={\sum }_{kl}{P}_{kl}(ij\mid kl)$; the 'ikjl,kl->ij' string permutes the index order of the tensor on the fly to give the exchange-type $K_{ij}={\sum }_{kl}{P}_{kl}(ik\mid jl)$. Index bookkeeping is going to be most of the work in the next several pages, and einsum is what keeps it readable.

Run it

Helium via the Fock-matrix machinery (STO-3G contraction held fixed):
  total energy E   = -2.807784   (matches step 4)
  orbital energy e = -0.876036
  check:  2 e - J  = -2.807784    ✓

Two things to confirm. (1) The total energy is bit-for-bit the same as page 4's $-2.807784$; we've changed the bookkeeping, not the physics. (2) The double-counting relation $E=2\epsilon -J$ closes. The orbital energy $\epsilon =-0.876036$ is the energy it costs to remove the electron — Koopmans' approximation to helium's first ionisation potential, which is experimentally about $0.904$ Hartree (about 24.6 eV). HF gets it within a few percent.

What we still don't have: a way to actually find the optimal $c$. We held it fixed at the STO-3G fit and asked the Fock machinery to evaluate the energy. The next page introduces the Roothaan-Hall equations — the matrix form of the HF eigenvalue problem — and turns the search for the optimal orbital into a generalized eigenvalue solve.