Hartree-Fock Theory and the Variational Principle

The Hartree-Fock method is one of the most fundamental approaches used in electronic structure calculations. Although it is an approximation, it serves as the foundation for more accurate quantum chemical methods. Understanding its principles is essential for exploring both qualitative and advanced quantitative models.

Hartree-Fock theory applies the variational principle to approximate solutions of the electronic Schrödinger equation for molecules. The method reduces the complex many-electron problem into a set of self-consistent one-electron equations.

Step 1: Define the Electronic Hamiltonian

The electronic Hamiltonian (Born-Oppenheimer approximation) is:

$$ \hat{H}_e = \sum_i \left( -\frac{1}{2} \nabla_i^2 - \sum_A \frac{Z_A}{r_{iA}} \right) + \sum_{i<j} \frac{1}{r_{ij}} $$

This includes kinetic energy, electron-nuclear attraction, and electron-electron repulsion.

Step 2: Choose a Trial Wavefunction (Slater Determinant)

To respect the Pauli principle, we use a Slater determinant of spin-orbitals:

$$ \Psi(1, 2, \ldots, N) = \frac{1}{\sqrt{N!}} \begin{vmatrix} \psi_1(1) & \psi_2(1) & \cdots & \psi_N(1) \\ \psi_1(2) & \psi_2(2) & \cdots & \psi_N(2) \\ \vdots & \vdots & \ddots & \vdots \\ \psi_1(N) & \psi_2(N) & \cdots & \psi_N(N) \end{vmatrix} $$

This represents a wavefunction for $N$ electrons using $N$ spin orbitals.

A spin orbital is the product of two functions: one depending on the electron's spatial coordinates $(x, y, z)$ and another depending on the spin coordinate. This is why we call it a spin orbital.

Spin orbitals are used to enforce the Pauli exclusion principle in the wavefunction. Note that the electronic Schrödinger equation does not include any spin terms, so the energy solution depends only on spatial coordinates.

Assumptions of the Slater Determinant

Each electron moves in its own orbital: $ \psi_a(i) $ is a one-electron wavefunction.
All orbitals are orthogonal.
The molecular wavefunction vanishes as $ (x, y, z) \rightarrow \infty $.
Electrons with parallel spins experience correlated motion due to exchange terms in the Slater determinant.

The probability of finding electron 1 at position 1, electron 2 at position 2, ..., and electron $ n $ at position $ n $ is:

$$ P(1,2,\ldots,n) = \psi^*(1,2,\ldots,n)\psi(1,2,\ldots,n) $$

Using a Slater determinant, this probability is not just a product of independent terms. The determinant expansion introduces exchange terms:

$$ \sum_{i=1}^{n} \left| \int \frac{d\mathbf{r}_j \psi_k^*(j) h_{v.e.}(i,j)}{r_{ij}} \psi_k(i) \right| $$

This structure enforces the Pauli exclusion principle: no two electrons with the same spin can occupy the same orbital. If they do, the Slater determinant vanishes. Since the Schrödinger equation itself does not include any explicit term enforcing this principle, the use of the Slater determinant is essential.

Each one-electron orbital is expanded as:

$$ \psi_i = \sum_{\mu=1}^{K} C_{\mu i} \phi_\mu $$

Where the $ \phi_\mu $ are atomic orbitals centered on the nuclei. In the limit $ K \rightarrow \infty $, the expansion becomes exact. Most basis functions are composed of Gaussian primitives:

$$ N x^a y^b z^c e^{-\alpha r^2} $$

The exponents $ a, b, c $ determine the angular momentum ($ L = a + b + c $), $ \zeta $ controls the size of the orbital, and $ N $ is a normalization constant.

Choosing a basis set involves a trade-off: larger basis sets give higher accuracy but require more computational resources. Advances in hardware have made high-accuracy calculations feasible for many molecules.

Spin orbitals can be expressed in different forms depending on the electron configuration:

Roothaan-Hall (Restricted)

$$ \psi_k(x) = \begin{pmatrix} \phi_k(r)\alpha(\omega) \\ \phi_k(r)\beta(\omega) \end{pmatrix} $$

