Chapter 6. A Brief Overview of Density-Functional Theory.

	Development note
	Changes to this chapter were last made on 16 January 2005.

6.1. The Hohenberg-Kohn theorems.

In this section, a brief overview of the formal foundations of density-functional theory will be given [56] [58][59][60][2][3] [Leach,p.127]. The seminal work on DFT appeared in two key papers in 1964 and 1965. In the first paper, it was proved by Hohenberg and Kohn (HK) [56] that for a N-particle system all of the electronic ground state properties are determined by the electron density ρ ⁢ ( r ) , and that if ρ ⁢ ( r ) is the true ground state density, then the total energy functional, E ⁡ [ ρ ⁢ ( r ) ] , is a minimum. To prove the Hohenberg-Kohn theorems, we will start by considering the same Hamiltonian as in Eq. (5.3), rewritten as

(6.1)

H ^ ≡ T ^ + V ^ ee + ∑ i = 1 N v ⁢ ( r i )

where the kinetic energy and electron-electron repulsion energy operators are

(6.2)

T ^ = ∑ i = 1 N ⁢ ( - 1 2 ⁢ ∇ i 2 ) V ^ ee = ∑ i < j N ⁢ 1 r ij

and v(r) is an arbitrary external potential that is often assumed simply to be the electron-nucleus attraction energy, in which case Eq. (6.1) would be the same as Eq. (5.3). The total energy, E ⁡ [ Ψ ] , written as a functional of the normalized N-electron wavefunction, Ψ ⁢ ( x 1 , … , x N ) , is given by

(6.3)

where we have used the fact that

(6.4)

ρ ⁢ ( r ) = N ⁢ ∫ | Ψ ⁢ ( x , x 2 , … , x N ) | 2 ⁢ ⅆ s ⅆ x 2 ⁢ ⋯ ⁢ ⅆ x N

and

(6.5)

N = ∫ ρ ⁢ ( r ) ⁢ ⅆ r

To show that v(r) is also an unique functional of ρ ⁢ ( r ) , we first assume that the opposite is true―that two different external potentials v(r) and v'(r) can have the same electron distribution ρ ⁢ ( r ) ―and then show that this leads to a contradiction. The first part of the Hamiltonian, T ^ + V ^ ee , is system independent and universally valid. It is the external potential v(r) which defines the difference between Hamiltonians. Thus, two different external potentials, v(r) and v'(r), correspond to two different Hamiltonians, wavefunctions, and ground state energies, H ^ , Ψ , E o , and H ^ ′ , Ψ ′ , E o ′ , respectively. If we let Ψ ′ serve as a trial wavefunction for the Hamiltonian H ^ , then the Rayleigh-Ritz minimal principle for the ground state,

(6.6)

E ⁡ [ Ψ ] ≥ E o = min Ψ ⁢ E ⁡ [ Ψ ]

tells us that

(6.7)

E o < 〈 Ψ ′ | H ^ | Ψ ′ 〉 = 〈 Ψ ′ | H ^ ′ | Ψ ′ 〉 + 〈 Ψ ′ | H ^ - H ^ ′ | Ψ ′ 〉 = E o ′ + 〈 Ψ ′ | H ^ - H ^ ′ | Ψ ′ 〉 = E o ′ + ∫ ρ ⁢ ( r ) ⁡ [ v ⁢ ( r ) - v ′ ⁢ ( r ) ]

We can repeat the derivation of Eq. (6.7) with the primed and unprimed quantities interchanged to obtain

(6.8)

E o ′ < E o + ∫ ρ ⁢ ( r ) ⁡ [ v ′ ⁢ ( r ) - v ⁢ ( r ) ]

By adding Eqs. (5.15) and (5.16) we find the following contradiction

(6.9)

E o + E o ′ < E o + E o ′

Thus, the ground state density ρ ⁢ ( r ) uniquely determines the external potential v(r), to within a trivial additive constant. Since v(r) fixes the Hamiltonian H ^ , it follows that the full N-particle ground state wavefunction, and all other properties determined by H ^ , are implicitly determined by the ground state density ρ ⁢ ( r ) .

For the second HK theorem (the HK variational principle) it was shown that if ρ ⁢ ( r ) is the true ground state density then the total energy functional E ⁡ [ ρ ⁢ ( r ) ] is a minimum. In other words, the HK variational principle states that

(6.10)

E ⁡ [ ρ ⁢ ( r ) ] ≥ E ⁡ [ ρ o ⁢ ( r ) ] ≡ E o

where ρ o ⁢ ( r ) and E o are the exact ground state density and energy, respectively, and the equality holds if ρ ⁢ ( r ) = ρ o ⁢ ( r ) . To prove this, HK defined for a given external potential, v(r), the following energy functional

(6.11)

E ⁡ [ ρ ⁢ ( r ) ] = ∫ ρ ⁢ ( r ) ⁢ v ⁢ ( r ) ⁢ ⅆ r + F ⁡ [ ρ ⁢ ( r ) ]

Equation (6.11) is justified by the first HK theorem, which established that H ^ , Ψ , and E are unique functionals of ρ(r). The functionals F[ρ] and E xc ⁡ [ ρ ] are defined by the following equations [Parr, p.56]

(6.12)

F ⁡ [ ρ ] = min Ψ → ρ ⁢ 〈 Ψ | T ^ + V ^ ee | Ψ 〉 = T ⁡ [ ρ ] + V ee ⁡ [ ρ ] = T s ⁡ [ ρ ] + J ⁡ [ ρ ] + E xc ⁡ [ ρ ]

where T s ⁡ [ ρ ] is the kinetic energy of a noninteracting system of electrons, and the classical part of V ee ⁡ [ ρ ] is given by

(6.13)

J ⁡ [ ρ ] = 1 2 ⁢ ∫ ∫ ⅆ r ⁢ ⅆ r′ ⁢ ρ ⁢ ( r ) ⁢ ρ ⁢ ( r′ ) | r - r′ |

The exchange-correlation functional E xc ⁡ [ ρ ] is defined to include the nonclassical part of V ee ⁡ [ ρ ] , as well as the difference between T s ⁡ [ ρ ] and the true kinetic energy T ⁡ [ ρ ] . The functionals in the last line of Eq. (6.12) will be important for our discussion of the Kohn-Sham theory.

