- 5.1. The many-electron Hamiltonian.
- 5.2. The Hartree-Fock equations.
- 5.3. The Roothaan-Hall equations: equations we can solve (numerically).
- 5.4. The basis sets in more detail.
- 5.5. The first example program.
- 5.6. The Hartree-Fock method applied to our first atom: helium.
- 5.7. The Hartree-Fock method applied to our first molecule: hydrogen.
- 5.8. Things to do.

Development note | |
---|---|

Changes to this chapter were last made on 27 February 2005. |

In Section 3.1 we described a fictitious system of noninteracting electrons in which a
*Slater determinant* of spin-orbitals would be an exact solution. A more realistic
molecular Hamiltonian operator representing a system of *N* electrons bound to
*M* nuclei could include more than thirty-two terms all together [Brown, p. 118].
As the next step toward improved accuracy, however, we will consider the Hamiltonian operator written as (in S.I. units)

and in atomic units (a.u.) it is given by

At this point, we will apply the Born-Oppenheimer approximation
and assume that the positions of the nuclei can be treated as though they are fixed
relative to the faster moving electrons. In this approximation, the nuclear repulsion energy can be calculated separately from the electronic energy. In the first demonstration program,
`freeMQM-1.cpp`, for example,
the last step of the energy calculation is the adding of the nuclear repulsion energy to the electronic energy to obtain the total energy.
The Hamiltonian that describes just the electronic degrees of freedom, without the kinetic and repulsion energy of the nuclei, is simply (in a.u.)

In principle, we are interested in solving an eigenvalue problem of the form

where, following convention, we have used the notation

In practice, $\Psi $ is expanded as a product of orthonormal one-electron spin-orbitals ${\chi}_{i}$, as described in Section 3.3, and an energy minimization procedure will lead to an approximate solution. For a noninteracting system of electrons, the total energy can be easily calculated by using the following equation for the matrix elements of an arbitrary single-particle operator $\hat{U}$ [Negele, p. 9]

and
${\hat{U}}_{i}$
operates on
${\mathcal{U}}_{i}$,
the subspace of the *i*th electron.
By using Eq. (5.10), it is easy to see that the energy expectation value of the Hamiltonian operator of Section 3.1, describing a system of *N* noninteracting electrons, is simply the sum of the one-electron energies

Equation (5.14) and Eq. (5.10) can be used to express the energy expectation value of the Hamiltonian in Eq. (5.3) as a sum of one-electron and two-electron integrals, as follows

and the Coulomb integrals are given by

and the exchange integrals are

where the fact that ${J}_{\mathrm{ii}}={K}_{\mathrm{ii}}$ has been used to include i = j in the sum. This is one point that must be kept in mind when trying to understand the physical meaning the energies being calculated. One may obtain, for example, one-electron atoms with a nonzero electron-electron repulsion energy and two electron atoms or molecules with a nonzero ground state exchange energy. Thus, caution must be exercised to be sure that no net self-interaction energy exists. In Hartree-Fock theory the self-interaction from the Coulomb contribution should be cancelled by exchange terms [Szabo,pp.86,135,210; Leach,p.51], but in density-function theory this kind of cancellation may require additional effort [Parr, p.181].

[Back to Table of Contents] [Top of Chapter 5] [scienceelearning.org]The Rayleigh-Ritz variational method tells us that a stationary value of the equation

will yield the solution to the eigenvalue problem in Eq. (5.5). In other words, we want to find a Ψ such that

In practice, the variational condition is recast in terms of the one-electron orbitals and the calculus of variations with constraints is applied [Szabo,pp.34,116; Leach,pp.18,52; Arfken,p.954]

with the constraint being the orthonormality of the orbitals

and the Lagrange multipliers, ${\epsilon}_{\mathrm{ij}}$, form a Hermitian matrix, ${\epsilon}_{\mathrm{ij}}={\epsilon}_{\mathrm{ji}}^{*}$, [Szabo, p.118; Leach, p. 52]. If a solution to the stationary value problem exists, then it must satisfy the Euler-Lagrange equation, expressed here as a functional derivative

Using the definition of the functional derivative

applied to an arbitrary functional in three-dimensions

one can obtain the following expression for the functional derivative

where *f* is a function of *y*(**r**) and its derivatives
$\nabla y\left(\mathbf{r}\right),{\nabla}^{2}y\left(\mathbf{r}\right),\dots ,{\nabla}^{n}y\left(\mathbf{r}\right)$.
If we now apply the functional derivative to the various terms in Eq. (5.25), rewritten here in its
expanded form,

we obtain from Eqs. (5.16) to (5.18),

Similarly, for the Coulomb and exchange terms in Eqs. (5.19) and (5.20), we find

respectively. Finally, for the last term we have

Collecting the terms together gives us the Hartree-Fock equations

where the Fock operator is given by

