Bretonnet Rev Basics

AIMS Materials Science, 4(6): 1372-1405.
DOI: 10.3934/matersci.2017.6.1372
Received: 11 September 2017
Accepted: 30 November 2017
http://www.aimspress.com/journal/Materials Published: 15 December 2017
Review
Basics of the density functional theory
Jean-Louis Bretonnet *
Université de Lorraine, LCP-A2MC, EA 3469, 1 Bd. François Arago, Metz, F-57078, France
∗
Correspondence: Email: [email protected].
Abstract: The density functional theory (DFT) established itself as a well reputed way to compute the
electronic structure in most branches of chemistry and materials science. In the formulation given by
Kohn, Hohenberg and Sham in the 1960’s, the many-electron wave function is replaced by the electron
density, so that the energy is just a functional of the latter. The DFT is applied, with low computational
cost and reasonable accuracy, to predict diverse properties as binding or atomization energies, shapes
and sizes of molecules, crystal structures of solids, energy barriers to various processes, etc. In the
mid 1980s, it became an attractive alternative to the well developed wave function techniques such as
Hartree-Fock, when crucial developments in exchange-correlation energy has been taken into account,
since the Hartree-Fock method treats exchange exactly but neglects correlation.
This article is an introduction to the conceptual basis of the DFT in a language accessible for readers
entering the field of quantum chemistry and condensed-matter physics. It begins with a presentation of
the Thomas-Fermi atomic model and follows by the essentials of the density functional theory based
on the works of Hohenberg, Kohn and Sham. With a discussion of the exchange and correlation effects,
possible improvements are then presented. Lastly, mention is made of the main hybrid functionals and
of the software packages successfully applied to diverse materials of chemical, physical and biological
interest.
Keywords: electronic structure; exchange/correlation; first-principles calculation methods
1. Introduction
One of the fundamental problems in condensed-matter physics and quantum chemistry is the
theoretical study of electronic properties of systems ranging from atoms and molecules to complex
materials. Since electrons are governed by the laws of quantum mechanics, all these systems are fully
described by the Schrödinger equation. Analytic solutions of the Schrödinger equation are obtainable
for very simple systems only. However, for systems with large number of atoms, the electrostatic
repulsion between the electrons makes its numerical resolution very difficult. In that case, it is natural
1373
to consider the many-electron wave function ψ(r1 , r2 , ..., ri , ...) with appropriate approximations. The
earliest and widely used approximation was that of Hartree [1], which expresses the wave function of
system as a product of one-electron wave functions, so that the problem reduces to a one-electron
Schrödinger equation. Then, considerable improvement of the energy computation was made by
incorporating the exchange effects with the so-called Hartree-Fock approximation [2], which replaces
the product of one-electron wave functions by a linear combination of orbitals [3].
When the adiabatic approximation is made that permits to separate the movements of ions and
electron, the electronic Schrödinger equation, describing the state of the electrons interacting in the
presence of ions, reads:
 
X ~2 ! 1 X X XX 
 − ∇2
i + U ee ri − r j + U ei (|R α − r |)
i 
 ψel = Eel ψel , (1)

i
2m 2 i j,i i α

or under the compact form:

[T + Uee + Uei ] ψel = Eel ψel , (2)
where the two first terms in the Hamiltonian are the kinetic and electron-electron interaction operators,
respectively, and the third term represents the potential energy operator of the electron i in the field of
all ions.
Within the Hartree approximation, the N-electron wave function ψel is the product of one-electron
wave functions: Y
ψel (r1 , r2 , ..., rN ) = ψi (ri ) = ψ1 (r1 )ψ2 (r2 )...ψN (rN ) . (3)
i
Consequently, the electronic energy Eel , as a sum of one-electron energies, is given by:
X
Eel = Ei = hψel |H| ψel i
i
Z Z
= ψ∗1 (r1 )...ψ∗N (rN ) [T + Uei ] ψ1 (r1 )...ψN (rN )dr1 ...drN
···
e2
Z Z " #
1 XX
+ · · · ψi (ri )ψ j (r j )
∗ ∗
ψi (ri )ψ j (r j )dri dr j . (4)
2 i j,i 4πε0 ri j
The last term in the previous equation is the so-called Coulomb energy. The simplification achieved in
this term comes from the fact that certain factors are equal to unity as a consequence of normalization
of the wave function. Now, the implementation of the Hartree approximation amounts to define the
Hartree potential W H (ri ) representing the potential energy of the electron i in the field of all other
electrons:
2
e2 ψ j (r j )
Z
1X
W (ri ) =
H
dr j , (5)
2 j,i 4πε0 ri j
in order to reduce the Schrödinger equation for multiple electrons (Eq. (1)) to the sum of one-electron
Schrödinger equations:
 
 ~2 2 X
− ∇i + Uei (|Rα − ri |) + W H (ri ) ψi = Ei ψi ,

(6)
2m α
AIMS Materials Science Volume 4, Issue 6, 1372-1405.

1374
whereby the wave function of the system ψel is equal to the product of one-electron wave functions
(ψel = ψi ) and the electronic energy of the system is equal to the sum of energies of each electron
Q
Pi
(Eel = i Ei ).
Within the Hartree-Fock approximation in place of the Hartree approximation, the N-electron wave
function ψel is constructed by taking a linear combination of products of one-electron wave functions
written in the determinantal form [4]:
1 X
ψel (r1 , r2 , ..., rN ) = √ (−1) p Pψ1 (r1 )ψ2 (r2 )...ψN (rN ) , (7)
N! P
where P is the permutation operator. As a result, an extra term known as the exchange energy must
be added onto equation (4), which comes from the exchange of electronic coordinates in the expanded
determinantal wave function. The exchange energy reads:
e2
Z Z " #
1 XX
− ψi (r j )ψ j (ri )
∗ ∗
ψi (r j )ψ j (ri )dri dr j , (8)
2 i j,i 4πε0 ri j
where the double sum is over the functions of same spin. Similarly to Eq. (5), its operator in the
one-electron Schrödinger equation for electron i in the field of all other electrons is:
X Z e2 ψ∗j (r j )ψi (r j )ψ j (ri )
− dr j . (9)
j,i
4πε0 ri j ψi (ri )
Beside these considerations, a different description of many-electron systems was made by Thomas [5]
and Fermi [6], who noted that the energy of electrons moving in the potential of nuclei can be directly
calculated from the electron density n(r) rather than from the wave function. If N is the number of
electrons in the system, the electron density is defined as:
Z Z
n(r) = N · · · ψ∗el (r, r2 , ..., rN )ψel (r, r2 , ..., rN )dr2 ...drN , (10)
with the normalization condition n(r)dr = N and the limit n(r → ∞) → 0. The term n(r)dr
R
represents the elemental probability of finding any electron in the volume element dr at position r. As
a result, the electronic energy Eel may be calculated in terms of the electron density as:
Z
Eel = hψel |H| ψel i = [T + Uee + Uei ] n(r)dr . (11)
Initially the authors [5, 6] assumed that the motions of electrons were uncorrelated and that the
corresponding kinetic energy could be described by a local approximation electron density. Shortly
later, Dirac [7] proposed that the exchange effects be included by incorporating a term obtained from
exchange energy density in the homogeneous system. It was demonstrated that the calculations based
on the homogeneous system could be satisfactory even if the electron density is far from
homogeneity. The Thomas-Fermi approach is a variational method providing an alternative to the
Schrödinger equation for calculating the ground state electronic energy associated to the electron
density n0 (r). For systems with many electrons, the advantage of the Thomas-Fermi (TF) method over

1375
the Hartree-Fock method is that, instead of solving the many-electron wave function ψ(r1 , r2 , ..., rN ) in
the Schrödinger equation (with 3N space coordinates and N spin coordinates), the TF method deals
with the electron density n(r) depending on 3 coordinates, whatever the number of electrons in the
system, plus one spin coordinate possibly.
Even if it was rigorously demonstrated by Hohenberg and Kohn [8] and Kohn and Sham [9] that
the ground state energy Eel (and related properties) of a system can be linked to the electron density
n(r), the Thomas-Fermi and Hartree-Fock methods can be regarded as forerunners of the
density-functional theory (DFT). Hohenberg and Kohn proved that the ground-state energy is
uniquely defined by a functional∗ E[n(r)] of the electron density profile, and Kohn and Sham
postulated that many-electron system can be equivalently represented by a non-interacting reference
system where each electron moves independently in the effective potential vKS (r) consisting of the
Coulomb repulsion between the electrons and the potential arising from the exchange and correlation
effects, plus the potential of nuclei whose the pseudopotential method is a representative application.
The main problem in the initial density-functional theory was to determine the functionals of the
kinetic energy and the exchange energy accounting for the electron-electron interaction. In the
modern version of the DFT, self-consistent equations are solved for a set of orbitals whose electron
density n(r) is exactly that of the real system. In practical calculations, the exchange and correlation
(xc) contributions are approximated. The simplest xc approximation is the LDA (local density
approximation) extensively used in 1970s. The GGA (generalized gradient approximations) was
developed later and, in the early 1990s, hybrid functionals were introduced by Becke [10], replacing a
fraction of the conventional GGA approximation with HF (Hartree-Fock) exchange. Among the
hybrid functionals, the most popular approximation in use in chemistry today, with about 80% of the
occurrences in the literature, is the B3LYP (exchange of Becke [11] and correlation of
Lee-Yang-Parr [12]).
This article provides a tutorial overview of the DFT. It starts with a presentation of the simple atomic
model of Thomas-Fermi—even if it has a limited success in describing the properties of real systems—
because it is a good introduction to the DFT formalism. Then, it sketches the outline of fundamentals
of the density-functional theory based on the works of Hohenberg, Kohn and Sham. As an application,
the Thomas-Fermi approximation is rederived properly from the DFT and possible improvements are
mentioned. A discussion of the exchange and correlation effects in the uniform electron gas is included
giving rise to the local density approximation (LDA), and the generalized-gradient approximations
(GGA), to treat the inhomogeneous electron gas, is also briefly presented. Lastly, mention is made
of the main hybrid functionals and the software packages successfully applied to diverse materials of
chemical, physical and biological interest.
2. Atomic Model of Thomas and Fermi
The Thomas-Fermi (TF) method has been developed independently by Thomas and Fermi
according to a semi-classical approach, which uses conventional functions instead of quantum
operators and borrows at the same time certain ideas from quantum mechanics like the electron
density. The TF method, which circumvents the concept of wave function, constitutes a good
∗
While a function is a relation that assigns a number f (x) to a number x, a functional is a relation associating a number F[ f (r)]
with a function f (r).

1376
introduction to the density-functional theory.

Contrary to the Bohr model, in the TF atomic model the electrons around the nucleus are assumed
to constitute a free electron gas at 0 K. It is well known that the average energy per electron in the
free electron gas is very large compared to kB T, so that the quantum states associated with each energy
level are so dense that it is possible to define the density of states:
!3/2
V 2m
g(E) = 2 2 E 1/2 . (12)
4π ~
As a consequence, the TF method is better suited for atoms with large atomic number Z (the number
of electrons in the system is N = Z). Taken into account that each quantum state can be occupied by
two electrons with opposite spins, the number of electrons in neutral atoms below the upper limit of
energy µ is:
Z µ !3/2
V 2m 2 3/2
Z=2 g(E)dE = 2 2 2 µ . (13)
0 4π ~ 3
This relation allows to write the chemical potential µ in terms of the electron density n = VZ :
2 ~2 kF2 p2
!
2 2/3 ~