The functional F[ρ] in Eq. (6.12) is defined by the constrained-search for the minimum value of 〈 Ψ | T ^ + V ^ ee | Ψ 〉 over all antisymmetric wavefunctions that yield the correct input density. This definition of the energy functional, together with the variational principle with respect to Ψ [Eq. (6.6)], are sufficient to prove the second HK theorem as follows

(6.14)

E o = min Ψ ⁢ 〈 Ψ | H ^ | Ψ 〉 = min ρ ⁢ { min Ψ → ρ ⁢ 〈 Ψ | H ^ | Ψ 〉 } = min ρ ⁢ E ⁡ [ ρ ]

This is a condensed version of the proof originally given by Levy [73][74] and Lieb [75]. It has a slightly broader range of validity than the one first put forth by HK [56]. Those interested in a more thorough discussion should read the reviews by Parr and Wang [2], Kohn [60], and Koch and Holthausen [3].

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

6.2. The Kohn-Sham method.

One year and three days after Hohenberg and Kohn presented their groundbreaking work setting the formal foundations of density-functional theory, Kohn and Sham (KS) provided a practical procedure to use DFT in calculations of electronic structure [58]. The derivation of the Kohn-Sham orbital equations and the self-consistent field (SCF) energy minimization procedure is similar to the derivation of the Hartree-Fock equations in Section 5.2. As in the Hartree-Fock case we can express the N-electron wavefunction in terms of of orthonormal one-electron spin states, | Φ 〉 = | ϕ 1 ⁢ ϕ 2 ⁢ ⋯ϕ N 〉 , where we use ϕ i to distinguish the Kohn-Sham spin-orbitals from the Hartree-Fock orbitals [e.g., Eq. (3.36) and Eq. (3.41)]. Φ can also be represented by a normalized N-orbital Slater determinant

(6.15)

Φ ⁢ ( x 1 , x 2 , … , x N ) = 1 N ! | ϕ 1 ⁢ ( x 1 ) ϕ 2 ⁢ ( x 1 ) ⋯ ϕ N ⁢ ( x 1 ) ϕ 1 ⁢ ( x 2 ) ϕ 2 ⁢ ( x 2 ) ⋯ ϕ N ⁢ ( x 2 ) ⋮ ⋮ ⋮ ϕ 1 ⁢ ( x N ) ϕ 2 ⁢ ( x N ) ⋯ ϕ N ⁢ ( x N ) |

For the general case of an arbitrary number of spin-orbitals, the kinetic energy functional T s ⁡ [ ρ ↑ , ρ ↓ ] of a system of noninteracting electrons [Eq. (6.12)] can be written in terms of the one-electron spin-orbitals as

(6.16)

T s ⁡ [ ρ ↑ , ρ ↓ ] = ∑ σ = ↑ , ↓ ⁢ ∑ i ⁢ f iσ ⁢ ∫ ⅆ rϕ iσ * ⁢ ( r ) ⁢ ( - 1 2 ⁢ ∇ 2 ) ⁢ ϕ iσ ⁢ ( r ) .

The Kohn-Sham orbitals can be written in a form similar to the spin-orbitals defined in Eq. (2.37)

(6.17)

ϕ i ⁢ ( x ) = 〈 x | ϕ i 〉 = 〈 r , σ ′ | ϕ iσ 〉 = ϕ iσ ⁢ ( r ) = ϕ i ⁢ ( r ) ⁢ σ ⁢ ( σ ′ )

where for the last term in Eq. (6.17) we have used the perhaps too subtle notation that without a spin label of some sort (i.e., x or σ) ϕ i , refers only to the spatial part of the spin-orbital (this notation may be changed later if it is found to result in confusion). The spin-orbitals in Eq. (6.17) are related to the orbitals written in the Slater determinant [Eq. (6.15)] as follows

(6.18)

ϕ 1 ⁢ ( x 1 ) = ϕ 1 ↑ ⁢ ( r 1 ) , ϕ 2 ⁢ ( x 2 ) = ϕ 1 ↓ ⁢ ( r 2 ) , ϕ 3 ⁢ ( x 3 ) = ϕ 2 ↑ ⁢ ( r 3 ) , ϕ 4 ⁢ ( x 4 ) = ϕ 2 ↓ ⁢ ( r 4 ) , ....etc.

The total density and the density of spin-up electrons ρ ↑ and spin-down electrons ρ ↓ are related by

(6.19)

ρ ⁢ ( r ) = ∑ σ = ↑ , ↓ ⁢ ρ σ ⁢ ( r ) = ∑ σ = ↑ , ↓ ⁢ ∑ i ⁢ f iσ ⁢ | ϕ iσ ⁢ ( r ) | 2

In this general case, the number of one-electron spin-orbitals ϕ iσ ⁢ ( r ) may be arbitrary, but there are restrictions on the values of their occupation numbers: the Pauli exclusion principle requires that 0 ≤ f iσ ≤ 1 and the occupation numbers f iσ are related to the total number of electrons, N, by

(6.20)

N = N ↑ + N ↓ = ∑ σ = ↑ , ↓ ⁢ N σ = ∑ σ = ↑ , ↓ ⁢ ∑ i ⁢ f iσ

For simplicity, we will consider here the special case that there are N occupied orbitals with f iσ = 1 , and f iσ = 0 for the rest. The kinetic energy functional T s ⁡ [ ρ ↑ , ρ ↓ ] is then given by

(6.21)

T s ⁡ [ ρ ↑ , ρ ↓ ] = ∑ σ = ↑ , ↓ ⁢ ∑ i = 1 N σ ⁢ ∫ ⅆ rϕ iσ * ⁢ ( r ) ⁢ ( - 1 2 ⁢ ∇ 2 ) ⁢ ϕ iσ ⁢ ( r )

and the density is

(6.22)

ρ ⁢ ( r ) = ∑ σ = ↑ , ↓ ⁢ ρ σ ⁢ ( r ) = ∑ σ = ↑ , ↓ ⁢ ∑ i = 1 N σ ⁢ | ϕ iσ ⁢ ( r ) | 2

The variational approach to the energy eigenvalue problem is analogous to the derivation of the Hartree-Fock equations. Thus, the stationary condition

