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

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

Table of Contents

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

(5.1)

H ^ = - ℏ 2 2 m ∑ i = 1 N ∇ i 2 ⁢ - ⁢ ℏ 2 2 ∑ A = 1 M ( 1 M A ⁢ ∇ A 2 ) ⁢ - ⁢ e 2 4 πε 0 ⁢ ∑ i = 1 N ( ∑ A ⁢ Z A | r i - R A | ) + ⁢ e 2 4 πε 0 ⁢ 1 2 ∑ i ≠ j N 1 r ij ⁢ + ⁢ e 2 4 πε 0 ⁢ 1 2 ∑ A ≠ B M Z A Z B R AB

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

(5.2)

H ^ = ∑ i = 1 N ⁢ ( - 1 2 ⁢ ∇ i 2 ) + ∑ A = 1 M ⁢ ( - 1 2 ⁢ M A ⁢ ∇ A 2 ) - ∑ i = 1 N ⁢ ( ∑ A ⁢ Z A | r i - R A | ) + 1 2 ⁢ ∑ i ≠ j N ⁢ 1 r ij + 1 2 ⁢ ∑ A ≠ B M ⁢ Z A ⁢ Z B R AB .

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

(5.3)

H ^ = ∑ i = 1 N ⁢ ( - 1 2 ⁢ ∇ i 2 ) - ∑ i = 1 N ⁢ ( ∑ A ⁢ Z A | r i - R A | ) + 1 2 ⁢ ∑ i ≠ j N ⁢ 1 r ij

The last term in Eq. (5.3),

(5.4)

V ^ ee = ∑ i < j N ⁢ 1 r ij ,

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

(5.5)

H ^ Ψ = E Ψ

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

(5.6)

E = E ⁡ [ Ψ ] = 〈 Ψ | H ^ | Ψ 〉 〈 Ψ | Ψ 〉 = ∫ Ψ * ⁢ H ^ ⁢ Ψ ⁢ ⅆ τ ∫ Ψ * ⁢ Ψ ⁢ ⅆ τ

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

(5.7)

ⅆ τ ≡ ⅆ x N = ⅆ x 1 ⁢ ⅆ x 2 ⁢ ⋯ ⅆ x N

where, following convention, we have used the notation

(5.8)

x i ≡ ( r i , σ i ) for i = 1 , 2 , … , N

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

(5.9)

∫ ⅆ x = ∫ ⅆ s ⁢ ⅆ r = ∑ σ = ↑ , ↓ ⁢ ∫ ⅆ r

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]

(5.10)

〈 φ 1 ⁢ φ 2 ⁢ ⋯φ N | U ^ | ϕ 1 ⁢ ϕ 2 ⁢ ⋯ϕ N 〉 = ∑ i = 1 N ⁢ 〈 φ 1 ⁢ φ 2 ⁢ ⋯φ N | U ^ i | ϕ 1 ⁢ ϕ 2 ⁢ ⋯ϕ N 〉   = ∑ i = 1 N ⁢ ∏ k ≠ i ⁢ 〈 φ k | ϕ k 〉 ⁢ 〈 φ i | U ^ 1 | ϕ i 〉

where

(5.11)

U ^ = ∑ i = 1 N ⁢ U ^ i

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

(5.12)

𝒰 1 ⊗ 2 ⊗ ⋯ ⊗ N = 𝒰 1 ⊗ 𝒰 2 ⊗ ⋯ ⊗ 𝒰 N

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

(5.13)

E = E ⁡ [ Ψ ] = 〈 Ψ | H ^ | Ψ 〉 = 〈 χ 1 ⁢ χ 2 ⁢ ⋯χ N | H ^ | χ 1 ⁢ χ 2 ⁢ ⋯χ N 〉   = ∑ i = 1 N ⁢ 〈 χ i | H ^ 1 | χ i 〉 = ∑ i = 1 N ⁢ ε i

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

(5.14)