µ = 3π n = = F , (14)
2m 2m 2m
1/3
where pF (= ~kF ) is the momentum and kF = 3π n 2
the wavevector at the Fermi level.
In the semi-classical atomic TF model, it is common to represent the kinetic energy, per electron, by
µ and to take in consideration its dependence on the position r in the atom through the electron density
2/3 2
n(r), according to the relation µ(r) = 3π2 n(r) ~
2m
. Now we consider the attractive potential energy
v(r) = −eϕ(r) of an electron at the position r from the nucleus, resulting from the electrostatic potential
ϕ(r). Thus, the total energy, E = µ(r) + v(r), of an electron at the position r, is:
~2 h 2 i2/3
E= 3π n(r) − eϕ(r) . (15)
2m
It should be stressed that the chemical potential must remain constant over the atom to avoid a
concentration of the electrons in a particular region, contrary to what happens in the Bohr model. In
fact, E is a constant negative or equal to zero since the electrons are not permitted to escape from the
atom. For a neutral atom, it is stated that E = 0 to guarantee that n(r) = 0 and ϕ(r) = 0 when r → ∞.
Eq. (15) is the central equation of the Thomas-Fermi model, which links the electron density n(r) to
the potential energy v(r), i.e., to the electrostatic potential ϕ(r), in place of the Schrödinger equation,
namely:
!3/2
1 2me
n(r) = ϕ3/2 (r) . (16)
3π2 ~2
This equation is the integral form of the TF one; it is possible to transform it into a differential
form by using the Poisson equation, which links the electrostatic potential ϕ(r) to the density of charge
ρ(r) = −en(r) as∗ :
ρ(r) en(r)
∇2 ϕ(r) = − = . (17)
ε0 ε0
∗
Here, the Poisson equation is written in the rationalized M.K.S. unit system (SI units) as ∇2 ϕ(r) = − ρ(r)
ε0
, rather than in the
unrationalized Gaussian unit system where ∇ ϕ(r) = −4πρ(r).
2

1377
Combining Eq. (16) and Eq. (17), we find:

!3/2
en(r) e 2me
∇ ϕ(r) =
2
= 2 ϕ3/2 (r). (18)
ε0 3π ε0 ~2
By using the radial part of the Laplacian in spherical coordinates, and considering the change of
variable u(r) = Zer ϕ(r), the TF equation reads:
!3/2 " #3/2
Ze d2 u(r)
!
1 d 2 dϕ(r) e 2me Zeu(r)
r = = 2 ,
r2 dr dr r dr2 3π ε0 ~2 r
!3/2 " 3 #1/2
d2 u(r) e (Ze)1/2 2me u (r)
= . (19)
dr 2 3π ε0
2 ~2 r
To solve this differential equation it is convenient to make the following change of variable x = αr,
1/3 2
with the parameter α = Z2 2/3 2me whose unit is C2/3 · s2/3 · kg−1/3 · m−2 . With this change of variable,
(3π ε0 ) ~2
the TF differential equation reduces to:
" 3 #1/2
d2 u(x) u (x)
2
= . (20)
dx x
The solution of this equation is obtained by numerical integration with the boundary conditions of
u(x), i.e., limu(x) = 1 since ϕ(r) → Zer near the nucleus, and lim u(x) → 0. The solution u(x) is a
x→0 x→∞
monotonous decreasing function of x. It has been shown that the curve starts from the point u(x = 0) =
1 with the slope a = −1.588 at the origin and approaches the x-axis asymptotically. For the numerical
integration it is advantageous to introduce the change of variable x = w2 /2. This makes the interval for
each step of the numerical integration small near the origin where u(x) changes appreciably, and larger
where u(x) has a slow variation. The detailed procedure of integration has been described by Feynman
et al. [13]. Besides, it is convenient to use the series u(x) = 1 − 1.588x + 4/3x3/2 valid for small values
of x to initiate the numerical integration.
It is interesting to see the effect of the atomic number Z on the electron density n(r). To do this,
consider the changes of variables x = αr and u(r) = Zer ϕ(r). At small values of r, it follows that
ϕ(r) ' Zer = α Zex = αϕ(x) and u(r) = Zex ϕ(x) = u(x), hence u(r) is independent of Z and ϕ(x) is
directly proportional to Z. In addition, since α is proportional to Z 1/3 , it results that ϕ(r) ∝ Z 4/3 and
n(r) ∝ ϕ3/2 (r) ∝ Z 2 .
Incidentally, note that
R ∞the substitution of the electron density (Eq. (16)) into the relation of electron
density normalization, 0 4πr2 n(r)dr = Z, allows us to find the following normalization:
Z ∞
x1/2 u3/2 (x)dx = 4πε0 , (21)
0
R∞
which becomes 0 x1/2 u3/2 (x)dx = 1 in the unrationazed Gaussian unit system.
The TF method has been found to give a rough description of the electron density and the
electrostatic potential in the atom, which is better when the number of electrons is large. Since the TF

1378
model does not agree well with the experiment, some authors tried to improve the model. For
instance, Dirac [7] introduced an exchange term that modifies Eq. (20) as follows:
#3
d2 u(x) u1/2 (x)
"
= x + 1/2 , (22)
dx2 x
where x has the same meaning as before and is defined as
!1/3
3
= Z −2/3 . (23)
32π2
Schwinger [14] provided some justification for the addition of two corrections to the ground-state
energy E of a neutral atom composed of Z electrons. One correction, Eex , lies in the exchange effects,
and the other, Equ , comes from the bulk electrons (the electrons near neither the nucleus nor the surface
of the atom) that experience a small varying potential. Both corrections are proportional to Z 5/3 , i.e.,
the total correction is Eex + Equ = 11 E = −0.2699 Z 5/3 Hartree. In addition, it has been argued [15]
9 ex
that E or its derivative is a discontinuous function of Z, and that E has a small correction [16, 17]
varying sinusoidally with Z, which have something to do with the filling of atomic shells.
3. Summary of the Density-Functional Theory
The rigorous developments of the density-functional theory were posed by Hohenberg, Kohn and
Sham [8, 9] that legitimized the model intuitively established by Thomas, Fermi and Dirac. For
practical reasons, they replace the term Uei with Vext to indicate that the electrons move in an arbitrary
external potential including the Coulomb energy due to point nuclei, so that the Schrödinger equation
(Eq. (2)) for N electrons reads:
[T + Uee + Vext ] ψel = Eel ψel , (24)
and the electronic energy (Eq. (11)) is a functional of the electron density profile n(r):
Z
Eel [n] = hψel |H| ψel i = F [n] + Vext (r)n(r)dr , (25)
where
F [n] = hψel |T | ψel i + hψel |Uee | ψel i = T [n] + Uee [n] (26)
is the so-called universal energy functional, in the
R sense that it does not contain the external potential
Vext (r) and can be determined independently of Vext (r)n(r)dr.
3.1. Formulation of Hohenberg and Kohn

The Hohenberg and Kohn (HK) approach may be understood in considering the ground state [18] of
a system of electrons at 0 K with the wave function ψ(1) , the Hamiltonian H (1) and the external potential
(1)
Vext (r). The HK approach asserts that the external potential can be uniquely determined by an electron
density, n(r), that minimizes the ground state energy:
D E Z
E = ψ H ψ
(1) (1) (1) (1)
= T [n] + Uee [n] + Vext
(1)
(r)n(r)dr . (27)

1379
Suppose that another Hamiltonian H (2) , only different from H (1) in the external potential
(2)
Vext (1)
(r),Vext (r), gives rise to the wave function ψ(2) that yields the same density n(r). Unless
(1)
Vext (r) − Vext(2)
(r) = const, ψ(2) can not be equal to ψ(1) since they satisfy different Schrödinger
equations. The wave function ψ(2) minimizes the ground-state energy for the Hamiltonian H (2) :
D E
E (2) = ψ(2) H (2) ψ(2) . (28)
By force of the variational principle, the expectation energy of the Hamiltonian is the lowest with its
“correct” wave function, hence
D E D E
E (1) = ψ(1) H (1) ψ(1) < ψ(2) H (1) ψ(2) . (29)
However, the last term of this inequality reads:

D E D E D E
ψ(2) H (1) ψ(2) = ψ(2) H (2) ψ(2) + ψ(2) H (1) − H (2) ψ(2) . (30)
The
D operators T EandR Uee being the same in the two Hamiltonians, it follows that
ψ H (1) − H (2) ψ(2) = [Vext
(2) (1) (2)
(r) − Vext (r)] n(r)dr, hence Eq. (29) becomes:
Z h i
E <E +
(1) (2) (1)
Vext (2)
(r) − Vext (r) n(r)dr . (31)
The same reasoning applied to Eq. (28) leads to:

Z h i
E <E +
(2) (1) (2)
Vext (1)
(r) − Vext (r) n(r)dr , (32)
and the sum of the two previous relations leads to the inconsistency:
E (1) + E (2) < E (1) + E (2) , (33)
indicating that the starting hypothesis is wrong. Thus the HK approach, based upon the variational
principle and reductio ad absurdum, implies that it is sufficient to use the electron density (and not the
full wave function) as the variational property, from which, once it is found, all ground-state properties
can be deduced. The considerations above are restricted to nondegenerate ground states; however, a
generalization towards degenerate states has been suggested [19] within a different and more systematic
scheme.
While the universal energy functional F [n] (composed of kinetic and potential energies of the
electron gas) remains unchanged whatever the external potential Vext (r) (see Eq. (26)), two different
external potentials can not yield the same ground state energy. The HK approach asserts that there
is a one-to-one mapping between external potential Vext (r), electron density n(r) and wave function
ψ of the system in its ground state. It implies that, given the electron density, only a unique external
potential and, consequently, its unique wave function can be determined. Conversely, given the external
potential, the corresponding
R electron density and wave function are uniquely defined. It should also
be stressed that the term Vext (r) n(r)dr is relatively easy to calculate for a given external potential,
whereas the expressions of the functionals T [n] and Uee [n] are more difficult to evaluate and some
approximations have to be made. In the initial density-functional theory, the main difficulty was to
determine the functionals of the kinetic and potential energies for the interacting electron gas.

1380
3.2. Formulation of Kohn and Sham

The Kohn and Sham (KS) approach consists in prescribing a convenient recipe to calculate the
universal energy functional F [n] = T [n] + Uee [n] for interacting electronic systems. The authors
assumed that it was reasonable to evaluate the kinetic energy T [n] of any interacting electronic system
by T 0 [n], which is the kinetic energy of the corresponding non homogeneous and non interacting
electron gas of the same density, plus the exchange-correlation energy E xc [n], to be further specified
(see Eq. (36)). The potential energy Uee [n] which, within the Hartree approximation, is that of electrons
in the field of all other electrons (Eq. (5)), is approximated by the electron density interacting with
itself, including self-interaction:
1 e2 n(r)n(r0 )
ZZ
Uee [n] ∼ W [n] =
H
drdr0 . (34)
2 4πε0 |r − r0 |
With the KS prescription, the total energy functional (Eq. (25)) reads:
Z
Eel [n] = T 0 [n] + E xc [n] + W [n] +
H
Vext (r)n(r)dr . (35)
In the KS approach, the kinetic energy T 0 [n] is not an explicit functional of density, but of orbitals of
which this density is constructed. The exchange-correlation energy functional E xc [n] incorporates the
difference between real and approximated universal energy functionals:

E xc [n] = (T [n] + Uee [n]) − T 0 [n] + W H [n] . (36)
It should be noted that E xc [n] is the only unknown functional in the right-hand side of Eq. (35), since
the other terms can be calculated exactly. In practical terms, E xc [n] is a small part of the total energy
because the magnitudes of T [n] and T 0 [n] are similar, and the same is true of Uee [n] and W H [n].
Note also that, if E xc [n] vanishes, the physical content of the KS procedure becomes identical to that
of the Hartree approximation given by Eq. (4). The simplest and reasonable approximation for E xc [n]
proposed by Kohn and Sham is the so-called local density approximation (LDA):
Z
LDA
E xc [n(r)] = e xc (n(r)) n(r)dr , (37)
where e xc (n(r)) is the exchange-correlation energy, per electron, of the uniform electron gas of local
LDA
electron density n(r). A fairly good result for E xc [n(r)] cannot be overstated when n(r) has a very
low spatial variation.
Considering the characteristics of the functionals in Eq. (35), it is not surprising that the
many-electron Schrödinger equation in the KS procedure reduces to a set of one-electron equations.
Indeed, the electron density would minimize the total energy, i.e., δEel [n] = 0, under the
normalization condition n(r)dr = N. Therefore, the method of Lagrange with the multiplier µ may
R
be used to obtain the equation requiring that the variation of the functional {Eel [n] − µN} with respect
to n(r) is equal to zero, i.e., δ {Eel [n] − µN} = 0. This variation allows us to define the functional

