See discussions, stats, and author profiles for this publication at: https://www.researchgate.

net/publication/30405238

Density Functional Theory: Basics, New Trends and Applications

Article · January 2003

Source: OAI

CITATIONS READS

52 6,567

2 authors, including:

Nikitas I. Gidopoulos
Durham University
78 PUBLICATIONS 1,095 CITATIONS

SEE PROFILE

All content following this page was uploaded by Nikitas I. Gidopoulos on 21 October 2014.

The user has requested enhancement of the downloaded file.

Density Functional Theory: Basics, New Trends and
Applications
J. Kohanoff and N.I. Gidopoulos
Volume 2, Part 5, Chapter 26, pp 532–568

in

Handbook of Molecular Physics and Quantum Chemistry

(ISBN 0 471 62374 1)

Edited by

Stephen Wilson

 John Wiley & Sons, Ltd, Chichester, 2003

Chapter 26
Density Functional Theory: Basics, New Trends and
Applications

J. Kohanoff1 and N.I. Gidopoulos2

1
Queen’s University Belfast, Belfast, Northern Ireland
2
Rutherford Appleton Laboratory, Oxfordshire, UK

describe the system by a number of nuclei and electrons

1 The Problem of the Structure of Matter 1 interacting through coulombic (electrostatic) forces. For-
2 The Electronic Problem 3 mally, we can write the Hamiltonian of such a system in
the following general form:
3 Density Functional Theory 4
4 Exchange and Correlation 10  h̄2 2  h̄2 2
P N
5 Exact Exchange: The Optimized Potential
=−
H ∇I − ∇
I =1
2MI 2m i
Method 19 i=1

6 Towards an Accurate Correlation Functional 23 e2   ZI ZJ e2  

P P N N
1
7 Comparison and Salient Features of the + +
2 I =1 J =I |RI − RJ | 2 i=1 j =i |ri − rj |
Different Approximations 27
Notes 35 
P 
N
ZI
References 35 − e2 (1)
I =1 i=1
|RI − ri |

where R = {RI }, I = 1, . . . , P , is a set of P nuclear

1 THE PROBLEM OF THE STRUCTURE coordinates and r = {ri }, i = 1, . . . , N , is a set of N elec-
tronic coordinates. ZI and MI are the P nuclear charges
OF MATTER and masses, respectively. Electrons are fermions, so that
the total electronic wave function must be antisymmetric
The microscopic description of the physical and chemical with respect to exchange of two electrons. Nuclei can be
properties of matter is a complex problem. In general, we fermions, bosons or distinguishable particles, according to
deal with a collection of interacting atoms, which may also the particular problem under examination. All the ingredi-
be affected by some external field. This ensemble of par-
ents are perfectly known and, in principle, all the properties
ticles may be in the gas phase (molecules and clusters) or
can be derived by solving the many-body Schrödinger
in a condensed phase (solids, surfaces, wires), they could
equation:
be solids, liquids or amorphous, homogeneous or hetero-
geneous (molecules in solution, interfaces, adsorbates on

(r, R) = E
(r, R)
H (2)
i i i
surfaces). However, in all cases we can unambiguously
In practice, this problem is almost impossible to treat in a
Handbook of Molecular Physics and Quantum Chemistry, full quantum-mechanical framework. Only in a few cases a
Edited by Stephen Wilson. Volume 2: Molecular Electronic Struc- complete analytic solution is available, and numerical solu-
ture.  2003 John Wiley & Sons, Ltd. ISBN: 0-471-62374-1. tions are also limited to a very small number of particles.
2 Electronic structure of large molecules

There are several features that contribute to this difficulty. where the electronic wave function m (R, r) [m (R, r) is
First, this is a multicomponent many-body system, where normalized for every R] is the mth stationary state of the
each component (each nuclear species and the electrons) electronic Hamiltonian
obeys a particular statistics. Second, the complete wave
function cannot be easily factorized because of coulombic ĥe = T
e + U

+V
ee

=H
ne

− T
− U
n


nn (4)
correlations. In other words, the full Schrödinger equation
cannot be easily decoupled into a set of independent equa- T
n and U
are the kinetic and potential nuclear oper-
nn
tions so that, in general, we have to deal with (3P + 3N ) ators, T
e and U
the same for electrons, and V
ee

the
ne
coupled degrees of freedom. The dynamics is an even more electron–nuclear interaction. The corresponding eigenvalue
difficult problem, and very few and limited numerical tech- is noted εm (R). In the electronic (stationary) Schrödinger
niques have been devised to solve it. The usual choice is to equation, the nuclear coordinates R enter as parameters,
resort to some sensible approximations. The large majority while the nuclear wave function m (R, t) obeys the time-
of the calculations presented in the literature are based on dependent Schrödinger equation
(i) the adiabatic separation of nuclear and electronic degrees ∂m (R, t) 

of freedom (adiabatic approximation) and (ii) the classical i h̄
+ ε (R)  (R, t)
= Tn + U (5)
nn m m
∂t
treatment of the nuclei.
or the stationary version
 
1.1 Adiabatic approximation
+ ε (R)  (R) = E  (R)
T
n + U nn m m m m (6)
(Born–Oppenheimer)
In principle, m can be any electronic eigenstate. In practice,
however, most of the applications in the literature are
The first observation is that the timescale associated to focused on the ground state (m = 0).
the motion of the nuclei is usually much slower than
that associated to electrons. In fact, the small mass of
the electrons as compared to that of the protons (the 1.2 Classical nuclei approximation
most unfavourable case) is about 1 in 1836, meaning
that their velocity is much larger. In this spirit, it was Solving any of the two last equations (5) or (6) is a
proposed in the early times of quantum mechanics that formidable task for two reasons: First, it is a many-body
the electrons can be adequately described as following equation in the 3P nuclear coordinates, the interaction
instantaneously the motion of the nuclei, staying always in potential being given in an implicit form. Second, the deter-
the same stationary state of the electronic Hamiltonian.(1) mination of the potential energy surface εn (R) for every
This stationary state will vary in time because of the possible nuclear configuration R involves solving M 3P
coulombic coupling of the two sets of degrees of freedom times the electronic equation, where M is, for example, a
but if the electrons were, for example, in the ground state, typical number of grid points. The largest size achieved up
they will remain there forever. This means that as the to date using non-stochastic methods is six nuclear degrees
nuclei follow their dynamics, the electrons instantaneously of freedom.
adjust their wave function according to the nuclear wave In a large variety of cases of interest, however, the solu-
function. tion of the quantum nuclear equation is not necessary. This
This approximation ignores the possibility of having is based on two observations: (i) The thermal wavelength
√
non-radiative transitions between different electronic eigen- for a particle of mass M is λ̄T = h̄/ MkB T , so that regions
states. Transitions can only arise through coupling with an of space separated by more than λT do not exhibit quan-
external electromagnetic field and involve the solution of tum phase coherence. The least favourable case is that of
the time-dependent Schrödinger equation. This has been hydrogen, and even so, at room temperature λ̄T ≈ 0.4 Å,
achieved, especially in the linear response regime, but also while inter-atomic distances are normally of the order of
in a non-perturbative framework in the case of molecules 1 Å. (ii) Potential energy surfaces in typical bonding envi-
in strong laser fields. However, this is not the scope of this ronments are normally stiff enough to localize the nuclear
section, and electronic transitions will not be addressed in wave functions to a large extent. For instance, a proton in
the following. a hydroxyl group has a width of about 0.25 Å.
Under the above conditions, the full wave function fac- This does not mean that quantum nuclear effects can be
torizes in the following way: neglected altogether. In fact, there is a variety of questions
in condensed matter and molecular physics that require a

(R, r, t) = m (R, t)m (R, r) (3) quantum-mechanical treatment of the nuclei. Well-known
 Density functional theory: Basics, new trends and applications 3

examples are the solid phases of hydrogen, hydrogen- viscosity, ionic conductivity and so forth. Excited electronic
bonded systems such as water and ice, fluxional molecules states (or the explicit time dependence) also give access to
and even active sites of enzymes. There is, however, an another wealth of measurable phenomena such as electronic
enormous number of systems where the nuclear wave pack- transport and optical properties.
ets are sufficiently localized to be replaced by Dirac’s
δ-functions. The centres of these δ-functions are, by def-
inition, the classical positions Rcl . 2.1 Quantum many-body theory: chemical
The connection between quantum and classical mechan- approaches
ics is achieved through Ehrenfest’s theorem for the mean
values of the position and momentum operators.(2) The The first approximation may be considered the one pro-
quantum-mechanical analog of Newton’s equations is posed by Hartree (as early as in 1928, in the very beginning
of the age of quantum mechanics).(3) It consists of postu-
d2 RI  lating that the many-electron wave function can be written
MI = −∇RI εn (R) (7) as a simple product of one-electron wave functions. Each
dt 2
of these verifies a one-particle Schrödinger equation in an
where the brackets indicate quantum expectation values. effective potential that takes into account the interaction
The classical nuclei approximation consists of identifying with the other electrons in a mean-field way (we omit the
RI  with Rcl
I . In this case, the nuclear wave function is dependence of the orbitals on R):
represented by a product of δ-functions, then ∇εm (R) =
∇εm (Rcl ). The latter is strictly valid only for δ-functions (R, r) = i ϕi (ri ) (8)
or for harmonic potentials. In the general case, the leading  
error of this approximation is proportional to the anhar- h̄2 2 (i)
− ∇ + Veff (R, r) ϕi (r) = εi ϕi (r) (9)
monicity of the potential and to the spatial extension of the 2m
wave function.
Assuming these two approximations, we are then left with
with the problem of solving the many-body electronic 
N
ρj (r )
Schrödinger equation for a set of fixed nuclear positions. 
j  =i
This is a major issue in quantum mechanics, and we shall (i)
Veff (R, r) = V (R, r) + dr (10)
devote the remainder of this chapter to it. |r − r |

where
2 THE ELECTRONIC PROBLEM ρj (r) = |ϕj (r)|2 (11)

The key problem in the structure of matter is to solve is the electronic density associated with particle j . The
the Schrödinger equation for a system of N interacting second term in the right-hand side (rhs) of equation (10) is
electrons in the external coulombic field created by a the classicalelectrostatic potential generated by the charge
collection of atomic nuclei (and may be some other external distribution jN=i ρj (r). Notice that this charge density does
field). It is a very difficult problem in many-body theory not include the charge associated with particle i, so that the
and, in fact, the exact solution is known only in the case Hartree approximation is (correctly) self-interaction-free. In
of the uniform electron gas, for atoms with a small number this approximation, the energy of the many-body system is
of electrons and for a few small molecules. These exact not just the sum of the eigenvalues of equation (9) because
solutions are always numerical. At the analytic level, one the formulation in terms of an effective potential makes
always has to resort to approximations. the electron–electron interaction to be counted twice. The
However, the effort of devising schemes to solve this correct expression for the energy is
problem is really worthwhile because the knowledge of 

N
1
N
ρi (r)ρj (r )
the electronic ground state of a system gives access to EH = εn − dr dr (12)
many of its properties, for example, relative stability n=1
2 i=j |r − r |
of different structures/isomers, equilibrium structural
information, mechanical stability and elastic properties, The set of N coupled partial differential equations (9)
pressure–temperature (P-T) phase diagrams, dielectric can be solved by minimizing the energy with respect to
properties, dynamical (molecular or lattice) properties a set of variational parameters in a trial wave function
such as vibrational frequencies and spectral functions, or, alternatively, by recalculating the electronic densities
(non-electronic) transport properties such as diffusivity, in equation (11) using the solutions of equation (9), then
4 Electronic structure of large molecules

casting them back into the expression for the effective Møller–Plesset perturbation theory of second (MP2) or
potential (equation 10), and solving again the Schrödinger fourth (MP4) order,(6) or by configuration interaction (CI)
equation. This procedure can be repeated several times, methods using a many-body wave function made of a linear
until self-consistency in the input and output wave function combination of Slater determinants, as a means for intro-
or potential is achieved. This procedure is called self- ducing electronic correlations. Several CI schemes have
consistent Hartree approximation. been devised during the past 40 years, and this is still an
The Hartree approximation treats the electrons as dis- active area of research. Coupled clusters (CC) and complete
tinguishable particles. A step forward is to introduce active space (CAS) methods are currently two of the most
Pauli exclusion principle (Fermi statistics for electrons) by popular ones.(7,8)
proposing an antisymmetrized many-electron wave function Parallel to the development of this line in electronic
in the form of a Slater determinant: structure theory, Thomas and Fermi proposed, at about the
same time as Hartree (1927–1928), that the full electronic
(R, r) density was the fundamental variable of the many-body
= SD{ϕj (ri , σi )} problem and derived a differential equation for the den-
sity without resorting to one-electron orbitals.(9,10) The
ϕ1 (r1 , σ1 ) ϕ1 (r2 , σ2 ) ··· ϕ1 (rN , σN ) Thomas–Fermi (TF) approximation was actually too crude
1 ϕ2 (r1 , σ1 ) ϕ2 (r2 , σ2 ) ··· ϕ2 (rN , σN ) because it did not include exchange and correlation effects
=√ .. .. .. .. and was also unable to sustain bound states because of the
N! . . . .
ϕN (r1 , σ1 ) ϕj (r2 , σ2 ) · · · ϕN (rN , σN ) approximation used for the kinetic energy of the electrons.
However, it set up the basis for the later development of
(13) density functional theory (DFT), which has been the way
This wave function introduces particle exchange in an exact of choice in electronic structure calculations in condensed
manner.(4,5) The approximation is called Hartree–Fock matter physics during the past 20 years and recently, it
(HF) or self-consistent field (SCF) approximation and has also became accepted by the quantum chemistry commu-
been for a long time the way of choice of chemists for nity because of its computational advantages compared to
calculating the electronic structure of molecules. In fact, HF-based methods [1].
it provides a very reasonable picture for atomic systems
and, although many-body correlations (arising from the
fact that, owing to the two-body Coulomb interactions, the
total wave function cannot necessarily be written as an 3 DENSITY FUNCTIONAL THEORY
antisymmetrized product of single-particle wave functions)
are completely absent, it also provides a reasonably good The total ground state energy of an inhomogeneous system
description of inter-atomic bonding. HF equations look the composed by N interacting electrons is given by
same as Hartree equations, except for the fact that the
exchange integrals introduce additional coupling terms in E = |T
+ V

+U

|
ee
the differential equations:
= |T
| + |V

| + |U

|
ee
 

N
 ρj (r , σ )  where | is the N -electron ground state wave function,
 h̄2  
 σ ,j =1  which has neither the form given by the Hartree approxi-
− ∇ + V (R, r) +
2
dr  ϕi (r, σ)
 2m |r − r |  mation (8) nor the HF form (13). In fact, this wave func-
 
tion has to include correlations amongst electrons, and its
general form is unknown. T
is the kinetic energy, V
is
 

N   ϕj∗ (r , σ )ϕi (r , σ )

the interaction with external fields, and Uee is the elec-
− dr ϕj (r, σ) tron–electron interaction. We are going to concentrate now
j =1 σ
|r − r |
on the latter, which is the one that introduces many-body

N effects.
= λij ϕj (r, σ) (14)

N 
N
j =1

= |U
U
| = | 1 1
|
ee ee
2 i=1 j  =i
|ri − rj |
Notice that also in HF the self-interaction cancels exactly.

Nowadays, the HF approximation is routinely used as ρ2 (r, r )
= dr dr (15)
a starting point for more elaborated calculations like |r − r |
 Density functional theory: Basics, new trends and applications 5

with In HF theory (equation 13) we can rewrite the elec-

tron–electron interaction as
1
ρ2 (r, r ) = |
σ† (r)
σ† (r )
σ (r )
σ (r)| (16) 
2 σ,σ 1 ρHF (r)ρHF (r )
HF
Uee = dr dr
2 |r − r |
the two-body density matrix expressed in real space,   
being
and
† the creation and annihilation operators  HF |ρHF
σ (r, r )|2
for electrons, which obey the anti-commutation relations 1 ρ (r)ρHF (r )  − σ

 dr dr
+ 
{
σ (r),
σ† (r )} = δσ,σ δ(r − r ). We define now the two- 2 
|r − r | ρ (r)ρ (r ) 
HF HF 

body direct correlation function g(r, r ) in the following

way: (22)
meaning that the exact expression for the exchange deple-
ρ2 (r, r ) = 12 ρ(r, r)ρ(r , r ) g(r, r ) (17)
tion (also called exchange hole) is
where ρ(r, r ) is the one-body density matrix (in real 
 2
space), whose diagonal elements ρ(r) = ρ(r, r) correspond |ρHF
σ (r, r )|
to the electronic density. The one-body density matrix is σ
gX (r, r ) = 1 − (23)
defined as ρHF (r)ρHF (r )

ρ(r, r ) = ρσ (r, r ) (18) The density and density matrix are calculated from the HF
σ ground state Slater determinant.
ρσ (r, r ) = |
σ† (r)
σ (r )| (19) The calculation of the correlation hole – gC (r, r ) – is a
major problem in many-body theory and, up to the present,
With this definition, the electron–electron interaction is it is an open problem in the general case of an inhomoge-
written as neous electron gas. The exact solution for the homogeneous
 electron gas is known numerically(11,12) and also in a num-
1 ρ(r)ρ(r )
Uee = dr dr ber of different analytic approximations (see below). There
2 |r − r | are several approximations that go beyond the homoge-

1 ρ(r)ρ(r ) neous limit by including slowly varying densities through
+ [g(r, r ) − 1] dr dr (20)
2 |r − r | its spatial gradients (gradient corrections) and also expres-
sions for the exchange-correlation energy that aim at taking
The first term is the classical electrostatic interaction energy into account very weak, non-local interactions of the van
corresponding to a charge distribution ρ(r). The second der Waals type (dispersion interactions).(13)
term includes correlation effects of both classical and quan- The energy of the many-body electronic system can, then,
tum origin. Basically, g(r, r ) takes into account the fact be written in the following way:
that the presence of an electron at r discourages a second
electron to be located at a position r very close to r because 
1 ρ(r)ρ(r )
of the Coulomb repulsion. In other words, it says that E =T +V + dr dr + EXC (24)
the probability of finding two electrons (two particles with 2 |r − r |
where
charges of the same sign, in the general case) is reduced
with respect to the probability of finding them at infinite

P 
N P 

distance. This is true already at the classical level and it V = | v(ri − RI )| = ρ(r)v(r − RI ) dr
is further modified at the quantum level. Exchange further I =1 i=1 I =1
diminishes this probability in the case of electrons having (25)
the same spin projection, owing to the Pauli exclusion.
2 
N 
To understand the effect of exchange, let us imagine h̄ h̄2  2 
T = | − ∇i2 | = − ∇r ρ(r, r ) r =r dr
that we stand on an electron with spin ↑ and we look at 2m 2m
i=1
the density of the other (N − 1) electrons. Pauli principle
forbids the presence of electrons with spin ↑ at the origin, (26)
but it says nothing about electrons with spin ↓, which can and EXC is the exchange and correlation energy
perfectly be located at the origin. Therefore,

1 ρ(r)ρ(r )
  EXC = [g(r, r ) − 1] dr dr (27)
gX (r, r ) −−−→ 1
2 for r −−−→ r (21) 2 |r − r |
6 Electronic structure of large molecules