〈 φ 1 ⁢ φ 2 ⁢ ⋯φ N | V ^ | ϕ 1 ⁢ ϕ 2 ⁢ ⋯ϕ N 〉 = 1 2 ⁢ ∑ i ≠ j N ⁢ ∏ k ≠ j k ≠ i ⁢ 〈 φ k | ϕ k 〉 ⁡ [ 〈 φ i ⁢ φ j | V ^ 12 | ϕ i ⁢ ϕ j 〉 - 〈 φ i ⁢ φ j | V ^ 12 | ϕ j ⁢ ϕ i 〉 ]

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

(5.15)

E = E ⁡ [ Ψ ] = 〈 Ψ | H ^ | Ψ 〉 = 〈 χ 1 ⁢ χ 2 ⁢ ⋯ χ N | H ^ | χ 1 ⁢ χ 2 ⁢ ⋯ χ N 〉   = ∑ i = 1 N ⁢ H ii + 1 2 ⁢ ∑ i , j = 1 N ⁢ ( J ij - K ij )

where

(5.16)

H ii = T ii + V ii nucl

(5.17)

T ii = ∫ ⅆ x χ i * ⁢ ( x ) ⁢ ( - 1 2 ⁢ ∇ 2 ) ⁢ χ i ⁢ ( x )

(5.18)

V ii nucl = ∫ ⅆ x χ i * ⁢ ( x ) ⁢ ( - ∑ A ⁢ Z A | r - R A | ) ⁢ χ i ⁢ ( x )

and the Coulomb integrals are given by

(5.19)

J ij = ∫ ∫ ⅆ x 1 ⁢ ⅆ x 2 ⁢ χ i * ⁢ ( x 1 ) ⁢ χ i ⁢ ( x 1 ) ⁢ 1 r 12 ⁢ χ j * ⁢ ( x 2 ) ⁢ χ j ⁢ ( x 2 )

and the exchange integrals are

(5.20)

K ij = ∫ ∫ ⅆ x 1 ⁢ ⅆ x 2 ⁢ χ i * ⁢ ( x 1 ) ⁢ χ j ⁢ ( x 1 ) ⁢ 1 r 12 ⁢ χ i ⁢ ( x 2 ) ⁢ χ j * ⁢ ( x 2 )

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] [scienceelearning.org]

The Rayleigh-Ritz variational method tells us that a stationary value of the equation

(5.21)

〈 Ψ | H ^ | Ψ 〉 - E 〈 Ψ | Ψ 〉 = 0

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

(5.22)

δ ⁡ { E ⁡ [ Ψ ] - E   〈 Ψ | Ψ 〉 } = 0

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]

(5.23)

δE ⁡ [ { χ i } ] - δ ⁢ ∑ i ⁢ ∑ j ⁢ ε ij ⁢ S ij = 0

with the constraint being the orthonormality of the orbitals

(5.24)

S ij ≡ ∫ χ i * ⁢ ( x ) ⁢ χ j ⁢ ( x ) ⁢ ⅆ x = 〈 χ i | χ j 〉= δ ij

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

(5.25)

δ δχ i * ⁢ ( x ) ⁡ [ E ⁡ [ { χ i } ] - ∑ i ⁢ ∑ j ⁢ ε ij ⁢ S ij ] = 0

Using the definition of the functional derivative

(5.26)

δF ⁡ [ y ⁢ ( x ) ] δy ⁢ ( z ) = lim ε -> 0 ⁢ F ⁡ [ y ⁢ ( x ) + ε δ ⁢ ( x - z ) ] - F ⁡ [ y ⁢ ( x ) ] ε

applied to an arbitrary functional in three-dimensions

(5.27)

F ⁡ [ y ] = ∫ f ⁢ ( y ⁢ ( r ) , ∇ y ⁢ ( r ) , ∇ 2 ⁢ y ⁢ ( r ) , … , ∇ n ⁢ y ⁢ ( r ) ) ⁢ ⅆ r

one can obtain the following expression for the functional derivative

(5.28)

δF ⁡ [ y ] δy ⁢ ( r ) = ∂ f ∂ y - ∇ · ( ∂ f ∂ ∇ y ) + ∇ 2 ⁢ ( ∂ f ∂ ∇ 2 ⁢ y ) - ⋯ + ( - 1 ) n ⁢ ∇ n ⁢ ( ∂ f ∂ ∇ n ⁢ y )

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,