(6.23)

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

leads to the Euler-Lagrange equation

(6.24)

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

with the constraint being the orthonormality of the Kohn-Sham spin-orbitals

(6.25)

S ij ≡ ∫ ϕ iσ * ⁢ ( r ) ⁢ ϕ jσ ⁢ ( r ) ⁢ ⅆ r = 〈 ϕ iσ | ϕ jσ 〉= δ ij

and the Lagrange multipliers, ε ij σ , form a Hermitian matrix, ε ij σ = ε ji σ* [Szabo, p.118; Leach, p. 52]. E[ρ] is defined in Eq. (6.11), and F[ρ] is defined in Eq. (6.12), the last line of which, we repeat here for clarity [98][Parr, p.172],

(6.26)

F ⁡ [ ρ ↑ , ρ ↓ ] = T s ⁡ [ ρ ↑ , ρ ↓ ] + J ⁡ [ ρ ↑ + ρ ↓ ] + E xc ⁡ [ ρ ↑ , ρ ↓ ]

In Eq. (6.26), we have explicitly shown the spin labels, and one must keep in mind that KS spin-orbitals themselves are functionals of the density, i.e., ϕ iσ ⁡ [ ρ ] ( r ) , as was established by the first HK theorem.

To evaluate the Euler-Lagrange equation [Eq. (6.24)], we will use the chain rule for functional derivatives. If, for example, we consider some arbitrary functional F[y], where y [ w ] ( r ) itself is a functional of the function w(r), then we can evaluate the functional derivative with respect to y(r) by using

(6.27)

δF δy ⁢ ( r ) = ∫ δF δw ⁢ ( r ′ ) ⁢ δw ( r ′ ) δy ⁢ ( r ) ⁢ ⅆ r ′

Equation (6.27) will be valid if w(r) is a function of y(r). It can also be obtained from the chain rule by using [Parr, p.250]

(6.28)

δF δw ⁢ ( r ) = ∫ δF δy ⁢ ( r ′ ) ⁢ δy ⁢ ( r ′ ) δw ⁢ ( r ) ⁢ ⅆ r ′

and

(6.29)

∫ δy ( r ) δw ⁢ ( r ′ ) ⁢ δw ⁢ ( r ′ ) δy ⁢ ( r ⁢ ″ ) ⁢ ⅆ r ′ = δ ⁢ ( r - r ⁢ ″ )

If the preceding discussion was a bit confusing, then we can illustrate with an example, which happens to be the case were interested in. In other words, let y = ϕ iσ * and w = ρ σ , then Eq. (6.27) gives us

(6.30)

δF δϕ iσ * ⁢ ( r ) = ∫ δF δρ σ ⁢ ( r ′ ) ⁢ δρ σ ⁢ ( r ′ ) δϕ iσ * ⁢ ( r ) ⁢ ⅆ r ′ = ∫ δF δρ σ ⁢ ( r ′ ) ⁢ ϕ iσ ⁢ ( r ′ ) ⁢ δ ⁢ ( r ′ - r ) ⁢ ⅆ r ′ = δF δρ σ ⁢ ( r ) ⁢ ϕ iσ ⁢ ( r )

We can use Eq. (6.30) to evaluate the various contributions to the Euler-Lagrange equation from Eqs. (6.11) and (6.12). Starting with the system dependent term containing the external potential, we find

(6.31)

δ δϕ iσ * ⁢ ( r ) ⁢ ∫ ρ ⁢ ( r ) ⁢ v ⁢ ( r ) ⁢ ⅆ r = ϕ iσ ⁢ ( r ) ⁢ δ δρ σ ⁢ ( r ) ⁢ ∫ ρ ⁢ ( r ) ⁢ v ⁢ ( r ) ⁢ ⅆ r = ϕ iσ ⁢ ( r ) ⁢ v ⁢ ( r )

From the system independent functional, F ⁡ [ ρ ↑ , ρ ↓ ] , we have a contribution from the kinetic energy T s ⁡ [ ρ ↑ , ρ ↓ ] , which is given by

(6.32)

δT s ⁡ [ ρ ↑ , ρ ↓ ] δϕ iσ * ⁢ ( r ) = - 1 2 ⁢ ∇ 2 ⁢ ϕ iσ ⁢ ( r ) ,

and the contribution from the classical Coulomb energy, J ⁡ [ ρ ] , is

(6.33)

δJ ⁡ [ ρ ] δϕ iσ * ⁢ ( r ) = ϕ iσ ⁢ ( r ) ⁢ δJ ⁡ [ ρ ] δρ σ ⁢ ( r ) = ϕ iσ ⁢ ( r ) ⁢ ∫ ⅆ r′ ⁢ ρ ⁢ ( r′ ) | r - r′ |

Finally, the contribution from exchange-correlation functional is given by

(6.34)

δE xc ⁡ [ ρ ↑ , ρ ↓ ] δϕ iσ * ⁢ ( r ) = ϕ iσ ⁢ ( r ) ⁢ δE xc ⁡ [ ρ ↑ , ρ ↓ ] δρ σ ⁢ ( r ) = ϕ iσ ⁢ ( r ) ⁢ v xc σ ⁢ ( r )

where we have defined the exchange-correlation potential as

(6.35)

v xc σ ⁢ ( r ) ≡ δE xc ⁡ [ ρ ↑ , ρ ↓ ] δρ σ ⁢ ( r )

and, as before [i.e., Eq. (5.33)], the contribution from the constraint is

(6.36)

δ ⁢ S ij δϕ iσ * ⁢ ( r ) = ϕ jσ ⁢ ( r )

Collecting the terms together gives us the Kohn-Sham orbital equations

(6.37)

h ^ σ ⁢ ϕ iσ ⁢ ( r ) = ∑ j = 1 N σ ⁢ ε ij σ ⁢ ϕ jσ ⁢ ( r )

where the Kohn-Sham operator is defined as

(6.38)

h ^ σ = - 1 2 ⁢ ∇ 2 ⁢ + v eff σ ⁢ ( r )

and the definition of the effective potential is

(6.39)

v eff σ ⁢ ( r ) ≡ v ⁢ ( r ) + ∫ ⅆ r′ ρ ( r′ ) | r - r′ | + v xc σ ⁢ ( r )

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

