← Part I/Hartree-Fock Method

3. The Hartree-Fock Method

Reading time: ~55 minutes | 5 full derivations | 2 interactive simulations

1. Introduction: The Many-Electron Problem

The Hartree-Fock method is the cornerstone of ab initio quantum chemistry, providing the best single-determinant approximation to the exact many-electron wavefunction.

For an atom or molecule with $N$ electrons and $M$ nuclei, the non-relativistic electronic Hamiltonian (in atomic units, $\hbar = m_e = e = 1$) is:

$$\hat{H} = -\sum_{i=1}^{N}\frac{1}{2}\nabla_i^2 - \sum_{i=1}^{N}\sum_{A=1}^{M}\frac{Z_A}{r_{iA}} + \sum_{i<j}^{N}\frac{1}{r_{ij}}$$

The first term is the electron kinetic energy, the second is the nuclear-electron attraction, and the third is the electron-electron repulsion. It is the final term that makes the Schrödinger equation unsolvable in closed form for $N \geq 2$.

The electron-electron repulsion $\sum_{i<j} 1/r_{ij}$ couples the coordinates of every pair of electrons, preventing separation of variables. For helium alone, we face a six-dimensional partial differential equation. For a molecule like benzene with 42 electrons, exact solution is utterly intractable. The central idea of the Hartree-Fock (HF) method is to replace this intractable many-body problem with an effective one-electron problem: each electron moves in the average field created by all the other electrons.

The Independent Particle Model

As a zeroth approximation, we might try writing the $N$-electron wavefunction as a simple product of one-electron functions (orbitals):

$$\Psi(\mathbf{x}_1, \mathbf{x}_2, \ldots, \mathbf{x}_N) \stackrel{?}{=} \chi_1(\mathbf{x}_1)\chi_2(\mathbf{x}_2)\cdots\chi_N(\mathbf{x}_N)$$

where $\mathbf{x}_i = (\mathbf{r}_i, \sigma_i)$ denotes both spatial and spin coordinates. This is the Hartree product. However, it fails to satisfy the antisymmetry requirement for fermions. A proper fermionic wavefunction must change sign under exchange of any two electrons: $\Psi(\ldots, \mathbf{x}_i, \ldots, \mathbf{x}_j, \ldots) = -\Psi(\ldots, \mathbf{x}_j, \ldots, \mathbf{x}_i, \ldots)$.

2. Derivation 1: The Slater Determinant and Antisymmetry

Constructing an Antisymmetric Wavefunction

Given $N$ orthonormal spin orbitals $\{\chi_1, \chi_2, \ldots, \chi_N\}$, we construct the antisymmetrized product as a determinant. For $N = 2$ electrons, antisymmetrization of $\chi_1(\mathbf{x}_1)\chi_2(\mathbf{x}_2)$ gives:

$$\Psi(\mathbf{x}_1, \mathbf{x}_2) = \frac{1}{\sqrt{2}}\begin{vmatrix} \chi_1(\mathbf{x}_1) & \chi_2(\mathbf{x}_1) \\ \chi_1(\mathbf{x}_2) & \chi_2(\mathbf{x}_2) \end{vmatrix}$$

This is a $2 \times 2$ determinant. Expanding:

$$\Psi(\mathbf{x}_1, \mathbf{x}_2) = \frac{1}{\sqrt{2}}\left[\chi_1(\mathbf{x}_1)\chi_2(\mathbf{x}_2) - \chi_1(\mathbf{x}_2)\chi_2(\mathbf{x}_1)\right]$$

Swapping $\mathbf{x}_1 \leftrightarrow \mathbf{x}_2$ flips the sign — antisymmetry is built in.

General N-Electron Slater Determinant

For $N$ electrons in $N$ spin orbitals, the Slater determinant is:

$$\Psi(\mathbf{x}_1, \ldots, \mathbf{x}_N) = \frac{1}{\sqrt{N!}}\begin{vmatrix} \chi_1(\mathbf{x}_1) & \chi_2(\mathbf{x}_1) & \cdots & \chi_N(\mathbf{x}_1) \\ \chi_1(\mathbf{x}_2) & \chi_2(\mathbf{x}_2) & \cdots & \chi_N(\mathbf{x}_2) \\ \vdots & \vdots & \ddots & \vdots \\ \chi_1(\mathbf{x}_N) & \chi_2(\mathbf{x}_N) & \cdots & \chi_N(\mathbf{x}_N) \end{vmatrix}$$

In compact notation: $|\Psi\rangle = |\chi_1 \chi_2 \cdots \chi_N\rangle$.

Derivation of the Normalization Factor

The determinant expands as a sum over all $N!$ permutations $\hat{P}$of the symmetric group $S_N$:

$$\Psi = \frac{1}{\sqrt{N!}}\sum_{P}(-1)^p \hat{P}\left[\chi_1(\mathbf{x}_1)\chi_2(\mathbf{x}_2)\cdots\chi_N(\mathbf{x}_N)\right]$$

where $(-1)^p$ is the parity of permutation $P$. To verify normalization, compute $\langle\Psi|\Psi\rangle$:

$$\langle\Psi|\Psi\rangle = \frac{1}{N!}\sum_{P}\sum_{Q}(-1)^{p+q}\langle \hat{P}[\chi_1\cdots\chi_N] | \hat{Q}[\chi_1\cdots\chi_N]\rangle$$

Since the spin orbitals are orthonormal, $\langle\chi_i|\chi_j\rangle = \delta_{ij}$, the only non-vanishing terms occur when $P = Q$, giving:

$$\langle\Psi|\Psi\rangle = \frac{1}{N!}\sum_{P}(-1)^{2p} \cdot 1 = \frac{1}{N!}\cdot N! = 1$$

Since $(-1)^{2p} = 1$ always. This confirms the $1/\sqrt{N!}$ prefactor ensures proper normalization.