(5.29)

δ δχ i * ⁢ ( x 1 ) ⁡ [ ∑ i = 1 N ⁢ H ii + 1 2 ⁢ ∑ i , j = 1 N ⁢ ( J ij - K ij ) - ∑ i , j = 1 N ⁢ ε ij ⁢ S ij ] = 0

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

(5.30)

δH ii δχ i * ⁢ ( x 1 ) = δT ii δχ i * ⁢ ( x 1 ) + δV ii nucl ⁢ δχ i * ⁢ ( x 1 )   = { - 1 2 ⁢ ∇ 1 2 ⁢ - ∑ A ⁢ Z A | r 1 - R A | } ⁢ χ i ⁢ ( x 1 )   = ℋ core ⁢ ( r 1 ) ⁢ χ i ⁢ ( x 1 )

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

(5.31)

1 2 ⁢ δJ ij ⁢ δχ i * ⁢ ( x 1 ) = ∫ ⅆ x 2 ⁢ χ i ⁢ ( x 1 ) ⁢ 1 r 12 ⁢ χ j * ⁢ ( x 2 ) ⁢ χ j ⁢ ( x 2 )   = 𝒥 j ⁢ ( x 1 ) ⁢ χ i ⁢ ( x 1 )

and

(5.32)

1 2 ⁢ δK ij δχ i * ⁢ ( x 1 ) = ∫ ⅆ x 2 ⁢ χ j ⁢ ( x 1 ) ⁢ 1 r 12 ⁢ χ i ⁢ ( x 2 ) ⁢ χ j * ⁢ ( x 2 )   = 𝒦 j ⁢ ( x 1 ) ⁢ χ i ⁢ ( x 1 )

respectively. Finally, for the last term we have

(5.33)

δS ij δχ i * ⁢ ( x 1 ) = χ j ⁢ ( x 1 )

Collecting the terms together gives us the Hartree-Fock equations

(5.34)

ℱ ^ χ j ⁢ ( x 1 ) = ∑ j ⁢ ε ij χ j ⁢ ( x 1 )

where the Fock operator is given by

(5.35)

ℱ ^ ⁢ ( x 1 ) = ℋ core ⁢ ( r 1 ) + ∑ j = 1 N ⁢ { 𝒥 j ⁢ ( x 1 ) - 𝒦 j ⁢ ( x 1 ) }

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

(5.36)

ℱ ^ ⁢ ( x 1 ) = ℋ core ⁢ ( r 1 ) + ∑ j = 1 N / 2 ⁢ { 2 ⁢ 𝒥 j ⁢ ( x 1 ) - 𝒦 j ⁢ ( x 1 ) }

and the "core" Hamiltonian operator is defined as

(5.37)

ℋ core ⁢ ( r 1 ) = - 1 2 ⁢ ∇ 1 2 ⁢ - ∑ A M ⁢ Z A | r 1 - R A |

with

(5.38)

𝒥 j ⁢ ( x 1 ) = ∫ ⅆ x 2 ⁢ χ j * ⁢ ( x 2 ) ⁢ 1 r 12 ⁢ χ j ⁢ ( x 2 )

being the Coulomb operator, and the exchange operator is

(5.39)

𝒦 j ⁢ ( x 1 ) ⁢ χ i ⁢ ( x 1 ) = [ ∫ ⅆ x 2 ⁢ χ j * ⁢ ( x 2 ) ⁢ 1 r 12 ⁢ χ i ⁢ ( x 2 ) ] ⁢ χ j ⁢ ( x 1 )

If we define the Hartree-Fock potential operator of the ith electron as

(5.40)

v ^ HF ⁢ ( x i ) = ∑ j = 1 N ⁢ [ 𝒥 j ⁢ ( x i ) - 𝒦 j ⁢ ( x i ) ]

then the Fock operator can then be written as

(5.41)

ℱ ^ ⁢ ( i ) = ℋ core ⁢ ( i ) + v ^ HF ⁢ ( i )

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

(5.42)

ℱ ^ χ i ⁢ ( x ) = ε i χ i ⁢ ( x )