(6.40)

h ^ σ ⁢ ϕ iσ ⁢ ( r ) = ε iσ ⁢ ϕ iσ ⁢ ( r )

The Kohn-Sham orbital energies can be obtained from the expectation value of the single-particle Hamiltonian operator

(6.41)

ε iσ ≡ 〈 ϕ iσ | h ^ σ | ϕ iσ 〉 = 〈 ϕ iσ | - 1 2 ⁢ ∇ 2 ⁢ + v ⁢ ( r ) + ∫ ⅆ r′ ρ ( r′ ) | r - r′ | + v xc σ ⁢ ( r ) | ϕ iσ 〉

To understand the relationship between the Kohn-Sham orbital energies and the total energy we note that the sum of the orbital energies gives us

(6.42)

∑ σ = ↑ , ↓ ⁢ ∑ i = 1 N σ ε iσ ⁢ = T s ⁡ [ ρ ↑ , ρ ↓ ] + ∫ ρ ⁢ ( r ) ⁢ v ⁢ ( r ) ⁢ ⅆ r + 2 ⁢ J ⁡ [ ρ ] + ∑ σ = ↑ , ↓ ⁢ ∫ ρ σ ⁢ ( r ) ⁢ v xc σ ⁢ ( r ) ⁢ ⅆ r

By comparing Eq. (6.42) with Eqs. (6.11) and (6.12), it can be seen that in a calculation of the total energy we need to correct for the overcounting of the classical Coulomb energy, J ⁡ [ ρ ] , and the additional contribution from the exchange-correlation potential v xc σ . The total energy [Eq. (6.11)] is then related to the orbital energies by

(6.43)

E = ∑ σ = ↑ , ↓ ⁢ ∑ i = 1 N σ ε iσ ⁢ - 1 2 ⁢ ∫ ∫ ⅆ r ⁢ ⅆ r′ ⁢ ρ ⁢ ( r ) ⁢ ρ ⁢ ( r′ ) | r - r′ | + E xc ⁡ [ ρ ↑ , ρ ↓ ] - ∑ σ = ↑ , ↓ ⁢ ∫ ρ σ ⁢ ( r ) ⁢ v xc σ ⁢ ( r ) ⁢ ⅆ r

or for closed-shell restricted spin-orbitals we have

(6.44)

E = ⁢ 2 ∑ i = 1 N / 2 ⁢ ε i - 1 2 ⁢ ∫ ∫ ⅆ r ⁢ ⅆ r′ ⁢ ρ ⁢ ( r ) ⁢ ρ ⁢ ( r′ ) | r - r′ | + E xc ⁡ [ ρ ] - ∫ ρ ⁢ ( r ) ⁢ v xc ⁢ ( r ) ⁢ ⅆ r

which is the form used in our first example program, freeMQM-1.cpp.

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

6.3. A few practical details in implementing the Kohn-Sham SCF procedure.

One method of solving the Kohn-Sham equations is to expand the spatial part of the orbitals ϕ i using a complete set of basis functions { φ μ ⁢ ( r ) }

(6.45)

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

where the function φ μ ( r ) represents a member of the basis. Substituting the expansion for the orbitals into the Kohn-Sham equations will lead to a matrix representation, analogous to the Pople-Nesbet equations,

(6.46)

F α ⁢ C α = SC α ⁢ ε α

(6.47)

F β ⁢ C β = SC β ⁢ ε β

where σ = α = ↑ ≐ + 1 2 labels the spin-up electrons and σ = β = ↓ ≐ - 1 2 labels the spin-down electrons.

Equations (6.46) and (6.47) represent the unrestricted Kohn-Sham (UKS) matrix formalism, analogous to the Pople-Nesbet equations. As before, in the case of our Hartree-Fock examples, we will consider closed shell systems and solve the Kohn-Sham equivalent of the Roothaan-Hall equations in the restricted Kohn-Sham (RKS) matrix formalism

(6.48)

F ⁢ C = SC ⁢ ε

where F is the Fock matrix (or Kohn-Sham matrix), C is the matrix of molecular orbital coefficients, S is the overlap matrix, and ε is the diagonal matrix of orbital energies. In the UKS formalism, the Fock matrix F α becomes

(6.49)

F α = H core + J + F XCα

and in the RKS formalism it can be written as

(6.50)

F = H core + J + F XC

where the exchange-correlation part of the Fock matrix in the RKS formalism, F XC , is equal to exchange-correlation part of the Fock matrix in the UKS formalism, F XCα , with ρ α = ρ / 2 .

The exchange-correlation contribution to the Fock matrix is related to the exchange-correlation potential by

(6.51)

F μν xcα = 〈 μ | v xc α | ν 〉

where v xc α is obtained from the functional derivative of the exchange-correlation energy functional

(6.52)

v xc α = δ ⁢ E xc δ ⁢ ρ α

In the local-density approximation (LDA) or the generalized gradient approximation (GGA), a general expression for the exchange-correlation energy functional can be written as

(6.53)

E xc ⁡ [ ρ ] = ∫ f ⁢ ( ρ α , ρ β , γ αα , γ αβ , γ ββ ) ⁢ ⅆ r

where γ αα ≡ | ∇ ρ α ⁢ | 2 , γ ββ ≡ | ∇ ρ β ⁢ | 2 , and γ αβ ≡ ∇ ρ α · ∇ ρ β . Applying the functional derivative from Eq. (5.28), together with the chain rule,

(6.54)

∂ f ∂ ∇ ρ α = ∂ γ αα ∂ ∇ ρ α ⁢ ∂ f ∂ γ αα + ∂ γ αβ ∂ ∇ ρ α ⁢ ∂ f ∂ γ αβ

leads to

(6.55)

v xc α = ∂ f ∂ ρ α - 2 ⁢ ∇ · ( ∂ f ∂ γ αα ⁢ ∇ ρ α ) - ∇ · ( ∂ f ∂ γ αβ ⁢ ∇ ρ β )

At this point, we can use Eq. (6.55) together with Green's Theorem in the following form

(6.56)

