Hartree–Fock theory#

Hartree–Fock equation#

In the Hartree–Fock approximation, the many-electron wave function takes the form of a Slater determinant

\[\begin{split} | \Psi \rangle = \frac{1}{\sqrt{N!}} \begin{vmatrix} \psi_{1}(\mathbf{r}_1) & \cdots & \psi_{N}(\mathbf{r}_1) \\ \vdots & \ddots & \vdots \\ \psi_{1}(\mathbf{r}_N) & \cdots & \psi_{N}(\mathbf{r}_N) \\ \end{vmatrix} \end{split}\]

where \(\psi_i\) are the single-electron wave functions known as spin orbitals. The Hartree–Fock energy and the associated state is found by minimizing the energy functional

\[ E_\mathrm{HF} = \min_{\psi} E[\psi] \]

under the constraint that the spin orbitals remain orthonormal. Here, \(\psi\) collectively refers to the entire set of \(N\) spin orbitals. Such a constrained minimization is conveniently performed by means of the technique of Lagrange multipliers.


In Hartree–Fock theory, we introduce the real-valued Lagrangian

\[ L[\psi] = E[\psi] - \sum_{i,j=1}^N \varepsilon_{ji} \big( \langle \psi_i | \psi_j \rangle - \delta_{ij} \big) \]

and search for the set of spin orbitals, \(\psi\), that results in a first variation that vanishes

\[ \delta L = 0 \]

Expressing the energy as the expectation value of the electronic Hamiltonian with respect to a Slater determinant and using the general expressions for matrix elements, we arrive at

\[\begin{align*} \delta L & = \sum_{i=1}^N \langle \delta \psi_i | \hat{h} | \psi_i \rangle + \sum_{i,j=1}^N \big( \langle \delta \psi_i \psi_j | \hat{g} | \psi_i \psi_j\rangle - \langle \delta \psi_i \psi_j | \hat{g} | \psi_j \psi_i\rangle - \varepsilon_{ji} \langle \delta \psi_i | \psi_j \rangle \big) + \mbox{complex conjugate} \\ &= \sum_{i=1}^N \langle \delta \psi_i | \big( \hat{f} | \psi_i \rangle - \sum_{j=1}^N \varepsilon_{ji} | \psi_j \rangle \big) + \mbox{complex conjugate} \end{align*}\]

where we have introduced the one-electron Fock operator

\[ \hat{f} = \hat{h} + \sum_{j=1}^N \big( \hat{J}_j - \hat{K}_j \big) \]