where the Hartree-Fock orbital energies   ε i   are given by

(5.43)

ε i ≡ 〈 χ i | ℱ ^ | χ i 〉 = H ii + ∑ j = 1 N ⁢ ( J ij - K ij )

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

(5.44)

E = ∑ i = 1 N ⁢ ε i - 1 2 ⁢ ∑ i = 1 N ⁢ ∑ j = 1 N ⁢ ( J ij - K ij )

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

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 } , 

(5.45)

ψ i ⁢ ( r ) = ∑ μ = 1 K ⁢ C μi ⁢ φ μ ⁢ ( r )

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

(5.46)

φ = N ⁢ x l ⁢ y m ⁢ z n ⁢ e - αr 2

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

(5.47)

F ⁢ C = SC ⁢ ε

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

(5.48)

F α ⁢ C α = SC α ⁢ ε α

(5.49)

F β ⁢ C β = SC β ⁢ ε β

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

(5.50)

F α = H core + J - K

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

(5.51)

E corr ≡ E exact - E HF .

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

(5.52)

F μν = H μν core + G μν

and the overlap matrix S is

(5.53)

S μν = ∫ ⅆ r φ μ * ⁢ ( r ) ⁢ φ ν ⁢ ( r )

Here we have defined a "core" Hamiltonian as

(5.54)

H μν core = T μν + V μν nucl

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

(5.55)

G μν = 2 J μν - K μν

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]

(5.56)

J μν = 1 2 ∑ λσ ⁢ P λσ ⁢ ( μν | λσ ) ,     and     K μν = 1 2 ∑ λσ ⁢ P λσ ⁢ ( μλ | νσ )

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

(5.57)

P λσ = 2 ⁢ ∑ i = 1 N / 2 ⁢ C λi * ⁢ C σi

The two-electron integrals are given by

(5.58)

( μν | λσ ) = ∫ ∫ ⅆ r 1 ⁢ ⅆ r 2 ⁢ φ μ * ⁢ ( r 1 ) ⁢ φ ν ⁢ ( r 1 ) ⁢ 1 r 12 ⁢ φ λ * ⁢ ( r 2 ) ⁢ φ σ ⁢ ( r 2 )

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:

  1. First, calculate all the matrix elements that remain constant throughout the calculation and are not changed by the SCF energy minimization procedure. These include the overlap matrix S,   H core ,   ( μν | λσ ) , and V. The matrix V diagonalizes the overlap matrix, i.e.,   V † SV = I .  

  2. Make an initial guess for the coefficients, e.g.,   C μi = 1 .  

  3. Using the coefficients   C μi ,  calculate the matrix G and add it to   H core   to obtain the Fock matrix F.

  4. Calculate the matrix   F ' = V † FV ,  and solve   F ' C ' = εC '   for the eigenvalues and eigenvectors. ε is a diagonal matrix containing the orbital energies, and C' is the matrix of new coefficients containing the eigenvectors in rows or columns, depending on whether one is programming in C/C++ or Fortran.

  5. Perform the transformation back to the original basis   C = VC ' .   At this point, we should be sure that C†SC=I,   F ⁢ C = SC ⁢ ε ,   and C†FC should be diagonal.

  6. Compare the total energy to the previous value, and if the difference is significant, repeat the procedures starting from step 3. At this point, one may also calculate the total energy using the new   C μi .

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] [scienceelearning.org]

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

(5.59)

φ = N ⁢ x l ⁢ y m ⁢ z n ⁢ e - αr 2 ,

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]

(5.60)

φ = N ⁢ r l ⁢ e - αr 2 2 ⁢ ( Y l - m ± Y l m )

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

(5.61)

φ = N ⁢ r n - 1 ⁢ e - αr 2 ⁢ Y l m ,

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

(5.62)

φ = N ⁢ r n - 1 ⁢ e - αr 2 2 ⁢ ( Y l - m ± Y l m ) ,

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

(5.63)

φ = N ⁢ r n - 1 ⁢ e - ζr 2 ⁢ ( Y l - m ± Y l m ) .

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],

(5.64)