Pauli Exclusion Principle

If two electrons occupy the same spin orbital, say $\chi_i = \chi_j$ with$i \neq j$, then the determinant has two identical columns and vanishes identically:

$$\chi_i = \chi_j \implies \Psi = 0$$

This is precisely the Pauli exclusion principle: no two electrons can occupy the same quantum state. The antisymmetry of the Slater determinant automatically encodes this fundamental requirement.

3. Derivation 2: Hartree-Fock Equations from the Variational Principle

Energy of a Slater Determinant

Using Slater-Condon rules, the expectation value of the Hamiltonian with a single Slater determinant $|\Psi_0\rangle$ is:

$$E_0 = \langle\Psi_0|\hat{H}|\Psi_0\rangle = \sum_{i=1}^{N}\langle i|\hat{h}|i\rangle + \frac{1}{2}\sum_{i=1}^{N}\sum_{j=1}^{N}\left[\langle ij|ij\rangle - \langle ij|ji\rangle\right]$$

where the one-electron integrals are:

$$\langle i|\hat{h}|i\rangle = \int \chi_i^*(\mathbf{x})\left[-\frac{1}{2}\nabla^2 - \sum_A\frac{Z_A}{r_A}\right]\chi_i(\mathbf{x})\,d\mathbf{x}$$

and the two-electron integrals in physicist's notation are:

$$\langle ij|ij\rangle = \iint \frac{|\chi_i(\mathbf{x}_1)|^2|\chi_j(\mathbf{x}_2)|^2}{r_{12}}\,d\mathbf{x}_1\,d\mathbf{x}_2$$

$$\langle ij|ji\rangle = \iint \frac{\chi_i^*(\mathbf{x}_1)\chi_j(\mathbf{x}_1)\chi_j^*(\mathbf{x}_2)\chi_i(\mathbf{x}_2)}{r_{12}}\,d\mathbf{x}_1\,d\mathbf{x}_2$$

The first is the Coulomb integral and the second is the exchange integral. Note that the exchange integral has no classical analogue — it arises purely from antisymmetry.

Variational Minimization with Constraints

We seek the spin orbitals $\{\chi_i\}$ that minimize $E_0$ subject to orthonormality constraints $\langle\chi_i|\chi_j\rangle = \delta_{ij}$. Using Lagrange multipliers $\varepsilon_{ij}$, we define the functional:

$$\mathcal{L}[\{\chi_i\}] = E_0[\{\chi_i\}] - \sum_{i=1}^{N}\sum_{j=1}^{N}\varepsilon_{ji}\left(\langle\chi_i|\chi_j\rangle - \delta_{ij}\right)$$

Setting the functional derivative $\delta\mathcal{L}/\delta\chi_i^* = 0$ requires:

$$\hat{h}(\mathbf{x}_1)\chi_i(\mathbf{x}_1) + \sum_{j=1}^{N}\left[\hat{J}_j(\mathbf{x}_1) - \hat{K}_j(\mathbf{x}_1)\right]\chi_i(\mathbf{x}_1) = \sum_{j=1}^{N}\varepsilon_{ji}\chi_j(\mathbf{x}_1)$$

A unitary transformation among the occupied orbitals diagonalizes the Lagrange multiplier matrix, giving the canonical Hartree-Fock equations.

The Coulomb and Exchange Operators

The Coulomb operator $\hat{J}_j$ represents the classical electrostatic repulsion from the charge density of electron $j$:

$$\hat{J}_j(\mathbf{x}_1)\chi_i(\mathbf{x}_1) = \left[\int \frac{|\chi_j(\mathbf{x}_2)|^2}{r_{12}}\,d\mathbf{x}_2\right]\chi_i(\mathbf{x}_1)$$

The exchange operator $\hat{K}_j$ has no classical analogue and arises from the antisymmetry requirement:

$$\hat{K}_j(\mathbf{x}_1)\chi_i(\mathbf{x}_1) = \left[\int \frac{\chi_j^*(\mathbf{x}_2)\chi_i(\mathbf{x}_2)}{r_{12}}\,d\mathbf{x}_2\right]\chi_j(\mathbf{x}_1)$$

Note how $\hat{K}_j$ exchanges the orbital labels: $\chi_i$ enters under the integral but $\chi_j$ appears outside. This operator is non-local — it depends on $\chi_i$ everywhere in space, not just at the point $\mathbf{x}_1$.

The Fock Operator and Canonical HF Equations

Defining the Fock operator:

$$\hat{f}(\mathbf{x}_1) = \hat{h}(\mathbf{x}_1) + \sum_{j=1}^{N}\left[\hat{J}_j(\mathbf{x}_1) - \hat{K}_j(\mathbf{x}_1)\right]$$

For a closed-shell system with $N/2$ doubly occupied spatial orbitals, summing over spatial orbitals with the factor of 2 for each occupation gives:

$$\hat{f} = \hat{h} + \sum_{j=1}^{N/2}\left(2\hat{J}_j - \hat{K}_j\right)$$

The canonical Hartree-Fock equations are the pseudo-eigenvalue problem:

$$\boxed{\hat{f}|\chi_i\rangle = \varepsilon_i|\chi_i\rangle}$$

These are pseudo-eigenvalue equations because $\hat{f}$ itself depends on its eigenfunctions through $\hat{J}_j$ and $\hat{K}_j$. The solutions must therefore be obtained iteratively — this is the self-consistent field (SCF) procedure.

Roothaan-Hall Equations

In practice, we expand each molecular orbital (MO) as a linear combination of $K$known basis functions $\{\phi_\mu\}$:

$$\psi_i(\mathbf{r}) = \sum_{\mu=1}^{K} C_{\mu i}\,\phi_\mu(\mathbf{r})$$