or for a closed-shell system [Szabo,p.84;Leach,p.53]

and the "core" Hamiltonian operator is defined as

being the Coulomb operator, and the exchange operator is

If we define the Hartree-Fock potential operator of the *i*th electron as

then the Fock operator can then be written as

where the Hartree-Fock orbital energies ${\epsilon}_{i}$ are given by

By comparing Eq. (5.15) and Eq. (5.43), it can be seen that the Coulomb and exchange terms picked up an additional factor of two in the derivation of the Hartree-Fock equations. This additional contribution is because the sum of the orbital energies counts the electron-electron interactions twice and must be corrected for in the calculation of the total energy [Szabo, p. 125]. Thus, the total energy is

[Back to Table of Contents] [Top of Chapter 5] [scienceelearning.org]

where the
$\left\{{\phi}_{\mu}\left(\mathbf{r}\right)\right\}$
could represent any complete set of basis functions, but preferably they will
be good approximations to the atomic orbitals so that a small subset will
give reasonable results. A finite set of *Slater-type orbitals* (STO) or
*Gaussian type functions* (GTF) [Szabo,p.153] are two
common choices. For the moment, we will restrict ourselves to GTF's, which take the following form

where *N* is a normalization constant and *L = l + m + n*
is used to identify Gaussian function as being s-, p-, d-, or f-type, etc., for *L* = 0, 1, 2, 3...(s, p, d, f,...), respectively.
We will also use the notation that Greek letters
μ, ν, σ, λ, etc., will be used to label the basis functions and Latin letters
*i, j, k,*, etc., will label the orbitals, such as in Eq. (5.45) where
the orbitals are expanded as a linear combination of the basis functions.

Substituting our expansion of the spatial orbitals in terms of known basis functions [Eq. (5.45)] into the Hartree-Fock equation [Eq. (5.42)] leads to the Roothaan-Hall equations

valid for open-shell unrestricted spin-orbitals, and in this case the Fock matrix can be written as

Here we have defined a "core" Hamiltonian as

The "core" part of the Fock matrix ${\mathbf{H}}^{\mathrm{core}}$ does not depend on the values of the expansion coefficients, ${C}_{\mathrm{\mu i}}$, and remains constant throughout SCF energy minimization procedure. The matrix ${G}_{\mathrm{\mu \nu}}$ [Szabo,p.141; Thijssen,p.75] does depend on the expansion coefficients, and it includes contributions from both the Coulomb repulsion matrix and the exchange matrix

where the factor of two results from the assumption that we are dealing with a closed-shell system [Szabo,p.84]. The Coulomb and exchange matrices are related to the two-electron integrals as follows [Leach,p.58; Koch,p.96]

The two-electron integrals are given by

The results of Clementi and Davis [113] and
Saunders [114] were used in the program to calculate
the overlap matrix
${S}_{\mathrm{\mu \nu}}$, the kinetic energy matrix
${T}_{\mathrm{\mu \nu}}$, the electron-nucleus attraction energy matrix
${V}_{\mathrm{\mu \nu}}^{\mathrm{nucl}}$, and the two-electron repulsion energy matrix
$\left(\mathrm{\mu \nu}|\mathrm{\lambda \sigma}\right)$.
${S}_{\mathrm{\mu \nu}}$
is calculated by the function `Smatrix`,
`KE` calculates the kinetic energy matrix
${T}_{\mathrm{\mu \nu}}$,
${V}_{\mathrm{\mu \nu}}^{\mathrm{nucl}}$ is calculated by the function `e_nucl`,
$\left(\mathrm{\mu \nu}|\mathrm{\lambda \sigma}\right)$ is calculated by the function `ERI`,
and the matrix elements
${H}_{\mathrm{\mu \nu}}^{\mathrm{core}}$
are calculated by `CalcHMatrix`.

Basis functions constructed from Gaussian type functions (GTF), written previously as

will not, in general, be eigenstates of the orbital angular momentum operators ${\hat{\mathbf{L}}}^{2}\text{\hspace{0.17em}and\hspace{0.17em}}{\hat{L}}_{z}$; however, linear combinations of the Gaussians can be used to create functions of the form [Levine,p.463]

and, of course, a *spherical* form of the GTF can be written as

a form that closely resembles the STOs used to describe nonlinear molecules, i.e.,

For most three-dimensional problems, Cartesian GTFs [113] will be used [114][Szabo, p. 181;Levine,p.462],

The Cartesian GTFs can be further generalized to be used as part of a contracted basis set; these contracted Gaussian functions (CGF)[Szabo,p.155] are written as