g ⁢ ( A , α μ , l μ , m μ , n μ ) = N ( x - x A ) l μ ⁢ ( y - y A ) m μ ⁢ ( z - z A ) n μ ⁢ e - α μ ⁢ | r - r A | 2 ,

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

(5.65)

φ μ ⁢ ( r ) = ∑ p = 1 L ⁢ d pμ ⁢ φ p ⁢ ( r )

where the dpμ 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

(5.66)

g ⁢ ( A , α pμ , l pμ , m pμ , n pμ ) = N ( x - x A ) l pμ ⁢ ( y - y A ) m pμ ⁢ ( z - z A ) n pμ ⁢ e - α pμ ⁢ | r - r A | 2 .

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 pμ , α pμ , l pμ , m pμ , and n pμ   carry the double indices of matrices, but in the example program freeMQM-1.cpp they will be stored as one-dimensional arrays.   L = l pμ + m pμ + n pμ   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

(5.67)

ψ i ⁢ ( r ) = ∑ μ = 1 K ⁢ ∑ p = 1 L ⁢ C μi ⁢ d pμ ⁢ φ p ⁢ ( r )

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] [scienceelearning.org]

Program details.

A compressed bzip2 file TutorialC.tar.bz2 can be downloaded, which includes the XHTML version of the text that you are reading right now. Those using Internet Explorer 5.5 with MathPlayer 1.0 can download the tutorial with the HTML version of the text TutorialC.tar.gz. The bzip2 file can be decompressed by the compression utility PowerArchiver for Windows users; Linux/Unix users can decompress the file using file-roller, or by typing the command:

prompt>tar -xjvf TutorialC.tar.bz2

This should create the directory Tutorial containing the subdirectories program-1, Basis_sets, and routines. The subdirectory program-1 contains the makefile GNUmakefile, the source code for the first example program, and the precompiled binary freeMQM-1_32bit. The binary freeMQM-1_32bit was compiled using gcc on Redhat Linux 8.0, but the main program itself is very simple and should compile on any platform (assuming LAPACK and GSL are installed). To recompile the program, simply type prompt>make from inside the directory containing GNUmakefile.

The program makes a few function calls to the open source GNU Scientific Library (GSL) and the LAPACK libraries. The GSL libraries are written in C and are a part of most Linux distributions. They can also be downloaded at http://www.gnu.org/software/gsl, and binaries compiled for Windows using Visual C++ used to be available for download at http://www.network-theory.co.uk/gsl/. LAPACK is written in Fortran, but the LAPACK routines can still be used in C/C++ programs by including the cfortran.h header in the C/C++ code and by including prototype definitions of the kind declared the beginning of the eigen.cpp file. The header file cfortran.h can be found in the routines folder, and documentation and other details can be found at http://www-zeus.desy.de/~burow/cfortran/. The public domain Numerical Recipes header file nrtypes_nr.h has been included in freeMQM-1.h to enable array dimensions to be defined at run time. Similarly, Mat4d.h has been included to do the same for the four dimensional arrays that are required for the two-electron integrals.

Windows users can also download the Cygwin package from http://www.cygwin.com. This package allows the user to create a basic Linux-like environment, and a range of programs, such as the GNU C/C++ and FORTRAN compilers, have been ported to the Cygwin API. Apparently, the LAPACK and BLAS numerical libraries are not currently included in the Cygwin package. These are easy to install, however. Simply download the file lapack.tgz from from http://www.netlib.org/lapack and follow the installation instructions. On a system running Cygwin, this should generate the libraries lapack_LINUX.a and blas_LINUX.a. The GNUmakefile included with freeMQM-1.cpp wants these libraries to be renamed as liblapack.a and libblas.a, and they should be copied to a directory were the compiler can find them, such as /YOUR_CYGWIN_INSTALLATION_DIRECTORY/cygwin/lib/.

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

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 dpμ 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] [scienceelearning.org]

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

Basis SetGround State EnergyBond 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.
[Back to Table of Contents] [Top of Chapter 5] [scienceelearning.org]

  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] [scienceelearning.org]
Prev Up Next
Chapter 4. Angular momentum in many-electron atoms and molecules. Home Chapter 6. A Brief Overview of Density-Functional Theory.