Chapter 5. Including Electron-Electron Interactions and the Hartree-Fock Approximation.

[Note]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

The first term is the kinetic energy of the electrons, the second term is the kinetic energy of the nuclei, the third term is the electron-nucleus attraction energy, the fourth term is the electron-electron repulsion energy, and the last term is the nuclear repulsion energy. The double sums in the last two potential energy terms run over both sets of indices (i,j and A,B, respectively) excluding the self interaction terms and the factor of ½ is to prevent double counting (This can also be expressed with relational operator < in the sum instead of    and no factor of ½).

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.)

The last term in Eq. (5.3),

describes the electron-electron interactions; it is a nonlocal two-electron operator which prevents a single Slater determinant from being an exact solution to the many-body Schrödinger equation. Both the Hartree-Fock method and density-functional theory are approaches to solving the many-body problem which involve a transformation of the two-body operator   V ^ ee   into an effective one-body operator. In the Hartree-Fock theory this effective one-body operator will remain nonlocal, but in Kohn-Sham theory   V ^ ee   will be transformed into a local one-body operator.

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

where the total energy E can be written as a functional of the N-electron wavefunction   Ψ ( x 1 , , x N )

and in what follows we will assume that Ψ is normalized,   Ψ | Ψ = 1 . The symbol   τ  represents the integration metric over all spatial and spin coordinates

where, following convention, we have used the notation

to denote the combined spatial and spin coordinates. The notation for integration over the combined coordinate is given by

In practice, Ψ  is expanded as a product of orthonormal one-electron spin-orbitals χ 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   U ^   [Negele, p. 9]


and   { | ϕ i }  and  { | φ i }   are sets of single-particle states used to construct the many-body wavefunction. U ^   operates on a direct product space

and   U ^ i   operates on   𝒰 i ,  the subspace of the ith 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

If the Hamiltonian includes a two-electron operator,   V ^ ,  then we can use the following expression to calculate its matrix elements

where   { | ϕ i }  and  { | φ i }   are sets of single-particle states used to construct the many-body wavefunction.

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 ii = K 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] []

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,   ε ij , form a Hermitian matrix,   ε ij = ε 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   y ( r ) , 2 y ( r ) , , n y ( r ) . 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 ith electron as

then the Fock operator can then be written as

where we have used the short-hand notation   i = x i   to label the coordinate of the ith electron. We can see that the electron-electron repulsion operator   1 / r ij   has been transformed into an effective one-electron operator in which   v ^ HF ( i )   is the average repulsive field experienced by the ith electron due to all the other electrons in the system.

A unitary transformation (i.e.   UεU , with  UU = 1  ) allows us to diagonalize the matrix of Lagrange multipliers and write the Hartree-Fock equations in their canonical form

where the Hartree-Fock orbital energies   ε 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] []

One method of solving the Hartree-Fock equations is to expand the spatial part of the orbitals   ψ i   with a set of K basis functions   { φ μ ( r ) | μ = 1 , 2 , , K }

where the   { φ μ ( r ) }   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

in the closed-shell restricted Hartree-Fock (RHF) formalism, where F is the Fock matrix, C is the matrix of molecular orbital coefficients, S is the overlap matrix, and   ε   is the diagonal matrix of orbital energies.

We will only mention at this point that the Roothaan-Hall equations are a special case valid for closed-shell systems in which   C μi α = C μi β = C μi , where,   σ = α = + 1 2 ,  labels the spin-up electrons, and   σ = β = - 1 2   labels the spin-down electrons. The generalization of the Roothaan-Hall equations are the Pople-Nesbet equations

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

for spin-up electrons ( σ = α = ), with a similar expression for spin-down electrons ( σ = β = ). J-K is the two-electron repulsion contribution to Fock matrix, with J being the Coulomb part and K is the exchange contribution. It should be noted that in density-functional theory (DFT), K will not be included in the Fock matrix; instead, the exchange terms will be accounted for by contributions from an exchange functional. Similarly, electron correlations, that is, the tendency of electrons to avoid each other, are not included in the Hartree-Fock approximation. In fact, the electron correlation energy is defined as the difference between the exact energy and the Hartree-Fock-limit energy

In DFT the electron correlation energy will be accounted for by a contribution from a correlation functional.

In the rest of this chapter, we will begin to work toward our first example problem, and for the sake of simplicity we will consider only closed-shell systems in the restricted Hartree-Fock formalism and not show the explicit spin dependence. In later chapters, however, spin degrees of freedom will need to be explicitly included in our discussion of excited states and optical processes. The Fock matrix for a closed-shell system expressed in terms of our basis functions is

and the overlap matrix S is

Here we have defined a "core" Hamiltonian as

The "core" part of the Fock matrix   H core   does not depend on the values of the expansion coefficients,   C μi ,  and remains constant throughout SCF energy minimization procedure. The matrix   G μν   [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]

and again the factor of one-half is a result of the factor of two in our definition of the density matrix in the closed-shell restricted Hartree-Fock formalism

The two-electron integrals are given by

The two-electron integrals   ( μν | λσ )   should not be confused with an overlap integral and following convention Greek letters   μ,  ν,  σ λ , etc., will be used to label the basis functions and Latin letters i, j, k, etc., will label orbitals.

We are now ready to solve the Roothaan-Hall equations and calculate the orbital energies   ε i ,  the expansion coefficients   C μi ,  and the total energy E. To do this we will implement the self-consistent field (SCF) procedure as follows:

The results of Clementi and Davis [113] and Saunders [114] were used in the program to calculate the overlap matrix   S μν , the kinetic energy matrix   T μν , the electron-nucleus attraction energy matrix   V μν nucl , and the two-electron repulsion energy matrix   ( μν | λσ ) .   S μν   is calculated by the function Smatrix, KE calculates the kinetic energy matrix   T μν ,   V μν nucl   is calculated by the function e_nucl,   ( μν | λσ )   is calculated by the function ERI, and the matrix elements   H μν core   are calculated by CalcHMatrix.

[Back to Table of Contents] [Top of Chapter 5] []

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

will not, in general, be eigenstates of the orbital angular momentum operators   L ^ 2  and  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

which are simultaneous eigenfunctions of   L ^ 2  and  L ^ z . The real part of the spherical Gaussians can be obtained from

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

The goal, however, is not to construct an eigenstate of the orbital angular momentum operators, which will not be possible for most molecules in the first place; rather it is to find an accurate solution to the Schrödinger equation using the minimum number of basis functions. For molecules, this will often require Gaussians with L-values that account for the properties of the constituent atoms.

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

where the d are the contraction coefficients for a given CGF, and   φ p ( r )   are the primitive Gaussian functions, which, in the case of CGFs, have the more general form

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 , α , l , m , and n   carry the double indices of matrices, but in the example program freeMQM-1.cpp they will be stored as one-dimensional arrays.   L = l + m + n   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   (i.e., L = 1 ) .  (Please see Revision Note 0.112 and Revision Note 0.116 about the bug involving the use of the 6-31G** basis set.)

[Back to Table of Contents] [Top of Chapter 5] []
[Back to Table of Contents] [Top of Chapter 5] []

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 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 μi .

[Back to Table of Contents] [Top of Chapter 5] []

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.

[Back to Table of Contents] [Top of Chapter 5] []
  1. 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.

  2. 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.

  3. The main program, freeMQM-1.cpp, uses the function BuildG or the function BuildG4 to build the matrix   G μν , 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.

[Back to Table of Contents] [Top of Chapter 5] []