∫ S ⁢ u ⁢ ∇ w · n ⁢ ⅆ a = ∫ V ⁢ u ⁢ ∇ 2 ⁢ w ⁢ ⅆ v + ∫ V ⁢ ∇ u · ∇ w ⁢ ⅆ v ,

to calculate the Fock matrix Eq. (6.51), where n is an unit normal vector outward pointing from the surface S. Then, with the following substitutions

(6.57)

u = φ μ ⁢ φ ν ⁢ ∂ f ∂ γ αα ,and w = ρ α

we can obtain the following expression [84][86]

(6.58)

F μν XCα = ∫ { ∂ f ∂ ρ α ⁢ φ μ ⁢ φ ν + [ 2 ⁢ ( ∂ f ∂ γ αα ⁢ ∇ ρ α ) + ( ∂ f ∂ γ αβ ⁢ ∇ ρ β ) ] · ∇ ( φ μ ⁢ φ ν ) } ⁢ ⅆ r

where the surface integral contributes nothing when integrated over all space. Equation (6.58) has a form which readily lends itself to numerical integration; it is used by the integration routines to calculate the exchange-correlation (XC) functional contribution to the Fock matrix, which has the variable name FMatrix[I][J] in freeMQM-1.cpp.

To perform a density-functional calculation we will need to make a couple changes to the procedure used previously for the Hartree-Fock method. The first change that needs to be made is to replace the exchange contribution to the Hartree-Fock equations, K μν , with a contribution from the density functional. The second change is that the exchange-correlation energy must be calculated during each iteration of the SCF procedure and added to the total energy. As pointed out before in Eq. (6.43), it is also necessary to subtract the additional contribution from the exchange-correlation potential v xc , to calculate the total energy. With these changes, the rest of the SCF procedure used in the Hartree-Fock method will be the same.

For a Hartree-Fock calculation, the Fock matrix is calculated by the function CalcFockMatrix, which is called from inside the self-consistent loop of the example program freeMQM-1.cpp. In our DFT calculation, CalcFockMatrix still calculates

(6.59)

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

but the contribution to the Fock matrix from the density-functional, F μν XCα , is calculated by xyz_Integrate, or by R_Integrate in cases of spherical symmetry, and then added to F. Both the exchange-correlation (XC) functional, E xc ⁡ [ ρ ] , and the XC contribution to the Fock matrix, F μν xcα ( F μν xcβ ) , are calculated within the integration routines, xyz_Integrate and R_Integrate, which in turn make function calls to routines in functionals.cpp. The program file functionals.cpp contains formulas for various exchange and correlation functionals and their derivatives[50]. By default, the program freeMQM-1.cpp uses the functional BLYP, but there are several other choices, and functionals.cpp can easily be altered to include other functionals. BLYP uses the Becke exchange functional B88 [87] and the LYP [96] correlation functional.

In general, it can be seen from the form of the function f ⁢ ( ρ α , ρ β , γ αα , γ αβ , γ ββ ) that the values of ρ α , ρ β , γ αα , γ ββ , and γ αβ , must be passed to functionals.cpp to receive back the value of the function f and its derivatives. For spherically symmetric problems, such as the ground state of the helium atom, the values of ρ α , ρ β , γ αα , γ ββ , and γ αβ , are calculated by the function R_Integrate, and in all other cases they are calculated by xyz_Integrate. These values are then passed to the program file functionals.cpp for each set of spatial coordinates, where the variables in the program are defined as follows:

(6.60)

gaa ≡ γ αα = | ∇ ρ α ⁢ | 2 ⁢ , gbb ≡ γ ββ = | ∇ ρ β ⁢ | 2 ⁢ , and gab ≡ γ αβ = ∇ ρ α · ∇ ρ β

The routines in functionals.cpp return the values of the function f and the derivatives,

(6.61)

(6.62)

(6.63)

(6.64)

(6.65)

dfdgab ≡ ∂ f ∂ γ αβ ,

from which we can construct the functional derivative and the exchange-correlation part of the Fock matrix [Eq. (6.58)]. Further simplifications can be made if we assume that we are dealing with a closed shell system, in which case the densities of spin-up ( σ = α = ↑ ≐ + 1 2 ) and spin-down electrons ( σ = β = ↓ ≐ - 1 2 )

(6.66)

ρ α ⁢ ( r ) = ∑ i N α ⁢ | ϕ iα ⁢ ( r ) | 2 ,and ρ β ⁢ ( r ) = ∑ i N β ⁢ | ϕ iβ ⁢ ( r ) | 2

are both equal to ρ / 2 , where

(6.67)

ρ ⁢ ( r ) = ρ α ⁢ ( r ) + ρ β ⁢ ( r )

is the total electron density. This simplification allows us to drop the spin labels and write the total density as a sum over doubly occupied spatial orbitals

(6.68)

ρ ⁢ ( r ) = 2 ⁢ ∑ i N / 2 ⁢ | ϕ i ⁢ ( r ) | 2 .

As before in the Hartree-Fock problem, the spatial part of the Kohn-Sham orbitals can be written as a sum over our K Gaussian basis functions [Eq. (5.45)]

(6.69)

ϕ i ⁢ ( r ) = ∑ μ K ⁢ C μi ⁢ φ μ ⁢ ( r )

When applied to density-functional theory (DFT) this is referred to as the restricted Kohn-Sham (RKS) model, as opposed to the unrestricted Kohn-Sham (UKS) formalism, which we might discuss in more detail in a future revision of this document. The file functionals.cpp contains formulas for both the RKS and UKS version of the functionals and their derivatives.

In the closed-shell RKS formalism the gradient of the total electron density is given by

(6.70)

∇ ρ ⁢ ( r ) = 4 ⁢ ∑ N / 2 i ⁢ ϕ i ⁢ ( r ) ⁢ ∇ ϕ i ⁢ ( r )

and ∇ ϕ i ⁢ ( r ) is equal to the linear superposition of the gradients of the basis functions

(6.71)

∇ ϕ i ⁢ ( r ) = ∑ μ K ⁢ C μi ⁢ ∇ φ μ ⁢ ( r )

To include ∇ ρ in the integration routines, R_Integrate and xyz_Integrate, we will need to take the gradients of the Gaussian basis functions by hand and add the necessary code to the program. If we recall that the basis sets are constructed from contracted Gaussian functions (CGF) [1]