1381
δ
derivative, denoted by δn(r)
{Eel [n] − µN}, by the relation∗ :
δ
Z
δ {Eel [n] − µN} = {Eel [n] − µN} δn(r)dr, (38)
δn(r)
where, according to Eqs. (34, 35 and 37), the functional {Eel [n] − µN} reads:
1 e2 n(r)n(r0 )
( Z Z Z Z Z )
T 0 [n] + drdr + e xc (n) n(r)dr + Vext (r)n(r)dr−µ n(r)dr .
0
(39)
2 4πε0 |r − r0 |
Now, with Eq. (38), the variational principle δ {Eel [n] − µN} = 0 yields:
δ 1 e2 δ n(r)n(r0 ) δ
Z Z Z ( ) Z
0 = {T 0 [n]} δndr + δndrdr +
0
{e xc (n) n(r)} δndr
δn 2 4πε0 δn |r − r | 0 δn
δ δ
Z Z
+ {Vext (r)n(r)} δndr−µ {n(r)} δndr, (40)
δn δn
δ e2 n(r0 ) 0 δ
Z ( Z )
{T 0 [n]} + dr + e xc [n] + n(r) {e xc [n]} + Vext (r) − µ δn(r)dr = 0.
δn 4πε0 |r − r0 | δn
Kohn and Sham postulated that the many-electron system described by the functional derivative in
brackets of the previous equation can be equivalently represented by a non interacting reference system
described by the following one-electron Schrödinger equation:
~2 2
!
− ∇i + vKS ψi = Ei ψi , (41)
2m
where each electron moves independently in the effective potential vKS (r) consisting of the Hartree
potential wH (r), the exchange-correlation potential v xc (r) and the external potential Vext (r), namely:
vKS (r) = wH (r) + v xc (r) + Vext (r), (42)
by setting:
~2 2 δ
− ∇i = {T 0 [n(r)]} ,
2m δn Z
e2 n(r0 ) 0
wH (r) = dr , (43)
4πε0 |r − r0 |
δ
v xc (r) = e xc [n(r)] + n(r) {e xc [n(r)]} ,
δn
Ei = X
µ,
and n(r) = |ψi |2 . (44)
i
∗
The functional derivative is the natural extension of the standard derivative. If F f (r) is the functional associated with the
δF f (r)
function f (r), its variation δF f (r) generated by the variation δ f (r) defines the functional derivative, denoted by δ[f (x) ] , with the

general relation:
δF f (r)
Z
δF f (r) = δ f (r)dr.

δ f (r)

1382
Thus, the electron density in atoms, molecules and solids can be regarded as a sum of densities due
to non-interacting quasiparticles moving in the effective potential vKS (r). In the KS equation (Eq. (41)),
the potential vKS (r) is a function of the electron density n(r), which can be calculated from the wave
functions (Eq. (44)) of independent electrons (Kohn-Sham orbitals). Provided that e xc [n(r)] is known,
the KS equation may be solved in a self-consistent procedure as follows: from an initial arbitrary
electron density n1 (r), the effective potential vKS (r) is determined and Eq. (41) is solved. Then, a new
electron density n2 (r) is deduced with Eq. (44), and the procedure is repeated until convergence of the
potential vKS (r) and the electron density n(r).
It should be noted that the KS orbitals are calculated at each iteration, in terms of a set of basis
functions, and that the coefficients of linear combination of basis functions are determined as in the
Hartree-Fock calculations, but the computational time is shorter. The choice of the set of basis
functions is of great importance in the KS calculations.
The KS approach gives useful results for most physical and chemical applications [20], but a large
number of subsequent approximations must be used for correcting some of the defects of the local
density approximation E xc LDA
[n(r)] .
4. The TF Equation from the KS Approach
In this section, we rederive the Thomas-Fermi model as an example, to illustrate the use of the
KS approach. All the components of Eq. (35) are written as a function of the electron density for
the interacting electron gas enclosed in a box and, at the end, the integral form of the generalized
Thomas-Fermi equation is recovered.
We start by determining a functional form for the kinetic energy of a non interacting electron gas,
which is exact for the homogeneous electron gas (n = const.). In that case, the calculation of the
kinetic energy of the electron gas, T 0 , may be performed in k-space by doing the sum of the energy
states contained in the Fermi sphere of radius kF :
Z kF 2 2
X 2V }k
T0 = E(k) = 3
dk ,
k<kF (2π) 0 2m
Z kF 2 4
2V }k V }2 5
= 4π dk = k . (45)
(2π)3 0 2m π2 10m F
The total number of electrons in the Fermi sphere is obtained similarly as:
Z kF
2V V 3
N= 4πk2 dk = k , (46)
3
(2π) 0 3π2 F
k3
so that the electron density is n = NV = 3πF2 . Then, the kinetic energy density t (kinetic energy per unit
volume) of the homogeneous electron gas may be expressed as a function of the electron density as:
T0 3}2 2 2/3 5/3
t= = 3π n . (47)
V 10m
Incidentally, note that the kinetic energy density may be written as a function of the Seitz radius,
1/3
r s = 4πn
3
, instead of the electron density, so that the kinetic energy per electron becomes:
T0 T0 V t
= = ,
N V N n

1383
T0 3}2 2 2/3 2/3 1.105

= 3π n = , (48)
N 10m r2s
where r s is expressed in atomic unit of length and TN0 in atomic unit of energy (Hartree)∗ .
For a non homogeneous electron gas where the electron density n(r) depends on position, the same
functional form is assumed for theR kinetic energy density, hence the kinetic energy functional is
determined by the relation T 0 [n] = t [n(r)] dr, which becomes:
Z
T 0 [n] = C K {n(r)}5/3 dr, (49)
2/3
3}2
where C K = 10m 3π2 .
For an interacting electron gas, it is necessary to evaluate the other terms of Eq. (35). In the case
of the kinetic energy functional, one must add the self-interaction of electrons W H [n] given by Eq.
(34) in the RHartree approximation and, if the electron gas is subject to the external potential Vext (r), the
functional Vext (r)n(r)dr must be also taken into account.
As far as the unknown functional E xc [n] is concerned, it can be approximated by the LDA (Eq.
(37)). It is generally divided into two separate terms,
E xc [n(r)] = E x [n(r)] + Ec [n(r)] , (50)
where the exchange functional E x [n(r)] follows from the Pauli exclusion principle and affects the
electrons of same spin, whereas the correlation functional Ec [n(r)] manifests the mutual avoiding of
electrons as charged particles and, consequently, does not depend on spin. While Ec [n(r)] can be
obtained through computer simulations, E x [n(r)] is exactly known as:
Z
E x [n] = C x {n(r)}4/3 dr, (51)
1/3
e2
where C x = − 34 4πε 0
3
π
has been specified by Dirac [7]. This is the LDA for E x [n]. A dimensional
argument may be used to know the power of the electron h density: ithe dimensional homogeneity of
Eq. (35) says that each term has dimension of energy M·L2 ·T−2 , so that the energy density has
h i h i
dimension M·L−1 ·T−2 , and the electron density n(r) has dimension L−3 . On the other hand, the
term E x [n] arising from the self-interaction of electronsh W H [n]i (Eq. (34)) requires that the energy
dE x [n] 2
density scales dr ∼ e rnn r3 , where e2 has dimension M·L3 ·T−2 . Then, if the power law assumed for
the energy density is dEdrx [n] ∝ e2 n(r)ν , the dimensional analysis of this relation yields M·L−1 ·T−2 =

M·L3 ·T−2 ·L−3ν , leading to the value of the exponent ν = 4/3. Thus, in the absence of correlation, the
total energy functional for a non homogeneous and interacting electron gas reads explicitly:
Z Z Z Z
Eel [n] = C K {n(r)} dr + w (r)n(r)dr + C x {n(r)} dr + Vext (r)n(r)dr.
5/3 H 4/3
(52)
In order to recover the generalized Thomas-Fermi equation, let’s seek the distribution of the electron
density that would minimize the total energy, i.e., δEel [n] = 0, under the normalization condition
∗
Atomic units are used throughout this article. With m = }2 = e2 = 4πε0 = 1, it is found that 1 a.u. of length = 0.529×10−10 m and
1 Hartree = 2 Rydberg = 27.21 eV = 4.36×10−18 J.

1384
n(r)dr = N. Using the method of Lagrange with the multiplier µ, write the equation requiring that
R
δ
δn(r)
{Eel [n] − µN} = 0, as in Eq. (38):
e2 n(r0 ) 0 4
Z ( Z )
5
C K {n(r)}2/3 + dr + C x {n(r)} + Vext (r) − µ δn(r)dr = 0.
1/3
(53)
3 4πε0 |r − r0 | 3
In brackets is the functional derivative that must be equal to zero to fulfill the minimization condition.
Note that this functional derivative corresponds to the generalized TF equation:
e2 n(r0 ) 0 4
Z
5
C K {n(r)} +
2/3
dr + C x {n(r)}1/3 + Vext (r) = µ. (54)
3 4πε0 |r − r0 | 3
Comparing this equation to Eq. (15) shows that the Lagrange multiplier µ represents the classical
energy of an electron on the Fermi sphere. It could be identified with the Fermi energy and it really
has the meaning of the chemical potential since it is defined by ∂E∂N
el [n]
= µ.
5. Improved Functional Approximations in DFT
With the KS equation, the DFT is an exact theory to the extent that the electron density minimizing
the total energy is the true electron density of the interacting electron gas. However, the theory is
unworkable as long as the exchange-correlation energy functional E xc [n] is unknown.
In this section, we describe some improvements of the DFT for an electron gas with slow spatial
variations of the electron density such as:
n(r) = n0 + n1 (r), (55)
subject to the conditions n1 (r) n0 and n1 (r)dr = 0. This assumption is appropriate for the
R
intermediate region in atoms or molecules where the electron density is almost constant, but it is no
longer valid near the atomic radius where the electron density is very dense and far away where the
electron density is very low.
The merit of Hohenberg and Kohn [8] has been to show that the universal energy functional
containing the kinetic energy T 0 [n] , the Hartree potential W H [n] and the exchange-correlation energy
E xc [n] in Eq. (35), i.e., F [n] = T 0 [n] + W H [n] + E xc [n] , has an expansion of the following form∗ :
ZZ
1
F [n] = F [n0 ] − K(|r − r0 |)n1 (r)n1 (r0 )drdr0 + ... (56)
2
As expected, this expansion has no linear term in n1 (r) because of the condition n1 (r)dr = 0.
R
5.1. Without Exchange-correlation Effects

In the absence of exchange-correlation effects, HK suggest to expand the total energy (Eq. (35)),
free from the exchange-correlation energy, according to the expression:
ZZ Z h
1 i
Eel [n] = T 0 [n0 ] − K(|r − r |)n1 (r)n1 (r )drdr +
0 0 0
wH (r) + Vext (r) n(r)dr, (57)
2
∗
see Appendix II of Ref. [9].

1385
where the first term T 0 [n0 ] is given by Eq. (49), i.e., T 0 [n0 ] = C K (n0 )5/3 dr, whereas the unknown
R
function K(|r − r0 |) remains to be determined. To define it, it is convenient to express the total energy
in the q-space owing to the Fourier transform∗ , whereby the functional becomes:
Z Z h
1 1 1 i
Eel [n] = T 0 [n0 ] + K(q) |n(q)|2
dq + w H
(q) + Vext (q) n(q)dq . (58)
2 (2π)3 (2π)3
As usual, the electron density would minimize the total energy, i.e., δEel [n] = 0. In q-space, the
minimization of Eel [n] leads to the functional derivative:
δ
Z h
1 i
{Eel [n]} = K(q)n(q) + w H
(q) + V ext (q) dq = 0, (59)
δn(q) (2π)3
so that the FT of the electron density reads:
wH (q) + Vext (q)

n(q) = − . (60)
K(q)
Now, to determine the function K(q), we proceed to define the linear response function of the
uniform interacting electron gas. Recall briefly the implementation of the method. According to the
Poisson equation, the Hartree potential wH (r) is related to the varying part of the electron density n1 (r)
2 2
by the relation ∇2 wH (r) = − εe0 n1 (r) in real space and by wH (q) = ε0eq2 n(q) in reciprocal space. But the
h i
Hartree potential wH (q) is also related to the total potential wH (q) + Vext (q) by the relation:
e2 h H i
wH (q) = w (q) + V ext (q) X0 (q) , (61)
ε0 q2
where X0 (q) is the linear response function of the uniform electron gas defined by the expression [21]:
4 − η2 η+2
" ! #
mkF
X0 (q) = − 2 2 1 + ln < 0, (62)
2π ~ 4η η−2
with η = q/kF . Hence, the FT of the electron density n(q) is connected to the total potential by the
following relation: h i
n(q) = wH (q) + Vext (q) X0 (q). (63)
From the comparison of Eqs. (60) and (63), it is found that K(q) identifies with − 1/X0 (q) . As