Here two electrons of opposite spin occupy the same spacial orbital. When we use electron orbitals like these we talk about the Rothaan-Hall equations. These equations can only be used to describe closed shell electron configurations because two electron are forced to reside in the same spatial orbital.

Pople-Nesbet (Unrestricted)

$$ \psi_k(x) = \begin{pmatrix} \phi_k^\alpha(r)\alpha_k(\omega) \\ \phi_k^\beta(r)\beta_k(\omega) \end{pmatrix} $$

Here each electron has its own spacial orbital. When we use electron orbitals like these we talk about the Pople-Nestbeth equations. Since each electron has his own spacial orbital these equations can be used to describe open shell electron configurations. The Pople-Nesbeth equations result in one set of equations for the alpha eletrons and a second set of equations for the beta electrons. These equations can be used to represent open shell electron configurations.

Step 3: Variational Principle with a Slater Determinant

When using a Slater determinant as the trial wavefunction, the variational principle can be expressed as follows. The total electronic energy is calculated as the expectation value of the Hamiltonian over the trial wavefunction:

$$ E[\Psi^{\text{SD}}] = \frac{ \int \Psi^{\text{SD}*}(\mathbf{x}_1, \ldots, \mathbf{x}_N)\, \hat{H}\, \Psi^{\text{SD}}(\mathbf{x}_1, \ldots, \mathbf{x}_N) \, d\mathbf{x}_1 \ldots d\mathbf{x}_N }{ \int |\Psi^{\text{SD}}(\mathbf{x}_1, \ldots, \mathbf{x}_N)|^2 \, d\mathbf{x}_1 \ldots d\mathbf{x}_N } $$

Since the Slater determinant $ \Psi^{\text{SD}} $ is normalized, the denominator equals 1, and the expression simplifies to:

$$ E[\Psi^{\text{SD}}] = \langle \Psi^{\text{SD}} | \hat{H} | \Psi^{\text{SD}} \rangle $$

This expectation value gives an upper bound to the exact ground state energy. By optimizing the spin-orbitals that define the Slater determinant, the Hartree-Fock method seeks the minimum of this functional.

HF Spin Treatment and Koopmans’ Caveat

UHF vs. ROHF: UHF allows different $\alpha/\beta$ spatial orbitals but can suffer spin contamination ($\langle \hat S^2\rangle$ deviates from the correct spin value). ROHF enforces exact spin but uses more involved equations.
Koopmans’ theorem (HF): $-\varepsilon_\text{HOMO}$ approximates the first ionization energy only in frozen‑orbital HF and neglects correlation; it is a guideline, not an exact result.

Step 4: Derive the Hartree-Fock Equations

Minimization leads to the Fock equations:

Operator Form

The Hartree-Fock equation in operator form is:

$$ \hat{F}(i) \psi_a(i) = \varepsilon_a \psi_a(i) $$

Where the Fock operator $ \hat{F}(i) $ consists of:

$$ \hat{F}(i) = -\frac{1}{2} \nabla_i^2 - \sum_{A=1}^{K} \frac{Z_A}{r_{iA}} + \sum_{j=1}^{n} \left[ \hat{J}_j(i) - \hat{K}_j(i) \right] $$

$ \nabla_i^2 $: the Laplacian (kinetic energy) of electron $ i $
$ Z_A $: the nuclear charge of nucleus $ A $
$ r_{iA} $: distance between electron $ i $ and nucleus $ A $
$ \hat{J}_j(i) $: Coulomb operator (electron repulsion)
$ \hat{K}_j(i) $: Exchange operator (Pauli antisymmetry)

Integral Form

In integral form, the action of the Fock operator on a spin orbital $ \psi_a(i) $ is:

$$ \hat{F}(i)\psi_a(i) = - \frac{1}{2} \nabla_i^2 \psi_a(i) + \sum_{A=1}^{K} \frac{Z_A}{r_{iA}} \psi_a(i) + \sum_{j=1}^{n} \left[ \left( \int \frac{|\psi_j(j)|^2}{r_{ij}} d\mathbf{r}_j \right)\psi_a(i) - \left( \int \frac{\psi_j^*(j)\psi_a(j)}{r_{ij}} d\mathbf{r}_j \right)\psi_j(i) \right] $$

Definitions