Substituting into the HF equation $\hat{f}\psi_i = \varepsilon_i\psi_i$, multiplying on the left by $\phi_\nu^*$, and integrating:

$$\sum_{\mu} F_{\nu\mu}C_{\mu i} = \varepsilon_i\sum_{\mu} S_{\nu\mu}C_{\mu i}$$

where $F_{\nu\mu} = \langle\phi_\nu|\hat{f}|\phi_\mu\rangle$ is the Fock matrix and$S_{\nu\mu} = \langle\phi_\nu|\phi_\mu\rangle$ is the overlap matrix. In matrix form:

$$\boxed{\mathbf{FC} = \mathbf{SC}\boldsymbol{\varepsilon}}$$

This is the Roothaan-Hall equation (Roothaan, 1951; Hall, 1951). It converts the integro-differential HF equation into a matrix eigenvalue problem. The Fock matrix elements are:

$$F_{\mu\nu} = H_{\mu\nu}^{\text{core}} + \sum_{\lambda\sigma}P_{\lambda\sigma}\left[(\mu\nu|\lambda\sigma) - \frac{1}{2}(\mu\lambda|\nu\sigma)\right]$$

where $P_{\lambda\sigma} = 2\sum_{i}^{N/2} C_{\lambda i}C_{\sigma i}$ is the density matrix and $(\mu\nu|\lambda\sigma)$ are two-electron integrals in chemist's notation.

4. Derivation 3: The Self-Consistent Field (SCF) Procedure

Why Iteration Is Necessary

The Fock matrix $\mathbf{F}$ depends on the density matrix $\mathbf{P}$, which in turn depends on the MO coefficients $\mathbf{C}$, which are the eigenvectors of $\mathbf{F}$. This circular dependency means we must solve the problem self-consistently: guess an initial $\mathbf{P}$, build $\mathbf{F}$, solve for new $\mathbf{C}$ and $\mathbf{P}$, and repeat until convergence.

Step-by-Step SCF Algorithm

Specify the molecule and basis set. Choose nuclear coordinates $\{R_A, Z_A\}$ and a basis set $\{\phi_\mu\}$ of$K$ functions.
Compute all one- and two-electron integrals. Calculate$S_{\mu\nu}$, $H_{\mu\nu}^{\text{core}}$, and all $K^4/8$unique two-electron integrals $(\mu\nu|\lambda\sigma)$.
Diagonalize the overlap matrix. Compute$\mathbf{S}^{-1/2}$ for the orthogonalization transformation. Using the eigendecomposition$\mathbf{S} = \mathbf{U}\mathbf{s}\mathbf{U}^\dagger$:
$$\mathbf{S}^{-1/2} = \mathbf{U}\mathbf{s}^{-1/2}\mathbf{U}^\dagger$$
Initial guess for density matrix. Common choices: diagonalize $\mathbf{H}^{\text{core}}$ (core Hamiltonian guess), use extended Hückel theory, or use a superposition of atomic densities (SAD).
Build the Fock matrix. From the current density matrix $\mathbf{P}$:
$$F_{\mu\nu} = H_{\mu\nu}^{\text{core}} + \sum_{\lambda\sigma}P_{\lambda\sigma}\left[(\mu\nu|\lambda\sigma) - \tfrac{1}{2}(\mu\lambda|\nu\sigma)\right]$$
Transform to orthogonal basis. Compute$\mathbf{F}' = \mathbf{S}^{-1/2\dagger}\mathbf{F}\mathbf{S}^{-1/2}$.
Diagonalize. Solve$\mathbf{F}'\mathbf{C}' = \mathbf{C}'\boldsymbol{\varepsilon}$.
Back-transform. Obtain$\mathbf{C} = \mathbf{S}^{-1/2}\mathbf{C}'$.
Form new density matrix. Using the $N/2$lowest eigenvectors:
$$P_{\mu\nu} = 2\sum_{i=1}^{N/2}C_{\mu i}C_{\nu i}$$
Check convergence. If$|\Delta E| < \epsilon_E$ and $||\mathbf{P}^{(n)} - \mathbf{P}^{(n-1)}|| < \epsilon_P$, the SCF has converged. Otherwise, return to step 5.

Convergence Criteria and Acceleration

Typical convergence thresholds are $\epsilon_E \sim 10^{-8}$ Hartree for energy and $\epsilon_P \sim 10^{-6}$ for the density matrix RMS change. The total electronic energy at convergence is:

$$E_{\text{elec}} = \frac{1}{2}\sum_{\mu\nu}P_{\nu\mu}\left(H_{\mu\nu}^{\text{core}} + F_{\mu\nu}\right)$$

Plain SCF iteration often oscillates or diverges. Common acceleration techniques include:

DIIS (Direct Inversion in the Iterative Subspace): Pulay's method (1980) extrapolates the Fock matrix from previous iterations to minimize the error vector $\mathbf{e} = \mathbf{FPS} - \mathbf{SPF}$
Damping: Mix new and old density matrices: $\mathbf{P}^{(n+1)} = \alpha\mathbf{P}^{\text{new}} + (1-\alpha)\mathbf{P}^{(n)}$
Level shifting: Add a constant shift to virtual orbital energies to increase the HOMO-LUMO gap

5. Derivation 4: Koopmans's Theorem

Statement and Proof

Koopmans's theorem (1934): Within the frozen orbital approximation, the ionization energy for removing an electron from occupied orbital $\chi_k$ equals the negative of its orbital energy:

$$\boxed{IP_k = E(N-1, k) - E(N) = -\varepsilon_k}$$

Proof: The $N$-electron HF energy can be written:

$$E(N) = \sum_{i=1}^{N}\langle i|\hat{h}|i\rangle + \frac{1}{2}\sum_{i=1}^{N}\sum_{j=1}^{N}\left(\langle ij||ij\rangle\right)$$