(6.72)

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

and, as before, φ p ⁢ ( r ) are the primitive Gaussian functions, then the gradient of the spatial part of the Kohn-Sham orbitals in terms of primitive Gaussians can be written as

(6.73)

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

Further simplifications can be made for both the helium atom and H₂, since they have a single ground-state orbital occupied by one spin-up electron and one spin-down electron. Therefore the total electron density is simply ρ ⁢ ( r ) = 2 ⁢ | ϕ ⁢ ( r ) | 2 , and the gradient of the total charge density is ∇ ρ = 4 ⁢ ϕ ⁢ ∇ ϕ . In the integration routines, R_Integrate and xyz_Integrate, the C/C++ code representation of these these quantities are rho ≡ ρ , psi ≡ ϕ , delpsi ≡ ∇ ϕ , and phi ≡ φ p . In the RKS approach to the helium atom and H₂, the values of gaa, gbb, gab, are equal to g/4, where

(6.74)

g ≡ γ ⁢ ( r ) = | ∇ ρ ⁢ ( r ) | 2 = 8 ⁢ ρ ⁢ ( r ) ⁢ | ∇ ϕ ⁢ ( r ) | 2

As before [i.e., Eq. (5.66)], the most general form of the primitive Gaussian functions [113] [114][Szabo, p. 181;Levine,p.462], φ p ⁢ ( r ) , we will consider are described by

(6.75)

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

It follows that the gradient of the primitive Gaussian functions in the spherically symmetric case ( L = l pμ + m pμ + n pμ = 0 ) is simply given by

(6.76)

∇ φ p = ∇ g = r ^ ∂ g ∂ r = - 2 ⁢ α pμ ⁢ r g .

In the more general non-spherically symmetric case we will need to evaluate

(6.77)

∇ φ p = ∇ g = ( ∂ g ∂ x , ∂ g ∂ y , ∂ g ∂ z )

where the partial derivatives of the Gaussian primitive functions take the more general form

(6.78)

∂ φ p ∂ x = ∂ g ∂ x = l pμ ⁢ ( x - x A ) l pμ - 1 ⁢ ( y - y A ) m pμ ⁢ ( z - z A ) n pμ ⁢ e - α pμ ⁢ | r - r A | 2 - 2 ⁢ α pμ ⁢ ( x - x A ) ⁢ g

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

6.4. Numerically integrating the exchange-correlation terms.

The integration methods used by the functions R_Integrate and xyz_Integrate are summarized in [85]. The function R_Integrate uses the Euler-Maclaurin quadrature formula of the form

(6.79)

∫ 0 ∞ ⁢ r 2 ⁢ G ⁢ ( r ) ⁢ ⅆ r ≈ ∑ i = 1 N r ⁢ w i r ⁢ G ⁢ ( r i )

where N r is the number of radial points and the radial weights w i r and radial coordinates r i are given by

(6.80)

w i r = 2 ⁢ R 3 ⁢ ( N r + 1 ) ⁢ i 5 ⁢ ( N r + 1 - i ) - 7

and

(6.81)

r i = Ri 2 ⁢ ( N r + 1 - i ) - 2

respectively. R is referred to as the "atomic radius" [85] and takes the value of 1.0 for the hydrogen atom and 0.5882 for He.

The 3-dimensional quadrature formulas have the following form

(6.82)

I = ∫ F ⁢ ( r ) ⁢ ⅆ r = ∫ 0 ∞ ∫ 0 π ⁢ ∫ 0 2 π ⁢ r 2 ⁢ sinϑ ⁢ F ⁢ ( r , ϑ , φ ) ⁢ ⅆ r ⁢ ⅆ ϑ ⁢ ⅆ φ ≈ 4 π ∑ i = 1 N r ⁢ w i r ⁢ ∑ i = 1 N Ω ⁢ w j Ω ⁢ F ⁢ ( r i , ϑ j , φ j )

where for simplicity, the radial and angular integrations are considered separately. The angular integrations are performed using the using the method of Gauss quadratures on a sphere developed by Lebedev (see Refs. [105] to [111]), and as before the radial integrations use the Euler-Maclaurin formula [85]. The Euler-Maclaurin-Lebedev (EML) atomic grid is sometimes referred to simply as ( N r , N Ω ) [85] as opposed an alternative grid by Murray, Handy and Laming (MHL), which involves three separate integrations and has been referred to as ( N r , N θ , N φ ) , but is not discussed further here. The function xyz_Integrate implements the Lebedev quadratures by making function calls to Lebedev-Laikov.cpp, which returns the atomic grid coordinates and weights using the method of Gauss quadratures on a sphere [105]. The original source code for the quadrature formulas (written in C) was made available by Laikov [105] and has been translated into Fortran and can be downloaded at http://www.ccl.net/cca/software/SOURCES/FORTRAN/Lebedev-Laikov-Grids/. Lebedev-Laikov.cpp includes routines to produce Lebedev angular grids from 6-point grids to 5810-point grids. Although thirty-two of these routines are available in Lebedev-Laikov.cpp, currently freeMQM-1.cpp only uses the "standard" grid (50,194).

The quadrature formulas for the sphere Lebedev and Laikov presented in Ref. [105] are valid for all polynomials x^ky^lz^m, where k + l + m ≤ 131 . The grid parameters listed should give a relative accuracy on the order of 2 × 10 - 14 , which is limited by the number of significant digits used in the grid parameters and by computer accuracy. The integrals we are interested in, however, involve several Gaussian functions that will not all be centered at a single atom in the case of polyatomic molecules. The Gaussians can be Taylor expanded to arbitrary order and if the function cannot be approximated by an expansion that is spherically symmetric beyond a finite order then the Lebedev quadrature formula will break down unless an alternative method can be employed to ensure that the integrand satisfies k + l + m ≤ 131 .

Becke developed a multicenter integration scheme for polyatomic molecules, which we can apply to EML quadratures [88]. In this approach the integral over the polyatomic system is split into individual atomic centered integrals as follows

(6.83)

I = ∑ A ⁢ I A

where the atomic centered integrals I_A can be solved numerically using the EML integration formulas. If the original molecular integral is

(6.84)

I = ∫ F ⁢ ( r ) ⁢ ⅆ r