expected, one of the most significant feature of K(q) is the logarithm singularity at q = 2kF . There is a
slight inflection in the curve in the neighborhood of q = 2kF , not visible on any scale. Despite this
smallness, the singularity may have a large effect on properties depending upon the linear response
function X0 (q). The function K(q) is expanded well in a power series of q2 , but the following
∗
The Fourier transforms of the different functions look like that of the electron density:
Z
n(q) = n1 (r) exp (−iq · r) dr,
Z
1
n1 (r) = n(q) exp (iq · r) dq.
(2π)3

1386
approximation [22] may also be employed to approximate K(q) in the limit of low and large regions
of q:
!2 
π2 ~2 

q 
K(q) = γ + 3λ  > 0 . (64)
mkF 2kF
It is easy to check that this empirical relation fits the low-q region (slow density variations) with the
parameters γ = 1 and λ = 1/9, and the large-q region (rapid density variations) with the parameters
γ = −1 and λ = 1.
Lastly,
RR we calculate the inverse Fourier transform of K(q) in order to determine the correction
− 12 K(|r − r0 |)n1 (r)n1 (r0 )drdr0 to bring in the kinetic energy functional after the unknown function
K(|r − r0 |) is determined. In order to facilitate this, we use Eq. (64) to approximate the function K(q).
With the definitions:
Z
1
δ(r) = exp(iq · r)dq (65)
(2π)3
Z
1
and ∇r δ(r) =
2
(−q2 ) exp(iq · r)dq, (66)
(2π)3
2
π2 ~2 π2 ~2 π2 ~2
and − 34 mkπ 3~(r) λ∇2r δ(r), respectively, so that:
2 2
the inverse FT of mkF (r)
γ and mkF
3λ 2kqF are mkF (r)
γδ(r)
F
π2 ~2 3 π2 ~2
" #
K(r) = FT −1
K(q) = γ− λ∇ δ(r).
2

(67)
mkF (r) 4 mkF3 (r) r
Taking into account that 3π2 n(r) = kF3 (r), the correction to the kinetic energy becomes:
γ π2 ~2 n21 (r) λ ~2 |∇r n(r)|2

Z Z Z Z
1
− K(|r − r |)n1 (r)n1 (r )drdr = −
0 0 0
dr + dr. (68)
2 2 m kF (r) 8m n(r)
It should be mentioned that the gradient correction in the previous relation coincides with the von
Weizsäcker [23] correction based on intuitive and semiclassical arguments. With a positive
contribution, Rit improves the result of the kinetic energy, since the Thomas-Fermi contribution,
T 0 [n0 ] = C K (n0 )5/3 dr, is always smaller than the Hartree-Fock result. Consequently, the expansion
proposed by Hohenberg and Kohn for the kinetic energy is valid for both slow and rapid spatial
variations of the electron density [24].
5.2. With Exchange-correlation Effects

To treat the exchange-correlation effects, Hohenberg and Kohn proposed to use the same
approximation as in Eq. (56) to expand the exchange-correlation energy, namely:
ZZ
1
E xc [n] = E xc [n0 ] −
LDA
K xc (|r − r0 |)n1 (r)n1 (r0 )drdr0 , (69)
2
LDA
where the leading term E xc [n0 ] is the so-called local density approximation (LDA) founded on the
uniform electronR gas, whose the general form is given by Eq. (37), i.e.,
LDA
E xc [n(r)] = e xc (n(r)) n(r)dr. Usually, the exchange-correlation energy per electron of an

1387
interacting electron gas, e xc (n(r)) , splits into two components: the exchange contribution e x (n(r))
and the correlation contribution ec (n(r)) .
For estimating the exchange energy e x (n) , Slater [25] suggested the use of the uniform electron gas
at constant electron density n, known exactly. This contribution may be calculated with the Hartree-
Fock approximation via Eq. (9) representing the exchange energy of the electron i in the field of all
other electrons:
X Z e2 ψ∗j (r j )ψi (r j )ψ j (ri )
E x (k) = − dr j . (70)
j,i
4πε0 ri j ψi (ri )
From here the subscripts i and j will be replaced by k and k0 , so that E x (k) is the exchange energy for
the state k. Applied to the uniform electron gas with the wave function V −1/2 exp(ik · r), the exchange
energy for the state k becomes:
e2 X exp [i (k − k0 ) · (r2 − r1 )]
(Z )
E x (k) = − dr2 ,
4πε0 V k0 r12
e2 X exp [i (k − k0 ) · r2 ]
( Z )
= − 0
dr2 .

exp −i k − k · r1 (71)
4πε0 V k0 r12
To calculate the integral in brackets, the standard method consists in assuming, for the sake of
calculation, that the electrostatic potential ϕ(r) at point r is produced by the formal distribution of
charge exp [i (k − k0 ) · r] , even if it is not the real distribution. Hence, the electrostatic potential at
R exp[i(k−k0 )·r2 ]
point r1 is ϕ(r1 ) = r12
dr2 . Besides, ϕ(r) must satisfy the Poisson equation
∇ ϕ(r)
2
= 0
−4π exp [i (k − k ) · r] (see footnote in the next page) whose solution is
exp[i(k−k0 )·r]
ϕ(r) = 4π |k−k0 |2 . Consequently, Eq. (71) simplifies as:
e2 X 4π e2 1 dk0
Z
E x (k) = − = − , (72)
4πε0 V k0 |k − k0 |2 ε0 (2π)3 |k − k0 |2
V
where the use is made of the volume element dk0 containing (2π) 0
3 dk states, which allows one to
transform the sum k0 into the integral over the Fermi sphere. The previous integral on k0 may be
P
carried out with the spherical coordinates (k0 ∈ [0, kF ] ; θ ∈ [0, π] ; ϕ ∈ [0, 2π]), namely,∗
|k + k0 | 0
Z kF 0
dk0
Z
k
= 2π ln dk . (73)
|k − k0 |2 0 k |k − k0 |
Care must be taken to evaluate this integral because the integrand becomes infinite when k0 = k.
Fortunately, the two integrals from 0 to k and from k to kF converge, so that the integration by parts
∗
This integration is performed as follows:
kF 2π π
dk0 k0 dθ
Z Z Z Z
= dk0 k0 sin θdϕ ,
|k − k0 |2 0 0 0 k2 + k02 − 2kk0 cos θ
Z kF Z +1
dx
= 2π k02 dk0 ,
0 −1 k2 + k02 − 2kk0 x
kF
k0 k + k0
Z
= 2π ln dk0 .
0 k k − k0

1388
may be carried out. The final result is∗ :
dk0 kF2 − k2 k + kF
Z " #
= πkF 2 + ln , (74)
|k − k0 |2 kF k k − kF
and the exchange energy associated with the state k reads:
e2 kF kF2 − k2 kF + k
" #
E x (k) = − 2+ ln . (75)
4πε0 2π kF k kF − k
Now, to determine the total exchange energy of the uniform electron gas containing N electrons,
it is necessary to sum over the states k. As usual, the replacement of the sum by an integral over the
Fermi sphere yields the expression:
kF2 − k2 k + k0
Z kF "
e2 kF V
X #
E x (k) = − 3
2+ ln 0
4πk2 dk . (76)
k
4πε0 2π (2π) 0 kF k k−k
It is found that the integral above is equal to 4πkF3 . Hence, with kF3 = 3π2 n, the exchange energy per
electron of the uniform electron gas is:
1X e2 2V kF4
e x (n) = E x (k) = −
N k 4πε0 (2π)3 N
!1/3
e2 1 V 2 4/3 e2 3 3
= − 3π n =− n1/3 . (77)
4πε0 4π3 N 4πε0 4 π
Incidentally, note that the exchange energy per electron may be written as a function of the Seitz radius
1/3 1/3
r s = 4πn
3
= 9π
4
1
kF
, instead of the Fermi wavevector, as:
!1/3 !1/3
e2 3 3 3 1 0.458
e x (n) = − =− , (78)
4πε0 4 π 4π rs rs
2
e
where r s is expressed in atomic unit of length and e x (n) in Hartree, with 4πε 0
= 1. Note also that many
investigations [7, 25, 26] of the exchange energy have lead to the expression obtained by multiplying
the Slater exchange energy with a factor α:
!1/3
e2 3 3
e x (n) = −α n1/3 . (79)
4πε0 4 π
As an example, by calculating the effective exchange potential e x (n) just for a state at the Fermi level,
Kohn and Sham [9] have shown that α = 2/3 instead of 1 for Slater.
∗
Use is made of the integrals:
Z " !#
1 02 3 1 0
k0 ln |k + k0 | dk0 = k − k2 ln |k + k0 | + (k + k0 ) k− k ,
2 2 2
Z " !#
1 02 3 1 0
k0 ln |k − k0 | dk0 = k − k2 ln |k − k0 | + (k − k0 ) k+ k .
2 2 2

1389
Finally, the dependence of the exchange energy on r, aiming at taking account of the inhomogeneity
of the electron gas, may be introduced as in Eq. (51). As a result, the LDA exchange energy becomes:
!1/3 Z
e2 3 3
E xLDA [n(r)] = − n4/3 (r)dr. (80)
4πε0 4 π
For the correlation energy ec (n), which is the difference between the exact energy of the electron
gas and the kinetic and exchange energies, the situation is more complicated since ec (n) of the uniform
electron gas is not known exactly. Among the old estimations of ec (n) , that of Wigner [27], valuable
at intermediate densities (5 > r s ≥ 2) such as those in metals, is:
0.44
ec (n) = − , (81)
r s + 7.8
and that of Gell-Mann and Bruckner [28], for high densities (r s ≤ 2), is:
ec (n) = + 0.031 ln r s − 0.048 . (82)
By using quantum Monte Carlo (QMC) calculations for the uniform electron gas, Ceperley and
Alder [29,30] obtained precise numerical results for ec (n) , from which many approximate expressions
have been derived. For instance, Perdew and Wang [31] proposed the following expression, which fits
well the QMC results for a large domain of density:
 
 1 
ec (n) = −2c0 (1 + αr s ) ln 1 + 1/2  , (83)
2c0 β1 r s + β2 r s + β3 r3/2
s + β r
4 s
2 
with c0 = 0.03109, α = 0.21370, β1 = 7.5957, β2 = 3.5876, β3 = 1.6382, β4 = 0.49294.

It should be stressed that the second term of Eq. (69) resulting in a von Weizsäcker gradient
correction, as in Eq. (68), almost never improves on the local density approximation. Therefore, a
good level of approximation for the LDA is obtained just in doing the sum of the kinetic energy (Eq.
(48)), the exchange energy (Eq. (78)) and the correlation energy (Eq. (83)), per electron, of the
uniform electron gas at the electron density n, namely:
 
1.105 0.458  1 
e(n) = − − 2c0 (1 + αr s ) ln 
 1 +  . (84)
r2s rs 
2c0 β1 r s + β2 r s + β3 r s + β4 r s
1/2 3/2 2 
The LDA has been found to give quite good results for most applications with large negative values
for e x (n) and small negative values for ec (n). The reason of this relative success is in the
compensation of the errors, since the LDA underestimates the exchange energy e x (n) and
overestimates the correlation energy ec (n), and also in the fact that the exact boundary conditions
concerning the exchange-correlation are satisfied in the LDA. Levy et al. [32, 33] have addressed (on
the basis of mathematical requirements) the problems of analysis that arise in discussing Coulomb
systems through the density functional approach. Rigorous necessary requisites upon the
exchange-correlation functionals have been found in the form of inequalities [34], which are satisfied
within the LDA and are useful in building and improving functional approximations.