where $\langle ij||ij\rangle = \langle ij|ij\rangle - \langle ij|ji\rangle$ is the antisymmetrized two-electron integral. The orbital energy is:

$$\varepsilon_k = \langle k|\hat{h}|k\rangle + \sum_{j=1}^{N}\langle kj||kj\rangle$$

Now remove electron $k$ while keeping all other orbitals unchanged (frozen orbital approximation). The $(N-1)$-electron energy is:

$$E(N-1, k) = \sum_{i \neq k}\langle i|\hat{h}|i\rangle + \frac{1}{2}\sum_{\substack{i \neq k \\ j \neq k}}\langle ij||ij\rangle$$

Taking the difference:

$$E(N) - E(N-1, k) = \langle k|\hat{h}|k\rangle + \frac{1}{2}\sum_{j \neq k}\langle kj||kj\rangle + \frac{1}{2}\sum_{i \neq k}\langle ik||ik\rangle$$

Since $\langle kj||kj\rangle = \langle jk||jk\rangle$ (by symmetry of the antisymmetrized integrals), the two sums are identical:

$$E(N) - E(N-1, k) = \langle k|\hat{h}|k\rangle + \sum_{j \neq k}\langle kj||kj\rangle = \varepsilon_k$$

Since $\varepsilon_k < 0$ for bound orbitals,$IP_k = E(N-1,k) - E(N) = -\varepsilon_k > 0$. QED.

Limitations of Koopmans's Theorem

Orbital relaxation: In reality, when an electron is removed, the remaining$N-1$ electrons relax to lower the energy. This relaxation energy is neglected in the frozen orbital approximation, causing Koopmans's IP to overestimate the true IP.
Electron correlation: The correlation energy changes between the $N$- and$(N-1)$-electron systems are ignored. Correlation effects typically reduce the IP.
Fortuitous cancellation: For many molecules, the overestimation from neglecting relaxation is partially cancelled by the underestimation from neglecting differential correlation, making Koopmans's values surprisingly accurate (typically within 1–2 eV of experiment).
Electron affinities: Koopmans's theorem applied to virtual orbital energies gives very poor estimates of electron affinities because virtual orbitals are optimized for the$N$-electron system, not the $(N+1)$-electron system.

6. Derivation 5: Electron Correlation and Post-Hartree-Fock Methods

Definition of Correlation Energy

Löwdin (1959) defined the correlation energy as the difference between the exact non-relativistic energy and the Hartree-Fock limit:

$$\boxed{E_{\text{corr}} = E_{\text{exact}} - E_{\text{HF}}}$$

By the variational principle, $E_{\text{HF}} \geq E_{\text{exact}}$, so$E_{\text{corr}} \leq 0$ always. The correlation energy typically accounts for about 1% of the total electronic energy but is crucial for chemical accuracy ($\sim 1$ kcal/mol $\approx 1.6$ mHartree).

The HF method recovers ~99% of the total energy, but the missing 1% can represent the entire bond dissociation energy. For $\text{N}_2$, for example,$E_{\text{corr}} \approx -0.55$ Ha while the bond energy is only 0.36 Ha.

Configuration Interaction (CI)

The most conceptually straightforward post-HF method: expand the wavefunction in a basis of Slater determinants formed by exciting electrons from occupied to virtual orbitals:

$$|\Psi_{\text{CI}}\rangle = c_0|\Psi_0\rangle + \sum_{ia}c_i^a|\Psi_i^a\rangle + \sum_{ijab}c_{ij}^{ab}|\Psi_{ij}^{ab}\rangle + \cdots$$

Full CI (all possible excitations) gives the exact answer within the chosen basis set, but scales factorially with system size. Truncated CI (CISD = singles + doubles) is practical but not size-consistent: the energy of two non-interacting fragments$A + B$ does not equal $E(A) + E(B)$.

Møller-Plesset Perturbation Theory (MP2)

Treat electron correlation as a perturbation to the Fock operator. The zeroth-order Hamiltonian is $\hat{H}_0 = \sum_i \hat{f}(i)$ and the perturbation is $\hat{V} = \hat{H} - \hat{H}_0$. The second-order energy correction is:

$$E_{\text{MP2}} = \sum_{i<j}^{\text{occ}}\sum_{a<b}^{\text{virt}} \frac{|\langle ij||ab\rangle|^2}{\varepsilon_i + \varepsilon_j - \varepsilon_a - \varepsilon_b}$$

MP2 scales as $\mathcal{O}(N^5)$ and recovers 80–90% of the correlation energy for many systems. It is size-consistent but not variational.

Coupled Cluster Theory

The coupled cluster (CC) ansatz uses an exponential parameterization:

$$|\Psi_{\text{CC}}\rangle = e^{\hat{T}}|\Psi_0\rangle, \qquad \hat{T} = \hat{T}_1 + \hat{T}_2 + \cdots$$

where $\hat{T}_n$ generates all $n$-fold excitations. The exponential form ensures size-consistency at every truncation level. CCSD(T) — coupled cluster with singles, doubles, and perturbative triples — is often called the “gold standard” of quantum chemistry:

CCSD: $\mathcal{O}(N^6)$ — includes $\hat{T}_1 + \hat{T}_2$
CCSD(T): $\mathcal{O}(N^7)$ — adds perturbative triples correction
Chemical accuracy ($< 1$ kcal/mol) for most main-group chemistry

7. Applications: Basis Sets and Computational Considerations

Slater-Type vs. Gaussian-Type Orbitals

Slater-type orbitals (STOs) have the correct functional form:

$$\phi_{\text{STO}}(\mathbf{r}) = N r^{n-1} e^{-\zeta r} Y_l^m(\theta, \phi)$$

They have the correct cusp at the nucleus ($r = 0$) and correct exponential decay at large $r$, but multi-center two-electron integrals cannot be computed analytically.