3.1 Thomas–Fermi theory with µ the chemical potential. This equation can be inverted
to obtain the density as a unique function of the external
potential. This integral form in real space is inconvenient,
Thomas and Fermi (1927) gave a prescription for con-
but it can be easily inverted by Fourier transforming the
structing the total energy in terms only of the elec-
equation to obtain ρ(g).
tronic density.(14) They used the expression for the kinetic,
Exchange can be straightforwardly added to the expres-
exchange and correlation energies of the homogeneous
sion above by considering Slater’s expression for the
electron gas to construct the same quantities for the inhomo-
homogeneous electron gas: εX [ρ] = −CX ρ4/3 (r) dr, with
geneous system in the following way: Eα = εα [ρ(r)] dr,
CX = 3(3/π)1/3 /4. Expression (32) is modified by the
where εα [ρ(r)] is the energy density (corresponding to the
addition of the term −(4/3)CX ρ(r)1/3 . This level of approx-
piece α), calculated locally for the value of the density at
imation is called Thomas–Fermi–Dirac (TFD) theory.
that point in space. This was the first time that the local
Correlation can also be easily added by using any approx-
density approximation, or LDA, was used. For the homo-
imation to the homogeneous electron gas, for instance, the
geneous electron gas, the density is related to the Fermi
one proposed by Wigner: εC [ρ] = −0.056 ρ4/3 /[0.079 +
energy (εF ) by
ρ1/3 ].
 3/2 This is the best that can be done at the local level. Addi-
1 2m 3/2
ρ= εF (28) tional corrections to the kinetic, exchange and correlation
3π2 h̄2 energies due to non-locality have been postulated in the
form of gradient corrections, for example, as given by the
The kinetic energy of the homogeneous gas is T = 3ρεF /5, von Weiszäcker functional:(15)
so that the kinetic energy density is

1 |∇ρ|2
3 h̄2 TvW = dr (33)
t[ρ] = (3π2 )2/3 ρ5/3 (29) 8 ρ
5 2m
Also, terms that correct the linear response properties of

Then, the kinetic energy is written as TTF = Ck ρ(r)5/3 dr, the functional have been proposed,(16 – 18) and even the
with Ck = 3(3π2 )2/3 /10 = 2.871 atomic units. The inho- second-order response functions have been incorporated
mogeneous system is thought of as locally homogeneous. into this approach.(19) These have been developed in the
At variance with the usual approaches in modern DFT, here hopes that an explicit expression for the energy in terms of
the LDA is applied also to the kinetic energy. Neglecting the electronic density does really exist because an explicit-
exchange and correlation in expression (24), we arrive at density scheme requiring only the solution of the inverse
TF theory: problem is computationally much more efficient. But how
  do we know that the energy can be written as a functional
purely dependent on the density?
ETF [ρ] = Ck ρ(r)5/3 dr + v(r)ρ(r) dr

1 ρ(r)ρ(r )
+ dr dr (30) 3.2 Hohenberg–Kohn theorem
2 |r − r |

It can be seen that ETF depends only on the electronic den- In 1964, Hohenberg and Kohn(20) formulated and proved a
sity; it is a functional of the density. Assuming intuitively theorem that put on solid mathematical grounds the former
some variational principle, one can search for the density ideas, which were first proposed by Thomas and Fermi. The
ρ(r), which minimizes ETF [ρ], subjected to the constraint theorem is divided into two parts:
that the total
 integrated charge be equal to the number of
electrons: ρ(r) dr = N . This can be put in terms of func- THEOREM: The external potential is univocally
tional derivatives: determined by the electronic density, except for a trivial
additive constant.
  
δ
ETF [ρ] − µ ρ(r) dr = 0 (31) PROOF: We will suppose the opposite to hold, that
δρ(r)
the potential is not univocally determined by the density.
that is, Then one would be able to find two potentials v, v  such
that their ground state density ρ is the same. Let
and

|
 be the ground state and ground state energy
5 ρ(r ) E0 = 
|H
µ= C ρ(r)2/3 + v(r) + dr (32)
= T
+ U
+V
, and
 and E  = 
 |H
 |
  the
3 k |r − r | of H 0
 Density functional theory: Basics, new trends and applications 7

 = T
+ U
ground state and ground state energy of H
+V

 . The knowledge of F [ρ] implies that one has solved the full
Owing to the variational principle, we have many-body Schrödinger equation. It has to be remarked
that F [ρ] is a universal functional that does not depend

|
  = 
 |H
E0 < 
 |H
 |
  + 
 |H

−H

 |
  explicitly on the external potential. It depends only on the
 electronic density. In the Hohenberg–Kohn formulation,
= E0 + ρ(r)(v(r) − v  (r)) dr F [ρ] = 
|T
+ U
|
, where
is the ground state wave
function. These two theorems form the basis of DFT.
where we have also used the fact that different Hamiltoni- In Hohenberg–Kohn theorem the electronic density de-
ans have necessarily different ground states
=
 . This is termines the external potential, but it is also needed that
straightforward to show since the potential is a multiplica- the density corresponds to some ground state antisymmetric
tive operator. Now we can simply reverse the situation of wave function, and this is not always the case. However,

and
 (H and H  ) and readily obtain DFT can be reformulated in such a way that this is
not necessary, by appealing to the constrained search

 |
 = 
|H
E0 < 
|H
|
 + 
|H

 − H
|
 method.(21) By defining

= E0 − ρ(r)[v(r) − v  (r)] dr F [ρ] = min {
|T
+ U

|
}

→ρ

Adding these two inequalities, it turns out that E0 + E0 < for densities such that ρ(r) dr = N and
 non-negative
E0 + E0 , which is absurd. Therefore, there are no v(r) |∇ρ1/2 (r)|2 dr < ∞, which arise from an antisymmetric
= v  (r) that correspond to the same electronic density for wave function, the search is constrained to the subspace
the ground state. of all the antisymmetric
that give rise to the same
density ρ.
COROLLARY: Since ρ(r) univocally determines Using DFT, one can determine the electronic ground state
v(r), it also determines the ground state wave function
. density and energy exactly, provided that F [ρ] is known.
A common misleading statement is that DFT is a ground
THEOREM: Let ρ̃(r) be a non-negative density
state theory and that the question of excited states can-
normalized to N . Then E0 < Ev [ρ̃], for not be addressed within it. This is actually an incorrect
 statement because the density determines univocally the
Ev [ρ̃] = F [ρ̃] + ρ̃(r)v(r) dr (34) potential, and this in turn determines univocally the many-
body wave functions, ground and excited states, provided
with that the full many-body Schrödinger equation is solved.
For the ground state, such a scheme was devised by Kohn
F [ρ̃] = 
[ρ̃]|T
+ U

|
[ρ̃] (35)
and Sham and will be discussed in the next subsection. For
excited states there are a few extensions and generalizations
where
[ρ̃] is the ground state of a potential that has ρ̃ as
of Kohn–Sham (KS) theory, but only very recently these
its ground state density.
are beginning to be used with some degree of success. One
PROOF: We have such scheme, the ensemble DFT, proposed by Theophilou
 in 1979 and further developed by other authors,(22 – 25) is

|
[ρ̃] = F [ρ̃] +

[ρ̃]|H ρ̃(r)v(r) dr based on a generalization of Rayleigh–Ritz’s variational
principle applied to an ensemble of low-lying orthogo-

|

= Ev [ρ̃] ≥ Ev [ρ] = E0 = 
|H nal states. Another approach relies on an extension of
DFT to the time-dependent domain (time-dependent DFT,
The inequality follows from Rayleigh–Ritz’s variational or (TDDFT)).(26 – 28) Finally, a KS-like theory based on
principle for the wave function, but applied to the electronic the adiabatic connection between the eigenstates (not the
density. Therefore, the variational principle says ground state, but any eigenstate) of a non-interacting sys-
   tem with the same density as the fully interacting one was
recently proposed by Görling.(29,30)
δ Ev [ρ] − µ ρ(r) dr − N =0

and a generalized TF equation is obtained: 3.3 Kohn–Sham equations

δEv [ρ] δF [ρ] We have already briefly discussed the electron–electron

µ= = v(r) +
δρ δρ interaction potential U and we have seen that a reasonably
8 Electronic structure of large molecules

good description can be obtained by separating the electro- where we have chosen the occupation numbers to be 1 for
static (classical Coulomb energy), exchange and correlation i ≤ Ns (s = 1, 2) and 0 for i > Ns . This means that the
contributions. The biggest difficulty is to deal with corre- density is written as
lation. This is, in fact, an active field of research that has
produced significant improvements in the past decade. We 
2 
Ns

shall discuss this later on but for the moment let us men- ρ(r) = |ϕi,s (r)|2 (39)
tion that this issue is quite under control for most systems s=1 i=1

of interest. On the contrary, there is a problem with the

while the kinetic term is
expression of the kinetic energy 
|T
|
 in terms of the
electronic density. The only expression we have mentioned 
2 
Ns
∇2
up to now was the one proposed by Thomas and Fermi, TR [ρ] = ϕi,s | − |ϕ  (40)
2 i,s
which is local in the density. This is a severe shortcoming s=1 i=1
because this model does not hold bound states, and also the
electronic shell structure is absent. The main problem with The single-particle orbitals {ϕi,s (r)} are the Ns lowest
it is that the kinetic operator is inherently non-local, though eigenfunctions of ĥR = −(∇ 2 /2) + vR (r), that is,
short-ranged.  
In 1965, Kohn and Sham(31) proposed the idea of ∇2
− + vR (r) ϕi,s (r) = εi,s ϕi,s (r) (41)
replacing the kinetic energy of the interacting electrons with 2
that of an equivalent non-interacting system, because this
latter can be easily calculated. The density matrix ρ(r, r ) Using TR [ρ], the universal density functional can be rewrit-
that derives from the (interacting) ground state is the sum ten in the following form:
of the spin-up and spin-down density matrices, ρ(r, r ) =
 
 1 ρ(r)ρ(r )
s ρs (r, r )(s = 1, 2). The latter can be written as F [ρ] = TR [ρ] + dr dr + EXC [ρ] (42)
2 |r − r |
∞

ρs (r, r ) = ∗
ni,s ϕi,s (r)ϕi,s (r ) (36) where this equation defines the exchange and correlation
i=1 energy as a functional of the density.
The fact that TR [ρ] is the kinetic energy of the non-
where {ϕi,s (r)} are the natural spin orbitals and {ni,s } are interacting reference system implies that the correlation
the occupation numbers of these orbitals. The kinetic energy piece of the true kinetic energy has been ignored and has
can be written exactly as to be taken into account somewhere else. In practice, this
∞
is done by redefining the correlation energy functional in

2 
∇2 such a way as to include kinetic correlations.
T = ni,s ϕi,s | − |ϕ  (37)
s=1 i=1
2 i,s Upon substitution of this expression
 for F in the total
energy functional Ev [ρ] = F [ρ] + ρ(r)v(r) dr, the latter
In the following we shall assume that the equivalent non- is usually renamed the KS functional:
interacting system, that is, a system of non-interacting 
fermions whose ground state density coincides with that EKS [ρ] = TR [ρ] + ρ(r)v(r) dr
of the interacting system, does exist. We shall call this the

non-interacting reference system of density ρ(r), which is 1 ρ(r)ρ(r )
described by the Hamiltonian + dr dr + EXC [ρ] (43)
2 |r − r |
 

N
∇i2 In this way we have expressed the density functional in

=
H − + vR (ri ) (38)
R
2 terms of the N = N↑ + N↓ orbitals (KS orbitals), which
i=1
minimize the kinetic energy under the fixed density con-
where the potential vR (r) is such that the ground state den- straint. In principle, these orbitals are a mathematical object

equals ρ(r) and the ground state energy equals
sity of H constructed in order to render the problem more tractable
R
the energy of the interacting system. This Hamiltonian has and do not have a sense by themselves, but only in terms
no electron–electron interactions and, thus, its eigenstates of the density. In practice, however, it is customary to
can be expressed in the form of Slater determinants consider them as single-particle physical eigenstates. It is
usual to hear that the KS orbitals are meaningless and can-
1 not be identified as single-particle eigenstates, especially in

s (r) = √ SD[ϕ1,s (r1 )ϕ2,s (r2 ) · · · ϕNs ,s (rNs )]
N! the context of electronic excitations. A rigorous treatment,
 Density functional theory: Basics, new trends and applications 9

however, shows that KS eigenvalues differences are a well- counting terms have to be subtracted:
defined approximation to excitation energies.(29,30)

Ns  
ρ(r)ρ(r )
2
The KS orbitals satisfy equation (41) while the problem 1
EKS [ρ] = εi,s − dr dr
is to determine the effective potential vR , or veff as it is also
i=1 s=1
2 |r − r |
known. This can be done by minimizing the KS functional   
over all densities that integrate to N particles. For the + EXC [ρ] − ρ(r)µXC [ρ](r) dr (50)
minimizing (i.e., correct) density ρ, we have

δTR [ρ] ρ(r ) δE [ρ]
+ v(r) + dr + XC = µR (44) 3.3.1 Interpretation
δρ(r) |r − r | δρ(r)
By introducing the non-interacting reference system, we
The functional derivative δTR [ρ]/δρ(r) can be quickly were able to take into account the most important part of


found by considering the non-interacting Hamiltonian H the kinetic energy. The missing part (correlations) is due
R
(equation 38). Its ground state energy is E0 . We can to the fact that the full many-body wave function is not a
construct the functional single Slater determinant, otherwise HF theory would be
 exact. If we think of a true non-interacting system, then
EvR [ρ̃] = TR [ρ̃] + ρ̃(r) vR (r) dr (45) the KS scheme is exact, while TF theory is quite a poor
approximation that becomes reasonably good only when the
electronic density is very smooth, as in alkali metals.
Then, clearly EvR [ρ̃] ≥ E0 , and only for the correct den-
The price we have to pay for having a good description
sity ρ we will have EvR [ρ] = E0 . Hence, the functional
of the kinetic energy is that, instead of solving a single
derivative of EvR [ρ̃] must vanish for the correct density
equation for the density in terms of the potential, we have
leading to
to solve a system of N Euler equations. The main difference
δTR [ρ] between the KS and Hartree equations is that the effective
+ vR (r) = µR (46)
δρ(r) potential now includes exchange and correlation. Therefore,
the computational cost is of the same order as Hartree,
where µR is the chemical potential for the non-interacting but much less than HF, which includes the exact non-local
system. exchange. Now let us make some observations:
To summarize, the KS orbitals satisfy the well-known
self-consistent KS equations 1. The true wave function is not the Slater determinant of
  KS orbitals, although it is determined by the density,
∇2 and thus by the KS orbitals used to construct the
− + veff (r) ϕi,s (r) = εi,s ϕi,s (r) (47)
2 density.
2. The correlation functional has to be modified to account
where, by comparison of expressions 44 and 46, the for the missing part in the kinetic energy TR [ρ], which
effective potential vR or veff is given by corresponds to a non-interacting system. The exchange
 functional remains unchanged.
ρ(r )
veff (r) = v(r) + dr + µXC [ρ](r) (48) 3. Nothing ensures that the non-interacting reference sys-
|r − r | tem will always exist. In fact, there are examples like
the carbon dimer C2 , which do not satisfy this require-
and the electronic density is constructed with KS orbitals ment. In that case, a linear combination of Slater deter-
minants that include single-particle eigenstates ϕi,s (r)

Ns 
2
ρ(r) = |ϕi,s (r)|2 (49) with i > Ns can be considered. This is equivalent to
i=1 s=1 extending the domain of definition of the occupation
numbers ni,s from the integer values 0 and 1 to a
The exchange-correlation potential µXC [ρ](r) defined continuum between 0 and 1. In such a way we are
above is simply the functional derivative of the exchange- including excited single-particle states in the density.
correlation energy δEXC [ρ]/δρ. Notice the similitude At this point, some authors proposed to carry out the
between the KS and Hartree equations (equation 9). minimization of the energy functional not only with
The solution of the KS equations has to be obtained by respect to KS orbitals but also with respect to the occu-
an iterative procedure, in the same way as Hartree and HF pation numbers.(32) Although there is nothing wrong,
equations. As in these methods, the total energy cannot be in principle, with minimizing the functional constructed
written simply as the sum of the eigenvalues εi,s , but double with fractional occupation numbers, the minimization
10 Electronic structure of large molecules

with respect to them is not justified.(33) The introduc- the interacting ground state of the external potential that
tion of excited single-particle states does not mean that has ρ as the ground state density), were known, then we
the system is in a true excited state. This is only an could use the original definition of the exchange-correlation
artefact of the representation. The true wave function energy that does not contain kinetic contributions:
is the correlated ground state. 
Janak’s theorem is valid.(34) The ionization energy is 1 ρ(r)ρ(r )
4. 0
EXC [ρ] = [g(r, r ) − 1] dr dr (51)
given by I = −µ = −εmax (if the effective potential 2 |r − r |
vanishes at long distances), while the eigenvalues are
Since we are using the non-interacting expression for the
defined as the derivatives of the total energy with
kinetic energy TR [ρ], we have to redefine it in the following
respect to the occupation numbers: εi,s = ∂E/∂ni,s .
way:
5. In DFT there is no Koopman’s theorem that would
allows us to calculate electron removal energies as EXC [ρ] = EXC
0
[ρ] + T [ρ] − TR [ρ]
the difference between the ground state energy of an
(N + 1)-electron system and that of an N -electron sys- It can be shown that the kinetic contribution to the
tem. Excitations in DFT are still an open issue because, correlation energy (the kinetic contribution to exchange is
even if the density determines the whole spectrum just Pauli’s principle, which is already contained in TR [ρ]
via the many-body wave function, standard approxi- and in the density when adding up the contributions of
mations focus only on the ground state. Nevertheless, the N lowest eigenstates) can be taken into account by
extensions have been devised that made it possible to averaging the pair correlation function g(r, r ) over the
address the question of excited states within a DFT-like strength of the electron–electron interaction, that is,
framework, in addition to the traditional many-body 
scenarios.(22 – 30) 1 ρ(r)ρ(r )
EXC [ρ] = [g̃(r, r ) − 1] dr dr (52)
2 |r − r |
3.3.2 Summary
where
We have described a theory that is able to solve the  1

complicated many-body electronic ground state problem by g̃(r, r ) = gλ (r, r ) dλ (53)
mapping exactly the many-body Schrödinger equation into 0
a set of N coupled single-particle equations. Therefore,
and gλ (r, r ) is the pair correlation function corresponding
given an external potential, we are in a position to find
= T
+ V
+ λU
.(35) If we separate
to the Hamiltonian H ee
the electronic density, the energy and any desired ground
the exchange and correlation contributions, we have
state property (e.g., stress, phonons, etc.). The density
of the non-interacting reference system is equal to that 
of the true interacting system. Up to now the theory is |ρσ (r, r )|2
σ
exact. We have not introduced any approximation into the g̃(r, r ) = 1 − + g̃C (r, r ) (54)
ρ(r)ρ(r )
electronic problem. All the ignorance about the many-
fermion problem has been displaced to the EXC [ρ] term, with ρσ (r, r ) the spin-up and spin-down components of the
while the remaining terms in the energy are well known. one-body density matrix, which in general is a non-diagonal
In the next section we are going to discuss the exchange operator. For the homogeneous electron gas, the expression
and correlation functionals. But now, we would like to for the density matrix is well known, so that the exchange
know how far is TR [ρ] from T [ρ]. Both are the expectation contribution to g̃(r, r ) assumes an analytic closed form
values of the kinetic operator, but in different states. The
non-interacting one corresponds to the expectation value  2
  9 j1 (kF |r − r |)