1390
6. Beyond the Local Density Approximation
The LDA approach, which is based on the uniform electron gas at a given electron density n, has a
formal justification for systems with slow and weak variations. But the atomic and molecular systems,
spatially inhomogeneous, differ from the uniform electron gas, therefore more accurate approximations
have been introduced making use of the electron density and its gradients.
Since the exchange-correlation energy functional E xc [n] is not known exactly, the corrections to
the LDA approach should ideally contain terms (i) of well-defined physical origin, (ii) easily usable
for practical calculations in physics and chemistry and (iii) allowing transferability from one system to
another. Roughly speaking, these approximations are classified into two main groups. One is based on
an exact equation for E xc [n] derived from quantum mechanics including scaling relations, and the other
of a more empirical nature is based on the gradient expansion of E xc [n] in power of small deviation of
the density.
6.1. Average Exchange Hole

Before discussing these approximations, it is important to introduce the concept of average
exchange hole around a given point r. In any electronic system, the electrons do not move
independently but repel one another by the Coulomb repulsion. The coupling (correlation) among the
electronic motions are usually called the Coulomb correlations. It should be stressed that Coulomb
correlations are ignored in the Hartree and Hartree-Fock methods, because each electron is supposed
to move in the average charge distribution of electron system. However, in the Hartree-Fock method
where use is made of a determinantal wave function, there are correlations of another type among the
positions of electrons with parallel spins due to the Pauli principle. This type of correlation means
that (i) there is zero probability to find two such electrons at same position and (ii) there is a small
probability to find an electron in the immediate vicinity of each other. To take account of this aspect
of Pauli’s principle, it is usually assumed that each electron is surrounded by a hole in the distribution
of electrons with parallel spins, called the exchange hole of limited extent.
It can be recalled that the Hartree potential energy (Eq. (5)) of an electron at position r1 is due to a
charge distribution whose the density at r2 is :
X e 2
− ψ j (r2 ) . (85)
j
4πε0
In addition, the exchange energy (Eq. (9)) of an electron at position r1 is due to an exchange charge
distribution whose the density at r2 is:
X e ψ∗j (r2 )ψi (r2 )ψ j (r1 )
+ , (86)
j
4πε0 ψi (r1 )
so that the total exchange energy reads:
e ψ j (r2 )ψi (r2 )ψ j (r1 )

XZ ∗
e
+ dr2 = , (87)
j
4πε0 ψi (r1 ) 4πε0

1391
because of the orthogonality and normalization of wave functions, ψ∗j (r2 )ψi (r2 )dr j = δi j . Thus, this
R
result suggests that an electron may be viewed as surrounded by a positive cloud of total charge e, such
as the electron with its surrounding hole together form a neutral object.
Taking advantage of the definition of exchange charge density given by Eq. (86), we calculate the
exchange charge density for the uniform electron gas by using the free electron wave function. The
exchange charge density produced at position r2 by an electron in state k at r1 reads:
X e ψ∗ 0 (r2 )ψk (r2 )ψk0 (r1 )
k
. (88)
k 0 4πε 0 ψ k (r1 )
After introducing the free electron wave function and changing the sum by an integral over the Fermi
sphere, the exchange charge density at position r2 becomes:
Z Z
e 1 V e 1
exp i k−k ·(r2 −r1 ) dk =
0
exp [ik·(r2 −r1 )] exp −ik0 ·(r2 −r1 ) dk0 .
0
4πε0 V (2π) 3 4πε0 (2π) 3
(89)
As usual, the integral in this relation is calculated with spherical coordinates. The exchange charge
density produced at position r2 by an electron in state k at r1 reads∗ :
!
e 1 sin (kF r12 ) − (kF r12 ) cos (kF r12 )
exp [ik · (r2 − r1 )] 3
. (90)
4πε0 2π2 r12
To determine the average exchange charge density, Eq. (90) could be used for a state at the top
of Fermi distribution, with the value at k = kF . But it seems better to average Eq. (90) over all the
occupied states k, containing N/2 electrons with same spin [4]. In this case, the average exchange
charge density at position r2 for an electron at r1 is:
! Z
e 1 sin (kF r12 ) − (kF r12 ) cos (kF r12 ) 2 V
3
exp [ik · (r2 − r1 )] dk. (91)
4πε0 2π2 r12 N (2π)3
Owing to the result of the previous footnote, this expression reads:
!2
e kF6 sin (kF r12 ) − (kF r12 ) cos (kF r12 )
. (92)
4πε0 2π4 n 3
kF3 r12
In terms of the electron density (kF3 = 3π2 n), the average exchange charge density at position r2 for an
electron at r1 becomes: !2
e 9n sin (kF r12 ) − (kF r12 ) cos (kF r12 )
. (93)
4πε0 2 kF3 r12
3
∗
If θ is the angle between the vectors (r1 − r2 ) and k0 , the integral becomes:
Z Z kF Z 2π Z π
exp −ik0 · (r2 − r1 ) dk0 = dk0 k0 sin θdϕ exp (ik0 r12 cos θ) k0 dθ,

0 0 0
Z kF Z π
= 2π k02 dk0 exp (ik0 r12 cos θ) sin θdθ,
0 0
kF
k0 sin (k r12 ) 0 0
Z
= 4π dk ,
0 r12
!
sin (kF r12 ) − (kF r12 ) cos (kF r12 )
= 4π 3
.
r12

1392
This function decreases rapidly when the distance r12 increases. Besides, when r12 tends to zero its
e n
value is 4πε 0 2
. This value is opposite to the average density of electrons with same spin. It means that,
if the electronic motions were uncorrelated, the average density of electrons with same spin would be
e n
(− 4πε 0 2
), in accordance with the fact that an electron plus its surrounding hole is a neutral object. Since
the electronic motions are correlated, the average density charge distribution of electrons with same
e n
spin is given by adding Eq. (93) to (− 4πε 0 2
), i.e.:
 !2 
e n  sin (kF r12 ) − (kF r12 ) cos (kF r12 )  e n
− 1 − 9 3 3  = − g(r1 , r2 ). (94)
4πε0 2 kF r12 4πε0 2
It can be seen that the term in brackets vanishes for r12 = 0 and approaches 1 as r12 is increasing.
It has been shown [35] to be the pair correlation function g(r1 , r2 ) of the uniform electron gas with
density n. For uncorrelated electronic motions, g(r1 , r2 ) = 1. It describes the effect of the hole dug into
the electron density at r2 by the electron at r1 , and it allows us to determined the probability of finding
simultaneously two electrons of same spin at given points being a distance r12 apart.
6.2. Exchange-correlation Energy Functional

Now we are in position to determine the exact expression for the exchange-correlation energy
functional E xc [n(r)] in interacting systems resulting from the interaction between an electron and its
exchange-correlation hole n xc (r1 , r2 ). As already mentioned, the formal definition of the
exchange-correlation energy E xc [n(r)] is given by Eq. (36),
E xc [n] = (T [n] + Uee [n]) − T 0 [n] − W H [n] , (95)
but, in the DFT context, it has been proved useful to consider another expression, known as the
adiabatic connection formula, differing from the approximate form of Hohenberg and Kohn (Eq.
(69)) by the presence of the exchange-correlation hole n xc (r1 , r2 ) . To achieve this, it is assumed that a
system in a stationary state remains in its initial state when the parameter λ describing the system
changes very slowly. Using this parameter λ, the exchange-correlation energy functional E xc [n(r)]
containing the kinetic energy of non interacting electron gas may be written as:
E xc [n] = hψλ |T [n] + λUee [n]| ψλ iλ=1 − hψλ |T [n] + λUee [n]| ψλ iλ=0 − W H [n] ,
Z 1
d
= hψλ |T [n] + λUee [n]| ψλ i dλ − W H [n] . (96)
0 dλ
The adiabatic connection permits us to pass from the interacting system (λ = 1) to the noninteracting
system (λ = 0) by infinite number of possibilities. In the Hamiltonian of KS equation (Eq. (41)), the
e2 λ
external potential λVext (r) and the Coulomb interaction 4πε 0 r12
of partial strength λ are taken in order
that the density n(r) is the same at any λ between 0 and 1.
Eq. (96) may be transformed with the Hellman-Feynman theorem [36–38] under the following
form∗ : Z 1
E xc [n] = hψλ |Uee [n]| ψλ i dλ − W H [n] . (97)
0
∗
The Hellman-Feynman theorem allows us to know how theD energy Eλ E= hψλ |Hλ | ψλ i varies as a function of λ. Due to the properties
of the Hamiltonian and the wavefunction, it is found that dE
dλ
λ
= ψ∗λ ∂H
∂λ
λ
ψλ . The theorem has been first used to determine the forces in
a molecular system.

1393
Of course, W H [n] is given by Eq. (34).
1 e2
Z Z
n(r1 )n(r2 )
W [n] =
H
dr1 dr2 . (98)
2 4πε0 |r1 − r2 |
On the other hand, the first term of the right-hand side of Eq. (97) has been proved to be equal to [39]:
n(r1 )n2 (r1 , r2 )

Z 1
1 e2
Z Z
hψλ |Uee [n]| ψλ i dλ = dr1 dr2 , (99)
0 2 4πε0 |r1 − r2 |
with:
n2 (r1 , r2 ) = n(r2 ) + n xc (r1 , r2 ) , (100)
so that the exchange-correlation energy functional E xc [n(r)] simplifies as:
1 e2 n xc (r1 , r2 )
Z Z
E xc [n(r)] = n(r1 )dr1 dr2 . (101)
2 4πε0 |r1 − r2 |
As a generalization of Eq. (94), the exchange-correlation hole n xc (r1 , r2 ) is the following integral [40]
over the coupling parameter λ:
Z 1
n xc (r1 , r2 ) = n(r2 ) dλ g(r1 , r2 ; λ) − 1 .

(102)
0
Here, g(r1 , r2 ; λ) is the pair correlation function of the system with the physical density n(r) calculated
e λ
from the Coulomb interaction potential 4πε 0 r12
. By its definition, the exchange-correlation hole
n xc (r1 , r2 ) must satisfy the sum rule requiring that the hole contains one electron, for all r, i.e.:
Z
n xc (r1 , r2 ) dr2 = −1, (103)
reflecting total screening of the electron at r1 . The pair correlation function g(r1 , r2 ; λ) may be obtained
from different theoretical approaches and computer simulations. However, for inhomogeneous systems
little is known about it.
Several approximations based on Eq. (101) have been proposed to improve the LDA approach,
bypassing the pair correlation function, since R they only need the exchange-correlation energy per
electron e xc (n(r)) to determine E xc [n(r)] = e xc (n(r)) n(r)dr, where
1 e2
Z
n xc (r, r2 )
e xc (n(r)) = dr2 . (104)
2 4πε0 |r − r2 |
Gunnarsson et al. [41] proposed two nonlocal approximations: (i) the averaged-density (AD)
approximation, which replaces n(r2 ) in Eq. (102) by the averaged density n(r1 ) obtained with a
weight function to be determined, and (ii) the weighted-density (WD) approximation that uses the pair
correlation function of the homogeneous electron gas in Eq. (102). Langreth and Perdew [42, 43]
developed a similar theory for the exchange-correlation energy by investigating Eq. (101) in
reciprocal space. The AD and WD approximations have attracted much attention as a potential tool
for investigating solid state structures, metallic surfaces and quantum chemistry [44, 45]. During the

1394
same period of time, some simple parametrized expressions have been introduced to describe the
exchange-correlation energy [46–48].
Despite the use of the adiabatic connection formula making the explicit link between
exchange-correlation energy E xc [n(r)] and exchange-correlation hole n xc (r1 , r2 ), the results may
deviate substantially from the simulation ones. Therefore, the adiabatic connection formula is much
less used than the so-called gradient expansion approximations. In this class of approximations, one
tries to systematically calculate the gradient corrections by a generalization of the famous
lowest-order von Weizsäcker gradient correction.
Since any real system has a spatially varying density n(r), it would clearly be useful to include
information on the rate of variation in the exchange-correlation energy E xc [n(r)]. A first attempt at
doing this was to use Eq. (69) even if, as already mentioned, the lowest order of the expansion is often
disappointing compared with the LDA results. Then, higher order corrections have been included.
Using the new coordinates r = r0 − r00 and R = (r0 + r00 ) /2, for a slowly varying density, the electron
density reads [49]:
1 1 1
n(r0 ) = n(R) + ( r · ∇)n(R) + ( r · ∇)2 n(R) + ... (105)
2 2 2
and the exchange-correlation energy functional (Eq. (69)) may be written in terms of the following
gradient expansion:
π
Z Z Z ∞
E xc [n] = e xc (n(R)) n(R)dR − 2
|∇n(R)| dR r4 K xc (r; n(R))dr
3 0
π
Z h i Z ∞
− ∇n(R)·∇ ∇2 n(R) dR r6 K xc (r; n(R))dr + ..., (106)
60 0
LDA
where the leading term corresponds to the LDA contribution E xc [n(R)] . The other terms are very
difficult to calculate and little is known about them.
7. Separation of the Exchange and Correlation Functionals
The failures of the LDA for strongly inhomogeneous systems have lead to going beyond the lowest
order gradient correction with different sequences of the series given in Eq. (106). Important progress
has been made in deriving gradient expansion approximations for the exchange-correlation functional,
under the general form:
Z
E xc [n] = e xc (n(r)) f xc (n(r), ∇n(r)) n(r)dr. (107)
This relation is known as the generalized-gradient approximation (GGA); the function