Gaussian-type orbitals (GTOs) use a Gaussian radial part:

$$\phi_{\text{GTO}}(\mathbf{r}) = N x^{l_x} y^{l_y} z^{l_z} e^{-\alpha r^2}$$

The key advantage: the product of two Gaussians centered at different points is itself a Gaussian centered at a third point (Gaussian product theorem). This makes all two-electron integrals analytically tractable. Boys (1950) introduced GTOs for this reason. The trade-off is that more GTOs are needed to approximate an STO — hence the STO-$n$G contracted basis sets that fit $n$ Gaussians to each STO.

Common Basis Set Families

Basis Set	Type	Functions per atom	Use case
STO-3G	Minimal	1 per occupied AO	Qualitative, large systems
6-31G*	Split-valence + polarization	~15 (2nd row)	Routine geometry optimization
cc-pVDZ	Correlation-consistent DZ	~14 (2nd row)	Correlated calculations
cc-pVTZ	Correlation-consistent TZ	~30 (2nd row)	High accuracy
aug-cc-pVQZ	Augmented QZ	~80 (2nd row)	Near basis set limit

Computational Scaling

The computational cost of HF and post-HF methods scales as a power of the number of basis functions $K$:

Two-electron integrals: $\mathcal{O}(K^4)$ — the bottleneck for HF
HF (with screening): effectively $\mathcal{O}(K^{2-3})$ for large molecules
MP2: $\mathcal{O}(K^5)$
CCSD: $\mathcal{O}(K^6)$
CCSD(T): $\mathcal{O}(K^7)$
Full CI: $\mathcal{O}(K!)$ — exact but exponentially expensive

Modern HF implementations exploit integral screening, density fitting (RI approximation), and linear-scaling techniques to handle systems with thousands of basis functions.

Molecular Properties from HF

The converged HF wavefunction provides access to numerous molecular properties:

Total energy and equilibrium geometry (via gradient optimization)
Dipole moment: $\boldsymbol{\mu} = -\text{Tr}(\mathbf{P}\mathbf{D}) + \sum_A Z_A \mathbf{R}_A$ where $\mathbf{D}$ is the dipole integral matrix
Mulliken charges: $q_A = Z_A - \sum_{\mu \in A}(\mathbf{PS})_{\mu\mu}$
Harmonic frequencies from the Hessian (second derivatives of energy)
Ionization energies via Koopmans's theorem

8. Historical Context

1928 — Douglas Hartree

Proposed the self-consistent field method for atoms, treating each electron as moving in the average field of all others. His original formulation used simple product wavefunctions (no antisymmetry), and he performed calculations by hand with the help of his father, who was a skilled numerical analyst.

1930 — Vladimir Fock

Generalized Hartree's method to include antisymmetry (exchange), deriving the full Hartree-Fock equations using the variational principle with Slater determinants. This introduced the exchange operator $\hat{K}$ that has no classical analogue.

1930 — John C. Slater

Independently introduced the determinantal wavefunction (now called the Slater determinant) and developed Slater-type orbitals (STOs) with the correct radial behavior for atoms. His rules for estimating effective nuclear charges remain widely used.

1951 — Clemens C. J. Roothaan & George G. Hall

Independently formulated the matrix form of the HF equations ($\mathbf{FC} = \mathbf{SC}\boldsymbol{\varepsilon}$), transforming the integro-differential equations into an algebraic eigenvalue problem suitable for digital computers. This was the breakthrough that made HF calculations practical.

1950 — S. Francis Boys

Introduced Gaussian-type orbitals (GTOs) as basis functions, exploiting the Gaussian product theorem to make all molecular integrals analytically computable. This insight underlies every modern quantum chemistry program (Gaussian, ORCA, Q-Chem, etc.).

9. Python Simulation: SCF for Helium

The following code implements a complete restricted Hartree-Fock SCF calculation for the helium atom using a two-function STO-3G basis. It computes all one- and two-electron integrals analytically for s-type Gaussians centered at the origin, solves the Roothaan-Hall equations iteratively, and plots the energy convergence:

Hartree-Fock SCF for Helium Atom (STO-3G Basis)

Python

Complete RHF/STO-3G calculation for He. Computes overlap, kinetic, nuclear attraction, and two-electron integrals over contracted Gaussian basis functions. Iterates the Roothaan-Hall equations to self-consistency and plots energy convergence.

script.py198 lines

import numpy as np

# =============================================================
# Minimal Hartree-Fock SCF for the Helium atom
# using a 2-function STO-3G basis set (1s orbitals)
# =============================================================
# Each STO is approximated by 3 Gaussians (STO-3G contraction).
# We use TWO contracted 1s functions with different exponents
# to form a minimal double-zeta-like basis for He.

# STO-3G contraction coefficients and exponents for He 1s (zeta=1.6875)
# Standard STO-3G parameters for zeta=1.0
alpha_sto3g = np.array([0.109818, 0.405771, 2.22766])
coeff_sto3g = np.array([0.444635, 0.535328, 0.154329])

zeta1 = 1.6875   # optimized for He
zeta2 = 1.0      # second basis function (different exponent)

# Scale exponents for each basis function
alpha1 = alpha_sto3g * zeta1**2
alpha2 = alpha_sto3g * zeta2**2

# Lists of (exponents, coefficients) per basis function
basis = [
    (alpha1, coeff_sto3g),
    (alpha2, coeff_sto3g),
]
nbasis = len(basis)
Z = 2  # nuclear charge for helium

# --- Integral routines for Gaussian primitives ---

def overlap_prim(a, b):
    """Overlap integral between two s-type Gaussians."""
    return (np.pi / (a + b))**1.5