in the ground state of the kinetic operator, while the gX (r, r ) = gX (|r − r |) = 1 − (55)
2 kF |r − r |
interacting one corresponds to the ground state of the full
Hamiltonian. This means that TR [ρ] ≤ T [ρ], implying that where j1 (x) = [sin(x) − x cos(x)]/x 2 is the first-order
ẼC [ρ] contains a positive contribution arising from the spherical Bessel function.
kinetic correlations. In Figure 1, we reproduce from Perdew and Wang(36)
the shape of the non-oscillatory part of the pair-distribution
4 EXCHANGE AND CORRELATION function, g(r), and its coupling constant average, g̃(r), for
the unpolarized uniform electron gas of density parameter
If the exact expression for the kinetic energy including rs = 2. The same function within the Hartree approximation
correlation effects, T [ρ] = 
[ρ]|T
|
[ρ] (with
[ρ] being is the constant function 1, because the approximation
 Density functional theory: Basics, new trends and applications 11

Scaled distance R /rs a 0 used for ρ̃XC (r, r ). If we separate the exchange and
1.1 correlation contributions, it is easy to see that the displaced
1 electron comes exclusively from the exchange part, and it is

Pair distribution function

0.9 a consequence of the form in which the electron–electron
0.8 interaction has been separated. In the Hartree term we
0.7 have included the interaction of the electron with itself.
0.6 This unphysical contribution is exactly cancelled by the
0.5 exchange interaction of the full charge density with the
0.4 displaced density. However, exchange is more than that.
0.3 It is a non-local operator whose local component is less
0.2 the self-interaction. On the other hand, the correlation hole
0.1 integrates to zero, ρ̃C (r, r ) dr = 0, so that the correlation
0 energy corresponds to the interaction of the charge density
0 0.5 1 1.5 2 with a neutral charge distribution.
Figure 1. Pair correlation function (solid line) and its coupling
A general discussion on DFT and applications can be
constant average (dashed line) for the uniform electron gas found in Reference 37.
[Reproduced by permission of APS Journals from J.P. Perdew
and Y. Wang (1992) Phys. Rev. B, 46, 12 947.](36)
4.1 The local density approximation
completely neglects both, exchange and correlation, so that
one electron is insensitive to the location of the other The LDA has been for a long time the most widely
electron. Within the HF approximation, the exchange is used approximation to the exchange-correlation energy. It
treated exactly, but the correlation is ignored. Therefore, the has been proposed in the seminal paper by Kohn and
HF pair distribution only reveals the fact that the electrons Sham, but the philosophy was already present in TF
with like spins do not like to be at the same place, and hence theory. The main idea is to consider general inhomogeneous
the HF pair correlation function is given by formula (55), electronic systems as locally homogeneous, and then to
tending to the limit 1/2 for r → 0. use the exchange-correlation hole corresponding to the
We are now going to define the exchange-correlation hole homogeneous electron gas for which there are very good
ρ̃XC (r, r ) in the following form: approximations and also exact numerical (quantum Monte
 Carlo) results. This means that
1 ρ(r)ρ̃XC (r, r )
EXC [ρ] = dr dr (56)
2 |r − r | ρ̃LDA  
XC (r, r ) = ρ(r)(g̃ [|r − r |, ρ(r)] − 1)
h
(57)
or ρ̃XC (r, r ) = ρ(r )[g̃(r, r ) − 1]. with g̃ h [|r − r |, ρ(r)] the pair correlation function of the
Then, ẼXC [ρ] can be written as the interaction between homogeneous gas, which depends only on the distance
the electronic charge distribution and the charge distribution between r and r , evaluated at the density ρh , which locally
that has been displaced by exchange and correlation effects, equals ρ(r). Within this approximation, the exchange-
that is, by the fact that the presence of an electron at correlation energy density is defined as
r reduces the probability for a second electron to be
at r , in the vicinity of r. Actually, ρ̃XC (r, r ) is the  
1 ρ̃LDA
XC (r, r )
exchange-correlation hole averaged over the strength of the εLDA
XC [ρ] = dr (58)
2 |r − r |
interaction, which takes into account kinetic correlations.
The properties of g̃(r, r ) and ρ̃XC (r, r ) are very interesting and the exchange-correlation energy becomes
and instructive: 
  EXC [ρ] = ρ(r)εLDA
LDA
XC [ρ] dr (59)
1. g̃(r, r ) = g̃(r , r) (symmetry)

2. g̃(r, r )ρ(r ) dr = g̃(r, r )ρ(r) dr = N − 1


(normalization) In general, the exchange-correlation energy density is not a

 
3. ρ̃XC (r, r ) dr = ρ̃XC (r, r ) dr = −1. local functional of ρ. From its very definition it is clear that
it has to be a non-local object, because it reflects the fact
This means that the exchange-correlation hole contains that the probability of finding an electron at r depends on
exactly one displaced electron. This sum rule is very the presence of other electrons in the surroundings, through
important and it has to be verified by any approximation the exchange-correlation hole.
12 Electronic structure of large molecules

Looking at expression (57), it may seem that there is an asking for the lowest N = N↑ + N↓ to be occupied. This
inconsistency in the definition. The exact expression would defines a Fermi energy εF such that the occupied eigenstates
indicate to take ρ(r ) instead of ρ(r). However, this would have εi,s < εF .
render EXCLDA
[ρ] a non-local object that would depend on In the case of non-magnetic systems, ρ↑ (r) = ρ↓ (r), and
the densities at r and r , and we want to parametrize it with everything reduces to the simple case of double occupancy
the homogeneous gas, which is characterized by only one of the single-particle orbitals.
density, and not two. This is the essence of the LDA, and The equivalent of the LDA in spin-polarized systems
it is equivalent to postulate is the local spin density approximation (LSDA), which
  basically consists of replacing the XC energy density with
ρ(r) a spin-polarized expression:
g̃(r, r ) ≈ g̃ h [|r − r |, ρ(r)] (60)
ρ(r )
LSDA
EXC [ρ↑ (r), ρ↓ (r)]
Therefore, there are in fact two approximations embodied 
in the LDA: = [ρ↑ (r) + ρ↓ (r)]εhXC [ρ↑ (r), ρ↓ (r)] dr (62)
1. The LDA exchange-correlation hole is centred at r and
interacts with the electronic density at r. The true XC obtained, for instance, by interpolating between the fully
hole is actually centred at r instead of r. polarized and fully unpolarized XC energy densities using
2. The pair correlation function (g) is approximated by an appropriate expression that depends on ζ(r). The stan-
that of the homogeneous electron gas of density ρ(r) dard practice is to use von Barth and Hedin’s interpolation
corrected by the density ratio ρ(r)/ρ(r ) to compensate formula:(38)
the fact that the LDA XC hole is centred at r instead
of r . εhXC [ρ↑ , ρ↓ ] = f (ζ)εP [ρ] + [1 − f (ζ)]εU [ρ]
(1 + ζ)4/3 + (1 − ζ)4/3 − 2
f (ζ) = (63)
4.2 The local spin density approximation 24/3 − 2
or a more realistic formula based on the random
In magnetic systems or, in general, in systems where phase approximation (RPA), given by Vosko, Wilk and
open electronic shells are involved, better approximations Nussair.(39)
to the exchange-correlation functional can be obtained by A thorough discussion of the LDA and the LSDA can be
introducing the two spin densities, ρ↑ (r) and ρ↓ (r), such found in Reference 40. In the following we reproduce the
that ρ(r) = ρ↑ (r) + ρ↓ (r), and ζ(r) = (ρ↑ (r) − ρ↓ (r))/ρ(r) main aspects related to these approximations.
is the magnetization density. The non-interacting kinetic
energy (equation 40) splits trivially into spin-up and spin- 4.2.1 Why does the LDA work so well in many
down contributions, and the external and Hartree potential cases?
depend on the full density ρ(r), but the approximate XC
functional – even if the exact functional should depend only 1. It satisfies the sum rule that the XC hole contains
on ρ(r) – will depend on both spin densities independently, exactly one displaced electron:
EXC = EXC [ρ↑ (r), ρ↓ (r)]. KS equations then read exactly  
 
as in equation (47), but the effective potential veff (r) now ρ̃LDA
XC (r, r ) dr = ρ(r)g̃ h [|r − r |, ρ(r)] dr = −1
acquires a spin index: (64)
 because for each r, g̃ h [|r − r |, ρ(r)] is the pair cor-
↑ ρ(r )
veff (r) = v(r) + dr relation function of an existing system, that is, the
|r − r | homogeneous gas at density ρ(r). Therefore, the mid-
δEXC [ρ↑ (r), ρ↓ (r)] dle expression is just the integral of the XC hole of the
+ (61)
δρ↑ (r) homogeneous gas. For the latter, both approximations
 and numerical results carefully take into account that
ρ(r ) δEXC [ρ↑ (r), ρ↓ (r)]
↓
veff (r) = v(r) + dr 
+ the integral has to be −1.
|r − r | δρ↓ (r) 2. Even if the exact ρ̃XC has no spherical symmetry, in
the expression for the XC energy what really matters
The density given by expression (49) contains a double is the spherical average of the hole:
summation, over the spin states and over the number of
  
electrons in each spin state, Ns . The latter have to be 1 1
EXC [ρ] = − ρ(r) dr
determined according to the single-particle eigenvalues, by 2 R(r)
 Density functional theory: Basics, new trends and applications 13

with For rs ≤ 1, the expression arises from the RPA – calculated

  by Gell–Mann and Brückner(43) – which is valid in the
1 ρ̃XC (r, r )  ∞
= dr = 4π s ρ̃SA
XC (r, s) ds
limit of very dense electronic systems. For low densities,
R(r) |r − r | 0 Perdew and Zunger have fitted a Padé approximant to the
and Monte Carlo results. Interestingly, the second derivative
 of the above εC [ρ] is discontinuous at rs = 1. Another
1
XC (r, s) =
ρ̃SA ρ̃XC (r, r ) d popular parametrization is that proposed by Vosko, Wilk
4π 
and Nussair.(39)
The spherical average ρ̃SA
XC (r, s) is reproduced to a good
extent by the LDA, whose ρ̃XC is already spherical. 4.2.4 When does the LDA fail?
4.2.2 Trends within the LDA The LDA is very successful an approximation for many
systems of interest, especially those where the electronic
There are a number of features of the LDA that are
density is quite uniform such as bulk metals, but also for
rather general and well established by now. These are the
less uniform systems as semiconductors and ionic crystals.
following:
There are, however, a number of known features that the
1. It favours more homogeneous systems. LDA fails to reproduce:
2. It overbinds molecules and solids.
3. Chemical trends are usually correct. 1. In atomic systems, where the density has large varia-
4. For ‘good’ systems (covalent, ionic and metallic tions, and also the self-interaction is important.
bonds), geometries are good, bond lengths, bond angles 2. In weak molecular bonds, for example, hydrogen
and phonon frequencies are within a few percent, while bonds, because in the bonding region the density is very
dielectric properties are overestimated by about 10%. small and the binding is dominated by inhomogeneities.
5. For ‘bad’ systems (weakly bound), bond lengths are 3. In van der Waals (closed-shell) systems, because there
too short (overbinding). the binding is due to dynamical charge–charge corre-
6. In finite systems, the XC potential does not decay lations between two separated fragments, and this is an
as −e2 /r in the vacuum region, thus affecting the inherently non-local interaction.
dissociation limit and ionization energies. This is a 4. In metallic surfaces, because the XC potential decays
consequence of the fact that both the LDA and the exponentially, while it should follow a power law
LSDA fail at cancelling the self-interaction included in (image potential).
the Hartree term of the energy. This is one of the most 5. In negatively charged ions, because the LDA fails
severe limitations of these approximations. to cancel exactly the electronic self-interaction,
owing to the approximative character of the ex-
4.2.3 What parametrization of EXC is normally used change. Self-interaction-corrected functionals have
within the LDA? been proposed,(42) although they are not satisfactory
For the exchange energy density, the form deduced by Dirac from the theoretical point of view because the potential
is adopted:(41) depends on the electronic state, while it should be the
same for all states. The solution to this problem is the
 1/3  
3 3 3 9 1/3 1 exact treatment of exchange (see Section 5).
εX [ρ] = − ρ1/3 = −
4 π 4 4π2 rs 6. The energy band gap in semiconductors turns out to
be very small. The reason is that when one electron
0.458
=− au (65) is removed from the ground state, the exchange hole
rs becomes screened, and this is absent in the LDA. On
the other hand, HF also has the same limitation, but
where ρ−1 = 4πrs3 /3 and rs is the radius of the sphere that,
the band gap turns out to be too large.
on average, contains one electron.
For the correlation, a widely used approximation is
Perdew and Zunger’s parametrization(42) of Ceperley and 4.2.5 How can the LDA be improved?
Alder quantum Monte Carlo results, which are essentially
exact,(11,12) Once the extent of the approximations involved in the
LDA has been understood, one can start constructing better

A ln rs + B + Crs ln rs + Drs , rs ≤ 1 approximations. The amount of work done in that direction
εC [ρ] = √ (66)
γ/(1 + β1 rs + β2 rs ), rs > 1 is really overwhelming, and there are new developments in
14 Electronic structure of large molecules

many different directions because there is not a unique and known as generalized gradient approximation (GGA), is
obvious way of improving the LDA. easier to implement in practice, and computationally more
One of the key observations is that the true pair correla- convenient than full many-body approaches, and has been
tion function, g(r, r ), actually depends on the electronic quite successful in improving over some features of the
density at two different points, r and r . The homoge- LDA.
neous g(r, r ) is quite well known (see equation 55 for
the exchange part and Reference 36 for correlation), but
it corresponds to a density that is the same everywhere. 4.3 Generalized gradient approximations
Therefore, the question is which of the two densities are
to be used in an inhomogeneous environment. Early efforts The exchange-correlation energy has a gradient expansion
went into the direction of calculating the pair correlation of the type
function at an average density ρ̄(r), which in general is 
different from ρ(r), and incorporates information about the EXC [ρ] = AXC [ρ] ρ(r)4/3 dr
density at neighbouring points. Clearly, there is no unique 
recipe for the averaging procedure, but there is at least a + CXC [ρ] |∇ρ(r)|2 /ρ(r)4/3 dr + · · · (68)
crucial condition that this averaging has to verify, namely,
the normalization condition:(44 – 48)
which is asymptotically valid for densities that vary slowly
 
in space. The LDA retains only the leading term of
ρ̃XC (r, r ) dr = ρ(r ) g̃ h [|r − r |, ρ̄(r)] dr = −1
WDA  
equation (68). It is well known that a straightforward
(67) evaluation of this expansion is ill-behaved, in the sense
Approaches of this type receive the name of weighted that it is not monotonically convergent, and it exhibits
density approximations (WDA). There is still a lot of singularities that cancel out only when an infinite number
freedom in choosing the averaging procedure provided of terms is re-summed, as in the RPA. In fact, the first-
that normalization is verified and, indeed, several different order correction worsens the results and the second-order
approximations have been proposed.(44 – 51) One problem correction is plagued with divergences.(55,56) The largest
with this approach is that the r → r symmetry of g(r, r ) error of this approximation actually arises from the gradient
is now broken. Efforts in the direction of the WDA are contribution to the correlation term. Provided that the
intended to improve over the incorrect location of the problem of the correlation term can be cured in some way,
centre of the XC hole in the LDA. An exploration in as the real space cut-off method proposed by Langreth and
the context of realistic electronic structure calculations was Mehl,(57,58) the biggest problem remains with the exchange
carried out by Singh but the results reported were not energy.
particularly encouraging.(52) Nevertheless, this is a direction Many papers have been devoted to the improvement of
worth exploring in more depth. the exchange term within DFT. The early work of Gross
Another possibility is to employ either standard or and Dreizler(59) provided a derivation of the second-order
advanced many-body tools, for example, one could try to expansion of the exchange density matrix, which was later
solve Dyson’s equation for the electronic Green’s function, re-analysed and extended by Perdew.(60) This expansion
starting from the LDA solution for the bare Green’s contains not only the gradient but also the Laplacian of the
function.(53) In the context of strongly correlated systems, density. The same type of expansion was obtained, using
for example those exhibiting narrow d or f bands, where Wigner distribution – phase space – techniques, by Ghosh
the limitation of the LDA is at describing strong on- and Parr.(61)
site correlations of the Hubbard type, these features have One of the main lessons learnt from these works is that
been introduced a posteriori within the so-called LDA + U the gradient expansion has to be carried out very carefully
approach.(54) This theory considers the mean-field solution in order to retain all the relevant contributions to the desired
of the Hubbard model on top of the LDA solution, where order. The other important lesson is that these expansions
the Hubbard on-site interaction U are computed for the d easily violate one or more of the exact conditions required
or f orbitals by differentiating the LDA eigenvalues with for the exchange and the correlation holes. For instance,
respect to the occupation numbers. the normalization condition, the negativity of the exchange
Undoubtedly, and probably because of its computational density and the self-interaction cancellation (the diagonal of
efficiency and its similarity to the LDA, the most popular the exchange density matrix has to be minus a half of the
approach has been to introduce semi-locally the inhomo- density). Perdew has shown that imposing these conditions
geneities of the density, by expanding EXC [ρ] as a series to functionals that originally do not verify them results
in terms of the density and its gradients. This approach, in a remarkable improvement of the quality of exchange
 Density functional theory: Basics, new trends and applications 15

energies.(60) On the basis of this type of reasoning, a number where

of modified gradient expansions have been proposed along
Cc (∞) |∇ρ(r)|
the years, mainly between 1986 and 1996. These have  = 1.745f˜
received the name of GGA. Cc (ρ) ρ(r)7/6
GGAs are typically based either on theoretical develop- C2 + C3 rs + C4 rs2
ments that reproduce a number of exact results in some Cc (ρ) = C1 +
1 + C5 rs + C6 rs2 + C7 rs3
known limits, for example, 0 and ∞ density, or the corre-
lation potential in the He atom, or are generated by fitting a
being f˜ = 0.11, C1 = 0.001667, C2 = 0.002568, C3 =
number of parameters to a molecular database (training set).
0.023266, C4 = 7.389 × 10−6 , C5 = 8.723, C6 = 0.472,
Normally, these improve over some of the drawbacks of the
C7 = 7.389 × 10−2 .
LDA, although this is not always the case. These aspects
4. Perdew–Wang ’91 (PW91) exchange functional.(65)
will be discussed below, after presenting some popular
functionals.
εX
The basic idea of GGAs is to express the exchange-  
1 + a1 s sinh−1 (a2 s) + (a3 + a4 e−100s )s 2
2
correlation energy in the following form:
= εLDA
1 + a1 s sinh−1 (a2 s) + a5 s 4
X
 
EXC [ρ] = ρ(r) εXC [ρ(r)] dr + FXC [ρ(r), ∇ρ(r)] dr
(69) where a1 = 0.19645, a2 = 7.7956, a3 = 0.2743, a4 =
where the function FXC is asked to satisfy a number of −0.1508 and a5 = 0.004.
formal conditions for the exchange-correlation hole, such as 5. Perdew–Wang ’91 (PW91) correlation functional.(65)
sum rules, long-range decay and so on. This cannot be done
by considering directly the bare gradient expansion (68). εC = εLDA
C + ρH [ρ, s, t]
What is needed from the functional is a form that mimics a
re-summation to infinite order, and this is the main idea of with
the GGA, for which there is not a unique recipe. Naturally,  
not all the formal properties can be enforced at the same β 2α t 2 + At 4
H [ρ, s, t] = ln 1 +
time, and this differentiates one functional from another. 2α β 1 + At 2 + A2 t 4
A thorough comparison of different GGAs can be found in  
+ Cc0 Cc (ρ) − Cc1 t 2 e−100s
2

Reference 62. In the following we quote a number of them:

1. Langreth–Mehl (LM) exchange-correlation func- and

tional.(57) 2α  −2αεC [ρ]/β2 −1
A= e −1
  β
|∇ρ(r)|2 7
εX = εLDA −a + 18f 2
X
ρ(r)4/3 9 where α = 0.09, β = 0.0667263212, Cc0 = 15.7559,
|∇ρ(r)|  −F
2  Cc1 = 0.003521, t = |∇ρ(r)|/(2ks ρ) for ks = (4kF /
εC = εRPA +a 2 e + 18 f 2 π)1/2 , and ρεC [ρ] = εLDA
C [ρ].
C
ρ(r)4/3
6. Becke ’88 (B88) exchange functional.(66)
where F = b|∇ρ(r)|/ρ(r)7/6 , b = (9π)1/6 f , a = π/  
β x2
(16(3π2 )4/3 ) and f = 0.15. εX = εLDA 1 − 1/3
2. Perdew–Wang ’86 (PW86) exchange functional.(63)
X
2 Ax 1 + 6βx sinh−1 (x)

 m for x = 2(6π2 )1/3 s = 21/3 |∇ρ(r)|/ρ(r)4/3 , Ax = (3/4)

s2
εX = εLDA
X 1 + 0.0864 + bs 4
+ cs 6
(3/π)1/3 , and β = 0.0042.
m
7. Closed-shell, Lee–Yang–Parr (LYP) correlation func-
with m = 1/15, b = 14, c = 0.2 and s = |∇ρ(r)|/ tional.(67)
(2kF ρ) for kF = (3π2 ρ)1/3 .  
1 −2/3
3. Perdew–Wang ’86 (PW86) correlation functional.(64) εC = −a ρ + bρ CF ρ5/3 − 2tW
1 + dρ−1/3
  
|∇ρ(r)|2 1 1 −1/3
εC = εLDA + e− Cc (ρ) + tW + ∇ 2 ρ e−cρ
C
ρ(r)4/3 9 2
16 Electronic structure of large molecules

where spin-scaling factor. The quantity β is the same as for the

 
1 |∇ρ|
2 exchange term β = 0.066725, and γ = (1 − ln 2)/π2 =
tW = − ∇ 2ρ 0.031091. The function A has the following form:
8 ρ
β  −εLDA [ρ]/(γφ3 e2 /a0 ) −1
and CF = 3/10(3π ) , a = 0.04918, b = 0.132, c =
2 2/3 A= e C −1
γ
0.2533 and d = 0.349. This correlation functional is
not based on the LDA as the others, but it has So defined, the correlation correction term H satisfies
been derived as an extension of the Colle–Salvetti the following properties: (i) it tends to the correct
expression for the electronic correlation in Helium, to second-order gradient expansion in the slowly varying
other closed-shell systems. (high-density) limit (t → 0), (ii) it approaches minus
8. Perdew–Burke–Ernzerhof (PBE) exchange-correlation the uniform electron gas correlation −εLDA
C for rapidly
functional.(68,69) First, the enhancement factor FXC varying densities (t → ∞), thus making the correlation
over the local exchange is defined: energy to vanish (this results from the correlation hole
 sum rule), (iii) it cancels the logarithmic singularity
EXC [ρ] = ρ(r)εLDA of εLDA in the high-density limit, thus forcing the
X [ρ(r)]FXC (ρ, ζ, s) dr C
correlation energy to scale to a constant under uniform
scaling of the density.
where ρ is the local density, ζ is the relative
This GGA retains the correct features of LDA
spin polarization and s = |∇ρ(r)|/(2kF ρ) is the
(LSDA) and combines them with the inhomogeneity
dimensionless density gradient, as in PW86:
features that are supposed to be the most energetically
κ important. It sacrifices a few correct but less
FX (s) = 1 + κ− important features, like the correct second-order
1 + µs 2 /κ
gradient coefficients in the slowly varying limit, and
where µ = β(π2 /3) = 0.21951 and β = 0.066725 is the correct non-uniform scaling of the exchange energy
related to the second-order gradient expansion.(65) in the rapidly varying density region.
This form: (i) satisfies the uniform scaling condition,
In the beginning of the age of GGAs, the most popu-
(ii) recovers the correct uniform electron gas limit
lar recipe was to use B88 exchange complemented with
because Fx (0) = 1, (iii) obeys the spin-scaling rela-
Perdew ’86 correlation corrections (BP). For exchange, B88
tionship, (iv) recovers the LSDA linear response limit
remained preferred, while LYP correlation proved to be
for s → 0 (FX (s) → 1 + µs 2 ) and (v) satisfies the
more accurate than Perdew ’86, particularly for hydrogen-
local Lieb-Oxford bound,(70) εX (r) ≥ −1.679ρ(r)4/3 ,
bonded systems (BLYP). The most recent GGA in the mar-
that is, FX (s) ≤ 1.804, for all r, provided that κ ≤
ket is the PBE due to Perdew, Burke and Ernzerhof.(68,69)
0.804. PBE choose the largest allowed value κ =
This is very satisfactory from the theoretical point of view,
0.804. Other authors have proposed the same form, but
because it verifies many of the exact conditions for the XC
with values of κ and µ fitted empirically to a database
hole and it does not contain any fitting parameters. In addi-
of atomization energies.(71 – 73) The proposed values of
tion, its quality is equivalent or even better than BLYP.(74)
κ violate Lieb–Oxford inequality.
The different recipes for GGAs verify only some of
The correlation energy is written in a form similar to
the mathematical properties known for the exact-exchange-
PW91,(65) that is,
correlation hole. A compilation and comparison of different

  approximations can be found in the work of Levy and
ECGGA = ρ(r) εLDA C (ρ, ζ) + H [ρ, ζ, t] dr Perdew.(75)

with
4.3.1 Trends of the GGAs
H [ρ, ζ, t] The general trends of GGAs concerning improvements over
 2 !  " the LDA are the following:
e βγ 2 1 + At 2
= γφ ln 1 +
3
a0 t 1 + At 2 + A2 t 4 1. They improve binding energies and also atomic ener-
gies.
Here, t = |∇ρ(r)|/(2φks ρ) is a dimensionless density 2. They improve bond lengths and angles.
gradient, ks = (4kF /πa0 )1/2 is the TF screening wave 3. They improve energetics, geometries and dynamical
number and φ(ζ) = [(1 + ζ)2/3 + (1 − ζ)2/3 ]/2 is a properties of water, ice and water clusters. BLYP and
 Density functional theory: Basics, new trends and applications 17

PBE show the best agreement with experiment. In years is to treat the exchange term exactly. Some authors
general, they improve the description of hydrogen- call these third-generation XC functionals, in relation to
bonded systems, although this is not very clear for the the early TF-like, and successive LDA-like, functionals.(79)
case of the F · · ·H bond. Exact exchange methods are described in the next section,
4. Semiconductors are marginally better described within followed by methods that combine exact exchange (EXX)
the LDA than in GGA, except for the binding energies. with density functional perturbation theory for correlation.
5. For 4d–5d transition metals, the improvement of the The properties of this approach are very elegant, and the
GGA over LDA is not clear and will depend on how error cancellation property present in GGA, meta-GGA
well the LDA does in the particular case. and hybrid methods is very much reduced. The computa-
6. Lattice constants of noble metals (Ag, Au, Pt) are tional cost of these two approaches is, at present, very high
overestimated. The LDA values are very close to compared to standard GGA or meta-GGA-like functionals.
experiment, and thus any modification can only worsen Nevertheless, they are likely to become widespread in the
them. future.
7. There is some improvement for the gap problem (and Finally, another possibility is to abandon the use of the
consequently for the dielectric constant), but it is homogeneous electron gas as a reference system (used at
not substantial because this feature is related to the the LDA level) for some other reference state. A functional
description of the screening of the exchange hole when for ‘edge’ states, inspired in the behaviour of the density
one electron is removed, and this feature is usually not at the surface of a system, has been proposed by Kohn and
fully taken into account by GGA. Mattson,(80) and further developed by Vitos et al.(81,82)
8. They do not satisfy the known asymptotic behaviour,
for example, for isolated atoms:
4.4 Meta-GGA
LDA,GGA
(a) vXC (r) ∼ −e /r for r → ∞, while
2
vXC (r)
The second-order gradient expansion of the exchange
vanish exponentially.
energy introduces a term proportional to the squared
(b) vXC (r) → const. for r → 0, while vXC
LDA
(r) →
gradient of the density. If this expansion is further carried
const., but vXC (r) → −∞.
GGA
out to fourth order, as originally done by Gross and
Dreizler(59) and further developed by Perdew,(60) it also
4.3.2 Beyond the GGA introduces a term proportional to the square of the Laplacian
There seems, then, to exist a limit in the accuracy that of the density. The Laplacian term was also derived using
GGAs can reach. The main aspect responsible for this is a different route by Ghosh and Parr,(61) although it was
the exchange term, whose non-locality is not fully taken then dropped out when considering the gradient expansion
into account. A particularly problematic issue is that GGA only up to second order. More recently, a general derivation
functionals still do not compensate the self-interaction of the exchange gradient expansion up to sixth order,
completely. using second-order density response theory, was given by
This has motivated the development of approximations Svendsen and von Barth.(83) The fourth-order expansion of
that combine gradient-corrected functionals with exact, HF- that paper was then used by Perdew et al.(84) to construct a
type exchange. An example is the approximation known as practical meta-GGA that incorporates additional semi-local
B3LYP,(76 – 78) which reproduces very well the geometries information in terms of the Laplacian of the density. The
and binding energies of molecular systems, at the same philosophy for constructing the functional is the same as
level of correlated quantum chemistry approaches like MP2 that of PBE, namely, to retain the good formal properties
or even at the level of CI methods, but at a significantly of the lower-level approximation (the PBE GGA in this
lower computational cost. Even if the idea is appealing case), while adding others.
and physically sensible, at present there is no rigorous The gradient expansion of the exchange enhancement
derivation of it, and the functional involves a number of factor FX is
fitting parameters.
10 146 2 73
In the past few years there have been serious attempts to FX (p, q) = 1 + p+ q − qp
go beyond the GGA. Some are simple and rather successful, 81 2025 405
although not completely satisfactory from the theoretical + Dp 2 + O(∇ρ6 ) (70)
point of view, because they introduce some fitting parame-
ters for which there are no theoretical estimates. These are where
the meta-GGA described in the next subsection. A very |∇ρ|2
p=
interesting approach that became very popular in recent [4(3π2 )2/3 ρ8/3 ]
18 Electronic structure of large molecules

is the square of the reduced density gradient and all possible densities, is chosen in Reference 84 (exactly
as in References 68, 69). The coefficient D still remains
∇ 2ρ undetermined. In the absence of theoretical estimations,
q=
[4(3π2 )2/3 ρ5/3 ] PKZB proposed to fix D by minimizing the absolute error
in the atomization energies for a molecular data set. The
is the reduced Laplacian of the density. value so obtained is D = 0.113. The meta-GGA recovers
The first two coefficients of the expansion are exactly the exact linear response function up to fourth order in
known. The third one is the result of a difficult many-body k/2kF . This is not the case of PBE-GGA (and other
calculation, and has only been estimated numerically by GGA’s), which recovers only the LSDA linear response,
Svendsen and von Barth, to an accuracy better than 20%. and at the expense of sacrificing the correct second-order
The fourth coefficient D has not been explicitly calculated gradient expansion.
till date. The correlation part of the meta-GGA retains the correct
In the same spirit of PBE, Perdew, Kurth, Zupan and formal properties of PBE GGA correlation, such as the
Blaha (PKZB) proposed an exchange enhancement factor slowly varying limit and the finite limit under uniform
that verifies some of the formal relations and reduces to scaling. In addition, it is required that the correlation energy
the gradient expansion (70) in the slowly varying limit of be self-interaction-free, that is, to vanish for a one-electron
the density. The expression is formally identical to that of system. PKZB proposed the following form:
PBE:
κ ECMGGA [ρ↑ , ρ↓ ]
FXMGGA (p, q̄) = 1 + κ − (71)  
1 + x/κ  2 

 

 
τW
 
σ
  σ
where = ρεGGA (ρ↑ , ρ↓ , ∇ρ↑ , ∇ρ↓ ) 1 + C    
#  C
   
  $ 


τσ
10 146 2 73 1 10 2 2 σ
x= p+ q̄ − q̄p + D + p 
81 2025 405 κ 81


 W 2 

is a new inhomogeneity parameter that replaces µp in τσ
− (1 + C) ρσ εC (ρσ , 0, ∇ρσ , 0) dr (72)
GGA
τσ 

PBE. The variable q in the gradient expansion (the reduced σ 

Laplacian) is also replaced by a new variable q̄ defined as

3τ[ρ] 9 p where εGGA is the PBE-GGA correlation energy density

q̄ = − − C
[2(3π ) ρ ] 20 12
2 2/3 5/3 and τW σ is the von Weiszäcker kinetic energy density
given by expression (33) above, which is exact for a one-
which reduces to q in the slowly varying limit but remains electron density. Therefore, the correlation energy vanishes
finite at the position of the nucleus, while q diverges (an for any one-electron density, irrespectively of the value
unpleasant feature of most GGA). In the above expression, of the parameter C. For many-electron systems, the self-
τ[ρ] = τ↑ + τ↓ is the kinetic energy density for the non- interaction cancellation is not complete, but the error is
interacting system, with shifted to fourth order in the gradient, thus having little
effect on systems with slowly varying density. As in the
1 
occup
case of the exchange term, there is no theoretical estimate
τσ = |∇ψασ (r)|2
2 α available for the parameter C. Perdew et al. obtained
a value of C = 0.53 by fitting it to PBE-GGA surface
σ =↑, ↓. The connection between τ and the density is given correlation energies for jellium. Atomic correlation energies
by the second-order gradient expansion also agree, but are slightly less accurate. Just as an example,
the correlation energy for He is −0.84 H in LSDA, −0.68 H
3 1 |∇ρ|2 1 in PBE-GGA and −0.48 H in meta-GGA (MGGA), which
τGEA = (3π2 )2/3 ρ5/3 + + ∇ 2ρ basically coincides with the exact value.(86)
10 72 ρ 6
Unlike the PBE-GGA, the meta-GGA has two fitted
The formal conditions requested for this functional are (i) parameters, C and D. The reason for it is actually the
the spin-scaling relation, (ii) the uniform density-scaling unavailability of first-principles theoretical estimates for
relation(85) and the Lieb–Oxford inequality.(70) Actually, them. Other authors have proposed similar expansions con-
a value of κ = 0.804 in equation (71), corresponding to taining the kinetic energy density in addition to the density
the largest value ensuring that the inequality is verified for gradients. These, however, took the road of constructing the
 Density functional theory: Basics, new trends and applications 19

functional in a semiempirical way, by fitting a large num- same density, one can employ a continuous sequence
ber of parameters (of the order of 10 or 20) to chemical of partially interacting systems with the same density as
data, instead of using theoretically calculated values.(87,88) the fully interacting one. In this way, by starting from
The quality of the results of different meta-GGA function- the non-interacting system, the electron–electron Coulomb
als is quite similar. An assessment of the general quality interaction is gradually switched on and the system evolves
of the PKZB meta-GGA in comparison with GGA and continually towards the fully interacting system, always
hybrid EXX – GGA models of the B3LYP type – has maintaining the same electronic density. This procedure has
been published very recently.(89) The conclusion is that been named the adiabatic connection. Since the electronic
the kinetic energy density is a useful additional ingredi- density for both interacting and non-interacting systems is
ent. Atomization energies are quite improved in PKZB the same, and Hohenberg–Kohn theorem states that this
meta-GGA with respect to PBE-GGA, but unfortunately, density is univocally determined by the potential for any
geometries and frequencies are worsened. In particular, form of the electron–electron interaction (in particular,
bond lengths are far too long. Adamo et al.(89) argued that full Coulomb and no interaction at all), the electronic
a possible reason could be that in this functional the long- problem can be re-casted in the form of a non-interacting
range part of the exchange hole, which would help localize problem with the same density of the interacting problem.
the exchange hole, thus favouring shorter bond lengths, The potential, however, has to be different because the
is missing. Intriguingly enough, one of the semiempiri- interaction is different.
cal meta-GGA functionals(88) gives very good geometries The OPM is useful because it deals with the following
and frequencies. According to the preceding discussion, problem: having a general expression for the energy,
this effect on geometries should be due to the non-local which is a functional of the orbitals, it searches for the
properties of the exchange functional, a factor that the optimum potential whose eigenorbitals minimize the energy
kinetic energy density, being still a semi-local object, can- expression. The KS scheme can be viewed from the OPM
not account for. Therefore, this agreement must originate perspective, as a special case.
in error cancellations between exchange and correlation. Mathematically, this can be formulated in the following
way:
 
∇2
5 EXACT EXCHANGE: THE OPTIMIZED − + vR [ρ](r) ϕjσ (r) = εjσ ϕjσ (r) (73)
2
POTENTIAL METHOD
where the orbitals ϕjσ (r) = ϕjσ [ρ](r) are also functionals of
The one-to-one correspondence between electronic density the density, although implicitly through the potential vR [ρ].
and external potential embodied into Hohenberg–Kohn’s The energy of such a non-interacting electronic system can
theorem suggests that the variational problem of minimizing be written as
the energy functional could be also formulated for the 
potential, instead of the density. Historically, this idea EvR [ρ] = TR [ρ] + ρ(r)vR [ρ](r) dr (74)
was proposed in 1953 by Sharp and Horton,(90) well
before the formulation of DFT, and received the name of with
Optimized Potential Method (OPM). The formal proof of  
N 

σ
this equivalence was given later on by Perdew et al.(91,92) ∇2
TR [ρ] = ϕjσ∗ (r) − ϕjσ (r) dr (75)
This idea proved very useful in the context of DFT, 2
σ j =1
because one of the main limitations of KS theory is that
even though the exact exchange-correlation energy is a the exact kinetic energy of non-interacting electrons.
functional of the density, unfortunately this functional is not Coming back to the fully interacting system, the energy
explicitly known. This is the reason why approximations to functional can be written in terms of TR [ρ] by displacing all
this term are needed and have been proposed at different the ignorance about the electronic many-body problem into
levels of accuracy. the energy term EXC [ρ]. This contains the exchange con-
It is to be noticed that the same happens with the tribution and, in addition, all correlation effects including
kinetic energy functional, which is not explicitly known those omitted in the kinetic term:
in terms of the density. However, in the case of non- 
interacting electrons, the exact expression in terms of EKS [ρ] = TR [ρ] + ρ(r)v(r) dr
the orbitals is well known. This is actually the basis

for KS theory.(31) In order to visualize the mapping of 1 ρ(r)ρ(r )
+ dr dr + EXC [ρ] (76)
the interacting system to a non-interacting one with the 2 |r − r |
20 Electronic structure of large molecules

This last expression is simply the definition of EXC [ρ]. Notice that this construction of the mapping onto a non-
Now, by using the variational principle that ρ(r) minimizes interacting system is completely general and it relies only
EKS [ρ], we obtain on the assumption of v-representability of the interacting
 electronic density. In particular, if an explicit dependence of
δTR [ρ] ρ(r ) δE [ρ] EXC [ρ] on the density (or the density and its gradient as in
+ v(r) + 
dr + XC =0 (77)
δρ(r) |r − r | δρ(r) GGA or density, gradient and Laplacian, as in meta-GGA)
is assumed, the conventional KS scheme is recovered.
Using the non-interacting equation (73), and first-order per- The above equations are quite general and can be used
turbation theory for calculating the variation of the single- even when an approximate expression for EXC [ρ] is given
particle eigenvalues, it can be shown that the variation of as an implicit functional of the density, for example, in
TR [ρ] with respect to the electronic density is
terms of the orbitals. In order to deal with orbital-dependent
δTR [ρ] functionals, we have to calculate the density variation of
= −vR [ρ](r) (78) EXC [ρ] via its variation with respect to the orbitals. This
δρ(r)
can be done by applying the chain rule in the context of
namely, that density and effective potential are conjugated functional derivation:
fields. This, in conjunction with equation (77), gives rise Nν   
to the desired expression for the non-interacting reference δEXC [ρ]   δEXC [ρ]
VXCσ (r) = =
potential: δρσ (r) ν i=1
δϕiν (r )
  
ρ(r ) δϕiν (r )
vR [ρ](r) = v(r) + dr + VXC [ρ](r) + const. × dr + c.c. (82)
|r − r | δρσ (r)
(79)
where where we have included a spin index (σ) to be consistent
δEXC [ρ] with the spin-dependence of the exact exchange functional.
VXC [ρ](r) = (80) But the orbitals are connected only implicitly with the
δρ(r)
density, through the reference potential. Therefore, we have
is the definition of the exchange-correlation potential. to introduce another intermediate step of derivation with
Therefore, if the exact exchange-correlation energy func- respect to vR [ρ]:
tional is used, then the density obtained from equation (77)
 Nν     
is the exact interacting density. δEXC [ρ] δϕiν (r )
The potential vR [ρ] in equations (73) and (79) is chosen VXCσ (r) =
νµ i=1
δϕiν (r ) δvRµ (r )
so that the two energy functionals (74) and (76) have  
the same minimizing density ρ. Further, the constant in δvRµ (r )
equation (79) is chosen so that the two functionals at × dr dr + c.c. (83)
δρσ (r)
their common minimizing density have equal values. This
fact can be exploited to cast the variational problem in The second factor in the product is the variation of the
a tractable form in terms of the non-interacting reference non-interacting orbitals with respect to the potential, which
system. The solution can then be obtained by solving can be calculated using first-order perturbation theory:
equation (73) and constructing the density according to the
∞
# ∗  $
usual expression for non-interacting electrons, whose wave δϕiν (r )  ϕkµ (r )ϕkµ (r )
function is a single Slater determinant of the orbitals ϕj (r), = δν,µ ϕiµ (r )
δvRµ (r ) k=1,k =i
ε iµ − εkµ
that is,
= GRiµ (r , r )ϕiµ (r )
σ

N (84)
ρ(r) = |ϕjσ (r)|2 (81)
σ j =1
where GRiσ (r , r ) is the Green’s function of the non-
The price for this simplification from an interacting many- interacting system, given by
electron problem to an effective non-interacting problem is ∞
 ϕkσ (r )ϕkσ (r )
that the effective potential defined by equation (79) depends GRiσ (r , r ) = (85)
on the electronic density, which is constructed with the k=1,k =i
εiσ − εkσ
solutions of the single-particle equations. Therefore, this
problem has to be solved in a self-consistent way, by The third factor is the variation of the potential with respect
ensuring that the input and output densities do coincide. to the density, which is the inverse of the linear response
 Density functional theory: Basics, new trends and applications 21

function of the reference system χR , defined as that is,

δE[vRσ ]
δρσ (r) =0 (91)
χRσ,µ (r, r ) = δσ,µ (86) δvRσ (r)
δvRµ (r )

This is a well-known quantity for non-interacting systems, which, by applying again the functional chain rule, can
which is related to the Green’s function above by be shown to be strictly equivalent to the original Hohen-
berg–Kohn principle, stating that the energy functional is

Nσ a minimum at the ground state density.(91,92) The formu-
χRσ (r , r) = GRiσ (r , r)ϕiσ (r )ϕiσ
∗
(r) + c.c. (87) lation described above was originally proposed by Görling
i=1 and Levy (GL).(97,98)
It can be easily seen that if the XC energy functional
As, Rfrom equation (85) GRiσ is orthogonal to ϕiσ , we have depends explicitly on the density, and not on the orbitals,
χσ (r , r) dr = 0, and the linear response function is not XCiσ (r) = µXCσ [ρ](r) is also an orbital-independent
then uOEP
invertible. In a plane-wave representation, this means that functional (an explicit functional of the density), and it
the G = 0 component is zero and, therefore, it should be coincides with the usual XC potential in KS theory. In that
excluded from the basis set.(93,94) This is simple to do in case we can choose VXCσ (r) = µXCσ [ρ](r), and the OEP
plane waves, but somewhat more complicated when dealing equation is automatically satisfied. With this choice, the
with localized basis sets.(95) original definition of the reference potential (equation 79)
If the restricted χ̃Rσ (r , r) (no G = 0 component) is and the traditional KS scheme are recovered.
considered, then the expression for the local XC potential If this is not the case, then the OEP integral equa-
corresponding to orbital-dependent functionals assumes the tion (89) or, equivalently, equation (88) has to be solved
form: for the XC potential, which is then used to construct
Nσ    the reference potential vR (r). Orbital-dependent correlation
δEXC [ρ] R  
VXCσ (r) = )
Giσ (r , r )ϕiσ (r ) + c.c. functionals are not very common. Notable exceptions are
i=1
δϕiσ (r Colle–Salvetti’s functional(99,100) and the early Perdew and
 −1  Zunger’s attempt at correcting the self-interaction prob-
× χ̃Rσ (r , r) dr dr (88) lem of the LDA by considering orbital-dependent XC
functionals [self-interaction correction (SIC) approach].(42)
where the inversion step has to be carried out explicitly, and The exchange term, however, is perfectly well known
this is typically a rather expensive numerical operation. as an orbital-dependent functional, as given by the Fock
An equivalent formulation can be obtained by multiply- expression:
ing both sides of equation (88) with χRσ (r , r), integrating in
r, and replacing the expression (87) for the response func- Nµ  
tion. In this case, we obtain the following integral equation: 1  
EX [ρ] = −
2 µ=↑,↓ j,k=1
Nσ 
  OEP  
∗
ϕiσ (r ) VXCσ (r ) − uOEP 
XCiσ (r )
ϕj∗µ (r)ϕkµ
∗
(r )ϕkµ (r)ϕj µ (r )
× dr dr (92)
i=1 |r − r |
× GRiσ (r , r)ϕiσ
∗
(r) dr + c.c. = 0 (89)
so that its orbital functional derivative is
where we have defined
  ∗
δEX [{ϕkµ }] Nσ
ϕiσ (r)ϕj σ (r)
1 δEXC [{ϕj τ }] =− ϕj∗σ (r ) dr (93)
XCiσ (r) = ∗
uOEP (90) δϕiσ (r ) j =1
|r − r |
ϕiσ (r) δϕiσ (r)

The integral equation (89) is the so-called optimized effec- and uOEP
XCiσ (r) is obtained by using equation (90).
tive potential (OEP) equation, and was originally proposed As in the conventional KS theory, the OEP equations
by Sharp and Horton in 1953,(90) and re-derived and applied have to be solved self-consistently because the solution
to atomic calculations by Talman and Shadwick in 1976.(96) depends on the single-particle orbitals. This scheme can
However, in these works it was obtained as the solution to be implemented in its exact form, as it has been done for a
the problem of minimizing the HF energy functional (76) number of systems, or can be re-casted in an approximate,
with respect to the non-interacting reference potential vR [ρ], more easily solvable form.
22 Electronic structure of large molecules

5.1 The Krieger–Li–Iafrate approximation has been used in a number of different contexts, mainly
in atomic and molecular systems.(110) In order to distin-
The OEP formulation, although known for a long time, was guish the two approaches, we shall use EXX to refer to
not used in practical applications until the early nineties, the solution of the full integral equation and KLI to refer
except for the early work of Talman and Shadwick. The to the above approximation. Both are based on the OEP
main reason was that solving the full integral equation philosophy.
numerically was perceived as an extremely demanding task,
which could only be achieved for very symmetric (spheri-
cally symmetric) systems. In fact, even up to 1999, the exact 5.2 Properties of exact exchange in the OEP
exchange OEP method was used only to study spherically approach
symmetric atoms,(101 – 104) and subsequently solids within
the atomic sphere approximation (ASA).(105 – 108) The formal properties of both the EXX approach and the
In 1992, Krieger, Li and Iafrate(109) proposed an alterna- KLI approximation have been considered in detail by Grabo
tive, still exact expression for the OEP integral equation, et al.(110) The most important ones, concerning the EXX (no
using the differential equation that defines the Green’s correlation) functional are the following:
function of the reference system. After some algebraic
manipulation, the XC potential assumes the form: 1. Owing to the exact cancellation of the self-interaction,
the EXX and KLI exchange potentials decay (correctly)
1 
N σ
as −1/r at long distances, in vacuum regions. This is
OEP
VXCσ (r) = |ϕ (r)|2
2ρσ (r) i=1 iσ one of the most severe shortcomings of the LDA and
  OEP  GGA, which leads to a number of unpleasant features
× vXCiσ (r) + V̄XCiσ − ūXCiσ + c.c. (94) such as the incorrect dissociation limit for molecules
or the image potential at metallic surfaces showing an
where incorrect decay into the vacuum.
1   2. The principle of integer preference,(111) which states
vXCiσ (r) = uXCiσ (r) − ∇· ψ∗iσ (r)∇ϕiσ (r) (95) that the exact XC potential considered as a continuous
|ϕiσ (r)| 2
function of the number of electrons, is discontinuous
ψiσ (r) are the solutions of the inhomogeneous KS-like at integer values – so that integer numbers of elec-
equation: trons are preferred – is verified both in EXX and
in KLI. This is another great advantage over LDA
, -  OEP
ĥRσ − εiσ ψiσ (r) = − VXCiσ (r) − uXCiσ (r) and GGA, none of which verifies this property. A
consequence of the lack of integer preference is the
 OEP 
− V̄XCiσ − ūXCiσ ϕiσ (r) (96) well-known underestimation of the band gap in bulk
semiconductors.
and the bars indicate averages over |ϕiσ (r)|2 . 3. At variance with HF theory, which at first glance
This formulation is strictly equivalent to equation (89), may seem equivalent, the long-range decay of the
but it admits a reasonably well-controlled mean-field exchange potential into vacuum regions goes, correctly,
approximation, which is obtained by neglecting the sec- as −1/r for all states, irrespectively of whether they are
ond terms in equation (95), that is, vXCiσ (r) = uXCiσ (r). occupied or empty. In HF, the potential corresponding
Even if this is not generally true, it can be shown that to occupied orbitals decays as −1/r, but for empty
the average of the neglected term is zero.(110) This is orbitals it decays exponentially. Therefore, the HF
the Krieger–Li–Iafrate (KLI) approximation, and the KLI potential can support very few (if any) empty bound
expression for the XC potential is states. In the same way, the exchange potential in the
LDA and GGA decays exponentially for all states,
1 
N σ occupied and empty. Again, only a few bound excited
KLI
VXCσ (r) = |ϕ (r)|2 states are possible and, moreover, many negatively
2ρσ (r) i=1 iσ
  KLI  charged ions are not even bound. The OEP solves
× uXCiσ (r) + V̄XCiσ − ūXCiσ + c.c. (97) these problems and has been shown to support a whole
Rydberg molecular series.(95) In addition, in HF all the
where the averaged XC potential is obtained as the solu- occupied orbitals decay exponentially with the same
tion of a set of coupled linear equations. This equation is exponent, while in the OEP each orbital decays with
much simpler to solve than the original OEP equation and its own exponent, as it should be.
 Density functional theory: Basics, new trends and applications 23

4. In exchange-only calculations, that is, when neglect- η(r) = 1 + dρ−1/3 (r) (100)
ing the correlation term, the spin-unrestricted Hartree– −5/3 −1/3
ρ (r) exp[−cρ (r)]
Fock (SUHF) approach gives the variationally lowest ξ(r) = (101)
ground state energy. Therefore, ESUHF is a lower bound η(r)
for any other exchange-only scheme. It has been shown
with a = 0.04918, b = 0.132, c = 0.2533 and d = 0.349.
that the EXX approach gives energies EEXX , which are
A purely density-dependent functional – instead of orbit-
only marginally larger than ESUHF .(112) This fact might
al-dependent – inspired in Colle–Salvetti’s formula has
appear as an inconsistency because both the SUHF and
been proposed by Lee, Yang and Parr,(67) and became
the EXX approaches are exact. However, the nature
very popular in standard GGA calculations under the
of the HF (non-local X potential) and DFT (local X
name of LYP correlation functional (see Section 4.3). It
potential) approaches is different, in the sense that
is commonly used in conjunction with Becke exchange(66)
the partition between exchange and correlation ener-
(BLYP functional), and also in hybrid DFT-HF schemes
gies is different in both schemes. Therefore, this small
such as B3LYP,(76 – 78) which are nowadays the way of
difference is related to the fact that correlation has
choice in quantum-chemical (QC) applications.
been neglected and should disappear when the exact
The original CS functional is, however, orbital-dependent
correlation is considered. In the next subsection we dis-
and can be used in an OEP scheme in conjunction with the
cuss some attempts to combine EXX and/or KLI with
EXX or KLI exchange functional. This program has been
approximate orbital-dependent correlation functionals.
carried out by Grabo et al.(110) (KLICS scheme), and the
results compared with density-dependent GGA, correlated
QC calculations, and exact values, for atomic and molecular
5.3 Orbital-dependent correlation functionals systems.
The spirit of Colle–Salvetti functional is, however, the
Although the EXX and KLI formulations are general for same as that of GGA, in the sense that it is constructed in
orbital-dependent XC functionals, the exact correlation terms of the local density and the gradients of the orbitals
functional is not known. This approach has been normally and the density. Therefore, it is an approximation that is
implemented in its exchange-only form(102,103,113,114) or valid in the case of slowly varying orbital densities and
augmented with the usual LDA and/or GGA density-only cannot constitute a solution for the correlation problem
correlation functionals.(93,94,115) in the general case. A possible consistent approach to
An orbital-dependent correlation functional has been deal appropriately with the correlation functional will
proposed by Colle and Salvetti,(99,100) starting from a be discussed in the following section. Nevertheless, this
correlated Jastrow-type many-electron wave function, and is a very active area of research, and several different
performing a series of approximations. The expression of approaches are currently being explored. In the last section
the correlation energy, as given by Lee, Yang and Parr,(67) we shall review some applications and compare the results
is the following: obtained using different exchange functionals (LDA, GGA,
 M-GGA and exact), and also combined with the appropriate
EC ({ϕj σ }) = −ab dr γ(r)ξ(r)
CS correlation functionals.
# $
  1
× ρσ (r) |∇ϕiσ (r)| − |∇ρσ (r)|2
2

σ
4 6 TOWARDS AN ACCURATE
i
 CORRELATION FUNCTIONAL
− ab dr γ(r)ξ(r)
# $ While the exchange contribution to the energy and potential
1 1 is well known, and has been usually approximated in
× − ρ (r)∇ ρσ (r) + ρ(r)∇ ρ(r)
2 2
4 σ σ 4 DFT because of its computational cost, the correlation
 contribution, in the general case of an inhomogeneous
ρ(r)
− a dr γ(r) (98) electronic system, is still unknown in a closed form. A few
η(r)
simple cases, such as the homogeneous electron gas and
where some atomic systems (especially He), have been studied
numerically very accurately, so that nowadays there are a
ρ↑ (r)ρ↓ (r) number of benchmarks to compare the quality of different
γ(r) = 4 (99)
ρ2 (r) approximations to correlation.
24 Electronic structure of large molecules

In order for an approximation to exchange and correlation limit would be to make a connection with the RPA, which
to be reliable, it is needed that both terms are treated is known to treat properly long-range correlations.(117) The
consistently. This is one of the main achievements of the short-range behaviour of the RPA correlation is, however,
LDA, where both energy functionals are approximated in rather poor.(118) Therefore, a different approach is needed in
the same limit of a locally homogeneous system. Therefore, that region. One possibility is to connect with the standard
even if separately each term is not particularly accurate, GGA(119,120) – or new variants of the GGA(121) – at short
the sum of the two terms is rather precise, at least distances, where GGAs are rather accurate. Nevertheless,
from the energetic point of view. The same can be said other partitionings are also possible.(122,123)
about properly constructed GGAs, where exchange and These approaches for finding a good approximation to
correlation are treated consistently.(68,69) In this case, the the (orbital-dependent) correlation functional can be put
reference system is an electron gas whose density is slowly on a sounder basis by making a connection with quantum
varying in space. many-body theory. In quantum chemistry, the perturbative
Proper GGAs are, by construction, better approximations approach known as Møller-Plesset (MP) theory has been
than the LDA. They should include the LDA as a limiting used since the early days of quantum mechanics.(6) It
case when the density gradient terms in the functional starts from the HF solution and introduces electronic
are neglected. It is to be remarked that some popular correlations in a perturbative way on top of HF. The
GGA, like the BLYP functional, do not satisfy this limit. perturbation expansion is now customarily carried out to
In addition, GGAs are constructed to fulfil some other 2nd order (MP2), and sometimes also to 4th order (MP4). In
exactly known conditions that the LDA does not, like DFT, the analogous density functional perturbation theory
the correct linear response limit and the correct limit for (DFPT) has been developed and discussed in detail by
slowly varying densities, amongst others. Therefore, if a Görling and Levy(124,125) (GL theory). In the following, we
proper GGA performs worse than the LDA for some system sketch the salient features of GL theory and discuss some
when compared with experimental data, this means that applications.
the good performance of the LDA is actually fortuitous, Görling and Levy considered the many-body Hamilto-
and it is based mainly on cancellation of errors. This nian
happens, for example, for noble metals, where the LDA
lattice constant is virtually on top of the experimental
= T
+ V
H
+ λU

(102)
λ λ ee
value, while PBE gives a value that is expanded by a few
percent. where 0 ≤ λ ≤ 1 is a coupling constant representing the
This is, then, the main reason for the seemingly poor strength of the electron–electron interaction and V
is
λ
performance of exact exchange combined with local cor- constrained to be the local external potential that keeps
relation functionals, that is, EXX-LDA and EXX-GGA
, invariant
the ground state density ρλ , corresponding to Hλ
approaches. This is more dramatic in atoms and molecules for any value of λ and equal to the density of the fully
than in solids because, there, local correlation function- interacting system, ρ. The total energy for interaction
als based on the homogeneous electron gas are partic- strength λ is written as
ularly inappropriate. Even the Colle–Salvetti pair corre-
lation function, where the radius of the correlation hole Eλ [ρ] = E (0) [ρ] + λE (1) [ρ] + Ecλ [ρ] (103)
is parametrized to fit the correlation energy of the He
atom,(116) strongly departs from the correct long-range
where E (0) [ρ] is the energy associated with the non-
dependence. In fact, for the uniform electron gas it
interacting Hamiltonian, λE (1) [ρ] = Ex [ρ] is the exact
should decay as r −4 , but the assumption of a strong
exchange energy given by the Fock expression and Ecλ [ρ]
Gaussian damping for the pair Jastrow correlation fac-
is formally given by the expression
tor (the pair density matrix) prevents against such a
decay. . / . /
This, together with other studies, indicates that the Ecλ [ρ] = ψλ |T
+ λU

|ψ − ψ |T
+ λU
ee λ 0

|ψ
ee 0 (104)
main limitation of such a hybrid (exact exchange – local . / . /
correlation) approach is that the long-range tail of exchange, where ψλ and ψ0 minimize T
+ λU
and T
, respectively.
ee
which is treated exactly, is not properly compensated by a Owing to exact scaling relations derived by Levy
similar long-range tail of opposite sign in the correlation and Perdew,(126) the correlation energy at any interaction
term. Therefore, there is a clear need for improvement at strength λ can be written in terms of the correlation energy
the level of correlation functionals, which should now take of the fully interacting system (λ = 1), but at a uni-
into account this aspect. A possible route to include this formly scaled density ρ1/λ (x, y, z) = λ−3 ρ(x/λ, y/λ, z/λ).
 Density functional theory: Basics, new trends and applications 25

The scaling relation is The general expression for Ec(n) has been given by
Görling and Levy.(124,125) While, in general, the perturba-
1 λ tive terms in the above expansion are complicated and com-
Ec [ρ1/λ ] = E [ρ] (105)
λ2 c putationally very expensive, the second-order term assumes
the familiar form known from the usual second-order per-
so that turbation theory, although applied to KS states
∞

−V
−V

|ψ0 |2
1 . / |ψ00 |U
Ec [ρ1/λ ] = 2 ψλ |T
+ λU
|ψ
λ Ec(2) [ρ] =− ee H x k
(110)
λ ee
Ek0 − E00
. / k=1
− ψ0 |T
+ λU

|ψ
ee 0
1  . / with ψ0k the kth excited state of the unperturbed Hamil-
=
−H
Eλ − E0 − ψ0 |H
|ψ (106) tonian – which are KS determinantal states – and Ek0 the
λ
λ2 0 0
corresponding energies. Obviously, the ground state is
excluded from the summation. This can be put in terms
where the λ-interacting Hamiltonian Ĥλ is partitioned in
is a two-
of KS single-particle orbitals, recalling that U ee
the following way:
and V
particle operator, while V
are single-particle oper-
H x
ators. The expression obtained is(127)

=H
H

λ 0
! N  " 1   |ϕi ϕj |v̂ee |ϕα ϕβ |2
 δE [ρ ] EcGL2 ({ϕj }) = −

+
+λ U −u(ri ) − vx (ri ) − λ c λ 4 ij αβ εα + εβ − εi − εj
ee
i=1
δρλ (ri )
  |ϕ |v̂ − fˆ|ϕ |2
α
(107) − i x
(111)
and the term inside curly brackets is treated as a perturba- i α
εα − εi
tion. Here, u(r) is the Hartree (direct Coulomb) potential
and vx (r) is the (local) exchange potential (in the sense where fˆ is the Fock-like, non-local exchange operator,
of the OEP discussed above). By considering the Taylor but formed with the {ϕj }, which are KS single-particle
series in λ around the non-interacting limit (λ = 0), and orbitals with associated eigenvalues εj . The indices i and j
being consistent in the order in λ (care has to be taken with correspond to occupied single-particle orbitals, while α and
the last term in the perturbing potential in equation (107), β indicate empty orbitals.
which depends on λ and on the correlation potential itself), As in MP theory, the computational cost of higher-
the following expansion is obtained: order terms in the series is prohibitively expensive, so
that normally it would not be possible to afford more than
∞
 second-order DFPT. This is sometimes called second-order
Ec [ρ1/λ ] = λn−2 Ec(n) [ρ] (108) Görling–Levy theory (GL2). As such, GL2 theory has been
n=2 used by Ernzerhof to calculate the energetics of atomiza-
tion of molecular systems,(128) and compared with the more
where the terms Ec(n) [ρ] are now calculated perturbatively traditional DFT and QC approaches such as LSD, GGA,
on the non-interacting system, for which the orbitals are HF and MP2. One important conclusion of this work is
known (e.g., the KS orbitals). The total correlation energy that if the perturbative series for the correlation is simply
is obtained as a coupling constant integration, which is a cut at the GL2 level, then the resulting atomization ener-
result of the adiabatic connection formula.(117) The latter gies are particularly bad, even worse than the local spin
states that the correlation contribution to the kinetic energy, density (LSD) values. Ernzerhof suggested that the one
which is neglected in the non-interacting reference system important shortcoming of this type of approximation is that
(KS), can be absorbed into the correlation potential energy some known exact limits, for example, the limit of very
(Ec ), by averaging the correlation energy corresponding to strongly interacting systems (λ → ∞), are not fulfilled. He
interaction λ from the non-interacting case λ = 0 to the then proposed an empirical re-summation of the series, in
fully interacting case λ = 1: the spirit of the GGA, so that these exact limits are ver-
ified. The comparison of these results with experimental
 1 ∞

atomization energies is very favourable, even improving
1 over MP2 results in some difficult cases such as the F2 and
Ec [ρ] = dλ Ec [ρ1/λ ] = E (n) [ρ] (109)
0 n=2
n−1 c O3 molecules.
26 Electronic structure of large molecules

Another application of the bare GL2 theory, combined be expanded in Taylor series for all values of the coupling
with EXX, was carried out by Engel et al.,(129) who constant 0 ≤ λ ≤ 1.
analysed the problem of the van der Waals binding of This hypothesis was tested by Seidl et al.(135) on the
closed-shell-atom dimers. This is a difficult benchmark basis of the atomization energies calculated by Ernz-
problem in many-body theory, which, ultimately, any erhof in the GL2 approximation.(128) Surprisingly, they
correlation functional should address. It has attracted a lot found that, except for a few notable cases such as H2
of attention in the past, since the early work of Zaremba and and CH4 , the radius of convergence of the perturba-
Kohn,(13) up to recent proposals.(122,130 – 133) The origin of tive series is always λc < 1. In some cases it can be
the van der Waals interaction between two non-chemically very small indeed (0.06 for B2 ). As already remarked
bonded fragments is the coupling of the electric field by Ernzerhof, this is the reason for the poor values
generated by fluctuations in the electronic density of one obtained for the atomization energies. They have also
fragment, with the density of the other fragment. At long shown, in a model calculation, how this problem of the
distances, this interaction should approach the classical radius of convergence is manifested numerically, in a
dipole–dipole interaction, which decays as R −6 , with R real calculation. Basically, the truncated series in equa-
the distance between fragments. Special cases of van der tion (108) behaves well until λc , where it departs from
Waals systems are dimers of closed-shell atoms such as the expected behaviour, and shows its tendency to diverge
He2 . (diverges only if the infinite series is considered). Succes-
While most works concentrated precisely on this long- sive terms in the expansion have different signs, so that
range behaviour, the goal of a correlation functional is the series oscillates upon adding more and more terms.
to reproduce the correct behaviour of the whole poten- In the correlation energy, this is reflected as an oscilla-
tial energy surface, especially binding energies and bond tory behaviour as a function of the number of terms in the
expansion.
lengths. This analysis has been presented by Engel et al.(129)
Seidl et al. generalized, then, the re-summation ideas
for the case of He2 and Ne2 , where a comparison between
proposed by Ernzerhof, in the spirit of GGAs. They did
LDA, HF, KLI (x-only), MP2 and KLI-GL2 calculations
this by asking the correlation functional to verify the limit
and exact results is reported. It is very clear that the LDA
of very strong interaction (λ → ∞), which is the region
severely overbinds, while HF and x-only KLI do not bind
of interaction strengths where the perturbative expansion is
the dimers at all. Therefore, correlation is confirmed as a
likely to fail. In this limit, the electronic positions become
crucial ingredient. In fact, the two correlated approaches,
strictly correlated and give rise to Wigner crystallization.
MP2 and KLI-GL2, bind the dimer quite reasonably.
The correlation functional for such very strong interaction
Compared to exact results for Ne2 (De = 3.6 meV),(134)
was calculated in the point charge plus continuum (PC)
MP2 tends to underbind – giving a dissociation energy of model,(136) and given in terms of the electronic density and
2.3 meV, while KLI-GL2 overbinds (De = 8.3 meV). This its gradient. Then, an interpolation formula was proposed
is reflected at the level of geometry, where the MP2 bond (the interaction strength interpolation, or ISI), which has
length is somewhat long (6.06 bohr) against 5.48 bohr in the above limit for λ → ∞, and also reproduces the small
KLI-GL2, and 5.84 bohr in the exact calculation. Also, the λ limit, that is, EXC → Ex + EcGL2 , where Ex is the EXX.
dynamics follows the same trends: lower frequency in MP2 It was shown that, as a function of the number of terms
(23 cm−1 ) and higher in KLI-GL2 (46 cm−1 ), compared to in the expansion, the perturbative Taylor series oscillates
the exact value of 29 cm−1 . This indicates that correlations around the correlation energy given by the interpolation
to a level higher than second order are necessary to obtain formula.(135) Results presented for atomization energies of
a quantitative agreement. molecules, as in the case of Ernzerhof, compare extremely
These applications of GL2 theory show a general fea- well with experiment. The limit of weak interactions is
ture of perturbative expansions, namely, that unless the more dubious, because for the reference system, that is,
perturbation is really very weak, the simplest approach of the uniform electron gas, GL perturbation theory does not
cutting the expansion at some low order is not quite suc- work. The reason is that each term in the perturbative
cessful. There are two possible reasons for that: first, that series is divergent, and they have to be grouped together
the higher-order terms are not significantly smaller than into some re-summation scheme to give a finite correlation
the second-order one; second, that for some range of val- energy. Such a framework is provided, for example, by the
ues of the coupling parameter λ, the perturbative series RPA mentioned above. Then, some kind of interpolation
might even be divergent. This feature was already noticed scheme is needed in order to go from the weak to
by GL,(124,125) who clearly stated that their expansion was the strong interaction limit (from RPA to PC). Another
based upon the assumption that the correlation energy can interesting approach has been proposed by Casida,(137)
 Density functional theory: Basics, new trends and applications 27

who extended the OPM to include correlation, in addition 7 COMPARISON AND SALIENT
to exchange. He showed that the OEP method finds the FEATURES OF THE DIFFERENT
variationally best local potential for the non-interacting
reference system, which has the same density as the APPROXIMATIONS
interacting system. Actually, the exact self-energy function
(r, r , ω) is non-local both in space and time. The OEP 7.1 Atoms and molecules: exchange-only
method finds the best local approximation to this self-
energy. In order to do this, it is necessary to solve the In this subsection we summarize the properties of different
OEP integral equation, which involves calculating the exchange-only approaches (i.e., neglecting correlation), as
matrix elements of the self-energy operator. Casida has compared to exact HF calculations, xLDA and xGGA.
shown that this approach is closely connected with the 1. Total energies: EXX energies are marginally larger
Green’s function many-body theory developed by Sham than the SUHF ones. For instance, for atomic systems
and Schlüter.(138) Unfortunately, it seems that casting this the difference is only of a few ppm (2 ppm for Xe,
framework into a practical scheme for performing actual 40 ppm for Be, the worst case).(110)
calculations is not an easy task. In fact, no applications
of this approach are known, even if the formalism is very (a) The KLI approach, being an approximation to
appealing. the exact OEP (EXX), gives ground state ener-
As a summary, in Figure 2 we show a scheme of the gies that are higher. For atoms, however, the
different approaches to deal with the electronic kinetic, differences between EKLI and EEXX are also very
exchange and correlation contributions within DFT, includ- small – between 1 and 10 ppm – thus indicating
ing a number of approaches devised to deal also with that the KLI scheme is a rather good approxi-
electronic excitations within the DFT framework. mation for this class of systems.(110) The average

Local airy
gas (LAG)
Thomas
fermi Kohn Density-
Edge electron Sham
K, X, C dependent
gas X, C functionals

Orbital-dependent LDA
functionals DFT

LDA ... BP
PW
GGA
BLYP
PBE
Exact
exchange
C M-body PKZB
GW M-GGA
Adiabatic B
conn. VS
Homogeneous
electron gas
DFPT
Time-dep.
DFT
Ensemble
DFT
Excitations

Figure 2. Schematics of the different approximations devised to deal with the different energy terms in density functional theory. Also,
a few approaches to deal with electronic excitations are included: Ensemble DFT,(22 – 25) the adiabatic connection,(29,30) time-dependent
DFT(26 – 28) and many-body Dyson’s equation within the GW approximation.(54) The letters K, X and C indicate which parts of the
energy are approximated (kinetic, exchange and/or correlation).
28 Electronic structure of large molecules

error with respect to EXX, for 1 ≤ Z ≤ 54, is low-Z atoms (H and He are actually exact within
3.1 mH. the SIC approach). However, results worsen
(b) Usual gradient-corrected DFT local exchange notably for atoms with more electrons, leading
functionals (GGA) give much larger errors when to an average error of 300 mH.
compared to EXX results. Differences could 2. Eigenvalues and exchange potentials: The trends
be as large as several hundred ppm for low- illustrated for the ground state energies, that is, that
Z atoms, although they decrease significantly HF and EXX are virtually identical, KLI is a very
for high-Z species. The average error for the good approximation to EXX, xGGAs are one order
B88 functional(66) is 66 mH and for the PW91 of magnitude worse than KLI and xLDA is the worst
functional(65) is 53 mH. approximation of all, are preserved for most quantities
(c) The case of the x-only LDA functional (Slater’s of interest. For instance, the single-particle eigenval-
exchange) is somewhat worse, especially for low- ues, which are a measure of the quality of the exchange
Z atoms, owing to the severe lack of cancellation potential, show that KLI produces a high-quality VX (r),
of the self-interaction. It is well known that the while LDA and GGA are seriously in error. Interest-
x-LDA energy of the H-atom is 0.457 H, instead ingly, standard SIC results for eigenvalues are much
of the exact value of 0.5 H. For low-Z atoms, better than the LDA and GGA ones. This can also be
differences can be as large as 50 000 ppm. For seen by directly inspecting the shape of the exchange
higher-Z atoms, as in GGA, these differences potential. The EXX and KLI exchange potentials decay
decrease, leading to an average error of 160 mH. (correctly) as −1/r in the asymptotic region, as well
(d) For comparison, Grabo et al.(110) reported results as the SIC potential, while LDA and GGA’s poten-
obtained by using the SIC scheme of Perdew and tials decay too fast, either exponentially or with the
Zunger,(42) implemented as in the KLI approxi- incorrect power law. In the inner region, the inter-shell
mation. This approach solves the self-interaction maxima (peaks) are poorly reproduced in GGA, and
problem, thus adjusting the discrepancies for the are almost absent in SIC. Also, the GGA potentials

0.2

0.0

−0.2

−0.4

−0.6
νx (hartree)

−0.8
Exchange potentials of helium

−1.0
−1/r
Exact
−1.2 LDA
Langreth−Mehl
Perdew−Wang ’86
−1.4
Perdew−Wang ’91
Becke ’88
−1.6

−1.8
0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0
r (units of a 0)

Figure 3. The exchange potential for the He atom within different approximations [Reproduced by permission of APS Journals from
C. Umrigar and X. Gonze (1994) Phys. Rev. A, 50, 3827.](139)
 Density functional theory: Basics, new trends and applications 29

exhibit an unphysical divergence at the position of the functionals combined with appropriate correlation function-
nucleus, which is related to a pathology of the gradi- als. Following Grabo et al.,(110) we consider the KLICS
ent expansion. In fact, the LDA is not divergent there. approach (KLI exchange and Colle–Salvetti correlation),
Figure 3 shows the shape of the exchange potential in self-interaction corrected LDA (SIC-LDA), GGAs (BLYP
the He atom for the xLDA and for different xGGA, as and PW91) and QC approaches.
reported by Umrigar and Gonze.(139) Clearly, none of
these functionals is particularly accurate. 1. Total energies: For atomic systems in the KLICS
3. Magnetic properties: Spin polarization densities are approximation, these are of very high quality, compa-
well described by the EXX, KLI and SIC functionals, rable to QC calculations. LDA performs rather badly
although the approximate schemes miss the fine shell for these systems, which is significantly improved by
the SIC-LDA. Even better results for the energetics
structure. GGAs give radically different results, even
are obtained within GGA, almost as good as KLICS
predicting non-magnetic properties for magnetic atoms.
and QC. Average errors for first-row atoms are 380 mH
For the magnetization density at the nucleus, which is
(LDA), 130 mH (SIC-LDA), 10 mH (GGA) and 5 mH
an important quantity to interpret NMR experiments,
(KLICS and QC). Similar trends are obtained for
the KLI approach is somewhat erratic; it is in fair
second-row atoms.
agreement for some elements, but for others it even has
Individually, EXKLI and ECCS are much closer to the
the wrong sign. All other approaches are also erratic
exact values than their GGA counterparts, for example,
and of poorer quality than KLI.
EXB and ECLYP . The KLI (or EXX, which is almost
4. Diatomic and polyatomic systems: Diatomic mole-
identical) exchange energy is lower than the GGA one,
cules are quite well described in the KLI approxi-
while the CS correlation energy is higher. The sum
mation, as compared to HF results.(110) Total energies
of the two terms, however, is quite similar in both
are within ≈0.01%, highest orbital eigenvalues within
approaches, thus leading to the well-known, remarkable
≈0.1% and multipole moments within a few percent.
cancellation of errors between exchange and correlation
An interesting exception is the N2 molecule, for which
energies. In the LDA, this derives from the fact that
the ordering of the two highest occupied molecular
the exchange-correlation hole corresponds to a well-
orbitals (HOMO), 1πu and 3σg , is reversed in all
defined physical system (the homogeneous electron
DFT approaches with respect to HF. The xLDA results
gas), so that a number of crucial sum rules are
are of much poorer quality, giving much higher total
automatically verified. The SIC-LDA approach does
energies (a few percent off), and HOMO eigenvalues
rather poorly in both terms.
that are twice as large as the HF values, due to the
2. Ionization potentials: When calculated from ground
wrong exponential decay of the LDA exchange poten-
state energy differences, ionization potentials are
tial. Interestingly, the KLI approach also reproduces
very well described both in KLICS and GGA (see
some well-known failures of HF, such as the instabil-
Table 1). Average deviations from experimental values
ity of Be2 , the Beryllium dimer (a particularly difficult
are around 10 mH. QC approaches are one order
case, where correlation is a crucial ingredient).
of magnitude better, and the SIC-LDA somewhat
Exact exchange methods for polyatomic molecules
poorer. The quality of the XC potential is measured
have been devised by Görling,(95) Ivanov et al.(113) and
by the ionization potential calculated from the
van Gisbergen et al.(114) The conclusions are basically
HOMO of the neutral atom (Janak’s theorem). Here,
the same as for diatomics.
the KLICS approach is definitely superior to the
GGA. KLICS values systematically overestimate the
7.2 Atoms and diatomic molecules: correlation ionization potential by less than 100 mH (typically
effects 10%), while GGAs give results that are roughly half
of the experimental values, that is, it underestimates by
Correlation is the smallest contribution to the total energy. 50%. SIC-LDA is of much better quality in this respect.
Its magnitude for atoms is around 30 times smaller than the The reason is that both KLICS and SIC-LDA reproduce
exchange one. This does not mean that correlation is not the correct −1/r behaviour of the exchange potential
important. In fact, several physical and chemical properties at long distances from the nucleus, and this is a crucial
depend actually on the potential rather than the energy, and aspect for the correct determination of the HOMO.
very different potentials can correspond to similar energies. Therefore, the XC potential is much better described in
Therefore, even if the energetics is correct, it is important KLICS than in GGA. Values for the first 18 elements
to assess also the quality of the correlation potential. We of the periodic table are presented in Table 1, for two
now summarize the results obtained for different exchange different GGA and for the KLICS approach. LSDA
30 Electronic structure of large molecules

Table 1. Ionization energies from ground state energy differences for the first 18 elements in the periodic table, in different
approximations: LSDA,(140) two different GGA PW91 and BLYP),(110) a meta-GGA (MGGA),(87) HF(140) and KLICS.(110) The last three
columns present ionization energies calculated from the neutral atom eigenvalues, in the PW91, BLYP and KLICS approximations.(110)
Experimental values are from Reference 141. Units are eV.