\[\begin{align*} \hat{J}_j | \psi_i \rangle & = \Big[ \int \frac{e^2 |\psi_j(\mathbf{r}')|^2}{4\pi\varepsilon_0 |\mathbf{r} - \mathbf{r}'|} d^3\mathbf{r}' \Big] | \psi_i \rangle \\ % \hat{K}_j | \psi_i \rangle & = \Big[ \int \frac{e^2 \psi_j^\dagger(\mathbf{r}')\psi_i(\mathbf{r}')}{4\pi\varepsilon_0 |\mathbf{r} - \mathbf{r}'|} d^3\mathbf{r}' \Big] | \psi_j \rangle \end{align*}\]

Since the first-order variation in the Lagrangian is required to vanish for general variations in the spin orbitals, we have shown that the Hartree–Fock solution is given by

\[ \hat{f} | \psi_i \rangle - \sum_{j=1}^N \varepsilon_{ji} | \psi_j \rangle = 0 \]

This equation is known as the Hartree–Fock equation and it to be solved for the spin orbitals and the associated Lagrange multipliers. We note that the matrix elements of the Fock operator equal the multipliers

\[ f_{ki} = \langle \psi_k | \hat{f} | \psi_i \rangle = \sum_{j=1}^N \varepsilon_{ji} \langle \psi_k | \psi_j \rangle = \varepsilon_{ki} \]

Canonical form#

Apart from a trivial overall phase factor, unitary transformations among the occupied orbitals are shown to leave the Hartree–Fock wave function unchanged. We introduce a unitary transformation that diagonalizes the Hermitian Fock matrix

\[ \mathbf{f}' = \langle \overline{\psi}' | \hat{f} | \overline{\psi}' \rangle = \mathbf{U}^\dagger \langle \overline{\psi} | \hat{f} | \overline{\psi} \rangle \mathbf{U} = \mathbf{U}^\dagger \mathbf{f} \, \mathbf{U} \]

We have here adopted the compact overline notation of orbitals. In this basis of canonical spin orbitals, the Hartree–Fock equation takes the form

\[ \hat{f} | \psi_i \rangle = \varepsilon_{i} | \psi_i \rangle \]

which we recognize as an eigenvalue equation introducing the orbital energies, \(\varepsilon_{i}\), as the eigenvalues of the Fock operator. With an infinite number of solutions to the Hartree–Fock equation, the Hartree–Fock ground state is given by employing the \(N\) spin orbitals with lowest orbital energies in the Slater determinant.

In AO basis#

The spatial parts of the spin orbitals, or molecular orbitals (MOs), are expanded as linear combination of atomic orbitals (LCAO). In the basis of spin atomic orbitals, the Fock matrix becomes block diagonal

\[\begin{split} \mathbf{F} = \begin{pmatrix} \mathbf{F}^{\alpha\alpha} & \mathbf{0} \\ \mathbf{0} & \mathbf{F}^{\beta\beta} \end{pmatrix} \end{split}\]


Here, we adopt the convention of using \(\mathbf{F}\) for the Fock matrix in AO basis, as compared to above use of \(\mathbf{f}\) in the MO basis.

Using the bar notation to distinguish \(\alpha\)- and \(\beta\)-spin atomic orbitals, we get

\[\begin{align*} F_{\mu\nu} & = F^{\alpha\alpha}_{\mu\nu} = h_{\mu\nu} + \sum_{\gamma\delta} \Big( D_{\gamma\delta}(\mu\nu|\gamma\delta) - D^\alpha_{\gamma\delta}(\mu\delta|\gamma\nu) \Big) \\ F_{\bar{\mu}\bar{\nu}} & = F^{\beta\beta}_{\mu\nu} = h_{\mu\nu} + \sum_{\gamma\delta} \Big( D_{\gamma\delta}(\mu\nu|\gamma\delta) - D^\beta_{\gamma\delta}(\mu\delta|\gamma\nu) \Big) \\ F_{\mu\bar{\nu}} & = F_{\bar{\mu}\nu} = 0 \end{align*}\]


\[\begin{align*} D_{\gamma\delta} &= D^\alpha_{\gamma\delta} + D^\beta_{\gamma\delta} \\ D^\alpha_{\gamma\delta}& = \sum_{j=1}^{N_\alpha} \big[c_{\gamma j}^\alpha\big]^* c_{\delta j}^\alpha ; \quad D^\beta_{\gamma\delta} = \sum_{j=1}^{N_\beta} \big[c_{\gamma j}^\beta\big]^* c_{\delta j}^\beta \\ \end{align*}\]

The canonical Hartree–Fock equation thereby takes the form

\[ \mathbf{F C} = \mathbf{S C} \boldsymbol{\varepsilon} \, \]

where \(\mathbf{S}\) is the overlap matrix and \(\boldsymbol{\varepsilon}\) is a diagonal matrix collecting the orbital energies.

Hartree–Fock energy#

For a given density \(\mathbf{D}\), the Hartree–Fock energy becomes equal to

\[ E_\mathrm{HF} = \frac{1}{2} \mathrm{tr} \big[ (\mathbf{h} + \mathbf{F}) \mathbf{D} \big] + V^\mathrm{n-n} \, \]

where \(V^\mathrm{n-n}\) is the nuclear repulsion energy.

Koopmans theorem#

The orbital energies of occupied and unoccupied orbitals, respectively, equal

\[\begin{align*} \varepsilon_i & = \langle \psi_i |\hat{f} | \psi_i \rangle = \langle \psi_i |\hat{h} | \psi_i \rangle + \sum_{j\neq i}^N \big( \langle \psi_i | \hat{J}_j | \psi_i \rangle - \langle \psi_i | \hat{K}_j | \psi_i \rangle \big) \\ \varepsilon_a & = \langle \psi_a |\hat{f} | \psi_a \rangle = \langle \psi_a |\hat{h} | \psi_a \rangle + \sum_{j=1}^N \big( \langle \psi_a | \hat{J}_j | \psi_a \rangle - \langle \psi_a | \hat{K}_j | \psi_a \rangle \big) \end{align*}\]

where the cancellation between Coulomb and exchange terms for \(j=i\) has been used in the former case. It thus appears as if \(\varepsilon_i\) relates to the energy of an electron interacting with \((N-1)\) other electrons, whereas \(\varepsilon_a\) relates to the energy of an electron interacting with \(N\) other electrons. In accordance with these observations, it is readily shown from the expressions for matrix elements that the ionization energy (IE) and electron affinity (EA) become

\[\begin{align*} \mathrm{IE} &= E_i^{N-1} - E_\mathrm{HF}^N = - \varepsilon_i \\ \mathrm{EA} &= E_\mathrm{HF}^N - E_a^{N+1} = - \varepsilon_a \\ \end{align*}\]

where, in the frozen orbital approximation, \(E_i^{N-1}\) is the energy of the system after the removal of the electron in spin orbital \(i\) and \(E_a^{N+1}\) is the energy of the system after the addition of an electron in spin orbital \(a\).

Brillouin theorem#

Based on the expressions for matrix elements, we find

\[ \langle \Psi_\mathrm{HF} | \hat{H} | \Psi_i^a \rangle = \langle \psi_i | \hat{f} | \psi_a \rangle = 0 \]

which shows that there is no coupling between the Hartree–Fock ground state and single excited determinants. This result is known as the Brillouin theorem.

SCF procedure#

Due to the summation over occupied spin orbitals that expresses the effective electron interactions, the Fock operator depends on its eigenfunctions and the canonical Hartree–Fock equation is therefore solved iteratively by means of a self-consistent field (SCF) procedure, such as the Roothaan–Hall approach:


Roothaan–Hall scheme#

In the following, we will consider the spin-restricted formulation where \(\alpha\)- and \(\beta\)-spin orbitals have identical spatial parts. We also restrict the situation to the common case of a closed-shell system such that

\[\begin{align*} N_\alpha & = N_\beta = \frac{1}{2} N \\ D^\alpha_{\gamma\delta} & = D^\beta_{\gamma\delta} = \frac{1}{2} D_{\gamma\delta} = \sum_{j=1}^{N/2} c_{\gamma j}^* c_{\delta j} \end{align*}\]


When referring to closed-shell systems it is customary to refer to the density matrix as that for either of the spin components.