## 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   [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)  that for a N-particle system all of the electronic ground state properties are determined by the electron density  $\rho \left(\mathbf{r}\right)$, and that if  $\rho \left(\mathbf{r}\right)$  is the true ground state density, then the total energy functional,  $E\left[\rho \left(\mathbf{r}\right)\right]$,  is a minimum. To prove the Hohenberg-Kohn theorems, we will start by considering the same Hamiltonian as in Eq. (5.3), rewritten as

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

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\left[\Psi \right]\text{,}$  written as a functional of the normalized N-electron wavefunction,  $\Psi \left({\mathbf{x}}_{1},\dots ,{\mathbf{x}}_{N}\right)$,  is given by

where we have used the fact that

and

To show that v(r) is also an unique functional of  $\rho \left(\mathbf{r}\right)$, we first assume that the opposite is true―that two different external potentials v(r) and v'(r) can have the same electron distribution  $\rho \left(\mathbf{r}\right)$―and then show that this leads to a contradiction. The first part of the Hamiltonian,  $\stackrel{^}{T}+{\stackrel{^}{V}}_{\mathrm{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,  $\stackrel{^}{H}\text{,\hspace{0.17em}}\Psi \text{,\hspace{0.17em}}{E}_{o}$, and  $\stackrel{^}{H}\prime \text{,\hspace{0.17em}}\Psi \prime \text{,\hspace{0.17em}}{E}_{o}\prime$, respectively. If we let  $\Psi \prime$  serve as a trial wavefunction for the Hamiltonian  $\stackrel{^}{H}$, then the Rayleigh-Ritz minimal principle for the ground state,

tells us that

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

Thus, the ground state density  $\rho \left(\mathbf{r}\right)$  uniquely determines the external potential v(r), to within a trivial additive constant. Since v(r) fixes the Hamiltonian  $\stackrel{^}{H}$,  it follows that the full N-particle ground state wavefunction, and all other properties determined by  $\stackrel{^}{H}$, are implicitly determined by the ground state density  $\rho \left(\mathbf{r}\right)$.

For the second HK theorem (the HK variational principle) it was shown that if  $\rho \left(\mathbf{r}\right)$  is the true ground state density then the total energy functional  $E\left[\rho \left(\mathbf{r}\right)\right]$  is a minimum. In other words, the HK variational principle states that

where  ${\rho }_{o}\left(\mathbf{r}\right)$  and  ${E}_{o}$  are the exact ground state density and energy, respectively, and the equality holds if  $\rho \left(\mathbf{r}\right)={\rho }_{o}\left(\mathbf{r}\right)$. To prove this, HK defined for a given external potential, v(r), the following energy functional

Equation (6.11) is justified by the first HK theorem, which established that  $\stackrel{^}{H}\text{,\hspace{0.17em}}\Psi$, and E are unique functionals of  $\rho \left(\mathbf{r}\right)$. The functionals F[ρ] and  ${E}_{\mathrm{xc}}\left[\rho \right]$  are defined by the following equations [Parr, p.56]

where  ${T}_{s}\left[\rho \right]$  is the kinetic energy of a noninteracting system of electrons, and the classical part of  ${V}_{\mathrm{ee}}\left[\rho \right]$  is given by

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

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

## 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  ${\varphi }_{i}$  using a complete set of basis functions { ${\phi }_{\mu }\left(\mathbf{r}\right)$ }