Z Exact LSDA PW91 BLYP KLICS HF MGGA P-e B-e K-e

He 2 0.903 0.892 0.903 0.912 0.903 0.862 0.910 0.583 0.585 0.945
Li 3 0.198 0.200 0.207 0.203 0.203 0.196 0.202 0.119 0.111 0.200
Be 4 0.343 0.331 0.333 0.330 0.330 0.296 0.337 0.207 0.201 0.329
B 5 0.305 0.315 0.314 0.309 0.314 0.291 0.306 0.149 0.143 0.328
C 6 0.414 0.429 0.432 0.425 0.414 0.396 0.416 0.226 0.218 0.448
N 7 0.534 0.548 0.551 0.542 0.527 0.513 0.544 0.308 0.297 0.579
O 8 0.500 0.508 0.505 0.508 0.495 0.437 0.504 0.267 0.266 0.559
F 9 0.640 0.659 0.660 0.656 0.621 0.577 0.645 0.379 0.376 0.714
Ne 10 0.792 0.812 0.812 0.808 0.767 0.728 0.799 0.494 0.491 0.884
Na 11 0.189 0.195 0.198 0.197 0.191 0.182 0.192 0.113 0.106 0.189
Mg 12 0.281 0.283 0.281 0.280 0.275 0.243 0.282 0.174 0.168 0.273
Al 13 0.220 0.220 0.221 0.212 0.218 0.202 0.214 0.112 0.102 0.222
Si 14 0.300 0.302 0.305 0.294 0.294 0.281 0.294 0.171 0.160 0.306
P 15 0.385 0.386 0.389 0.376 0.379 0.369 0.382 0.233 0.219 0.399
S 16 0.381 0.385 0.379 0.379 0.380 0.331 0.381 0.222 0.219 0.404
Cl 17 0.477 0.484 0.482 0.476 0.471 0.433 0.478 0.301 0.295 0.506
Ne 18 0.579 0.585 0.583 0.576 0.575 0.542 0.582 0.380 0.373 0.619