The functions in Eq. (5.54) are the kind used in the first example program, `freeMQM-1.cpp`.
In Eq. (5.54), the constants
${d}_{\mathrm{p\mu}}\text{,}{\alpha}_{\mathrm{p\mu}}\text{,}{l}_{\mathrm{p\mu}}\text{,}{m}_{\mathrm{p\mu}}\text{, and}{n}_{\mathrm{p\mu}}$
carry the double indices of matrices, but in the example program
`freeMQM-1.cpp`
they will be stored as one-dimensional arrays.
$L={l}_{\mathrm{p\mu}}+{m}_{\mathrm{p\mu}}+{n}_{\mathrm{p\mu}}$
is the total-electronic-orbital-angular momentum quantum number of the
atomic orbital represented by the CGF and should not be confused with the
number of primitive Gaussian functions, also labeled *L*, used in the
expansion of the CGF. In terms of the primitive Gaussian functions, the molecular
orbitals now can be written as

The first example program `freeMQM-1.cpp` only offers three choices of basis sets: one non-contracted
basis, 6-31G*, and 6-31G**. The non-contracted basis has the same exponents as the
contracted s-type functions used in 6-31G* and 6-31G**.
For both the hydrogen and helium atoms, the difference between 6-31G* and 6-31G** is that 6-31G** includes
a single set of p-type polarization functions
$\text{(i.e.,}L=1\text{)}$.
(Please see Revision Note 0.112 and Revision Note 0.116 about the bug involving the use of the 6-31G** basis set.)

The first example program, `freeMQM-1.cpp`, calculates the ground state energy of the helium atom or the hydrogen molecule using the Hartree-Fock method, or
by performing a density-functional calculation. At the moment, the user has the choice of
three different basis sets found in the `Basis_sets` directory. The basis
set files are in the GAMESS format [52], with an additional line of code describing the contents of the
file (this requirement could easily be coded away). As can be seen from the
`GNUmakefile`, the function prototypes and global
variable definitions for `freeMQM-1.cpp` can be found in `freeMQM-1.h`.
The routines in `freeMQM-1.cpp` also make function calls to
functions in `eigen.cpp` that solve the eigenvalue problem. Choosing to perform the Hartree-Fock
calculation for helium gives a total energy of -2.85516 a.u. All three basis sets give the same
result. The non-contracted basis has the same exponents as the 6-31G* basis, and
the fact that all three basis sets give the same result illustrates two points. First, the contraction
coefficients
*
d _{pμ}
*
are optimized properly. Second, the polarization functions included in 6-31G** do not
contribute to the spherically symmetric ground state of the helium atom, as can be seen
from the calculated values of the expansion coefficients
${C}_{\mathrm{\mu i}}$.

The non-contracted basis gives a ground state energy of -1.12655 a.u. and a bond length of
1.389 a.u. The result using the basis 6-31G** is consistent
with that found on page 192 and page 201 of Szabo and Ostlund
[1], namely a total energy of -1.131 a.u. and a bond length of 1.385 a.u. As can be seen from the table below, the inclusion of the *L*=1 polarization
functions in 6-31G** does give a moderate improvement over the contracted basis set without polarization functions, 6-31G*.
In fact, in the program, the two hydrogen atoms are defined
to lie along the z-axis, and as can be seen from the values of the coefficients `C[I]`, it is only the
z-polarization function that contributes to
the molecular ground state wavefunction. In the program, the user may also vary the bond distance
to find the energy minimum.

**Table 5.1. Hartree-Fock results for H_{2}.**

Basis Set | Ground State Energy | Bond length |
---|---|---|

Exact/Experiment | -1.174474983 a.u. (Exact [Lowe,p.232]) | 1.401 a.u. (Experiment [1]) |

non-contracted | -1.12655 a.u. | 1.389 a.u. |

6-31G* | -1.12683 a.u. | 1.379 a.u. |

6-31G** | -1.13133 a.u. | 1.384 a.u. |

The program can read basis sets of the GAMESS format, but in most cases it does not know what to do with the information. The file

`freeMQM-1.cpp`needs to be generalized to process most p-type and d-type basis sets.The functions adapted from Clementi and Davis [113] and Saunders [114] are generally valid for any (i.e. s- ,p- , d-, and f-type ) GTF. This fact, however, is not taken advantage of in other parts of the program. The function that calculates the kinetic energy,

`KE`, for example, makes calls to the function`dSMatrix`, which calculates derivatives of the**S**matrix. The function`dSMatrix`could easily be generalized to perform derivatives on any**S**matrix involving GTFs.The main program,

`freeMQM-1.cpp`, uses the function`BuildG`or the function`BuildG4`to build the matrix ${G}_{\mathrm{\mu \nu}}$, depending on whether one is performing a Hartree-Fock calculation or a density-functional calculation. The function`BuildG`was adopted from [Thijssen,see Eq. (4.86) and chapter 9 example program`nucl_md.f`], and it takes full advantage of symmetries to reduce the number of integrals that need to be performed. The function`BuildG4`is used in the Hartree-Fock calculation and does not take advantage of any of the symmetry relations. It is a relic of the debugging process and should be eliminated.