then the I_A can be written as

(6.85)

I A = ∫ F A ⁢ ( r ) ⁢ ⅆ r

where

(6.86)

F A ⁢ ( r ) = w A ⁢ ( r ) ⁢ F ⁢ ( r )

The problem reduces to finding the weight functions w A ⁢ ( r ) that are unity near the nucleus A, but vanish near all other nuclei in a well-behaved manner. For the integrals to accurately reproduce Eq. (6.84) and give physically meaningful results

(6.87)

∑ A ⁢ w A ⁢ ( r ) = 1

so that

(6.88)

∑ A ⁢ F A ⁢ ( r ) = F ⁢ ( r )

Following the initial proposal of Becke [88] the weight functions w A ⁢ ( r ) are constructed from "cell functions" P A ⁢ ( r ) as follows

(6.89)

w A ⁢ ( r ) = P A ⁢ ( r ) / ∑ C ⁢ P C ⁢ ( r )

The "cell functions" are constructed to be unity inside a cell surrounding a given atom and decay rapidly outside the cell; Becke expanded P A ⁢ ( r ) as a product of step functions s(μ)

(6.90)

P A ⁢ ( r ) = ∏ C ≠ A ⁢ s ⁢ ( μ AC )

where μ AC is the elliptical coordinate

(6.91)

μ AC = r A - r C R AC

and r_A and r_C are the distances to nuclei A and C, respectively. R_AC is the internuclear distance between nuclei A and C. Becke proposed [88] a smoothly varying analog to the step function

(6.92)