values are similar to the GGA ones and HF results are this agreement is mainly due to the improvement in
similar to KLICS, although somewhat reduced due to the description of the exchange, as testified by x-only
the absence of correlation. calculations. Interestingly, including the CS correla-
3. Electron affinities: These are more stringent a test tion actually worsens the results obtained in x-only.
for VXC . The wrong asymptotic behaviour in LDA This is an indication that the CS correlation poten-
and GGA in many cases prevents the very existence tial has the wrong sign, and this has been confirmed
of bound states in negative ions, so that electron by analysing directly VC for the He atom. The exact
affinities cannot be calculated. Owing to the proper correlation potential is positive, while all the approx-
asymptotic decay, KLICS and SIC-LDA support such imations (LDA, PW, LYP and CS) are negative and
bound states. However, the calculated electron affinities longer ranged. Moreover, gradient-corrected function-
are rather poor. If calculated from ground state energy als exhibit spurious divergences at the origin. Recent
differences, they are underestimated by approximately work by Tao et al.(116) proved that the CS functional
10 mH, although in some cases such as O, it is actually recovers only 25% of the correlation energy of
underestimated by as much as 40 mH. In B, it has even the uniform electron gas. Previous calculations, which
the wrong sign. Calculating them from the HOMO of were the basis for adopting CS and derived func-
the negative ion, electron affinities are overestimated tionals in DFT calculations, misleadingly suggested
roughly by a factor of 2, the agreement worsening for a much better quality. This indicates that there is a
increasing number of valence electrons. QC approaches clear need for improvement on the correlation func-
are extremely accurate in this respect. This indicates tional. Short-range correlations appear to be described
that the CS correlation potential is rather poor and very well, but it misses in the long-range part. This
needs further improvement. is not extremely important in atoms, but in more
4. Correlation functional: The CS correlation functional extended systems such as molecules and solids it can
has been studied more in detail by analysing the case of be crucial.
two-electron systems (Helium-isoelectronic series),(110) 5. Bond lengths: KLICS bond lengths in diatomic
where the KLI is essentially exact and the only error molecules are shortened with respect to LDA and GGA
introduced is due to correlation. The correlation energy values. Actually, these bonds as well as HF ones are
is clearly superior in CS with respect to GGA and too short compared to experiment. For instance, in N2 ,
SIC-LDA. The average error is around 5 mH, com- it goes from 2.068 bohr in LDA to 2.079 bohr in GGA,
pared to 50 to 100 mH in other methods. Also, HOMO and to 1.998 bohr in KLICS, while the experimental
eigenvalues are in much better agreement. However, value is 2.074 bohr, remarkably similar to the GGA
 Density functional theory: Basics, new trends and applications 31