$ A = 1 \ldots K $: nuclei
$ i, j = 1 \ldots n $: electrons
$ a, b, \ldots, k $: spin orbitals
$ r_{ij} $: distance between electrons $ i $ and $ j $
$ r_{iA} $: distance between electron $ i $ and nucleus $ A $

This form of the Hartree-Fock equation is fundamental to computational chemistry and underlies the construction of molecular orbital approximations in quantum chemistry software.

Limitations of the Mean-Field Approximation

In Hartree-Fock theory, we assume that each electron moves in an average or mean field created by all other electrons. While this assumption simplifies the problem significantly, it is not an accurate reflection of how electrons actually behave.

In reality, electrons interact dynamically and instantaneously via the Coulomb force:

$$ V_{ee} = \sum_{i<j} \frac{1}{| \mathbf{r}_i - \mathbf{r}_j|} $$

This means that the position and motion of one electron immediately influences the motion of another. The Hartree-Fock approximation, however, only accounts for this through a time-independent averaged field, missing these real-time correlation effects.

More specifically:

Hartree-Fock includes exchange effects (due to antisymmetry) only for electrons with the same spin.
It completely neglects dynamic correlation between opposite-spin electrons.
The probability of finding two electrons close together is overestimated, especially for opposite spins.

Analogy: Imagine dancers on a crowded floor. The mean-field model assumes each dancer feels an average push from the crowd. But in reality, each dancer reacts dynamically to the exact location and movement of others to avoid collisions. This real-time adjustment is what Hartree-Fock ignores.

Step 5: Solve by Self-Consistent Field (SCF)

In Hartree-Fock theory, the electronic wavefunction is approximated by a single Slater determinant composed of molecular spin-orbitals. These orbitals are determined by minimizing the total energy using the variational principle. However, the Fock operator used to find these orbitals depends on the orbitals themselves, making the problem nonlinear. The solution is an iterative process known as the Self-Consistent Field (SCF) method.

1. Choose a Basis Set

Each orbital $ \psi_i $ is expressed as a linear combination of basis functions $ \chi_\mu $:

$$ \psi_i(\mathbf{r}) = \sum_{\mu} C_{\mu i} \chi_\mu(\mathbf{r}) $$

The coefficients $ C_{\mu i} $ are the parameters that are varied to minimize the total energy.

2. Initial Guess

Start with an initial guess for the orbital coefficients $ C_{\mu i} $. This is used to construct the density matrix:

$$ P_{\mu\nu} = \sum_i^{\text{occ}} C_{\mu i} C_{\nu i}^* $$

3. Build the Fock Matrix

Use the current density matrix to construct the Fock matrix $ F_{\mu\nu} $, which includes one-electron and two-electron contributions:

$$ F_{\mu\nu} = h_{\mu\nu} + \sum_{\lambda\sigma} P_{\lambda\sigma} \left[ \langle \mu\nu | \lambda\sigma \rangle - \frac{1}{2} \langle \mu\lambda | \nu\sigma \rangle \right] $$

4. Solve the Roothaan–Hall Equations

Solve the generalized eigenvalue problem:

$$ \mathbf{F} \mathbf{C} = \mathbf{S} \mathbf{C} \mathbf{\varepsilon} $$

This yields updated orbital coefficients $ C_{\mu i} $ and orbital energies $ \varepsilon_i $.

5. Update and Repeat

Recalculate the density matrix and Fock matrix using the new orbitals. Iterate this procedure until convergence is achieved — typically when the energy or density matrix changes less than a defined threshold.

6. Convergence

The procedure stops when the input and output orbitals are self-consistent, meaning the orbitals are compatible with the Fock operator they generate.

Summary of Varying Parameters

$ C_{\mu i} $: coefficients of molecular orbitals (varied)
$ P_{\mu\nu} $: density matrix (computed)
$ F_{\mu\nu} $: Fock matrix (depends on $ P_{\mu\nu} $)
$ \varepsilon_i $: orbital energies (from diagonalization)

The SCF method transforms a nonlinear many-electron problem into a loop of linear algebra operations, converging toward the optimal mean-field description of the electron system.

Step 6: Compute Properties

After convergence, compute:

Total electronic energy
Electron density and orbital energies
Inputs for post-Hartree-Fock methods