def kinetic_prim(a, b):
    """Kinetic energy integral between two s-type Gaussians."""
    return a * b / (a + b) * (3.0 * np.pi / (a + b)) * (np.pi / (a + b))**0.5

def nuclear_prim(a, b, Z_nuc):
    """Nuclear attraction integral for s-type Gaussians at origin."""
    return -Z_nuc * 2.0 * np.pi / (a + b)

def eri_prim(a, b, c, d):
    """Two-electron repulsion integral (ab|cd) for s-type Gaussians at origin."""
    return 2.0 * np.pi**2.5 / ((a + b) * (c + d) * np.sqrt(a + b + c + d))

# --- Build contracted integrals ---

def contracted_integral(bra, ket, prim_func, *args):
    """Compute contracted integral from primitive function."""
    alphas_i, coeffs_i = bra
    alphas_j, coeffs_j = ket
    val = 0.0
    for p in range(len(alphas_i)):
        for q in range(len(alphas_j)):
            val += coeffs_i[p] * coeffs_j[q] * prim_func(alphas_i[p], alphas_j[q], *args)
    return val

def contracted_eri(bas_a, bas_b, bas_c, bas_d):
    """Four-center two-electron integral over contracted functions."""
    aa, ca = bas_a
    ab, cb = bas_b
    ac, cc = bas_c
    ad, cd = bas_d
    val = 0.0
    for p in range(len(aa)):
        for q in range(len(ab)):
            for r in range(len(ac)):
                for s in range(len(ad)):
                    val += (ca[p] * cb[q] * cc[r] * cd[s]
                            * eri_prim(aa[p], ab[q], ac[r], ad[s]))
    return val

# Build one-electron matrices
S = np.zeros((nbasis, nbasis))
T = np.zeros((nbasis, nbasis))
V = np.zeros((nbasis, nbasis))

for i in range(nbasis):
    for j in range(nbasis):
        S[i, j] = contracted_integral(basis[i], basis[j], overlap_prim)
        T[i, j] = contracted_integral(basis[i], basis[j], kinetic_prim)
        V[i, j] = contracted_integral(basis[i], basis[j], nuclear_prim, Z)

H_core = T + V

# Build two-electron integral tensor
ERI = np.zeros((nbasis, nbasis, nbasis, nbasis))
for i in range(nbasis):
    for j in range(nbasis):
        for k in range(nbasis):
            for l in range(nbasis):
                ERI[i, j, k, l] = contracted_eri(basis[i], basis[j], basis[k], basis[l])

# --- SCF Procedure ---
max_iter = 50
convergence = 1e-10

# Solve generalized eigenvalue problem: H_core C = S C e (initial guess)
from numpy.linalg import eigh, inv

S_inv_half = np.zeros_like(S)
eigvals_s, eigvecs_s = eigh(S)
S_inv_half = eigvecs_s @ np.diag(1.0 / np.sqrt(eigvals_s)) @ eigvecs_s.T

# Initial guess: diagonalize core Hamiltonian
F_prime = S_inv_half.T @ H_core @ S_inv_half
eigvals, eigvecs_prime = eigh(F_prime)
C = S_inv_half @ eigvecs_prime

# Form initial density matrix (2 electrons in lowest orbital)
P = np.zeros((nbasis, nbasis))
P = 2.0 * np.outer(C[:, 0], C[:, 0])

energies = []
print("SCF Iteration    Electronic Energy (Ha)    Total Energy (Ha)      Delta E")
print("=" * 78)

E_old = 0.0
for iteration in range(1, max_iter + 1):
    # Build Fock matrix: F = H_core + G
    G = np.zeros((nbasis, nbasis))
    for i in range(nbasis):
        for j in range(nbasis):
            for k in range(nbasis):
                for l in range(nbasis):
                    # Coulomb - 0.5 * Exchange
                    G[i, j] += P[k, l] * (ERI[i, j, k, l] - 0.5 * ERI[i, l, k, j])
    F = H_core + G

# Electronic energy
    E_elec = 0.5 * np.sum(P * (H_core + F))
    E_total = E_elec + Z * Z / 1.0  # no internuclear repulsion for atom, but E_nuc = 0
    E_total = E_elec  # single atom, no nuclear repulsion

delta_E = abs(E_total - E_old)
    energies.append(E_total)
    print(f"    {iteration:3d}             {E_elec:16.10f}         {E_total:16.10f}       {delta_E:.2e}")

if delta_E < convergence and iteration > 1:
        print(f"\nSCF converged after {iteration} iterations!")
        break

E_old = E_total

# Solve Roothaan-Hall: F C = S C e
    F_prime = S_inv_half.T @ F @ S_inv_half
    eigvals_f, eigvecs_prime = eigh(F_prime)
    C = S_inv_half @ eigvecs_prime
    P = 2.0 * np.outer(C[:, 0], C[:, 0])

print(f"\nFinal HF Energy: {energies[-1]:.10f} Ha")
print(f"Exact He energy: -2.8616800 Ha")
print(f"HF limit:        -2.8616800 Ha")
print(f"Orbital energies: {eigvals_f}")

# Plot convergence
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 5.5))
fig.patch.set_facecolor('#0a0a0a')
ax1.set_facecolor('#0a0a0a')
ax2.set_facecolor('#0a0a0a')

iters = np.arange(1, len(energies) + 1)

# Panel 1: Energy convergence
ax1.plot(iters, energies, 'o-', color='#34d399', linewidth=2, markersize=6, label='SCF Energy')
ax1.axhline(y=energies[-1], color='#f97316', linestyle='--', alpha=0.7, label=f'Converged: {energies[-1]:.6f} Ha')
ax1.set_xlabel('SCF Iteration', color='white', fontsize=13)
ax1.set_ylabel('Total Energy (Hartree)', color='white', fontsize=13)
ax1.set_title('HF/STO-3G SCF Convergence for Helium', color='#34d399', fontsize=14, fontweight='bold')
ax1.legend(fontsize=10, facecolor='#1a1a2e', edgecolor='#34d399', labelcolor='white')
ax1.tick_params(colors='white')
for spine in ax1.spines.values():
    spine.set_color('#334155')