f xc (n(r), ∇n(r)) is the enhancement factor depending on the density and the density gradient. The
construction of f xc (n(r), ∇n(r)) has been made by considering dimensional analysis, sums rules as
well as asymptotic behavior of density and effective potentials around atoms or molecules. It is a
matter of fact that different GGAs differ greatly in the choice of the function f xc (n(r), ∇n(r)) .
Since the exchange energy E x is given by an exact expression, contrary to the correlation energy
Ec (Section 5), and since |E x | is greater than |Ec | , it seems natural to treat separately the exchange
energy and the correlation energy. In this section, two approximations for pure exchange functional and
correlation functional will be briefly presented, accounting for the most popular DFT investigations.

1395
7.1. Exchange Energy Functionals

For the exchange energy functionals, Becke [11] showed that E x can be constructed with one
parameter fitted to exchange energies of atoms. This corresponding semiempirical functional,
commonly denoted E xB88 , reproduces the exact asymptotic behavior of E x and the effective potential.
For spin polarized systems, its expression reads:
XZ xα2
E xB88 = E xLDA −β n4/3
α dr, (108)
α 1 + 6βxα sinh−1 xα
where xα = |∇nα |
n4/3
and nα stands for spin up and spin down densities. It should be noted that this
α
functional contains the single parameter β = 0.0042. If sinh−1 xα is expanded into a power series, the
E xB88 functional may be written under the gradient expansion:
!1/3 Z
e2 3 3 h i
E xB88 =− n4/3 (r) 1 + µxα2 + ... dr, (109)
4πε0 4 π
where the leading term is E xLDA [n(r)] of Eq. (80). The E xB88 functional is particularly well adapted to
the calculation of the atomization energies of molecules, i.e., the required energy to fully dissociate a
molecule to free atoms.
Other very convenient exchange energy functionals have been developed by the Perdew group.
These functionals are based on rational principles of quantum mechanics and theoretical constraints.
In particular, for spin unpolarized systems, the so-called PBE (Perdew-Burke-Ernzerhof [50, 51])
exchange is often used under the form:
!1/3 Z  
e2 3 3  κ 
E xPBE =− n4/3 (r) 1 + κ −  dr, (110)
4πε0 4 π (1 + µs2
)

κ
where s = 2k|∇n|
Fn
(with kF3 = 3π2 n) and the two universal parameters are µ = 0.2195 and κ = 0.804. The
PBE functional, bringing minor modifications to the PW (Perdew-Wang [31]) functional, increases the
precision of the results of total energies of atoms and atomization energies of molecules. In addition,
its recent version [52] improves equilibrium properties of packed solids and their surfaces.
7.2. Correlation Energy Functionals

Correlation energy functionals have been also constructed by the Perdew group. The first such
expression, known as PW (Perdew-Wang [31]) correlation, has been refined and accurately represented
by the expression: Z
Ec =
PW
ec (r s , ζ) + H(r s , ζ, t) n(r)dr.

(111)
1/3
In this relation, the variables are the Seitz radius r s = 4πn3
introduced in place of the Fermi
wavevector kF (with r s kF = 9π/4), the relative spin polarization ζ = nn11 +n
3 3 −n2
2
, which ranges from 0 for
unpolarized system to ±1 for a fully spin polarized system, and the reduced density gradient t