where  $\sigma =\alpha =↑\doteq +\frac{1}{2}$  labels the spin-up electrons and  $\sigma =\beta =↓\doteq -\frac{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

where F is the Fock matrix (or Kohn-Sham matrix), C is the matrix of molecular orbital coefficients, S is the overlap matrix, and  $\epsilon$  is the diagonal matrix of orbital energies. In the UKS formalism, the Fock matrix  ${\mathbf{F}}^{\alpha }$  becomes

and in the RKS formalism it can be written as

where the exchange-correlation part of the Fock matrix in the RKS formalism,  ${\mathbf{F}}^{\mathrm{XC}}$,  is equal to exchange-correlation part of the Fock matrix in the UKS formalism,  ${\mathbf{F}}^{\mathrm{XC\alpha }}$,  with  ${\rho }_{\alpha }=\rho /2\text{.}$

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

where  ${v}_{\mathrm{xc}}^{\alpha }$  is obtained from the functional derivative of the exchange-correlation energy functional

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

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

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

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}_{\mathrm{\mu \nu }}$,  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}_{\mathrm{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

but the contribution to the Fock matrix from the density-functional,  ${F}_{\mathrm{\mu \nu }}^{\mathrm{XC\alpha }}$,  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}_{\mathrm{xc}}\left[\rho \right]$,  and the XC contribution to the Fock matrix,  ${F}_{\mathrm{\mu \nu }}^{\mathrm{xc\alpha }}$   $\text{(}{F}_{\mathrm{\mu \nu }}^{\mathrm{xc\beta }}\text{)}$,  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. 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  and the LYP  correlation functional.

In general, it can be seen from the form of the function  $f\left({\rho }_{\alpha },{\rho }_{\beta },{\gamma }_{\mathrm{\alpha \alpha }},{\gamma }_{\mathrm{\alpha \beta }},{\gamma }_{\mathrm{\beta \beta }}\right)$  that the values of  ${\rho }_{\alpha }\text{,\hspace{0.17em}}{\rho }_{\beta }\text{,\hspace{0.17em}}{\gamma }_{\mathrm{\alpha \alpha }}\text{,\hspace{0.17em}}{\gamma }_{\mathrm{\beta \beta }}\text{,\hspace{0.17em}and\hspace{0.17em}}{\gamma }_{\mathrm{\alpha \beta }}$, 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  ${\rho }_{\alpha }\text{,\hspace{0.17em}}{\rho }_{\beta }\text{,\hspace{0.17em}}{\gamma }_{\mathrm{\alpha \alpha }}\text{,\hspace{0.17em}}{\gamma }_{\mathrm{\beta \beta }}\text{,\hspace{0.17em}and\hspace{0.17em}}{\gamma }_{\mathrm{\alpha \beta }}$, 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:

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

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 $\text{(}\sigma =\alpha =↑\doteq +\frac{1}{2}\text{)}$ and spin-down electrons $\text{(}\sigma =\beta =↓\doteq -\frac{1}{2}\text{)}$

are both equal to  $\rho /2\text{,}$  where

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

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

and  $\nabla {\varphi }_{i}\left(\mathbf{r}\right)$  is equal to the linear superposition of the gradients of the basis functions

and, as before,  ${\phi }_{p}\left(\mathbf{r}\right)$  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

As before [i.e., Eq. (5.66)], the most general form of the primitive Gaussian functions  [Szabo, p. 181;Levine,p.462],  ${\phi }_{p}\left(\mathbf{r}\right)$,  we will consider are described by

It follows that the gradient of the primitive Gaussian functions in the spherically symmetric case  $\text{(}L={l}_{\mathrm{p\mu }}+{m}_{\mathrm{p\mu }}+{n}_{\mathrm{p\mu }}=0\text{)}$  is simply given by

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

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

with similar expressions for the partial derivatives with respect to y and z.

## 6.4. Numerically integrating the exchange-correlation terms.

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

where  ${N}^{r}$  is the number of radial points and the radial weights  ${w}_{i}^{r}$  and radial coordinates  ${r}_{i}$  are given by

and

respectively. R is referred to as the "atomic radius"  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

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.  to ), and as before the radial integrations use the Euler-Maclaurin formula . The Euler-Maclaurin-Lebedev (EML) atomic grid is sometimes referred to simply as  $\text{(}{N}^{r},{N}^{\Omega }\text{)}$   as opposed an alternative grid by Murray, Handy and Laming (MHL), which involves three separate integrations and has been referred to as  $\left({N}^{r},{N}^{\theta },{N}^{\phi }\right)$, 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 . The original source code for the quadrature formulas (written in C) was made available by Laikov  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.  are valid for all polynomials xkylzm, where  $k+l+m\le 131$. The grid parameters listed should give a relative accuracy on the order of  $2\mathsf{×}{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\le 131$.

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

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

then the IA can be written as

where

so that

Following the initial proposal of Becke  the weight functions  ${w}_{A}\left(\mathbf{r}\right)$  are constructed from "cell functions"  ${P}_{A}\left(\mathbf{r}\right)$  as follows

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}\left(\mathbf{r}\right)$  as a product of step functions s(μ)

where  ${\mu }_{\mathrm{AC}}$  is the elliptical coordinate

and rA and rC are the distances to nuclei A and C, respectively. RAC is the internuclear distance between nuclei A and C. Becke proposed  a smoothly varying analog to the step function

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

where

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]. In upcoming sections we will discuss recent developments and alternative methods.

## 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. , which provides a discussion of some of the important "2nd generation" functionals, including the Becke (B88) exchange functional , the Vosko, Wilkin, and Nusair (VWN) correlation functional , and the Lee, Yang and Parr (LYP) correlation functional . This article also includes an appendix with the relevant functional derivatives. The CLRC Quantum Chemistry Group density-functional Repository  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  includes an appendix on how to perform functional derivatives and the article by Jemmer and Knowles  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

where

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

## 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}_{\mathrm{total}}$ ), Kohn-Sham orbital energies ( ${\epsilon }_{i}$ ), kinetic energy  $\text{(}{E}_{\mathrm{Kin}}\text{),}$ Hartree energy  $\text{(}{E}_{\mathrm{Hartree}}\text{),}$ electron-nucleus attraction energy ( ${E}_{\mathrm{e-n}}$ ), exchange energy ( ${E}_{x}$ ), correlation energy ( ${E}_{c}$ ), and the exchange-correlation energy ( ${E}_{\mathrm{xc}}$ ). Gill reported  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.  ${\epsilon }_{N}$ ). The negative of the highest occupied Kohn-Sham orbital energy ( $-{\epsilon }_{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}_{\mathrm{xc}}$ , calculated from most of standard functionals, is a poor approximation to the exact exchange-correlation potential   [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 .

 [a] From Umrigar and Gonze .[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). ______________________________________________________________________________________________________________________

Our first example of density-functional theory applied to a molecular system is again the hydrogen molecule. Gill reported  a total energy of E = -1.165182 a.u. for H2 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. ${\epsilon }_{N}$ ) and the experimental value due to the deficiencies in the approximate exchange-correlation potential  [Koch,p.88].
 [a] The experimental value of the bond distance can be found on page 201 of Szabo and the experimental value of  ${\epsilon }_{N}$  is on page 194.[b] The exact value of Etotal can be found in Lowe, p.232.[c] Computed using freeMQM-1 on a PII 266 MHz with 224 RAM (with ~ 50% of the resources available). __________________________________________________________________________________________________________________________________________