# Panel 2: Energy difference (log scale)
if len(energies) > 1:
    deltas = [abs(energies[i] - energies[i-1]) for i in range(1, len(energies))]
    ax2.semilogy(np.arange(2, len(energies) + 1), deltas, 's-', color='#a78bfa', linewidth=2, markersize=6)
    ax2.axhline(y=convergence, color='#ef4444', linestyle=':', alpha=0.7, label=f'Threshold: {convergence:.0e}')
    ax2.set_xlabel('SCF Iteration', color='white', fontsize=13)
    ax2.set_ylabel(r'$|\Delta E|$ (Hartree)', color='white', fontsize=13)
    ax2.set_title('Energy Convergence Rate', color='#a78bfa', fontsize=14, fontweight='bold')
    ax2.legend(fontsize=10, facecolor='#1a1a2e', edgecolor='#a78bfa', labelcolor='white')
    ax2.tick_params(colors='white')
    for spine in ax2.spines.values():
        spine.set_color('#334155')

plt.tight_layout()
plt.savefig('output.png', dpi=150, bbox_inches='tight', facecolor='#0a0a0a')
plt.close()
print("\nSCF convergence plot saved.")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

Key observations: The SCF converges rapidly (typically 5–10 iterations) to the HF/STO-3G energy for helium. The converged energy is higher than the exact value ($-2.8617$ Ha) because the minimal basis set is too inflexible. A larger basis would approach the Hartree-Fock limit. The difference between the HF limit and the exact energy is the correlation energy ($\approx -0.042$ Ha for He).

10. Fortran Simulation: HF Orbital Energies for H–Ne

This Fortran program uses tabulated Hartree-Fock orbital energies (from Clementi & Roetti, 1974) to demonstrate Koopmans's theorem across the first row of the periodic table. It compares the negative HOMO orbital energy with the experimental first ionization energy:

HF Orbital Energies and Koopmans's Theorem (H-Ne)

Fortran

Computes orbital energy sums and tests Koopmans's theorem (IP = -epsilon_HOMO) for atoms H through Ne using tabulated Clementi-Roetti HF orbital energies. Compares with experimental ionization energies.

hf_atomic_energies.f90134 lines

program hf_atomic_energies
  implicit none
  ! ============================================================
  ! Compute Hartree-Fock total energies for atoms H through Ne
  ! using tabulated HF orbital energies (in Hartree).
  ! Demonstrates Koopman's theorem: -epsilon_i ~ IP
  ! ============================================================

integer, parameter :: dp = selected_real_kind(15, 307)
  integer, parameter :: max_atoms = 10
  integer, parameter :: max_orbitals = 5

character(len=2) :: symbols(max_atoms)
  integer :: Z_vals(max_atoms)
  integer :: n_orb(max_atoms)                     ! number of occupied subshells
  integer :: occ(max_atoms, max_orbitals)          ! occupation of each subshell
  real(dp) :: eps(max_atoms, max_orbitals)         ! orbital energy (Hartree)
  character(len=4) :: orb_labels(max_atoms, max_orbitals)

real(dp) :: E_total, E_orb_sum, IP_koopman, E_hf_total
  real(dp) :: exp_IP(max_atoms)  ! experimental first IP in Hartree
  integer :: i, j

! --- Tabulated Hartree-Fock orbital energies (Clementi & Roetti, 1974) ---

symbols = (/ 'H ', 'He', 'Li', 'Be', 'B ', 'C ', 'N ', 'O ', 'F ', 'Ne' /)
  Z_vals  = (/ 1, 2, 3, 4, 5, 6, 7, 8, 9, 10 /)

! Experimental first ionization energies (Hartree)
  exp_IP = (/ 0.4998_dp, 0.9036_dp, 0.1981_dp, 0.3426_dp, 0.3049_dp, &
              0.4138_dp, 0.5341_dp, 0.5005_dp, 0.6404_dp, 0.7925_dp /)

! H: 1s
  n_orb(1) = 1
  occ(1,1) = 1;    eps(1,1) = -0.5000_dp;   orb_labels(1,1) = '1s  '

! He: 1s
  n_orb(2) = 1
  occ(2,1) = 2;    eps(2,1) = -0.9179_dp;   orb_labels(2,1) = '1s  '

! Li: 1s, 2s
  n_orb(3) = 2
  occ(3,1) = 2;    eps(3,1) = -2.4777_dp;   orb_labels(3,1) = '1s  '
  occ(3,2) = 1;    eps(3,2) = -0.1963_dp;   orb_labels(3,2) = '2s  '

! Be: 1s, 2s
  n_orb(4) = 2
  occ(4,1) = 2;    eps(4,1) = -4.7327_dp;   orb_labels(4,1) = '1s  '
  occ(4,2) = 2;    eps(4,2) = -0.3093_dp;   orb_labels(4,2) = '2s  '

! B: 1s, 2s, 2p
  n_orb(5) = 3
  occ(5,1) = 2;    eps(5,1) = -7.6953_dp;   orb_labels(5,1) = '1s  '
  occ(5,2) = 2;    eps(5,2) = -0.4947_dp;   orb_labels(5,2) = '2s  '
  occ(5,3) = 1;    eps(5,3) = -0.3099_dp;   orb_labels(5,3) = '2p  '

! C: 1s, 2s, 2p
  n_orb(6) = 3
  occ(6,1) = 2;    eps(6,1) = -11.3255_dp;  orb_labels(6,1) = '1s  '
  occ(6,2) = 2;    eps(6,2) = -0.7056_dp;   orb_labels(6,2) = '2s  '
  occ(6,3) = 2;    eps(6,3) = -0.4333_dp;   orb_labels(6,3) = '2p  '