1396
measuring the rapidity of variation of n(r) on the scale of the screening length 1/k s . The reduced
density gradient t is defined as:
|∇n|
t= , (112)
2φk s n
1/2 h i
where k s = 4kπF is the Thomas-Fermi wavevector and φ(ζ) = 21 (1 + ζ)2/3 + (1 − ζ)2/3 is the spin
interpolation formula taking account of the polarization dependence.
Accurate analytic representations of ec (r s , ζ) are available under the general form [31]:
ec (r s , ζ) = ec (r s , 0) + φ(ζ) [ec (r s , 1) − ec (r s , 0)] , (113)
where ec (r s , 0) is the correlation energy of the uniform electron gas given by Eq. (83). As far as
H(r s , ζ, t) in Eq. (111) is concerned, it reads:
β 2 1 + At2
( " #)
H(r s , ζ, t) = γφ ln 1 + t
3
, (114)
γ 1 + At2 + A2 t4
where the more recent values of the universal parameters [50] are γ = 0.031091, β = 0.066725 and
#−1
β ec (r s , ζ)
" !
A= exp − −1 . (115)
γ γφ3
The PW correlation energy functional (Eq. (111)) has been derived from the slowly varying limit
(t → 0) and the rapidly varying limit (t → ∞), as well as from the uniform scaling to the high
density limit (γ → ∞), which arise from numerical calculations of the local density and the gradient
contributions [50]. It is found that H(r s , ζ, t) ∼ βφ3 t2 when t → 0 and H(r s , ζ, t) ∼ −ec (r s , ζ) when
t → ∞, so that the correlation energy per electron vanishes in systems with rapid spatial variations of
the electron density. Since the second-order gradient contribution H(r s , ζ, t) is a growing function of t,
going from 0 to |ec (r s , ζ)| , EcPW tends to a small negative constant.
This is a much more difficult problem to test the correlation energy than the exchange energy, which
is accurately defined. Therefore, a useful way to think through gradient approximations is to compare
one GGA with another. There is an interesting correlation energy functional derived by Lee, Yang and
Parr [12] from the early work of Colle and Salvetti [53] and known as LYP correlation energy, whose
results compare favorably with those by EcPW . For the unpolarized system, the LYP correlation energy
may be written in terms of electron density as:
Z n h io n(r)1/3
Ec = −a
LY P
n(r) + bn(r)−2/3 [tHF (r) − 2tW (r)] exp −cn(r)−1/3 dr, (116)
n(r)1/3 + f
where the constants are a = 0.049, b = 0.132, c = 0.2533 and f = 0.349, tW (r) is the local Weizsacker
kinetic energy density, and tHF (r) the Hartree-Fock kinetic energy to second order defined as:
|∇n(r)|2 ∇2 n(r)
tW (r) = − , (117)
8n(r) 8
tW (r) ∇2 n(r)
" #
tHF (r) = C K n(r) +
5/3
+ , (118)
9 18

1397
2/3
with C K = 10
3
3π2 .
Note that the combination of E xB88 , for the exchange energy, with EcLY P , for the correlation energy,
forms the GAUSSIAN program package widely used by the chemistry community. Obviously, many
other combinations of the exchange energy and the correlation energy are possible. The most popular
exchange-correlation functionals are composed of E xB88 for exchange, and EcPW or EcLY P for correlation.
These functionals provide a good compromise between computational efficiency and numerical
accuracy, but are not sufficient for a correct description of some chemical and physical properties.
Particularly, they do not yield significantly better results than LDA for solid state properties [54–57]
and they fail for the van der Waals dispersion interactions [58–60]. However, promising exchange and
correlation functionals have been developed by including additional information depending explicitly
on higher order density gradients [61, 62], i.e., the second derivative of the density.
8. A Classification of the Approximations and Main DFT Programs
To rank the approximated functionals, Perdew [63] introduced a DFT Jacob’s ladder, as a guideline
from the most approximate to the more sophisticated approaches for the functionals. The ladder has
five rungs corresponding to different choices in which way exactly do the energy density n(r) and its
derivatives enter. In principle, when climbing the DFT ladder progressively, higher overall accuracy
can be expected until exact results are reached. However, as with the ladder, it does not necessarily
mean higher accuracy for every particular system and property but rather a smaller probability to fail.
The first rung stands for the local spin density (LSD) approximation, proposed initially by Kohn
and Sham [9] and using only the energy density nα (α being used to label the spin components).
The second rung corresponds to the GGA (generalized gradient approximation) using ∇nα . It is
an extension of the LSD that is widely used in quantum chemistry while the LSD is the most popular
approximation for calculating the electronic structure in solid state. Two different philosophies may
be adopted. The first one searches for exchange-correlation functionals that obey physical constraints,
scaling properties and sum rules. The second one seeks for a parametric form where the parameters
are fitted to a set of experiments. They fall into the category of semiempirical xc functionals.
The third rung is the meta-GGA including the higher order density gradients ∇2 nα . It requires
an explicit dependence on the kinetic energy density, and a combination of theoretical constraints
satisfying the exact function and functionals that use few parameters fitting properties measured or
calculated by other quantum chemistry methods. Several meta-GGA functionals have been constructed
using the orbital kinetic energy densities, giving greater accuracy over a wide range of systems and
properties [64–66]. Concise overview of these functionals are given by Perdew et al. [52].
The fourth rung is representative of the hyper-GGA making use of the exact exchange energy
E exact
x and of a correlation energy functional Ec . A notable improvement of the exchange-correlation
GGAs appeared when Becke [67, 68] proposed to mix some fraction of exchange and correlation
GGA functionals with the exact exchange energy E exact x . This approach allowed to generate a new
class of functionals, known as the hybrid ones, the simplest among which may read:
xc = aE x
E hyb + (1 − a)E GGA + EcGGA ,
exact
x (119)
where E exact
x is the exact exchange energy obtained from Hartree-Fock method, E GGA
x the exchange
GGA
GGA functional and Ec the correlation GGA functional. Note that a(∼ 0.25) differs from 1 because

1398
the exact exchange energy has a highly nonlocal character compared to the GGA exchange. It should
be mentioned that neither the Hartree method nor the Hartree-Fock method are able to calculate the
exact correlation energy. Expressions such as Eq. (119) are implemented in many DFT code packages
and give almost identical results in specific applications. These explicitly nonlocal, exact exchange
based functionals look very promising in use for chemical physics.
The fifth rung or fully nonlocal description refers to exact exchange combined with exact partial
correlation, making use of the occupied and unoccupied Kohn-Sham orbitals. The strategy amounts
to identify the part of the exchange-correlation energy for which the GGA is accurate and to treat
the remainder exactly [69, 70]. As an alternative to the standard Kohn-Sham procedure, other exact
realizations of density functional theory may be carried out. Such generalized Kohn-Sham schemes
minimize not only the kinetic energy for fixed density, as in the standard Kohn-Sham method, but also
part of the electron-electron energy [71, 72] with the result of gaining more flexibility.
While the performances of hybrid GGA and hybrid meta-GGA functionals are traditionally better
than just plain GGA, and particularly better than LDA, results vary from functional to functional. The
performances of each functional lined up on the DFT ladder depend on the type of system under
consideration and on the property under evaluation. The number of density functionals currently
available is very large and, among the plethora of density functionals, the difficulty lies in selecting
the best one for each specific application. There is a great number of density functionals implemented
in the DFT codes∗ . Limited comparisons of the performance for density functionals can be found in
the literature. For instance, the performances of a large number of DFT functionals have been
evaluated across about fifty different density functionals [73] in the calculation of several properties of
chemical interest against experimental data (structure, kinetics, thermochemistry, nonbounded
interactions). A rigorous extension of the theory to these classes of properties is possible, but it
requires the introduction of new functionals that depend not only on density but also on properties.
In this paragraph, we briefly present few ubiquitous hybrid functionals employed in the currently
available DFT program packages, for the calculation of a large variety of molecular and chemical
properties. The most popular hybrid functionals in computational chemistry are the so-called E xc B3PW91
,
B3LY P B3PBE
E xc and E xc functionals.
B3PW91
(i) In the E xc functional, Eq. (108) is used for exchange and Eq. (111) for correlation, both
mixed via the 3 parameters a, b and c, according to the relationship:
B3PW91
E xc = E xc
LDA
+ a(E exact
x − E xLDA ) + b(E xB88 − E xLDA ) + c(EcPW91 − EcLDA ), (120)
with a = 0.20, b = 0.72 and c = 0.81. The weights of the various contributions to the exchange and
correlation energies are obtained by a fit either of the atomization energy of reference data [74] or of
theoretical constraints [75].
B3LY P
(ii) The E xc functional employs the same expression as Eq. (120) with the same parameters a, b
and c, whereas the EcPW91 correlation functional is replaced by the EcLY P one. It should be noted that this
type of functional is less successful to predict the solid state properties than the chemical properties,
because of the difficulties in computing the exact (Hartree-Fock) exchange with a plane waves basis
set.
∗
Here are a few official sites of distributed code packages that the reader may consult for calculations of quantum chemistry:
ADF (http://www.scl.com/); CPMD (http://www.cpmd.org/); GAUSSIAN (http://www.gaussian.com/) and quantum physics: ABINIT
(http://www.abinit.org/); VASP (http://www.vasp.at/); WIEN (http://www.wien2k.at/); Quantum ESPRESSO (Giannozzi P, et al. (2009)
J Phys-Condens Mat 21: 395502.)

1399
B3BPE
(iii) The E xc functional is the combination of the E xB88 functional for the exchange with the EcPBE
one for the correlation, with a = 0.25, b = 0.75 and c = 1. The value of 25% of the exact exchange
energy, instead of 20%, is supported by a strong theoretical basis [76, 77]. These code packages are
flexible enough to generate new density functionals from any combination of exchange and correlation
functionals with a fraction of the Hartree-Fock exchange energy (HF/DFT).
Standard DFT is particularly efficient to study systems where the origin of interactions is
essentially electrostatic. But the interest of DFT is also to treat systems, such as biological systems,
complex materials, surfaces, ..., where the dispersion interactions play a crucial role. It is a matter of
fact, that the van der Waals dispersion interaction is highly problematic for conventional functionals,
because no known correlation functional could generate the London dispersion interactions (∼ −r−6 )
in contrast to the asymptotic interactions falling off exponentially. The treatment of the dispersion
interactions within DFT is in full expansion according two principal routes. The first one, generally
known as DFT-D, consists in adding a more or less empirical dispersion correction, −C6 /r−6 , that
would account for the missing long range attraction. In this approach, the dispersion coefficient C6 is
either fitted on experimental ionization potentials and polarizabilities of isolated atoms [78] or
obtained with a non empirical approach making C6 environment dependent [79–81]. The second
route, known as vdW-DF, describes the dispersion interaction through a DFT functional by
introducing a nonlocal correlation functional [82, 83]. It is a promising computational tool for weakly
bound systems [84]. Some versions of the vdW-DF are implemented in widely distributed DFT
software package. Compared to the first approach, the second one improves the precision of the
results but it is offset by the increase of computational time by about 50%.
DFT is a powerful tool to investigate the static properties of electronic systems (geometrical
structures, relative energies, ...). It is also a convenient tool to be used in conjunction with ab initio
molecular dynamics (AIMD) simulations for extended systems, as described by Car and
Parrinello [85]. Whereas classic MD describes trajectories of atoms as objects within an empirical
interaction potential, the AIMD computations are simulations that involve the motions of both nuclei
and electrons. For the latter, the DFT model is used, and the evolution of the kinetic energies of both
electrons and nuclei is observed in an identical way, through the use of fictitious masses. Plane waves
are usually employed for the description of valence orbitals, whereas pseudo-potential
approximations are used for taking into account the core electrons. The combination of Car-Parrinello
quantum mechanics/molecular mechanics (QM/MM) is certainly one of the most promising
theoretical tools available for theoretical chemistry. Its implementation for ab initio molecular
dynamics simulation of large systems has also been successfully applied to solid state and liquid state
physics, as well as to materials science.
B3LY P
For molecular properties, the semiempirical E xc functional is very successful in describing a
wide range of molecular properties. For periodic systems and crystalline calculations, nonempirical
HF/DFT functionals are necessary but they are computationally expensive, especially for systems with
metallic characteristics. The situation improves with the HSE functionals, developed by Heyd et al.
[86, 87], which use a screening technique to take advantage of the fast spatial decay of the short-
range Hartree-Fock exchange. Such functionals are based on the PBE exchange-correlation functional
[88]: the exchange energy term is split into short-range and long-range components and the HF long-
range is replaced by the PBE long-range. The results obtained with HSE applied to metals exhibit
B3LY P
significantly smaller errors than pure DFT calculations. It has been shown that the E xc functional

1400
yields unsatisfactory atomization energies of solids underestimated by approximately 17% and lattice
parameters overestimated by about 1%, when compared with those of PBE and HSE which perform
reasonably well for materials with localized electrons and free electron metals [89].
The theoretical description of matter as well as of many chemical, physical and biological processes
requires accurate methods for investigating atomic and molecular-scale interactions, whose origins and
perspectives have been widely discussed recently [20, 90, 91].
Acknowledgments
I am grateful to Andrei Postnikov for stimulating discussions and its very careful reading of the
manuscript.
Conflict of Interest
There is no conflict of interest regarding the publication of this article.
References
1. Hartree DR (1928) The Wave Mechanics of an Atom with a Non-Coulomb Central Field. Part
I. Theory and Methods. Mathematical Proceedings of the Cambridge Philosophical Society, 24:
89–110.
2. Fock V (1930) Näherungsmethode zur Lösung des quantenmechanischen Mehrkörperproblems. Z
Physik 61: 126–148.
3. Slater JC (1930) Note on Hartree’s Method. Phys Rev 35: 210–211.
4. Raimes S (1967) The Wave Mechanics of Electrons in Metals, Amsterdam: North-Holland.
5. Thomas LH (1927) The calculation of atomic fields. Mathematical Proceedings of the Cambridge
Philosophical Society, 23: 542–548.
6. Fermi E (1928) Eine Statistische Methode zur Bestimmung einiger Eigenschaften des Atoms und
ihre Anwendung auf die Theorie des periodischen Systems der Elemente. Z Physik 48: 73–79.
7. Dirac PAM (1930) Note on Exchange Phenomena in the Thomas Atom. Mathematical Proceedings
of the Cambridge Philosophical Society, 26: 376–385.
8. Hohenberg P, Kohn W (1964) Inhomogeneous Electron Gas. Phys Rev 136: B864–B871.
9. Kohn W, Sham LJ (1965) Self-Consistent Equations Including Exchange and Correlation Effects.
Phys Rev 140: A1133–A1138.
10. Becke AD (1996) Density-functional thermochemistry. IV. A new dynamical correlation
functional and implications for exact-exchange mixing. J Chem Phys 104: 1040–1046.
11. Becke AD (1988) Density-functional exchange-energy approximation with correct asymptotic
behavior. Phys Rev A 38: 3098–3100.
12. Lee C, Yang W, Parr RG (1988) Development of the Colle-Salvetti correlation-energy formula into
a functional of the electron density. Phys Rev B 37: 785–789.

1401
13. Feynman RP, Metropolis N, Teller E (1949) Equations of State of Elements Based on the
Generalized Fermi-Thomas Theory. Phys Rev 75: 1561–1573.
14. Schwinger J (1981) Thomas-Fermi model: the second correction. Phys Rev A 24: 2353–2361.
15. Shakeshalf R, Spruch L (1981) Remarks on the existence and accuracy of the Z −1/3 expansion of
the nonrelativistic ground-state energy of a neutral atom. Phys Rev A 23: 2118–2126.
16. Englert BG, Schwinger J (1985) Atomic-binding-energy oscollations. Phys Rev A 32: 47–63.
17. Spruch L (1991) Pedagogic notes on Thomas-Fermi theory (and on some improvements): atoms,
stars, and the stability of bulk matter. Rev Mod Phys 63: 151–209.
18. Mermin ND (1965) Thermal Properties of the Inhomogeneous Electron Gas. Phys Rev 137:
A1441–A1443. (Concerns the generalization to electrons at finite temperatures).
19. Levy M (1979) Universal variational functionals of electron densities, first-order density matrices,
and natural spin-orbitals and solution of the v-representability problem. Proc Natl Acad Sci U S A
76: 6062–6065.
20. Becke AD (2014) Perspective: Fifty years of density-functional theory in chemical physics. J
Chem Phys 140: 18A301.
21. Lindhard J (1954) On the properties of a gas of charged particles. Kgl Danske Videnskab Selskab
Mat Fys Medd 28: 8.
22. Jones RO, Young WH (1971) Density functional theory and the von Weizsacker method. J Phys C
4: 1322–1330.
23. von Weizsäcker CF (1935) Zur Theorie der Kernmassen. Z Physik 96: 431–458 .
24. Jones RO, Gunnarsson O (1989) The density functional formalism, its applications and prospects.
Rev Mod Phys 61: 689–746.
25. Slater JC (1951) A Simplification of the Hartree-Fock Method. Phys Rev 81: 385–390.
26. Robinson JE, Bassani F, Knox BS, et al. (1962) Screening Correction to the Slater Exchange
Potential. Phys Rev Lett 9: 215–217.
27. Wigner EP (1934) On the Interaction of Electrons in Metals. Phys Rev 46: 1002–1011.
28. Gell-Mann M, Brueckner K (1957) Correlation Energy of an Electron Gas at High Density. Phys
Rev 106: 364–368.
29. Ceperley DM (1978) Ground state of the fermion one-component plasma: A Monte Carlo study
in two and three dimensions. Phys Rev B 18: 3126–3138.
30. Ceperley DM, Alder BJ (1980) Ground State of the Electron Gas by a Stochastic Method. Phys
Rev Lett 45: 566–569.
31. Perdew JP, Wang Y (1992) Accurate and simple analytic representation of the electron-gas
correlation energy. Phys Rev B 45: 13244–13249.
32. Levy M (1982) Electron densities in search of Hamiltonians. Phys Rev A 26: 1200–1208.
33. Lieb EH (1983) Density functionals for coulomb-systems. Int J Quantum Chem 24: 243–277.
34. Levy M, Perdew JP (1985) Hellmann-Feynman, virial, and scaling requisites for the exact
universal density functionals. Shape of the correlation potential and diamagnetic susceptibility for
atoms. Phys Rev A 32: 2010–2021.

1402
35. Wigner E, Seitz F (1933) On the Constitution of Metallic Sodium. Phys Rev 43: 804–810.
36. Hellmann H (1933) Zur Rolle der kinetischen Elektronenenergie für die zwischenatomaren Kräfte.
Z Physik 85: 180–190.
37. Feynman RP (1939) Forces in Molecules. Phys Rev 56: 340–343.
38. Pupyshev VI (2000) The Nontriviality of the Hellmann-Feyman Theorem. Russ J Phys Chem 74:
S267–S278.
39. Perdew JP, Kurth S (2003) Density Functionals for Non-relativistic Coulomb Systems in the New
Century, In: Fiolhais C, Nogueira F, Marques MAL, A Primer in Density Functional Theory, Berlin
Heidelberg: Springer-Verlag, 1–55.
40. Gunnarsson O, Lundqvist BI (1976) Exchange and correlation in atoms, molecules, and solids by
the spin-density-functional formalism. Phys Rev B 13: 4274–4298; Gunnarsson O, Lundqvist BI
(1977) Erratum: Exchange and correlation in atoms, molecules, and solids by the spin-density-
functional formalism. Phys Rev B 15: 6006.
41. Gunnarsson O, Jonson M, Lundqvist BI (1979) Descriptions of exchange and correlation effects
in inhomogeneous electron systems. Phys Rev 20: 3136–3165.
42. Langreth DC, Perdew JP (1975) The Exchange-Correlation Energy of a Metallic Surface. Solid
State Commun 17: 1425–1429.
43. Langreth DC, Perdew JP (1977) Exchange-correlation energy of a metallic surface: Wave-vector
analysis. Phys Rev B 15: 2884–2901.
44. Jones RO, Gunnarsson O (1989) The density functional formalism, its applications and prospects.
Rev Mod Phys 61: 689–746.
45. Charlesworth JPA (1996) Weighted-density approximation in metals and semiconductors. Phys
Rev B 53: 12666–12673.
46. Gunnarsson O, Lundqvist BI (1976) Exchange and correlation in atoms, molecules, and solids by
the spin-density-functional formalism. Phys Rev B 13: 4274–4298.
47. Alvarellos JE, Chacon E, Tarazona P (1986) Nonlocal density functional for the exchange and
correlation energy of electrons. Phys Rev 33: 6579–6587.
48. Chacon E, Tarazona P (1988) Self-consistent weighted-density approximation for the electron gas.
I. Bulk properties. Phys Rev 37: 4020–4025.
49. Gupta AK, Singwi KS (1977) Gradient corrections to the exchange-correlation energy of electrons
at metal surfaces. Phys Rev B 15: 1801–1810.
50. Perdew JP, Burke K, Ernzerhof M (1996) Generalized Gradient Approximation Made Simple.
Phys Rev Lett 77: 3865–3868.
51. Perdew JP, Burke K, Wang Y (1996) Generalized gradient approximation for the exchange-
correlation hole of a many-electron system. Phys Rev B 54: 16533–16539.
52. Perdew JP, Ruzsinszky A, Csonka GI, et al. (2008) Restoring the Density-Gradient Expansion
for Exchange in Solids and Surfaces. Phys Rev Lett 100: 136406. (Modified PBE subroutines are
available from: http://dft.uci.edu./pubs/PBEsol.html.)

1403
53. Colle R, Salvetti O (1975) Approximate calculation of the correlation energy for the closed shells.
Theoret Chim Acta 37: 329–334.
54. Leung TC, Chan CT, Harmon BN (1991) Ground-state properties of Fe, Co, Ni, and their
monoxides: Results of the generalized gradient approximation. Phys Rev B 44: 2923–2927.
55. Singh DJ, Ashkenazi J (1992) Magnetism with generalized-gradient-approximation density
functionals. Phys Rev B 46: 11570–11577.
56. Körling M, Häglund J (1992) Cohesive and electronic properties of transition metals: The
generalized gradient approximation. Phys Rev B 45: 13293–13297.
57. Engel E, Vosko SH (1994) Fourth-order gradient corrections to the exchange-only energy
functional: Importance of ∇2 n contribution. Phys Rev B 50: 10498–10505.
58. Andersson Y, Langreth DC, Lundqvist BI (1996) van der Waals Interactions in Density-Functional
Theory. Phys Rev Lett 76: 102–105.
59. Dobson JF, Dinte BP (1996) Constraint Satisfaction in Local and Gradient Susceptibility
Approximations: Application to a van der Waals Density Functional. Phys Rev Lett 76: 1780–
1783.
60. Patton DC, Pederson MR (1997) Application of the generalized-gradient approximation to rare-gas
dimers. Phys Rev A 56: R2495–R2498.
61. Van Voorhis T, Scuseria GE (1998) A novel form for the exchange-correlation energy functional.
J Chem Phys 109: 400–410.
62. Tao J, Perdew JP, Staroverov VN, et al. (2003) Climbing the Density Functional Ladder:
Nonempirical Meta–Generalized Gradient Approximation Designed for Molecules and Solids.
Phys Rev Lett 91: 146401.
63. Van Doren VE, Van Alsenoy K, Greelings P (2001) Density Functional Theory and Its
Applications to Materials, Melville, NY: American Institute of Physics.
64. Gonis A, Kioussis N (1999) Electron Correlations and Materials Properties, New York: Plenum.
65. Perdew JP, Kurth S, Zupan A, et al. (1999) Accurate Density Functional with Correct Formal
Properties: A Step Beyond the Generalized Gradient Approximation. Phys Rev Lett 82: 2544–
2547; Perdew JP, Kurth S, Zupan A, et al. (1999) Erratum: Accurate Density Functional with
Correct Formal Properties: A Step Beyond the Generalized Gradient Approximation. Phys Rev
Lett 82: 5179.
66. Becke AD (2000) Simulation of delocalized exchange by local density functionals. J Chem Phys
112: 4020–4026.
67. Becke AD (1993) A new mixing of Hartree–Fock and local density-functional theories. J Chem
Phys 98: 1372–1377.
68. Becke AD (1993) Density-functional thermochemistry. III. The role of exact exchange. J Chem
Phys 98: 5648–5652.
69. Seidl M, Perdew JP, Kurth S (2000) Simulation of All-Order Density-Functional Perturbation
Theory, Using the Second Order and the Strong-Correlation Limit. Phys Rev Lett 84: 5070–5073.

1404
70. Fuchs M, Gonze X (2002) Accurate density functionals: Approaches using the adiabatic-
connection fluctuation-dissipation theorem. Phys Rev B 65: 235109.
71. Seidl A, Görling A, Vogl P, et al. (1996) Generalized Kohn-Sham schemes and the band-gap
problem. Phys Rev B 53: 3764–3774.
72. Savin A (1996) On degeneracy, near degeneracy and density functional theory, In: Seminario JM,
Recent Developments of Modern Density Functional Theory, Amsterdam: Elsevier, 327–357.
73. Sousa SF, Fernandes PA, Ramos MJ (2007) General Performance of Density Functionals. J Phys
Chem A 111: 10439–10452.
74. Curtiss LA, Raghavachari K, Redfern PC, et al. (1997) Assessment of Gaussian-2 and density
functional theories for the computation of enthalpies of formation. J Chem Phys 106: 1063–1079.
(This paper contains a set of 148 molecules having well-established enthalpies of formation at 298
K.)
75. Chan GKL, Handy NC (1999) Optimized Lieb-Oxford bound for the exchange-correlation energy.
Phys Rev A 59: 3075–3077. (As an example, the exchange-correlation energy must satisfy the
inegality |E xc | ≤ 2.27 E xLDA .)
76. Perdew JP, Ernzerhof M, Burke K (1996) Rationale for mixing exact exchange with density
functional approximations. J Chem Phys 105: 9982–9985.
77. Adamo C, Barone V (1998) Exchange functionals with improved long-range behavior and
adiabatic connection methods without adjustable parameters: The mPW and mPW1PW models.
J Chem Phys 108: 664–675.
78. Wu Q, Yang W (2002) Empirical correction to density functional theory for van der Waals
interactions. J Chem Phys 116: 515–524.
79. Becke AD, Johnson ER (2005) Exchange-hole dipole moment and the dispersion interaction. J
Chem Phys 122: 154104.
80. Tkatchenko A, Scheffler M (2009) Accurate Molecular Van Der Waals Interactions from Ground-
State Electron Density and Free-Atom Reference Data. Phys Rev Lett 102: 073005.
81. Grimme S, Antony J, Ehrlich S, et al. (2010) A consistent and accurate ab initio parametrization
of density functional dispersion correction (DFT-D) for the 94 elements H-Pu. J Chem Phys 132:
154104.
82. Dion M, Rydberg H, Schröder E, et al. (2004) Van der Waals Density Functional for General
Geometries. Phys Rev Lett 92: 246401.
83. Thonhauser T, Cooper VR, Li S, et al. (2007) Van der Waals density functional: Self-consistent
potential and the nature of the van der Waals bond. Phys Rev B 76: 125112 .
84. Klimes J, Michaelides A (2012) Perspective: Advances and challenges in treating van der Waals
dispersion forces in density functional theory. J Chem Phys 137: 120901.
85. Car R, Parrinello M (1985) Unified Approach for Molecular Dynamics and Density-Functional
Theory. Phys Rev Lett 55: 2471–2474.
86. Heyd J, Scuseria GE, Ernzerhof M (2003) Hybrid functionals based on a screened Coulomb
potential. J Chem Phys 118: 8207–8215.

1405
87. Heyd J, Scuseria GE (2004) Efficient hybrid density functional calculations in solids: Assessment
of the Heyd–Scuseria–Ernzerhof screened Coulomb hybrid functional. J Chem Phys 121: 1187–
1192.
88. Ernzerhof M, Scuseria GE (1999) Assessment of the Perdew–Burke–Ernzerhof exchange-
correlation functional. J Chem Phys 110: 5029–5036.
89. Paier J, Marsman M, Kresse G (2007) Why does the B3LYP hybrid functional fail for metals? J
Chem Phys 127: 024103.
90. Burke K (2012) Perspective on density functional theory. J Chem Phys 136: 150901.
91. Jones RO (2015) Density functional theory: Its origins, rise to prominence, and future. Rev Mod
Phys 87: 897–923.
c 2017, Jean-Louis Bretonnet, licensee AIMS Press.

This is an open access article distributed under the
terms of the Creative Commons Attribution License
(http://creativecommons.org/licenses/by/4.0)

Bretonnet Rev Basics

Uploaded by

Copyright:

Available Formats

Bretonnet Rev Basics

Uploaded by

Document Information

Original Description:

Original Title

Copyright

Available Formats

Share this document

Share or Embed Document

Sharing Options

Did you find this document useful?

Is this content inappropriate?

Copyright:

Available Formats

Bretonnet Rev Basics

Uploaded by

Copyright:

Available Formats

AIMS Materials Science, 4(6): 1372-1405.

or under the compact form:

AIMS Materials Science Volume 4, Issue 6, 1372-1405.

AIMS Materials Science Volume 4, Issue 6, 1372-1405.

2. Atomic Model of Thomas and Fermi

AIMS Materials Science Volume 4, Issue 6, 1372-1405.

introduction to the density-functional theory.

AIMS Materials Science Volume 4, Issue 6, 1372-1405.

Combining Eq. (16) and Eq. (17), we find:

AIMS Materials Science Volume 4, Issue 6, 1372-1405.

3. Summary of the Density-Functional Theory

3.1. Formulation of Hohenberg and Kohn

AIMS Materials Science Volume 4, Issue 6, 1372-1405.

However, the last term of this inequality reads:

The same reasoning applied to Eq. (28) leads to:

E (1) + E (2) < E (1) + E (2) , (33)

AIMS Materials Science Volume 4, Issue 6, 1372-1405.

3.2. Formulation of Kohn and Sham

AIMS Materials Science Volume 4, Issue 6, 1372-1405.

vKS (r) = wH (r) + v xc (r) + Vext (r), (42)

AIMS Materials Science Volume 4, Issue 6, 1372-1405.

4. The TF Equation from the KS Approach

AIMS Materials Science Volume 4, Issue 6, 1372-1405.

T0 3}2  2 2/3 2/3 1.105

E xc [n(r)] = E x [n(r)] + Ec [n(r)] , (50)

AIMS Materials Science Volume 4, Issue 6, 1372-1405.

5. Improved Functional Approximations in DFT

n(r) = n0 + n1 (r), (55)

5.1. Without Exchange-correlation Effects

AIMS Materials Science Volume 4, Issue 6, 1372-1405.

wH (q) + Vext (q)

AIMS Materials Science Volume 4, Issue 6, 1372-1405.

γ π2 ~2 n21 (r) λ ~2 |∇r n(r)|2

5.2. With Exchange-correlation Effects

AIMS Materials Science Volume 4, Issue 6, 1372-1405.

AIMS Materials Science Volume 4, Issue 6, 1372-1405.

may be carried out. The final result is∗ :

AIMS Materials Science Volume 4, Issue 6, 1372-1405.

ec (n) = + 0.031 ln r s − 0.048 . (82)

with c0 = 0.03109, α = 0.21370, β1 = 7.5957, β2 = 3.5876, β3 = 1.6382, β4 = 0.49294.

AIMS Materials Science Volume 4, Issue 6, 1372-1405.

6. Beyond the Local Density Approximation

6.1. Average Exchange Hole

so that the total exchange energy reads:

e ψ j (r2 )ψi (r2 )ψ j (r1 )

AIMS Materials Science Volume 4, Issue 6, 1372-1405.

AIMS Materials Science Volume 4, Issue 6, 1372-1405.

6.2. Exchange-correlation Energy Functional

AIMS Materials Science Volume 4, Issue 6, 1372-1405.

Of course, W H [n] is given by Eq. (34).

n(r1 )n2 (r1 , r2 )

AIMS Materials Science Volume 4, Issue 6, 1372-1405.

7. Separation of the Exchange and Correlation Functionals

This relation is known as the generalized-gradient approximation (GGA); the function

AIMS Materials Science Volume 4, Issue 6, 1372-1405.

7.1. Exchange Energy Functionals

7.2. Correlation Energy Functionals

AIMS Materials Science Volume 4, Issue 6, 1372-1405.

T0 3}2 2 2/3 2/3 1.105