Table 2. Atomization energies for a few selected molecules. The first column quotes experimental values (reproduced from Refer-
ence 84), the next three are for pure DFT approaches: LSD, a GGA (PBE) and a meta-GGA (PKZB),(84) the following two columns
report HF and correlated quantum-chemical MP2 results,(128) then orbital-dependent DFT KLICS,(110) a hybrid HF/DFT functional
(B3LYP)(88) and finally a DFT perturbative correlated approach (second-order Görling-Levy), bare and re-summed in the Interaction
Strength Interpolation (ISI).(135) Units are mH.

Mol. Exp. LSDA PBE MGGA UHF MP2 KLICS B3LYP GL2 ISI

H2 174.5 180.3 166.7 182.5 133.9 165.7 171.4 — 181.7 180.0

Li2 39.3 37.9 31.7 35.9 4.8 25.5 32.4 33.5 62.1 35.9
Be2 4.8 20.6 15.6 7.2 11.2 −1.6 −10.5 — 35.1 9.1
N2 364.0 427.1 387.6 365.2 183.3 368.1 287.3 365.6 545.0 373.9
F2 62.1 124.6 85.1 68.8 −15.9 111.6 −22.7 57.7 213.5 54.2
LiH 92.4 96.9 85.2 93.1 52.6 86.1 89.4 92.9 111.6 93.7
OH 169.6 197.9 175.0 171.8 108.4 165.7 — 172.3 204.0 173.1
HF 225.7 259.1 226.3 221.0 154.6 227.9 193.9 222.1 275.7 229.0
H2 O 370.0 424.9 373.2 366.7 245.4 366.5 — 368.1 436.6 375.6
NH3 473.9 537.5 480.8 476.2 318.7 462.1 — 478.4 541.8 479.5
CH4 668.2 737.2 669.0 671.1 522.7 661.3 — 670.4 723.5 674.7
CO 413.2 476.3 428.4 408.0 277.3 423.9 — 408.3 565.7 423.7
NO 243.7 316.2 273.9 252.6 84.5 242.2 — 248.0 422.3 251.6
Cl2 92.4 132.1 103.7 94.7 — — — 87.8 — —