s ⁢ ( μ AC ) = { 1 , - 1 ≤ μ AC ≤ 0 0 , 0 < μ AC ≤ + 1

and found that the simplest possible function that could satisfy the necessary constraints was

(6.93)

s ⁢ ( μ ) = 1 2 ⁡ [ 1 - p ⁢ ( μ ) ]

where

(6.94)

p ⁢ ( μ ) = 3 2 ⁢ μ - 1 2 ⁢ μ 3

This does not give a sufficiently rapid cutoff required of the step function; however, three iterations of the function p(μ) has been found to give excellent results [88]

(6.95)

s ⁢ ( μ ) = 1 2 ⁡ [ 1 - p ⁢ ( p ⁢ ( p ⁢ ( μ ) ) ) ]

This atomic-centered grid method was one of the first implemented for density-functional theory. Improved descendants are no doubt still used in modern ab initio software packages [Koch,p.107;Treutler and Ahlrichs, 1995][101]. In upcoming sections we will discuss recent developments and alternative methods.

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

6.5. Choosing an exchange/correlation functional.

Much time and effort has been spent trying to obtain exchange-correlation functionals which will give accurate results. Many useful review articles and Internet resources describe some of the most popular functionals. Mark Casida has a relatively exhaustive and up-to-date tabulated listing of most of the density functionals with links to further information. Earlier work includes that of Johnson et al. [86], which provides a discussion of some of the important "2nd generation" functionals, including the Becke (B88) exchange functional [87], the Vosko, Wilkin, and Nusair (VWN) correlation functional [72], and the Lee, Yang and Parr (LYP) correlation functional [96]. This article also includes an appendix with the relevant functional derivatives. The CLRC Quantum Chemistry Group density-functional Repository [50] has a nice summary of many of the functionals used in research and makes available the Fortran source code that calculates the values of functionals and the necessary derivatives. If one is interested in evaluating other functional derivatives, the book by Parr and Yang [86] includes an appendix on how to perform functional derivatives and the article by Jemmer and Knowles [102] describes a Mathematica program that calculates functional derivatives and gives the Becke88 functional and LYP functional as examples. The program can be obtained online at http://www.tc.bham.ac.uk/~peterk/papers/fderiv/ .

Currently the default functional used by freeMQM-1.cpp is BLYP, which consists of the Becke88 (B88) exchange functional and the LYP correlation functional. Just as an example, the function that contributes to the B88 exchange functional can be written as

(6.96)

f B88 X = f LDA X + Δf B88 X

where

(6.97)

f LDA X = - 3 4 ⁢ ( 6 π ) 1 3 ⁢ ( ρ α 4 3 + ρ β 4 3 )

is the Slater local-density approximation (LDA) exchange contribution and the Becke gradient correction is

(6.98)

Δf B88 X = - β ⁢ ∑ σ = ↑ , ↓ ⁢ ρ σ 4 3 ⁢ x σ 2 1 + 6 ⁢ β ⁢ x σ ⁢ sinh - 1 ⁢ x σ x σ = γ σ ρ σ 4 / 3 γ σ = γ σσ β = 0.0042

The relevant derivatives of f B88 X are given by [86][50]

(6.99)

∂ f B88 X ∂ ρ σ = 4 3 ⁢ ρ σ 1 3 ⁡ [ G ⁢ ( x σ ) - x σ ⁢ G ' ⁢ ( x σ ) ] ∂ f B88 X ∂ γ σσ = 1 2 ⁢ G ' ⁢ ( x σ ) γ σ G ⁢ ( x ) = A - β ⁢ x 2 1 + 6 ⁢ β ⁢ x ⁢ sinh - 1 ⁢ ( x ) G ' ⁢ ( x ) = 6 ⁢ β 2 ⁢ x 2 ⁢ ( x / 1 + x 2 - sinh - 1 ⁢ x ) - 2 ⁢ β ⁢ x ( 1 + 6 ⁢ β ⁢ x ⁢ sinh - 1 ⁢ x ) 2 where A = - 3 4 ⁢ ( 6 π ) 1 3

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

6.6. More program details.

The function xyz_Integrate only works with atoms and diatomic molecules at the moment, which is just a special case of the procedure outlined by Becke. Currently the step function s(μ) is calculated with the function StepS.

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

6.7. Density-functional theory applied to the helium atom.

From the self-consistent-field (SCF) method we obtain the Kohn-Sham (KS) orbitals, electron density, total ground state energy ( E total ), Kohn-Sham orbital energies ( ε i ), kinetic energy ( E Kin ), Hartree energy ( E Hartree ), electron-nucleus attraction energy ( E e-n ), exchange energy ( E x ), correlation energy ( E c ), and the exchange-correlation energy ( E xc ). Gill reported [85] a total energy of E = -2.897845 a.u. for the helium atom using BLYP/6-31G* and the "standard-grid", which is consistent with the results shown in Table 2. Note in Table 2 that the largest discrepancy between the results of the density-functional calculation and experiment are the Kohn-Sham orbital energies (e.g. ε N ). The negative of the highest occupied Kohn-Sham orbital energy ( - ε N ) should be exactly equal to the ionization energy, but as can be seen from Table 2 the calculated values are very inaccurate. The reason for the large discrepancy is that the exchange-correlation potential v xc , calculated from most of standard functionals, is a poor approximation to the exact exchange-correlation potential [104][174] [216][178][224] [Koch,p.88]. Furthermore, additional errors can be introduced in the linear combination of atomic orbitals (LCAO) method if the basis functions are not very accurate approximations to the exact atomic orbitals, as in the case of Gaussian-type orbitals (GTOs), or if too small of a basis is used in the expansion of the orbitals [185][216].

Table 6.1. Density-functional theory results for the helium atom. The energy values are in atomic units (hartrees). E total , ε N , E Kin , E e-n , E Hartree , E x , E c , and E xc , are the total ground state energy, highest occupied Kohn-Sham orbital energy, kinetic energy, Hartree energy, electron-nucleus attraction energy, exchange energy, correlation energy, and the exchange-correlation energy, respectively. The basis used was 6-31G* and the program required approximately 1 second to complete the calculation ^[c].

Exchange/Correlation Functional E total ε N E Kin E e - n E Hartree E x E c E xc Exact[a] - 2.903724377 - 0.903570389[b] 2.867082 - 6.753267 2.049137 - 1.024568 - 0.042107 - 1.066676 BLYP[c] - 2.89892 - 0.56913 2.85734 - 6.73979 2.05458 - 1.02612 - 0.0438572 - 1.06997 Slater-LDA/None[c] - 2.71464 - 0.499324 2.72349 - 6.57878 2.00367 - 0.863027 0 - 0.863027 B88/VWNRPA (BVWNRPA)[c] - 3.00461 - 0.613979 2.89172 - 6.78057 2.06763 - 1.03258 - 0.150815 - 1.1834 B88/PZ81[c] - 2.9531 - 0.588046 2.8841 - 6.77155 2.06474 - 1.03115 - 0.0992394 - 1.13039 PBE/PZ81[c] - 2.93909 - 0.539544 2.72941 - 6.58597 2.00592 - 0.990352 - 0.0980981 - 1.08845

^[a]From Umrigar and Gonze [104].

^[b]Experimental results from the NIST Atomic Spectra Database (ASD).

^[c]Computed using freeMQM-1 on a PII 266 MHz with 224 RAM (with ~ 50% of the resources available). ______________________________________________________________________________________________________________________

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

6.8. Density-functional theory applied to the hydrogen molecule.

Our first example of density-functional theory applied to a molecular system is again the hydrogen molecule. Gill reported [85] a total energy of E = -1.165182 a.u. for H₂ using BLYP/6-31G* and the "standard-grid". Note again in Table 3 the large discrepancy between the the Kohn-Sham orbital energy (e.g. ε N ) and the experimental value due to the deficiencies in the approximate exchange-correlation potential [104][174] [216][178][224][Koch,p.88].

Table 6.2. Density-functional results for the hydrogen molecule. The energy values are in atomic units (hartrees). E total , ε N , E Kin , E e-n , E Hartree , E x , E c , and E xc , are the total ground state energy, highest occupied Kohn-Sham orbital energy, kinetic energy, Hartree energy, electron-nucleus attraction energy, exchange energy, correlation energy, and the exchange-correlation energy, respectively. For the basis 6-31G* the program required approximately 3 minutes and 20 seconds to complete the calculation and the basis 6-31G** required an additional two minutes^[c].

Exchange/Correlation Functional Basis Set Bond distance E total ε N E Kin E e - n E Hartree E x E c E xc Exact/Experiment[a] - 1.401[a] - 1.174474983 [b] - 0.584[a] - - - - - - BLYP[c] 6-31G * 1.414 - 1.16526 - 0.387119 1.14369 - 3.63127 1.30805 - 0.654919 - 0.0380184 - 0.692938 BLYP[c] 6-31G ** 1.412 - 1.16791 - 0.384742 1.14338 - 3.64091 1.31719 - 0.657506 - 0.0382732 - 0.695779 Becke B88/None[c] 6-31G ** 1.424 - 1.12967 - 0.348525 1.13224 - 3.62006 1.30991 - 0.654014 0 - 0.654014 Slater-LDA/None[c] 6-31G * 1.476 - 1.03961 - 0.321777 1.05413 - 3.48348 1.25423 - 0.542 0 - 0.542 B88/VWNRPA (BVWNRPA)[c] 6-31G * 1.395 - 1.25786191 - 0.42053 1.1745 - 3.67992 1.32621 - 0.663969 - 0.131016 - 0.794986 B88/PZ81[c] 6-31G * 1.400 - 1.20988 - 0.395854 1.1676 - 3.66905 1.32214 - 0.661946 - 0.0829109 - 0.744857 PBE/PZ81[c] 6-31G * 1.404 - 1.19994 - 0.365907 1.09587 - 3.5768 1.28446 - 0.633946 - 0.081775 - 0.715721

^[a]The experimental value of the bond distance can be found on page 201 of Szabo and the experimental value of ε N is on page 194.

^[b]The exact value of E_total can be found in Lowe, p.232.

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

6.9. More things to do.

Having two integration functions is redundant. The function R_Integrate could be eliminated and be made a special case of xyz_Integrate.
Only LD0194 is currently used in xyz_Integrate, but all the sets of quadrature weights and coordinates described by Lebedev and Laikov in Ref. [105] are available in Lebedev-Laikov.cpp. The function xyz_Integrate could easily be generalized to enable it to access all the functions in Lebedev-Laikov.cpp.
Currently the weight functions w_A(r) are only calculated for diatomic molecules.

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


Chapter 5. Including Electron-Electron Interactions and the Hartree-Fock Approximation.		Chapter 7. Orbital-Dependent Functionals and Exact-Exchange Methods.