! N: 1s, 2s, 2p
  n_orb(7) = 3
  occ(7,1) = 2;    eps(7,1) = -15.6291_dp;  orb_labels(7,1) = '1s  '
  occ(7,2) = 2;    eps(7,2) = -0.9452_dp;   orb_labels(7,2) = '2s  '
  occ(7,3) = 3;    eps(7,3) = -0.5677_dp;   orb_labels(7,3) = '2p  '

! O: 1s, 2s, 2p
  n_orb(8) = 3
  occ(8,1) = 2;    eps(8,1) = -20.6686_dp;  orb_labels(8,1) = '1s  '
  occ(8,2) = 2;    eps(8,2) = -1.2443_dp;   orb_labels(8,2) = '2s  '
  occ(8,3) = 4;    eps(8,3) = -0.6319_dp;   orb_labels(8,3) = '2p  '

! F: 1s, 2s, 2p
  n_orb(9) = 3
  occ(9,1) = 2;    eps(9,1) = -26.3829_dp;  orb_labels(9,1) = '1s  '
  occ(9,2) = 2;    eps(9,2) = -1.5724_dp;   orb_labels(9,2) = '2s  '
  occ(9,3) = 5;    eps(9,3) = -0.7300_dp;   orb_labels(9,3) = '2p  '

! Ne: 1s, 2s, 2p
  n_orb(10) = 3
  occ(10,1) = 2;   eps(10,1) = -32.7725_dp; orb_labels(10,1) = '1s  '
  occ(10,2) = 2;   eps(10,2) = -1.9306_dp;  orb_labels(10,2) = '2s  '
  occ(10,3) = 6;   eps(10,3) = -0.8504_dp;  orb_labels(10,3) = '2p  '

! Tabulated total HF energies (Clementi & Roetti)
  ! We compute sum of orbital energies and compare

write(*,'(A)') '=================================================================='
  write(*,'(A)') ' Hartree-Fock Orbital Energies and Koopman''s Theorem'
  write(*,'(A)') '=================================================================='
  write(*,'(A)') ''
  write(*,'(A,A4,A6,A14,A16,A14,A14)') ' ', 'Atom', 'Z', 'HOMO (Ha)', &
       '-HOMO (Ha)', 'Exp IP (Ha)', 'Error (%)'
  write(*,'(A)') '------------------------------------------------------------------'

do i = 1, max_atoms
     ! HOMO energy is the last occupied orbital
     IP_koopman = -eps(i, n_orb(i))

write(*,'(A,A4,I5,F14.6,F16.6,F14.6,F12.2)') &
          ' ', symbols(i), Z_vals(i), eps(i, n_orb(i)), &
          IP_koopman, exp_IP(i), &
          100.0_dp * abs(IP_koopman - exp_IP(i)) / exp_IP(i)
  end do

write(*,'(A)') ''
  write(*,'(A)') '=================================================================='
  write(*,'(A)') ' Orbital Energy Decomposition'
  write(*,'(A)') '=================================================================='

do i = 1, max_atoms
     write(*,'(A)') ''
     write(*,'(A,A2,A,I2,A)') ' --- ', symbols(i), ' (Z=', Z_vals(i), ') ---'

E_orb_sum = 0.0_dp
     do j = 1, n_orb(i)
        write(*,'(A,A4,A,I1,A,F12.6,A)') '   ', orb_labels(i,j), &
             ' : occ=', occ(i,j), '  eps=', eps(i,j), ' Ha'
        E_orb_sum = E_orb_sum + occ(i,j) * eps(i,j)
     end do
     write(*,'(A,F14.6,A)') '   Sum(n_i * eps_i) = ', E_orb_sum, ' Ha'
     write(*,'(A)') '   (Note: E_HF /= Sum(n_i*eps_i) due to double-counting of e-e)'
  end do

write(*,'(A)') ''
  write(*,'(A)') '=================================================================='
  write(*,'(A)') ' Key: Koopman''s theorem states -eps_HOMO ~ IP'
  write(*,'(A)') ' Agreement is typically within 5-15% for light atoms.'
  write(*,'(A)') ' Deviations arise from orbital relaxation (frozen orbital approx).'
  write(*,'(A)') '=================================================================='

end program hf_atomic_energies

Click Run to execute the Fortran code

Code will be compiled with gfortran and executed on the server

Key observations: Koopmans's theorem predictions agree with experimental ionization energies to within 5–15% for light atoms. The systematic overestimation reflects the neglect of orbital relaxation, while the occasional underestimation for atoms like O and F reflects the role of differential correlation. The fortuitous partial cancellation of these two errors is a well-known feature of Koopmans's theorem.

Chapter Summary

The Hartree-Fock method provides the optimal single-determinant description of many-electron systems. Starting from the antisymmetric Slater determinant, application of the variational principle yields the Fock operator and the HF equations. The Roothaan-Hall matrix formulation converts these into an algebraic eigenvalue problem that is solved iteratively via the SCF procedure.

Key Results

The Slater determinant automatically satisfies the Pauli exclusion principle
The Fock operator: $\hat{f} = \hat{h} + \sum_j(2\hat{J}_j - \hat{K}_j)$
Roothaan-Hall equations: $\mathbf{FC} = \mathbf{SC}\boldsymbol{\varepsilon}$
Koopmans's theorem: $IP_k = -\varepsilon_k$ (frozen orbital approximation)
Correlation energy: $E_{\text{corr}} = E_{\text{exact}} - E_{\text{HF}} < 0$
Post-HF hierarchy: CI, MP2, CCSD, CCSD(T) for systematic improvement

← 2. Molecular Orbital Theory 4. Variational Method →