value. The reason for this disagreement is, partly, the different approximations, including many-body pertur-
correlation functional, but more importantly, the zero- bative approaches (see previous section) are presented
point motion of the nuclei. In fact, the anharmonicity in Table 2.
of the potential energy surface along the bond is 7. Exchange-correlation potential: The quality of the
such that the quantum average of the bond length is KLICS XC potential is, nevertheless, superior to that
displaced towards larger values [2]. In the case of the of GGA. This can be seen by inspecting the HOMO
H2 molecule, the bond length is increased by 3%, and eigenvalues of diatomics, which should coincide with
the vibrational frequency decreases by 400 cm−1 by the ionization potential. For instance, in N2 , the HOMO
considering the quantum mechanics of the protons.(142) eigenvalue goes from 0.3826 H in the LDA to 0.3804 in
GGA results are remarkably close to experimental GGA, and to 0.6643 H in KLICS. Experimental value
values, but for the wrong reason. When corrected for is 0.5726 H. As for atoms, this behaviour can be traced
zero-point motion, KLICS results are in much better back to the correct asymptotic decay of EXX.
agreement with experiment. 8. Polyatomic molecules: An EXX method for poly-
atomic molecules, where the correlation term is treated
6. Dissociation energies: KLICS dissociation energies of
within the usual LDA and GGA approaches, has been
diatomics are disappointingly far from experimental
developed by Görling.(95) Pure EXX atomization ener-
values, except for a few exceptions. The left–right cor-
gies are close to HF values. Inclusion of correla-
relation error in molecules is well known in HF theory,
tion improves the agreement with experiment, better
and it is clearly inherited by the KLICS approach. The
in LDA than in GGA. However, as in the case of
correlation functional has to compensate properly the diatomics, the pure GGA results (X and C treated
long-range part of the exchange potential, so that the within GGA) are in much better agreement with exper-
combined XC hole is shorter ranged than the X and iment, which is worsened when improving the descrip-
C holes separately. Evidently, the CS functional is not tion of exchange. This clearly indicates that the good
adequate to solve this problem. This feature is simi- performance of GGA heavily relies on error cancella-
lar to the one appearing when considering meta-GGA tion between exchange and correlation. Very interesting
functionals. Dissociation energies are much better in results have been obtained for eigenvalue spectra. In
the LDA and GGA because of the better compen- contrast to LDA, GGA and HF spectra, the EXX-GGA
sation of exchange and correlation in the long-range spectrum is physically meaningful: the HOMO lies in
region. LDA is known to overbind molecules, and the correct energetic region (coinciding with the ioniza-
GGA reduces this overbinding tendency. In particular, tion potential), and it also exhibits a Rydberg molecular
GGA values are remarkably close to the experimen- series, which is absent in all other methods. This is a
tal ones. Results for a few selected molecules within consequence of the correct asymptotic behaviour of the
32 Electronic structure of large molecules

very small, especially for GGA. This proves that most

0
of the error arises from the energy functional and not
6σ∗ from the self-consistent density. An analysis of the
2π∗
individual contributions to the total energy reveals
that the LDA underestimates the absolute values of
the kinetic, Hartree, exchange and electron–nuclear
interaction by 1 to 2%. The reason for this is that LDA
5σ
1π
densities are homogeneous in excess. GGA energetic
contributions are, individually, much closer than the
[eV]

−20
4σ LDA ones, thus indicating the high quality of the
GGA density.
Since EXX and HF densities are practically identical,
the density-dependent terms of the energy (Hartree
and electron–nuclear) are almost equal. Kinetic and
exchange energies, however, are somewhat different.
This is because even if both densities are constructed
−40 from a single-determinantal wave function, the HF
3σ
and EXX single-particle states are solutions to dif-
PBE EXX-PBE HF
ferent equations: the EXX determinant minimizes the
Figure 4. Eigenvalue spectrum of the CO molecule [Reproduced kinetic energy and corresponds to a local effective
by permission of APS Journals from A. Görling (1999) Phys. Rev. potential, while the HF determinant minimizes the
Lett., 83, 5459.](95) sum of kinetic and electrostatic energies, and the
effective potential is non-local. As a consequence,
exchange potential, not only for occupied states but the kinetic energy is lower in EXX and the exchange
also for empty states. Figure 4 shows the eigenvalue energy is lower in HF.
spectrum of the CO molecule within the PBE, HF and 2. Cohesive energies: Within the HF approximation,
EXX-PBE approaches, as calculated by Görling.(95) cohesive energies of solids are severely underesti-
mated, by more than 1 eV/atom (20–30%). Exchange-
only calculations give cohesive energies that are virtu-
7.3 Solids ally identical with HF results, as in the case of atoms.
On the other side, DFT calculations within the LDA
Exact exchange methods for solids have been developed are known to exhibit a significant overbinding, of the
by a number of authors. Kotani(105 – 108) implemented
order of 1–2 eV/atom (again 20–30%). Since EXX
a linear-muffin-tin-orbitals method within the atomic
contains the exact exchange, it is sensible to analyse
sphere approximation (LMTO-ASA), while Bylander and
the effect of combining EXX with usual LDA and
Kleinman(115) and Städele et al.(93,94) proposed, instead,
GGA correlation functionals. EXX-LDA cohesive
a pseudopotential plane-wave implementation. For the
energies are very much improved with respect to pure
correlation, they used the usual LDA and GGA functionals.
EXX, while the EXX-GGA ones are remarkably close
The accuracy of the KLI and another approximation –
to experimental values.(93,94) Similar results have
the average Fock approximation (AFA)(143,144) – in solids
been reported for HF-GGA calculations.(145) Figure 5
was investigated in Reference 115. A very thorough
shows the performance of different approaches for the
study of many properties of eight semiconductors in
cohesive energy of eight different semiconductors, as
several approximations was presented by Städele et al.(93,94)
reported by Städele et al.(93,94)
In the observations below, we mainly follow these
3. Lattice constants: Owing to the well-known
works.
overbinding effect, lattice constants in DFT-LDA are
1. Exchange energies: Exact-exchange energies have usually underestimated by 1 to 3%. GGAs normally
been compared against different approximations. over-correct this effect. EXX-LDA lattice constants
LDA values are small, roughly by 5%, while differ- are also in rather good agreement with experiment,
ent GGAs produce results of very high quality (lower somewhat better than the LDA. An important issue
by 1%) compared to EXX. The best agreement is in the case of pseudopotential calculations is how
obtained when the exchange energy is calculated at these pseudopotentials have been generated, that is,
the EXX converged density, although the difference is which approximation of exchange and correlation
 Density functional theory: Basics, new trends and applications 33

correlation are treated at different levels of accuracy.

C
Exchange is much better described than in the LDA,

SiC
8 but correlation does not compensate for this improve-
ment, and hence the discrepancies. Therefore, a better

AIN
treatment of correlation is needed in order to obtain
GaN

more accurate bulk moduli.

GaAs
Calculated energy [eV]

6 5. Energy band gaps: While the above properties probe

Si

the global energetics, the quality of the EXX poten-

Ge
AlAs

tial is probed by single-particle properties such as the

OEP eigenvalues. Band structures have been carefully
4
studied by several authors, and shown to improve
systematically in the EXX approach. The underes-
timation of the band gap – the smallest energy dif-
LDA
ference between valence and conduction bands in
2 EXX(GGA)
semiconductors – has been a long-standing problem
EXX
for density functionals. It was shown to be solved by
HF
performing the full many-body calculation, for exam-
0 ple, in the GW approximation.(53,147) This has been a
3 4 5 6 7 8 benchmark problem for approximate functionals, but
Experimental energy [eV] none of the usual GGA was able to improve signif-
icantly on it. The reason for this is again the lack
Figure 5. Cohesive energies for several semiconductors within
the LDA, pure EXX, EXX combined with PW86-GGA correlation
of self-interaction cancellation in LDA and GGA. In
and HF [Reproduced by permission of APS Journals from the EXX, occupied orbitals are self-interaction-free,
M. Städele, M. Moukara, J.A. Majevski, P. Vogl, and A. Görling and hence more localized and lower in energy than
(1999) Phys. Rev. B, 59, 10 031.](94) the LDA ones. At the same time, conduction states
are not affected by self-interaction because they are
has been used for the core electrons. EXX-LDA empty. Therefore, the gap is increased in the EXX,
calculations use EXX-LDA pseudopotentials. These and agreement with experiment is astonishingly good,
are less attractive than pure LDA ones, due to the as it is shown in Figure 6, extracted from the work of
better screening provided by the EXX. Therefore, Städele et al.(93,94) Hartree–Fock gaps are known to
when used in a pure LDA calculation, lattice constants be extremely large compared to experiment, the rea-
increase, the discrepancy being more important for son being that empty states ‘see’ a different potential
increasing number of core electrons. In a consistent than occupied states, which is not self-interaction-
EXX-LDA functional, which describes both core and free. The EXX potential is state-independent, treating
valence electrons at the same level, the previous occupied and empty states consistently.
effect is counteracted because the valence charge 6. XC discontinuity: The energy gap defined as the
density shrinks with respect to the LDA one. A orbital eigenvalue difference εgap is actually different
remarkably good agreement is obtained when LDA from the true band gap, which is defined as the energy
calculations are supplemented with pseudopotentials difference between states with N and N ± 1 electrons,
corrected for the non-additivity of the core and that is, Egap = E(N + 1) + E(N − 1) − 2E(N ):
valence charge (non-linear core corrections.)(146) This
effect is supposed to be smaller in EXX calculations, Egap = εgap + XC = εEXX
gap + εgap + XC
c
(112)
although it has not been assessed directly.
4. Bulk moduli: These are overestimated in the EXX- where εEXX
gap is the eigenvalue gap in the exchange-
LDA approach by approximately 20% with respect only EXX calculation, εcgap is a contribution to the
to experiment, while LDA values are in much bet- gap arising from correlation and XC is an energy
ter agreement. EXX pseudopotential LDA calcula- difference originated in the discontinuous jump of
tions underestimate them, so that the disagreement the exchange-correlation potential at integer num-
arises from the correlation of the valence electrons. In bers of electrons (the integral preference principle).
the pure LDA approach, the agreement is fortuitous, This quantity XC is usually called the discontinuity.
and mainly based on error cancellation. The prob- The magnitude of this discontinuity was controversial
lem of the EXX-LDA approach is that exchange and up to now, but EXX calculations(93,94) have shown
34 Electronic structure of large molecules

6 reproduced in comparison to experimental reflectiv-

ity data. The positions of the salient features, that is,

C
AIN
the absorption edge and peaks, and also the inten-
LDA sity of the features are very precise, except for some
Calculated band gap [eV]

EXX neglected effects like excitonic binding (a many-body

GaN
4 correlation effect) and spin-orbit coupling.
9. Exchange potential The quality of the exchange
potential in the LDA and GGA has been analysed
SiC
AIAs

by comparing them with the EXX (exact) potential,

GaAs

in a few semiconductors. The first observation is

2
that, in the EXX, the exchange potential is not a
Si

simple function of the density, so that there is a

Ge

range of values of VX corresponding to the same

