The Hartree-Fock Method

Molecular Electronic Structure - Lecture 4

The Hartree-Fock-Roothaan Method: Part I

In the previous lecture, we learned to write a simple approximate form for molecular wavefunctions, as Slater determinants Psi _SD, i.e. as antisymmetrised products of spin-orbitals:

$\begin{array}{cccc} \Psi_{SD}(\mbox{x}_1,\mbox{x}_2,\ldots \mbox{x}_2)= \displaystyle{ \mathcal{A} \{ \prod_{i=1}^{n_{elec}} \chi_i(\mbox{x}_i) \} } \\ \mbox{ } \\ \mbox{where:} \;\; \chi_i(\mbox{x}_i) = \psi_i(r_i)\alpha \; \mbox{or} \; \psi_i(r_i)\beta \\ \mbox{and:} \;\; \displaystyle{ \psi_i(r_i) = \sum_j^{n_{basis}} c_{ij} \phi_{ij}(r_i) } \\ \end{array}$

We now want to determine the value of the coefficients c_ij. To do this, we need to solve the electronic Schrödinger equation:

$\hat {H}_{elec} \Psi (\mbox{x})=E_{elec}\Psi (\mbox{x})$

In fact, because we have decided to use an approximate wavefunction, there is no exact solution to this equation (no values of c_ij which, inserted into the equation, make the two sides rigorously equal). Instead, we are looking for the best approximate solution. How do we define "best" ? The answer to this question is given by the Variational Principle: The best wavefunction is the one with the lowest energy.

The energy of an approximate wavefunction Psi is given by:

$E=\frac{ \int \Psi(\mbox{x}) \hat{H}_{elec} \Psi(\mbox{x}) \> \mbox{dx}} { \int \Psi(\mbox{x}) \Psi(\mbox{x}) \> \mbox{dx} }$

Or, more concisely, and provided the function is normalised:

$E= \langle \Psi | H | \Psi \rangle$

What we are trying to do, then, is to find the coefficients c_ij (which define the molecular orbitals psi in terms of the atomic orbitals phi ) that give the lowest possible energy in the expression above. Doing this is the object of the Hartree-Fock-Roothaan method.

Energy of a Slater Determinant

To calculate the energy of a Slater Determinant wavefunction, we "just" have to insert the mathematical expression for it into the expression above for the energy of a wavefunction.

$E_{SD}=\int \Psi_{SD}\hat{H}_{elec}\Psi_{SD} d\mbox{x}$

Where:

$\Psi(\mbox{x}_1,\mbox{x}_2,\mbox{x}_3,\ldots \mbox{x}_n)= \frac{1}{\sqrt{n_{elec}!}} \left | \begin{array}{ccccc} \mathbf{\chi_1(\mbox{x}_1)} & \chi_2(\mbox{x}_1) & \chi_3(\mbox{x}_1) & \cdots & \chi_n(\mbox{x}_1) \\ \chi_1(\mbox{x}_2) & \mathbf{\chi_2(\mbox{x}_2)} & \chi_3(\mbox{x}_2) & \cdots & \chi_n(\mbox{x}_2) \\ \chi_1(\mbox{x}_3) & \chi_2(\mbox{x}_3) & \mathbf{\chi_3(\mbox{x}_3)} & \cdots & \chi_n(\mbox{x}_3) \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ \chi_1(\mbox{x}_n) & \chi_2(\mbox{x}_n) & \chi_3(\mbox{x}_n) & \cdots & \mathbf{\chi_n(\mbox{x}_n)} \end{array} \right |$

And the Hamiltonian operator is given by:

$\hat{H}_{elec}= \sum_{i=1}^{n} {\frac{-\nabla_i^{2}}{2}} + \sum_{i}^{n}\sum_{A}^{N} \frac{-Z_{A}}{r_{iA}} + \sum_{i}^{n}\sum_{j>i}^{n} \frac{1}{r_{ij}} + \sum_{A}^{N}\sum_{B>A}^{N} \frac{Z_{A}Z_{B}}{R_{AB}}$

This is obviously going to be complicated!! Doing all the algebra required is, at the very least, tedious. We will more modestly try to get a qualitative idea of how the final expressions are reached. You are not expected to be able to reproduce any of the derivation!

First of all, it is useful to note that the Hamiltonian is a sum of three terms:

$\mbox{I.} \;\;\; \sum_{A}^{N}\sum_{B>A}^{N}\frac{Z_{A}Z_{B}}{R_{AB}}$

This first term does not depend on any electronic coordinates.

$\mbox{II.} \;\;\; \sum_{i=1}^{n} \hat{h}_i$

This is a sum of one-electron operators h, so called because each depends only on the coordinates of one electron, and which are given by:

$\hat{h}_i = \frac{-\nabla_i^{2}}{2} + \sum_{A}^{N} \frac{-Z_{A}}{r_{iA}}$

And:

$\mbox{III.} \;\;\; \sum_{i}^{n}\sum_{j>i}^{n} \frac{1}{r_{ij}}$

This is a sum of n(n-1)/2 two-electron contributions, which each depend on the coordinates of two electrons.

The Hamiltonian can thus be written more concisely as a sum of zero-, one- and two-electron terms:

$\hat{H}_{elec}= \sum_{A}^{N}\sum_{B>A}^{N}\frac{Z_{A}Z_{B}}{R_{AB}} + \sum_{i=1}^{n} \hat{h}_i + \sum_{i}^{n}\sum_{j>i}^{n} \frac{1}{r_{ij}}$

Inserting this into the expression for the energy leads to:

$\begin{array}{lr} E_{SD} = \displaystyle{ \int \Psi \left \{ \sum_{A}^{N}\sum_{B>A}^{N}\frac{Z_{A}Z_{B}}{R_{AB}} \right \} \Psi d\mbox{x} } \; + \\ \displaystyle{ \int \Psi \left \{ \sum_{i=1}^{n} \hat{h}_i \right \} \Psi d\mbox{x} \; + \int \Psi \left \{ \sum_{i}^{n}\sum_{j>i}^{n} \frac{1}{r_{ij}} \right \} \Psi d\mbox{x} } \\ \end{array}$

The first term in this sum is an integral over a constant, since the nuclear-nuclear repulsion energy does not depend on the electronic coordinates. This therefore simply gives:

$V_{NN}= \int \Psi \left \{ \sum_{A}^{N}\sum_{B>A}^{N}\frac{Z_{A}Z_{B}}{R_{AB}} \right \} \Psi d\mbox{x} = \sum_{A}^{N}\sum_{B>A}^{N}\frac{Z_{A}Z_{B}}{R_{AB}}$

V_NN is the potential energy due to nuclear-nuclear Coulombic repulsion.

The two following terms are both integrals of sums, which can be re-written as sums of integrals, leading to:

$E_{SD} = V_{NN} + \sum_{i=1}^n \left \{ \int \Psi \hat{h}_i \Psi d\mbox{x} \right \} + \sum_i^n\sum_{j>i}^n \left \{ \int \Psi \frac{1}{r_{ij}} \Psi d\mbox{x} \right \}$

Each one-electron operator h only acts on a very small part of the wavefunction, which leads to a considerable simplification of the first term. This term is also very easy to understand: it is a sum of the one-electron energies h_ii of each orbital.

$\sum_{i=1}^{n} \left \{ \int \Psi \hat{h}_i \Psi d\mbox{x} \right \} = \sum_{i=1}^{n} \left \{ \int \chi_i \hat{h} \chi_i dr d\omega \right \} = \sum_{i=1}^n h_{ii}$

Note that the one-electron operator h includes a part corresponding to the electron's kinetic energy, and another to the potential energy created by the attractive Coulombic interaction with the nuclei. Therefore, the one-electron energy of an orbital is the sum of a kinetic and a potential part:

$h_{ii}= \left \{ \begin{array}{lll} \displaystyle{\int \chi_i \hat{h} \chi_i dr d\omega} \\ \mbox{ } \\ \displaystyle{\int \chi_i \left ( \frac{-\nabla_i^{2}}{2} + \sum_{A}^{N} \frac{-Z_{A}}{r_{iA}} \right ) \chi_i dr d\omega} \\ \mbox{ } \\ \displaystyle {\int \chi_i \frac{-\nabla_i^{2}}{2} \chi_i dr d\omega + \int \chi_i \left ( \sum_{A}^{N} \frac{-Z_{A}}{r_{iA}} \right ) \chi_i dr d\omega} \\ \mbox{ } \\ T_{e,i} + V_{Ne,i} \end{array} \right.$

And:

$\sum_{i=1}^n h_{ii}= \sum_{i=1}^n ( T_{e,i} + V_{Ne,i}) = T_e + V_{Ne}$

T_e is the electronic kinetic energy, and V_Ne is the potential energy due to nuclear-electronic Coulombic attraction.

The second set of integrals, the two-electron terms, are more complicated, because there are more of them, and because of the antisymmetrisation of the wavefunction !! Each term in the summation looks like:

$\begin{array}{ll} {\displaystyle \iint} \left \{ (\chi_i(\mbox{x}_1)\chi_j(\mbox{x}_2)- \chi_i(\mbox{x}_2)\chi_j(\mbox{x}_1) \right \} {\displaystyle \frac {1}{r_{12}}}\;\;\;\;\; \\ \;\;\;\;\;\;\;\; \left \{ (\chi_i(\mbox{x}_1)\chi_j(\mbox{x}_2)- \chi_i(\mbox{x}_2)\chi_j(\mbox{x}_1) \right \} d\mbox{x}_1d\mbox{x}_2 \end{array}$

This rather daunting expression can be multiplied out, to give four terms:

$\begin{array}{llll} \iint \chi_i(\mbox{x}_1)\chi_j(\mbox{x}_2) {\displaystyle \frac {1}{r_{12}}} \chi_i(\mbox{x}_1)\chi_j(\mbox{x}_2) d\mbox{x}_1d\mbox{x}_2 \\ -\iint \chi_i(\mbox{x}_1)\chi_j(\mbox{x}_2) {\displaystyle \frac {1}{r_{12}}} \chi_i(\mbox{x}_2)\chi_j(\mbox{x}_1) d\mbox{x}_1d\mbox{x}_2 \\ -\iint \chi_i(\mbox{x}_2)\chi_j(\mbox{x}_1) {\displaystyle \frac {1}{r_{12}}} \chi_i(\mbox{x}_1)\chi_j(\mbox{x}_2) d\mbox{x}_1d\mbox{x}_2 \\ \iint \chi_i(\mbox{x}_2)\chi_j(\mbox{x}_1) {\displaystyle \frac {1}{r_{12}}} \chi_i(\mbox{x}_2)\chi_j(\mbox{x}_1) d\mbox{x}_1d\mbox{x}_2 \end{array}$

The first and fourth, and second and third, term are identical two by two. To understand what they both mean, it helps to rewrite them in terms of the spatial orbitals and coordinates (the spin terms are not very important here), and to rearrange them a bit. Taking the first term first (!), it can be written as:

$\iint \psi_i^2(r_1) \frac {1}{r_{12}} \psi_j^2(r_2) dr_1dr_2$

Now, psi ²(r) is the probability of finding an electron at a given point in space. So this first term is simply the energy of the Coulombic interaction between an electron in orbital i with an electron in orbital j. For this reason, this integral is often called the Coulomb Integral, written in shorthand as J_ij. Because 1/r is always positive, and psi ² obviously also, this term contributes a positive energy, i.e. a destabilisation - this is what you expect from a Coulombic repulsion between electrons.

The other integrals become:

$-\iint \psi_i(r_1)\psi_j(r_2) \frac {1}{r_{12}} \psi_i(r_2)\psi_j(r_1) dr_1dr_2$

What does this term, called an Exchange Integral, and written as K_ij, mean ? Unlike for the Coulomb integral, there is no immediate classical interpretation for the Exchange integral. The name Exchange Integral comes from the fact that the two electrons exchange their positions from the left to the right of the integrand. This suggests, correctly, that it has something to do with the Pauli principle.

Remember that the probability density for two electrons is very different in an antisymmetrised Slater Determinant than in a simple Hartree Product. The total density is not simply a sum of the densities associated with each molecular orbital. (See the graphs for the two-electrons-in-a-box model in the section concerning the Pauli Principle.) Yet the expression for the repulsion energy between electrons given by the Coulomb integrals suggests that the total density is a sum of densities - it neglects all the effect of the antisymmetrisation on the density!

However, the expressions we have been deriving are for the energy of a Slater Determinant - so the antisymmetrisation effect should be there, somewhere. In fact it is the the exchange integrals which "correct" the Coulomb integrals to take into account the antisymmetry of the wavefunction. We saw that the electrons (especially those of same spin) tend to avoid each other rather more in the Slater Determinant model than in the Hartree Product model, so the Coulomb integrals should exagerrate the Coulomb repulsion of the electrons. The Exchange integrals, which are negative, compensate for this exagerration, and are typically ca. 25% in magnitude of the Coulomb integrals.

The overall contribution to the total energy of the potential energy due to electronic-electronic Coulombic repulsion, V_ee is therefore given as a difference of two terms:

$V_{ee}=J_{ee} - K_{ee} = \sum_{i}^{n_{elec}}\sum_{j>i}^{n_{elec}}(J_{ij}-K_{ij})$

Overall, the energy of a Slater Determinant is given by adding up all the terms discussed above. For the general case with matrix elements expressed as spin orbitals, one reaches the following expression:

$E_{SD} = V_{NN} + \sum_i^{n_{elec}} h_{ii} + \sum_{i}^{n_{elec}}\sum_{j>i}^{n_{elec}}(J_{ij}-K_{ij})$

A slightly different expression, which is more commonly met in books, is reached if one considers matrix elements over the spatial orbitals, for a closed-shell system (= a spin singlet where all the occupied orbitals have two electrons in them),

$E_{SD}=V_{NN} + 2\sum_i^{n_{orb}} h_{ii} + \sum_{i}^{n_{orb}}\sum_{j}^{n_{orb}}(2J_{ij}-K_{ij})$

Or:

$E_{SD} = V_{NN} + T_{e} + V_{Ne} + V_{ee}$

When this expression is computed for the Hartree-Fock wavefunction (i.e. for the Slater determinant obtained by minimising the energy with respect to all the coefficients c_ij), it is called the Total Energy, V_tot. (In some textbooks, you will find that E_SD does not include V_NN, so that the total energy is given by E_SD + V_NN.)

An Example

To get a feel for how all of the above quantities vary upon forming a molecule, let us consider the bonding in the water molecule.

(A note on units. Throughout this course, we are implicitly using a set of units called Atomic Units. These are units such that many of the constants involved in the various equations have values of 1.000 and are thus very convenient. The atomic unit of energy is also called a Hartree and is equal to 4.35 10^-18 J - rather small !! - or equivalently to 2625 kJ/mol. The values in the Table are in kJ/mol). Notice also that the zero of energy is the state where all the nuclei and electrons are infinitely separated and at rest.

It is instructive to think about the origins of chemical bonding in terms of the four contributions shown here. In previous lecture courses, the term which has usually been used to account for bonding is V_Ne - by accumulating electron density between the nuclei, you lower the overall potential energy. As you can see here, V_Ne does indeed drop dramatically upon bond formation. But T_e and V_ee go up (the electrons are confined to a smaller region of space), and V_NN does so too. The overall effect is very finely balanced! In fact, even in antibonding situations, V_Ne will decrease - but not enough to compensate for the increase in the other terms.

The overall energetic situation is summarized on the following picture:

(Note also that the experimental difference in energy is of ca. -930 kJ/mol. The difference between this value and the HF computed number is due to the approximations involved in HF theory).

The Hartree-Fock Roothan Equations

The previous paragraphs show how to calculate the energy of a Slater Determinant. The Hartree-Fock method is then very easy to understand: It simply involves optimizing the spin-orbitals chi (or in other words the coefficients c_ij defining the MOs in terms of the AOs) to give the lowest energy possible. This is done using a modified version of the variational method discussed for Hückel theory in Dr. Mulholland's second year lecture course. Like in that case, a differential equation is reduced to a set of homogeneous equations which can be solved using linear algebra. It turns out that this needs to be done in an iterative way, because to solve what are called the Fock equations (which give the coefficients of the moleculars orbitals), one already needs to know the form of all the occupied orbitals (basically because one needs to be able to calculate J_ij and K_ij)!!

In practice, a self-consistent method is used:

Choose initial orbitals. This is done using an approximate version of MO theory, e.g. some version of Hückel theory.
Solve the Fock equations, to give an improved set of orbitals.
Compare these orbitals to the previous set; If they are identical within some or other numerical convergence threshold, stop.
Repeat step 2., using the improved orbitals as input.

With a reasonable set of "guess" orbitals to start off with, 10-20 cycles are usually enough to converge the procedure to numerical accuracy.

Click here to return to the list of online lectures

Click here to return to the previous lecture.

Click here to go on to the next lecture.

This page created by Jeremy Harvey, Bristol, 2001.