density. The LDA value mimics the average density
0 dependence of EXX at low electronic densities, but
0 2 4 6 at higher densities it departs, becoming much less
Experimental band gap [eV] attractive than it should be, owing to the residual
self-interaction. Looking specifically at the VX (r)
Figure 6. Energy band gaps for several semiconductors within the
LDA and EXX [Reproduced by permission of APS Journals from profile along the covalent bonds, it was observed
M. Städele, M. Moukara, J.A. Majevski, P. Vogl, and A. Görling that the LDA is not attractive enough in the bonding
(1999) Phys. Rev. B, 59, 10 031.](94) region, again due to the self-repulsion. In the core
region, however, the LDA does reasonably well.
that Egap and εEXXgap are actually very close, so that GGAs, which attempt at correcting the problem of
XC ≈ −εcgap . Typical values for semiconductors are self-interaction, effectively do much better than LDA
of the order of −0.1 eV in the LDA and 0.2 eV in in the bonding region, which is the main aspect
the GGAs. Discriminating the exchange and corre- responsible for chemical properties. However, they
lation parts of the discontinuity was also possible, exhibit unphysical peaks in the region approaching
and indicated values of the order of 5 eV for x . the nuclear sites, because of the artificial divergence
This implies values of c of the order of −5 eV, and of the approximation (see above). These features can
a massive cancellation between exchange and corre-
lation discontinuities. The discontinuity, that is, the
[111] direction
difference between real quasiparticle gaps (Egap ) and 20
EXX eigenvalue gaps (εEXX gap ), is actually the excitonic Si
binding energy, and is very much dependent on the
particular system and also on external conditions like GGA
10
pressure.(148)
Vx [eV]

7. Band structures: The full k-dependence of the band

structure across the Brillouin zone is very well
LDA
reproduced by the EXX approach, in particular, the 0
ordering of the conduction-band minima, which is a
well-known LDA problem in some semiconductors EXX(X)
like Ge (negative LDA direct gap at , where an
−10
indirect gap at L is experimentally observed). Part of
this tendency, especially the direct gaps, is already Si Si
corrected at the level of pseudopotentials, when the
core electrons are treated within the EXX. However, Ga As
this is not the case at every k-point. Bandwidths, due
to the absence of self-interaction, are decreased in the
EXX relative to the LDA values. Figure 7. Exchange potential for bulk Si along the Si–Si
bond, for the LDA, GGA and EXX functionals [Reproduced
8. Optical properties: As a consequence, optical prop- by permission of APS Journals from M. Städele, M. Moukara,
erties like the dielectric function, which depend on J.A. Majevski, P. Vogl, and A. Görling (1999) Phys. Rev. B, 59,
the details of the band structure, are remarkably well 10 031.](94)
 Density functional theory: Basics, new trends and applications 35

be observed in Figure 7, reproduced from the work 15. von Weiszäcker, C.F. (1935) Z. Phys., 96, 431.
of Städele et al.(93,94) 16. Perrot, F. (1994) J. Phys.: Condens. Matter, 6, 431.
10. KLI approximation: The KLI approximation to 17. Wang, L.-W. and Teter, M.P. (1992) Phys. Rev. B, 45,
EXX proved very good for atomic systems and was 13 397.
assumed to be also good for solids, but without 18. Smargiassi, E. and Madden, P.A. (1994) Phys. Rev. B, 49,
proof.(115) Städele et al.(93,94) have actually shown 5220.
that total energies are a few tenth of an eV/atom 19. Foley, M. and Madden, P.A. (1996) Phys. Rev. B, 53,
higher than EXX ones, while energy gaps are under- 10 589.
estimated by about 0.5 eV. These conclusions do not 20. Hohenberg, P. and Kohn, W. (1964) Phys. Rev., 136, B864.
depend much on the pseudopotentials used (KLI or
21. Levy, M. (1982) Phys. Rev. A, 26, 1200.
EXX), but mainly on the description of exchange for
22. Theophilou, A. (1979) J. Phys. C, 12, 5419.
the valence electrons. Probably the reason lies on the
averaging of the denominator in the Green’s function 23. Hadjisavvas, N. and Theophilou, A. (1985) Phys. Rev. A, 32,
720.
in the KLI approximation, which could be too crude
in solids because of the k-vector dependence of the 24. Kohn, W. (1986) Phys. Rev. A, 34, 737.
eigenvalues. 25. Gidopoulos, N.I., Papaconstantinou, P.G., and Gross, E.K.U.
(2002) Phys. Rev. Lett., 88, 033003.
26. Gross, E.K.U., Dobson, J.F., and Petersilka, M. (1996) Den-
NOTES sity functional theory, in Springer Series Topics in Current
Chemistry, ed. R.F. Nalewajski, Springer, Berlin.
[1] The mentor of modern density functional theory, Prof. 27. Runge, E. and Gross, E.K.U. (1984) Phys. Rev. Lett., 52,
Walter Kohn, has been awarded the 1998 Nobel prize 997.
for chemistry together with Prof. John Pople, who 28. Petersilka, M., Gossmann, U.J., and Gross, E.K.U. (1996)
popularized quantum chemical calculations by means Phys. Rev. Lett., 76, 1212.
of the computational package GAUSSIAN. 29. Görling, A. (1996) Phys. Rev. A, 54, 3912.
[2] J. Kohanoff (unpublished).
30. Görling, A. (2000) Phys. Rev. Lett., 85, 4229.
31. Kohn, W. and Sham, L.J. (1965) Phys. Rev., 140, A1133.
REFERENCES 32. Pederson, M.R. and Jackson, K.A. (1991) Phys. Rev. B, 43,
7312.
1. Born, M. and Oppenheimer, J.R. (1927) Ann. der Phys., 84, 33. Valiev, M.M. and Fernando, G.W. (1995) Phys. Rev. B, 52,
457. 10 697.
2. Messiah, A. (1961) Quantum Mechanics, Vol. 1, North- 34. Janak, J.F. (1978) Phys. Rev. B, 18, 7165.
Holland, Amsterdam. 35. Langreth, D.C. and Perdew, J.P. (1977) Phys. Rev. B, 15,
3. Hartree, D.R. (1928) Proc. Cambridge. Philos. Soc., 24, 89. 2884.
4. Fock, V. (1930) Z. Phys., 61, 126. 36. Perdew, J.P. and Wang, Y. (1992) Phys. Rev. B, 46, 12 947.
5. Slater, J.C. (1930) Phys. Rev., 35, 210. 37. Parr, R.G. and Yang, W. (1989) Density Functional Theory
6. Møller, C. and Plesset, M.S. (1934) Phys. Rev., 46, 618. of Atoms and Molecules, Oxford University Press, New
York.
7. Szabo, A. and Ostlund, N.S. (1989) Modern Quantum
Chemistry: Introduction to Advanced Electronic Structure 38. von Barth, U. and Hedin, L. (1979) J. Phys. C, 12, 5419.
Theory, McGraw-Hill, New York. 39. Vosko, S.H., Wilk, L., and Nussair, M. (1980) Can. J. Phys.,
8. Jensen, F. (1999) Introduction to Computational Chemistry, 58, 1200.
John Wiley & Sons, Chichester. 40. Jones, R.O. and Gunnarsson, O. (1989) Rev. Mod. Phys., 61,
9. Thomas, L.H. (1927) Proc. Cambridge. Philos. Soc., 23, 689.
542. 41. Mahan, G.D. (1990) Many Particle Physics, Plenum Pub-
10. Fermi, E. (1928) Z. Phys., 48, 73. lishing, New York.
11. Ceperley, D.M. (1978) Phys. Rev. B, 18, 3126. 42. Perdew, J.P. and Zunger, A. (1981) Phys. Rev. B, 23, 5048.
12. Ceperley, D.M. and Alder, B.J. (1980) Phys. Rev. Lett., 45, 43. Gell-Mann, M. and Brückner, K.A. (1957) Phys. Rev., 106,
566. 364.
13. Zaremba E. and Kohn, W. (1976) Phys. Rev. B, 13, 2270. 44. Alonso, J.A. and Girifalco, L.A. (1977) Solid State Com-
14. March, N.H. (1983) in Theory of the Inhomogeneous Elec- mun., 24, 135.
tron Gas, eds S. Lundqvist and N.H. March, Plenum Pub- 45. Alonso, J.A. and Girifalco, L.A. (1978) Phys. Rev. B, 17,
lishing, New York. 3735.
36 Electronic structure of large molecules

46. Gunnarsson, O. and Jones, R.O. (1980) Phys. Scripta, 21, 79. Goerling, A. (1999) Phys. Rev. Lett., 83, 5459.
394. 80. Kohn, W. and Mattson, A.E. (1998) Phys. Rev. Lett., 81,
47. Gunnarsson, O. and Jones, R.O. (1980) J. Chem. Phys., 72, 3487.
5357.
81. Vitos, L., Johansson, B., Kollár, J., and Skriver, H.K. (2000)
48. Gunnarsson, O., Jonson, M., and Lundqvist, B.I. (1979) Phys. Rev. B, 62, 10 046.
Phys. Rev. B, 20, 3136.
82. Vitos, L., Johansson, B., Kollár, J., and Skriver, H.K. (2000)
49. Denton, A.R. and Ashcroft, N.W. (1989) Phys. Rev. A, 39, Phys. Rev. A, 61, 052511.
4701.
83. Svendsen, P.S. and von Barth, U. (1996) Phys. Rev. B, 54,
50. Denton, A.R., Nielaba, P., Runge, K.J., and Ashcroft, N.W. 17 402.
(1990) Phys. Rev. Lett., 64, 1529.
84. Perdew, J.P., Kurth, S., Zupan, A., and Blaha, P. (1999)
51. Lutsko, J.F. and Baus, M. (1990) Phys. Rev. Lett., 64, 761. Phys. Rev. Lett., 82, 2544.
52. Singh, D.J. (1993) Phys. Rev. B, 48, 14 099. 85. Levy, M. and Perdew, J.P. (1985) Phys. Rev. A, 32, 2010.
53. Hedin, L. (1965) Phys. Rev., 139, A796. 86. Seidl, M., Perdew, J.P., and Levy, M. (1999) Phys. Rev. A,
54. Anisimov, V.I., Zaanen, J., and Andersen, O.K. (1991) 59, 51.
Phys. Rev. B, 44, 943. 87. Becke, A.D. (1998) J. Chem. Phys., 109, 2092.
55. Ma, S.-K. and Brueckner, K.A. (1968) Phys. Rev., 165, 18. 88. van Voorhis, T. and Scuseria, G.E. (1998) J. Chem. Phys.,
56. Fetter, A.L. and Walecka, J.D. (1971) Quantum Theory of 109, 400.
Many-Particle Systems, McGraw-Hill, New York.
89. Adamo, C., Ernzerhof, M., and Scuseria, G.E. (2000) J.
57. Langreth, D.C. and Mehl, M.J. (1981) Phys. Rev. Lett., 47, Chem. Phys., 112, 2643.
446.
90. Sharp, R.T. and Horton, G.K. (1953) Phys. Rev., 90, 317.
58. Langreth, D.C. and Mehl, M.J. (1983) Phys. Rev. B, 28,
91. Sahni, V., Gruenebaum, J., and Perdew, J.P. (1982) Phys.
1809.
Rev. B, 26, 4371.
59. Gross, E.K.U. and Dreizler, R.M. (1981) Z. Phys. A, 302,
103. 92. Perdew, J.P. and Norman, M.R. (1982) Phys. Rev. B, 26,
5445.
60. Perdew, J.P. (1985) Phys. Rev. Lett., 55, 1665.
93. Städele, M., Majevski, J.A., Vogl, P., and Görling, A.
61. Ghosh, S.K. and Parr, R.G. (1986) Phys. Rev. A, 34, 785. (1997) Phys. Rev. Lett., 79, 2089.
62. Filippi, C., Umrigar, C.J., and Taut, M. (1994) J. Chem. 94. Städele, M., Moukara, M., Majevski, J.A., Vogl, P., and
Phys., 100, 1295. Görling, A. (1999) Phys. Rev. B, 59, 10 031.
63. Perdew, J.P. and Wang, Y. (1986) Phys. Rev. B, 33, 8800. 95. Görling, A. (1999) Phys. Rev. Lett., 83, 5459.
64. Perdew, J.P. (1986) Phys. Rev. B, 33, 8822. 96. Talman, J.D. and Shadwick, W.F. (1976) Phys. Rev. A, 14,
65. Perdew, J.P. and Wang, Y. (1991) Phys. Rev. B, 45, 13 244. 36.
66. Becke, A.D. (1988) Phys. Rev. A, 38, 3098. 97. Görling, A. and Levy, M. (1994) Phys. Rev. A, 50, 196.
67. Lee, C., Yang, W., and Parr, R.G. (1988) Phys. Rev. B, 37, 98. Görling, A. (1996) Phys. Rev. B, 53, 7024.
785.
99. Colle, R. and Salvetti, D. (1975) Theor. Chim. Acta, 37, 329.
68. Perdew, J.P., Burke, K., and Ernzerhof, M. (1996) Phys.
Rev. Lett., 77, 3865. 100. Colle, R. and Salvetti, D. (1979) Theor. Chim. Acta, 53, 55.

69. Perdew, J.P., Burke, K., and Ernzerhof, M. (1997) Phys. 101. Li, Y., Krieger, J.B., and Iafrate, G.J. (1992) Chem. Phys.
Rev. Lett., 78, 1396 (E). Lett., 191, 38.
70. Lieb, E.H. and Oxford, S. (1981) Int. J. Quantum Chem., 102. Krieger, J.B., Li, Y., and Iafrate, G.J. (1992) Phys. Rev. A,
19, 427. 45, 101.
71. Becke, A.D. (1986) J. Chem. Phys., 84, 4524. 103. Li, Y., Krieger, J.B., and Iafrate, G.J. (1993) Phys. Rev. A,
47, 165.
72. Zhang, Y. and Yang, W. (1998) Phys. Rev. Lett., 80, 890.
104. Engel, E. and Vosko, S.H. (1993) Phys. Rev. A, 47, 2800.
73. Perdew, J.P., Burke, K., and Ernzerhof, M. (1998) Phys.
Rev. Lett., 80, 891. 105. Kotani, T. (1994) Phys. Rev. B, 50, 14 816.
74. Ernzerhof, M. and Scuseria, G.E. (1999) J. Chem. Phys., 106. Kotani, T. (1995) Phys. Rev. B, 51, 13 903 (E).
110, 5029. 107. Kotani, T. (1995) Phys. Rev. Lett., 74, 2989.
75. Levy, M. and Perdew, J.P. (1994) Int. J. Quantum Chem., 108. Kotani, T. and Akai, H. (1995) Phys. Rev. B, 52, 17 153.
49, 539.
109. Krieger, J.B., Li, Y., and Iafrate, G.J. (1992) Phys. Rev. A,
76. Becke, A.D. (1993) J. Chem. Phys., 98, 5648. 46, 5453.
77. Becke, A.D. (1996) J. Chem. Phys., 104, 1040. 110. Grabo, T., Kreibich, T., Kurth, S., and Gross, E.K.U. (1998)
78. Becke, A.D. (1997) J. Chem. Phys., 107, 8554. in Strong Coulomb Correlations in Electronic Structure:
 Density functional theory: Basics, new trends and applications 37

Beyond the Local Density Approximation, ed. V.I. Anisimov, 130. Rapcewicz, K. and Ashcroft, N.W. (1991) Phys. Rev. B, 44,
Gordon & Breach, Tokyo. 4032.
111. Perdew, J.P., Parr, R.G., Levy, M., and Balduz, J.L. (1982) 131. Andersson, Y., Langreth, D.C., and Lundqvist, B.I. (1996)
Phys. Rev. Lett., 49, 1691. Phys. Rev. Lett., 76, 102.
112. Görling, A. and Ernzerhof, M. (1995) Phys. Rev. A, 51, 132. Dobson, J.F. and Dinte, B.P. (1996) Phys. Rev. Lett., 76,
4501. 1780.
113. Ivanov, S., Hirata, S., and Bartlett, R.J. (1999) Phys. Rev. 133. Lein, M., Dobson, J.F., and Gross, E.K.U. (1999) J. Com-
Lett., 83, 5455. put. Chem., 20, 12.
114. van Gisbergen, S.J.A., Schipper, P.R.T., Grischenko, O.V., 134. Aziz, R.A. and Slaman, M.J. (1991) J. Chem. Phys., 94,
Baerends, E.J., Snijders, J.G., Champagne, B., and Kirt- 8047.
man, B. (1999) Phys. Rev. Lett., 83, 694.
135. Seidl, M., Perdew, J.P., and Kurth, S. (2000) Phys. Rev.
115. Bylander, D.M. and Kleinman, L. (1995) Phys. Rev. Lett., Lett., 84, 5070.
74, 3660.
136. Seidl, M., Perdew, J.P., and Kurth, S. (2000) Phys. Rev. A,
116. Tao, J., Gori-Giorgi, P., Perdew, J.P., and McWeeny, R. 62, 012502.
(2001) Phys. Rev. A, 63, 032513.
137. Casida, M.E. (1995) Phys. Rev. A, 51, 2005.
117. Langreth, D.C. and Perdew, J.P. (1977) Phys. Rev. B, 15,
138. Sham, L.J. and Schlüter, M. (1983) Phys. Rev. Lett., 51,
2884.
1888.
118. Singwi, K.S., Tosi, M.P., Land, R.H., and Sjölander, A.
(1963) Phys. Rev., 176, 589. 139. Umrigar, C. and Gonze, X. (1994) Phys. Rev. A, 50, 3827.

119. Kurth, S. and Perdew, J.P. (1999) Phys. Rev. B, 59, 10 461. 140. Perdew, J.P., Chevary, J.A., Vosko, S.H., Jackson, K.A.,
Pederson, M.R., Singh, D.J., and Fiolhais, C. (1992) Phys.
120. Kurth, S. and Perdew, J.P. (1999) Phys. Rev. B, 60, 11 212 Rev. B, 46, 6671.
(E).
141. Radzig, A.A. and Smirnov, M.M. (1985) Reference Data on
121. Yan, Z., Perdew, J.P., and Kurth, S. (2000) Phys. Rev. B, 61, Atoms and Molecules, Springer, Berlin.
16 430.
142. Kohanoff, J. and Scandolo, S. (1998) Mater. Res. Soc. Symp.
122. Kohn, W., Meir, Y., and Makarov, D.E. (1998) Phys. Rev. Proc., 499, 329.
Lett., 80, 4153.
143. Kleinman, L. (1994) Phys. Rev. B, 49, 14 197.
123. Leininger, T., Stoll, H., Werner, H.-J. and Savin, A. (1997)
Chem. Phys. Lett., 275, 151. 144. Gu, Y.M., Bylander, D.M., and Kleinman, L. (1994) Phys.
Rev. B, 50, 2227.
124. Görling, A. and Levy, M. (1993) Phys. Rev. B, 47, 13 105.
145. Causa, M., Dovesi, R., and Roetti, C. (1991) Phys. Rev. B,
125. Görling, A. and Levy, M. (1995) Phys. Rev. A, 52, 4493. 43, 11 937.
126. Levy, M. and Perdew, J.P. (1985) Phys. Rev. A, 32, 2010. 146. Louie, S.B., Froyen, S., and Cohen, M.L. (1982) Phys. Rev.
127. Görling, A. and Levy, M. (1994) Phys. Rev. A, 50, 196. B, 26, 1738.
128. Ernzerhof, M. (1996) Chem. Phys. Lett., 263, 499. 147. Gygi, F. and Baldereschi, A. (1986) Phys. Rev. B, 34, 4405.
129. Engel, E., Höck, A., and Dreizler, R.M. (2000) Phys. Rev. 148. Städele, M. and Martin, R.M. (2000) Phys. Rev. Lett., 84,
A, 61, 032502. 6070.

View publication stats