Review of Approximations For The Exchange-Correlation Energy in Density-Functional Theory

Chapter in the book Density Functional Theory edited by Éric Cancès and Gero Friesecke.
Review of approximations for the exchange-correlation energy in

density-functional theory
Julien Toulouse
Laboratoire de Chimie Théorique (LCT), Sorbonne Université and CNRS, F-75005 Paris, France
arXiv:2103.02645v2 [physics.chem-ph] 9 Sep 2022
Institut Universitaire de France, F-75005 Paris, France
31 August, 2022
In this chapter, we provide a review of ground-state Kohn–Sham density-functional theory

of electronic systems and some of its extensions, we present exact expressions and constraints
for the exchange and correlation density functionals, and we discuss the main families of ap-
proximations for the exchange-correlation energy: semilocal approximations, single-determinant
hybrid approximations, multideterminant hybrid approximations, dispersion-corrected approx-
imations, as well as orbital-dependent exchange-correlation density functionals. The chapter
aims at providing both a consistent bird’s-eye view of the field and a detailed description of
some of the most used approximations. It is intended to be readable by chemists/physicists and
applied mathematicians. For more coverage of the subject, the reader may consult for example
the books of Refs. [1–5] and the review articles of Refs. [6–15].
Contents
1 Basics of density-functional theory 2
1.1 The many-body problem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2
1.2 The universal density functional . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
1.3 The Kohn–Sham scheme . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
1.3.1 Decomposition of the universal functional . . . . . . . . . . . . . . . . . . 4
1.3.2 The Kohn–Sham equations . . . . . . . . . . . . . . . . . . . . . . . . . . 6
1.3.3 Extension to spin density-functional theory . . . . . . . . . . . . . . . . . 7
1.4 The generalized Kohn–Sham scheme . . . . . . . . . . . . . . . . . . . . . . . . . 9
2 Exact expressions and constraints for the Kohn–Sham exchange and correlation func-
tionals 9
2.1 The exchange and correlation holes . . . . . . . . . . . . . . . . . . . . . . . . . . 9
2.2 The adiabatic connection . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
2.3 One-orbital and one-electron spatial regions . . . . . . . . . . . . . . . . . . . . . 13
2.4 Coordinate scaling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
2.4.1 Uniform coordinate scaling . . . . . . . . . . . . . . . . . . . . . . . . . . 15
2.4.2 Non-uniform coordinate scaling . . . . . . . . . . . . . . . . . . . . . . . . 17
2.5 Atoms in the limit of large nuclear charge . . . . . . . . . . . . . . . . . . . . . . 17
2.6 Lieb–Oxford lower bound . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
1
3 Semilocal approximations for the exchange-correlation energy 19
3.1 The local-density approximation . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
3.2 The gradient-expansion approximation . . . . . . . . . . . . . . . . . . . . . . . . 22
3.3 Generalized-gradient approximations . . . . . . . . . . . . . . . . . . . . . . . . . 23
3.4 Meta-generalized-gradient approximations . . . . . . . . . . . . . . . . . . . . . . 28
4 Single-determinant hybrid approximations 34

4.1 Hybrid approximations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34
4.2 Range-separated hybrid approximations . . . . . . . . . . . . . . . . . . . . . . . 36
5 Multideterminant hybrid approximations 40

5.1 Double-hybrid approximations . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40
5.2 Range-separated double-hybrid approximations . . . . . . . . . . . . . . . . . . . 43
5.2.1 Range-separated one-parameter double-hybrid approximations . . . . . . 43
5.2.2 Range-separated two-parameter double-hybrid approximations . . . . . . 47
6 Semiempirical dispersion corrections and nonlocal van der Waals density functionals 49
6.1 Semiempirical dispersion corrections . . . . . . . . . . . . . . . . . . . . . . . . . 49
6.2 Nonlocal van der Waals density functionals . . . . . . . . . . . . . . . . . . . . . 50
7 Orbital-dependent exchange-correlation density functionals 51

7.1 Exact exchange . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51
7.2 Second-order Görling–Levy perturbation theory . . . . . . . . . . . . . . . . . . . 52
7.3 Adiabatic-connection fluctuation-dissipation approach . . . . . . . . . . . . . . . 54
7.3.1 Exact adiabatic-connection fluctuation-dissipation expression . . . . . . . 54
7.3.2 Random-phase approximations . . . . . . . . . . . . . . . . . . . . . . . . 55
References 56
1 Basics of density-functional theory

1.1 The many-body problem
We consider an N -electron system (atom or molecule) in the Born–Oppenheimer and non-
relativistic approximations. The electronic Hamiltonian in the position representation is, in
atomic units,
N N N N
1X 2 1 XX 1 X
H=− ∇r i + + vne (ri ), (1.1)
2 2 |ri − rj |
i=1 i=1 j=1 i=1
i6=j
where ∇2ri = ∆ri is the Laplacian with respect to the electron coordinate ri and vne (ri ) =
− N
P n
α=1 Zα /|ri − Rα | is the nuclei-electron interaction depending on the positions {Rα } and
charges {Zα } of the Nn nuclei. The stationary electronic states are determined by the time-
independent Schrödinger equation,
HΨ(x1 , x2 , ..., xN ) = EΨ(x1 , x2 , ..., xN ), (1.2)
2
where Ψ(x1 , x2 , ..., xN ) is a wave function written with space-spin coordinates xi = (ri , σi ) ∈
R3 × {↑, ↓} (with {↑, ↓} ∼ = Z2 being the set of spin coordinates) which is antisymmetric with
respect to the exchange of two coordinates, and E is the associated energy.
Using Dirac notation, the Schrödinger equation (1.2) can be rewritten in a representation-
independent formalism,
Ĥ|Ψi = E|Ψi, (1.3)
where the Hamiltonian is formally written as
Ĥ = T̂ + Ŵee + V̂ne , (1.4)

with the kinetic-energy operator T̂ , the electron-electron interaction operator Ŵee , and the
nuclei-electron interaction operator V̂ne .
The quantity of primary interest is the ground-state energy E0 . The variational theorem
establishes that E0 can be expressed as an infimum,
E0 = inf hΨ|Ĥ|Ψi, (1.5)
Ψ∈W N
where the search is over the set of N -electron antisymmetric normalized wave functions Ψ having
a finite kinetic energy,
n N
^ o
N
W = Ψ∈ L2 (R3 × {↑, ↓}; C), Ψ ∈ H 1 ((R3 × {↑, ↓})N ; C), hΨ|Ψi = 1 , (1.6)
where N is the N -fold antisymmetrized tensor product, L2 and H 1 are the standard Lebesgue
V
and Sobolev spaces (i.e., respectively, the space of functions that are square integrable and the
space of functions that are square integrable together with their first-order derivatives), and h·|·i
designates the L2 inner product. Density-functional theory (DFT) is based on a reformulation
of the variational theorem in terms of the one-electron density defined as1
Z
ρΨ (r) = N |Ψ(x, x2 , ..., xN )|2 dσdx2 ...dxN , (1.7)
{↑,↓}×(R3 ×{↑,↓})N−1
R
which is normalized to the electron number, R3 ρΨ (r)dr = N .
1.2 The universal density functional

Building on the work of Hohenberg and Kohn [16], Levy [17] and Lieb [18] proposed to define
the following universal density functional F [ρ] using a constrained-search approach,
F [ρ] = min hΨ|T̂ + Ŵee |Ψi = hΨ[ρ]|T̂ + Ŵee |Ψ[ρ]i, (1.8)
Ψ∈WρN
where the minimization is done over the set of N -electron wave functions Ψ yielding the fixed
density ρ [via Eq. (1.7)],
WρN = {Ψ ∈ W N , ρΨ = ρ}. (1.9)
In Eq. (1.8), for a given density ρ, Ψ[ρ] denotes a minimizing wave function, which is known to
exist [18] but is possibly not unique. This so-called Levy–Lieb functional F [ρ] is defined on the
set of N -representable densities [18]:
DN = {ρ | ∃Ψ ∈ W N s.t. ρΨ = ρ}
√
Z
1 3
= {ρ ∈ L (R ) | ρ ≥ 0, ρ(r)dr = N, ρ ∈ H 1 (R3 )}. (1.10)
R3
1
In Eq. (1.7), an integration over a spin coordinate σ just means a sum over the two values σ ∈ {↑, ↓}.
3
We note that an alternative universal density functional can be defined by a Legendre–
Fenchel transformation, or equivalently by a constrained-search over N -electron ensemble density
matrices [18]. This so-called Lieb functional has the advantage of being convex but in this chapter
we will simply use the Levy–Lieb functional of Eq. (1.8).
The exact ground-state energy can then be expressed as
Z
E0 = inf F [ρ] + vne (r)ρ(r)dr , (1.11)
ρ∈D N R3
and if a minimizer exists then it is a ground-state density ρ0 (r) for the potential vne (r). Hence,
the ground-state energy can be in principle obtained by minimizing over the density ρ, i.e.
a simple function of 3 real variables, which is a tremendous simplification compared to the
minimization over a complicated many-body wave function Ψ. However, the explicit expression
of F [ρ] in terms of the density is not known, and the direct approximations for F [ρ] that have
been tried so far turn out not to be accurate enough.
If there is a unique wave function Ψ[ρ] (up to a phase factor) in Eq. (1.8), we can define
kinetic and potential contributions to F [ρ],
F [ρ] = T [ρ] + Wee [ρ], (1.12)
where T [ρ] = hΨ[ρ]|T̂ |Ψ[ρ]i and Wee [ρ] = hΨ[ρ]|Ŵee |Ψ[ρ]i. The kinetic-energy functional T [ρ]
is the contribution which is particularly difficult to approximate as an explicit functional of the
density.
1.3 The Kohn–Sham scheme

1.3.1 Decomposition of the universal functional
Following the idea of Kohn and Sham (KS) [19], the difficulty of approximating F [ρ] directly
can be circumvented by decomposing F [ρ] as
F [ρ] = Ts [ρ] + EHxc [ρ], (1.13)
where Ts [ρ] is the non-interacting kinetic-energy functional which can be defined with a con-
strained search2 ,
Ts [ρ] = min hΦ|T̂ |Φi = hΦ[ρ]|T̂ |Φ[ρ]i, (1.14)
Φ∈SρN
where the minimization is over the set of N -electron single-determinant wave functions Φ yielding
the fixed density ρ:
SρN = {Φ ∈ S N , ρΦ = ρ}. (1.15)
Here, S N is the set of N -electron single-determinant wave functions built from orthonormal spin
orbitals
S N = {Φ = φ1 ∧ φ2 ∧ ... ∧ φN | ∀i φi ∈ H 1 (R3 × {↑, ↓}; C), ∀i, j hφi |φj i = δi,j }, (1.16)
where φ1 ∧φ2 ∧...∧φN designates the normalized N -fold antisymmetrized tensor product of N spin
orbitals. The functional Ts [ρ] is defined over the entire set of N -representable densities D N since
2
It is also possible to define the non-interacting kinetic-energy functional analogously to the Levy–Lieb func-
tional in Eq. (1.8) by minimizing over wave functions Ψ ∈ WρN , i.e. Ts,LL [ρ] = minΨ∈WρN hΨ|T̂ |Ψi [18]. In this
case, the corresponding minimizing KS wave function can generally be a linear combination of Slater determi-
nants. However, we often have Ts,LL [ρ] = Ts [ρ], in particular for densities ρ that come from a non-interacting
ground-state wave function which is not degenerate. In this chapter, we will usually assume this nondegeneracy
condition.
4
any N -representable density can be obtained from a single-determinant wave function [18,20,21].
In Eq. (1.14), for a given density ρ, Φ[ρ] denotes a minimizing single-determinant wave function
(again known to exist [18] but possibly not unique), also referred to as KS wave function. The
remaining functional EHxc [ρ] that Eq. (1.13) defines is called the Hartree-exchange-correlation
functional. The idea of the KS scheme is then to use the exact expression of Ts [ρ] by reformulating
the minimization over densities in Eq. (1.11) as a minimization over single-determinant wave
functions Φ, n o
E0 = inf hΦ|T̂ + V̂ne |Φi + EHxc [ρΦ ] , (1.17)
Φ∈S N
and if a minimum exists then any minimizing single-determinant wave function in Eq. (1.17)
gives a ground-state density ρ0 (r). Thus, the exact ground-state energy can in principle be ob-
tained by minimizing over single-determinant wave functions only. Even though a wave function
has been reintroduced compared to Eq. (1.11), it is only a single-determinant wave function
Φ and therefore it still represents a tremendous simplification over the usual variational theo-
rem involving a correlated (multideterminant) wave function Ψ. The advantage of Eq. (1.17)
over Eq. (1.11) is that a major part of the kinetic energy can be treated exactly with the
single-determinant wave function Φ, and only EHxc [ρ] needs to be approximated as an explicit
functional of the density.
In practice, EHxc [ρ] is decomposed as
EHxc [ρ] = EH [ρ] + Exc [ρ], (1.18)
where EH [ρ] is the Hartree energy functional,

1 ρ(r1 )ρ(r2 )
Z
EH [ρ] = dr1 dr2 , (1.19)
2 R3 ×R3 |r1 − r2 |
representing the classical electrostatic repulsion energy for the charge distribution ρ(r) and which
is calculated exactly, and Exc [ρ] is the exchange-correlation energy functional that remains to be
approximated. If there is a unique KS wave function Φ[ρ] (up to a phase factor), we can further
decompose Exc [ρ] as
Exc [ρ] = Ex [ρ] + Ec [ρ], (1.20)
where Ex [ρ] is the exchange energy functional,
Ex [ρ] = hΦ[ρ]|Ŵee |Φ[ρ]i − EH [ρ], (1.21)
and Ec [ρ] is the correlation energy functional,
Ec [ρ] = hΨ[ρ]|T̂ + Ŵee |Ψ[ρ]i − hΦ[ρ]|T̂ + Ŵee |Φ[ρ]i = Tc [ρ] + Uc [ρ], (1.22)
which contains a kinetic contribution Tc [ρ] = hΨ[ρ]|T̂ |Ψ[ρ]i − hΦ[ρ]|T̂ |Φ[ρ]i and a potential
contribution Uc [ρ] = hΨ[ρ]|Ŵee |Ψ[ρ]i − hΦ[ρ]|Ŵee |Φ[ρ]i. Using the fact that Φ[ρ] is a single-
determinant wave function, it can be shown that the exchange functional can be expressed as
1 X |γσ (r1 , r2 )|2
Z
Ex [ρ] = − dr1 dr2 , (1.23)
2 R3 ×R3 |r1 − r2 |
σ∈{↑,↓}
where γσ , for σ ∈ {↑, ↓}, is the spin-dependent one-particle KS density matrix,

Z
′
γσ (r, r ) = N Φ[ρ](r′ , σ, x2 , ..., xN )∗ Φ[ρ](r, σ, x2 , ..., xN )dx2 ...dxN , (1.24)
(R3 ×{↑,↓})N−1
which shows that Ex [ρ] ≤ 0. Moreover, from the variational definition of F [ρ], we see that
Ec [ρ] ≤ 0.
5
1.3.2 The Kohn–Sham equations
The single-determinant wave function Φ in Eq. (1.17) is constructed from a set of N or-
thonormal occupied spin-orbitals {φi }i=1,...,N . To enforce Sz spin symmetry, each spin-orbital
is factorized as φi (x) = ϕi (r)χσi (σ) where ϕi ∈ H 1 (R3 , C) is a spatial orbital and χσi is a spin
function from {↑, ↓} to {0, 1} such that ∀σi , σ ∈ {↑, ↓}, χσi (σ) = δσi ,σ (σi is the spin of the spin-
orbital i). Alternatively, when this is convenient, we will sometimes reindex the spatial orbitals,
{ϕi } −→ {ϕiσ }, including explicitly the spin σ in the index. Writing the total electronic energy
in Eq. (1.17) in terms of spin-orbitals and integrating over the spin variables, we obtain:
N
1X
Z Z
E[{ϕi }] = |∇ϕi (r)|2 dr + vne (r)ρ(r)dr + EHxc [ρ], (1.25)
2 R3 R3
i=1
where the density is expressed in terms of the orbitals as

N
|ϕi (r)|2 .
X
ρ(r) = (1.26)
i=1
The minimization over Φ can then be recast into a minimization of E[{ϕi }] with respect to the
spatial orbitals {ϕi } with the constraint of keeping the orbitals orthonormalized. The stationary
condition with respect to variations of ϕi (r) leads to the KS equations [19],

1 2
− ∇ + vne (r) + vHxc (r) ϕi (r) = εi ϕi (r), (1.27)
2
where εi is the Lagrange multiplier associated to the normalization condition of ϕi and vHxc (r)
is the Hartree-exchange-correlation potential defined as the functional derivative of EHxc [ρ] with
respect to ρ(r),
δEHxc [ρ]
vHxc (r) = , (1.28)
δρ(r)
which is itself a functional of the density. The orbitals satisfying Eq. (1.27) are called the KS
orbitals. They are the eigenfunctions of the KS one-electron Hamiltonian,
1
hs (r) = − ∇2 + vs (r), (1.29)
2
where
vs (r) = vne (r) + vHxc (r) (1.30)
is the KS potential, and εi are then the KS orbital energies. Note that Eq. (1.27) constitutes a
set of coupled self-consistent equations since the potential depends on all the occupied orbitals
{ϕi }i=1,...,N through the density [Eq. (1.26)]. The operator hs (r) defines the KS system which
is a system of N non-interacting electrons in an effective external potential vs (r) ensuring that
the density ρ(r) in Eq. (1.26) is the same as the exact ground-state density ρ0 (r) of the physical
system of N interacting electrons. The exact ground-state energy E0 is then obtained by injecting
the KS orbitals in Eq. (1.25). The other (unoccupied) eigenfunctions in Eq. (1.27) define virtual
KS orbitals {ϕa }a≥N +1 .
Note that to define the potential vHxc (r) in Eq. (1.28) a form of differentiability of the func-
tional EHxc [ρ], also referred to as v-representability of the density, has been assumed. Justifying
this is in fact subtle and has been debated [22–28] (see also the chapter by Kvaal in the vol-
ume). Here, we will simply assume that a form of differentiability of EHxc [ρ] holds on at least
6
a restricted set of densities that allows one to define the potential vHxc (r) up to an additive
constant. For a further restricted set of densities that should include ground-state densities of
electronic Hamiltonians of molecular systems [Eq. (1.1)], it is expected that the KS potential
vs (r) tends to a constant as |r| → ∞ and we choose this constant to be zero. Note also that the
assumption of the existence of the KS potential vs (r) in Eq. (1.29), which does not depend on
spin coordinates, implies that each spin-orbital must indeed have a definite Sz spin value.
Following the decomposition of EHxc [ρ] in Eq. (1.18), the potential vHxc (r) is written as
vHxc (r) = vH (r) + vxc (r), (1.31)
where vH (r) = δEH [ρ]/δρ(r) = R3 ρ(r′ )/|r − r′ |dr′ is the Hartree potential and vxc (r) =
R
δExc [ρ]/δρ(r) is the exchange-correlation potential. Likewise, following the decomposition of

Exc [ρ] in Eq. (1.20), and assuming that both Ex [ρ] and Ec [ρ] are differentiable with respect to
ρ, the potential vxc (r) can be further decomposed as
vxc (r) = vx (r) + vc (r), (1.32)
where vx (r) = δEx [ρ]/δρ(r) is the exchange potential and vc (r) = δEc [ρ]/δρ(r) is the correlation
potential. Thus, the KS equations are similar to the Hartree-Fock (HF) equations, with the
difference that they involve a local exchange potential vx (r) instead of the nonlocal HF exchange
potential, and an additional correlation potential. At least for ground-state densities of finite
molecular systems, the exchange potential has the long-range asymptotic behavior (see, e.g.,
Ref. [29]),
1
vx (r) ∼ − , (1.33)
|r|→∞ |r|
whereas the correlation potential decays faster [30].
1.3.3 Extension to spin density-functional theory
To deal with an external magnetic field, DFT has been extended from the total density to
spin-resolved densities [31, 32]. Without external magnetic fields, this spin density-functional
theory is in principle not necessary, even for open-shell systems. However, the dependence on
the spin densities allows one to construct approximate exchange-correlation functionals that are
more accurate, and is therefore almost always used in practice for open-shell systems.
The spin density ρσ,Ψ with σ ∈ {↑, ↓} associated to a wave function Ψ is defined as
Z
ρσ,Ψ (r) = N |Ψ(x, x2 , ..., xN )|2 dx2 ...dxN , (1.34)
(R3 ×{↑,↓})N−1
and integrates to the number of σ-spin electrons Nσ , i.e. R3 ρσ,Ψ (r)dr = Nσ . For ρ↑ ∈ D N↑ and
R
ρ↓ ∈ D N↓ , the universal density functional is now defined as [33],
F [ρ↑ , ρ↓ ] = min hΨ|T̂ + Ŵee |Ψi = hΨ[ρ↑ , ρ↓ ]|T̂ + Ŵee |Ψ[ρ↑ , ρ↓ ]i, (1.35)
Ψ∈WρN ,ρ
↑ ↓
where the search is over the set of normalized antisymmetric wave functions Ψ with N = N↑ +N↓
electrons and yielding the fixed spin densities ρ↑ and ρ↓
WρN↑ ,ρ↓ = {Ψ ∈ W N , ρ↑,Ψ = ρ↑ , ρ↓,Ψ = ρ↓ }. (1.36)
In Eq. (1.35), Ψ[ρ↑ , ρ↓ ] designates a minimizing wave function.
7
A spin-dependent KS scheme is obtained by decomposing F [ρ↑ , ρ↓ ] as
F [ρ↑ , ρ↓ ] = Ts [ρ↑ , ρ↓ ] + EH [ρ] + Exc [ρ↑ , ρ↓ ], (1.37)
where Ts [ρ↑ , ρ↓ ] is defined as
Ts [ρ↑ , ρ↓ ] = min hΦ|T̂ |Φi = hΦ[ρ↑ , ρ↓ ]|T̂ |Φ[ρ↑ , ρ↓ ]i, (1.38)

Φ∈SρN ,ρ
↑ ↓
with a constrained search over the set of single-determinant wave functions Φ yielding the fixed
spin densities ρ↑ and ρ↓
SρN↑ ,ρ↓ = {Φ ∈ S N , ρ↑,Φ = ρ↑ , ρ↓,Φ = ρ↓ }. (1.39)
Here, Φ[ρ↑ , ρ↓ ] denotes a minimizing KS single-determinant wave function, EH [ρ] is the Hartree
energy which is a functional of the total density ρ = ρ↑ + ρ↓ only [Eq. (1.19)], and Exc [ρ↑ , ρ↓ ]
is the spin-dependent exchange-correlation energy functional. The ground-state energy is then
obtained as n o
E0 = inf hΦ|T̂ + V̂ne |Φi + EH [ρΦ ] + Exc [ρ↑,Φ , ρ↓,Φ ] . (1.40)
Φ∈S N
Writing the spatial orbitals of the spin-unrestricted determinant as {ϕiσ }i=1,...,N (with the
index explicitly including the spin σ now for clarity), we arrive at the spin-dependent KS equa-
tions,
1 2
− ∇ + vne (r) + vH (r) + vxc,σ (r) ϕiσ (r) = εiσ ϕiσ (r), (1.41)
2
with the spin-dependent exchange-correlation potential,
δExc [ρ↑ , ρ↓ ]
vxc,σ (r) = , (1.42)
δρσ (r)
and the spin density,

Nσ
|ϕiσ (r)|2 .
X
ρσ (r) = (1.43)
i=1
As before, if there is a unique KS wave function Φ[ρ↑ , ρ↓ ] (up to a phase factor), we can
decompose Exc [ρ↑ , ρ↓ ] into exchange and correlation contributions,
Exc [ρ↑ , ρ↓ ] = Ex [ρ↑ , ρ↓ ] + Ec [ρ↑ , ρ↓ ], (1.44)
with Ex [ρ↑ , ρ↓ ] = hΦ[ρ↑ , ρ↓ ]|Ŵee |Φ[ρ↑ , ρ↓ ]i − EH [ρ]. It turns out that the spin-dependent ex-
change functional Ex [ρ↑ , ρ↓ ] can be exactly expressed in terms of the spin-independent exchange
functional Ex [ρ] [34],
1
Ex [ρ↑ , ρ↓ ] = (Ex [2ρ↑ ] + Ex [2ρ↓ ]) , (1.45)
2
which is known as the spin-scaling relation and stems directly from the fact the ↑- and ↓-spin
electrons are uncoupled in the exchange energy [see Eq. (1.23)]. Therefore, any approximation
for Ex [ρ] can be easily extended to an approximation for Ex [ρ↑ , ρ↓ ]. Unfortunately, there is no
such relation for the spin-dependent correlation functional Ec [ρ↑ , ρ↓ ].
Obviously, in the spin-unpolarized case, i.e. ρ↑ = ρ↓ = ρ/2, this spin-dependent formalism
reduces to the spin-independent one.
8
1.4 The generalized Kohn–Sham scheme
An important extension of the KS scheme is the so-called generalized Kohn–Sham (GKS)
scheme [35], which recognizes that the universal density functional F [ρ] of Eq. (1.8) can be
decomposed in other ways than the KS decomposition of Eq. (1.13). In particular, we can
decompose F [ρ] as
n o
F [ρ] = min hΦ|T̂ |Φi + EH [ρΦ ] + S[Φ] + S̄[ρ], (1.46)
Φ∈SρN
where S[Φ] is any functional of a single-determinant wave function Φ ∈ S N leading to a mini-

mum in Eq. (1.46), and S̄[ρ] is the corresponding complementary density functional that makes
Eq. (1.46) exact. Defining the S-dependent GKS exchange-correlation functional as
S
Exc [Φ] = S[Φ] + S̄[ρΦ ], (1.47)
we can express the exact ground-state energy as

n o
S
E0 = inf hΦ|T̂ + V̂ne |Φi + EH [ρΦ ] + Exc [Φ] , (1.48)
Φ∈S N
and if a minimum exists then any minimizing single-determinant wave function in Eq. (1.48)
gives a ground-state density ρ0 (r). Similarly to the KS equations [Eq. (1.27)], Eq. (1.48) leads
to the one-electron GKS equations,

1 2 δS[Φ]
− ∇ + vne (r) + vH (r) + vS̄ (r) ϕiσ (r) + ∗ = εiσ ϕiσ (r), (1.49)
2 δϕiσ (r)
where vS̄ (r) = δS̄[ρ]/δρ(r) is a local potential and δS[Φ]/δϕ∗iσ (r) generates a one-electron (pos-
sibly nonlocal) operator.
In the special case S[Φ] = 0, we recover the KS exchange-correlation density functional:
S=0
Exc [Φ] = Exc [ρΦ ]. (1.50)
Due to the freedom in the choice of S[Φ], there is an infinity of GKS exchange-correlation
S [Φ] giving the exact ground-state energy via Eq. (1.48). This freedom and the
functionals Exc
fact that Φ carries more information than ρΦ gives the possibility to design more accurate
approximations for the exchange-correlation energy.
Of course, by starting from the density functional F [ρ↑ , ρ↓ ] in Eq. (1.35), this GKS scheme
can be extended to the spin-dependent case, leading to GKS exchange-correlation functionals of
S [Φ] = S[Φ] + S̄[ρ
the form Exc ↑,Φ , ρ↓,Φ ].
2 Exact expressions and constraints for the Kohn–Sham exchange and

correlation functionals
2.1 The exchange and correlation holes
Let us consider the pair density associated with the wave function Ψ[ρ] defined in Eq. (1.8),
Z
ρ2 (r1 , r2 ) = N (N − 1) |Ψ[ρ](x1 , x2 , ..., xN )|2 dσ1 dσ2 dx3 ...dxN , (2.1)
{↑,↓}2 ×(R3 ×{↑,↓})N−2
which
R is a functional of the density, and normalized to the number of electron pairs,
R3 ×R3 ρ2 (r1 , r2 )dr1 dr2 = N (N − 1). The pair density is proportional to the probability density
9
of finding two electrons at positions (r1 , r2 ) with all the other electrons being anywhere. The pair
density is useful to express the expectation value of the electron-electron interaction operator,
1 ρ2 (r1 , r2 )
Z
hΨ[ρ]|Ŵee |Ψ[ρ]i = dr1 dr2 . (2.2)
2 R3 ×R3 |r1 − r2 |
Mirroring the decomposition of the Hartree-exchange-correlation energy performed in the KS

scheme [Eq. (1.18)], the pair density can be decomposed as
ρ2 (r1 , r2 ) = ρ(r1 )ρ(r2 ) + ρ2,xc (r1 , r2 ). (2.3)
The product of the densities ρ(r1 )ρ(r

R 2 ) corresponds to the case of independent electrons [up
to a change of normalization, i.e. R3 ×R3 ρ(r1 )ρ(r2 )dr1 dr2 = N 2 instead of N (N − 1)] and the
exchange-correlation pair density ρ2,xc (r1 , r2 ) represents the modification of the pair density due
to exchange and correlation effects between the electrons. It can be further written as
ρ2,xc (r1 , r2 ) = ρ(r1 )hxc (r1 , r2 ), (2.4)
where hxc (r1 , r2 ) is the exchange-correlation hole. Introducing the conditional density ρcond (r1 , r2 ) =
ρ2 (r1 , r2 )/ρ(r1 ) of the remaining N − 1 electrons at r2 given that one electron has been found
at r1 , the exchange-correlation hole can be interpreted as the modification of ρcond (r1 , r2 ) due
to exchange and correlation effects:
ρcond (r1 , r2 ) = ρ(r2 ) + hxc (r1 , r2 ). (2.5)
The positivity of ρ2 (r1 , r2 ) implies that
hxc (r1 , r2 ) ≥ −ρ(r2 ). (2.6)
Moreover, from Eq. (2.5), we have the following sum rule:

Z
∀r1 ∈ R3 , hxc (r1 , r2 )dr2 = −1. (2.7)
R3
We can separate the exchange and correlation contributions in the exchange-correlation hole.
For this, consider the pair density ρ2,KS (r1 , r2 ) associated with the KS single-determinant wave
function Φ[ρ] defined in Eq. (1.14). It can be decomposed as
ρ2,KS (r1 , r2 ) = ρ(r1 )ρ(r2 ) + ρ2,x (r1 , r2 ), (2.8)
where ρ2,x (r1 , r2 ) is the exchange pair density, which is further written as
ρ2,x (r1 , r2 ) = ρ(r1 )hx (r1 , r2 ), (2.9)
where hx (r1 , r2 ) is the exchange hole. Just like the exchange-correlation hole, the exchange hole
satisfies the conditions
hx (r1 , r2 ) ≥ −ρ(r2 ), (2.10)
and Z
∀r1 ∈ R3 , hx (r1 , r2 )dr2 = −1. (2.11)
R3
Moreover, since the exchange hole can be written as [compare with Eq. (1.23)]
1 X
hx (r1 , r2 ) = − |γσ (r1 , r2 )|2 , (2.12)
ρ(r1 )
σ∈{↑,↓}
10
where γσ (r1 , r2 ) = N
P σ ∗
i=1 ϕiσ (r2 )ϕiσ (r1 ) is the spin-dependent one-particle KS density matrix, it
thus appears that the exchange hole is always non-positive,
hx (r1 , r2 ) ≤ 0. (2.13)
From Eqs. (1.21), (2.2), (2.8), and (2.9), it can be seen that the exchange energy functional can
be written in terms of the exchange hole,
1 ρ(r1 )hx (r1 , r2 )

Z
Ex [ρ] = dr1 dr2 , (2.14)
2 R3 ×R3 |r1 − r2 |
leading to the interpretation of Ex as the electrostatic interaction energy of an electron and its
exchange hole. It is useful to write the exchange energy functional as
Z
Ex [ρ] = ρ(r1 )εx [ρ](r1 )dr1 , (2.15)
R3
where εx [ρ](r1 ) is the exchange energy density per particle
1 hx (r1 , r2 )
Z
εx [ρ](r1 ) = dr2 , (2.16)
2 R3 |r1 − r2 |
which is itself a functional of the density. It is also convenient to define the exchange energy
density ex [ρ](r) = ρ(r)εx [ρ](r). For finite systems, we have the exact asymptotic behavior [36,37]
1
εx [ρ](r) ∼ − . (2.17)
|r|→+∞ 2|r|
The correlation hole is defined as the difference
hc (r1 , r2 ) = hxc (r1 , r2 ) − hx (r1 , r2 ), (2.18)
and, from Eqs. (2.7) and (2.11), satisfies the sum rule
Z
3
∀r1 ∈ R , hc (r1 , r2 )dr2 = 0, (2.19)
R3
which implies that the correlation hole has negative and positive contributions3 . In contrast with
the exchange hole which is a smooth function of the interelectronic coordinate r12 = r2 − r1 , the
correlation hole satisfies the electron-electron cusp condition (i.e., it has a derivative discontinuity
in r12 ) [38, 39],
∀r1 ∈ R3 , h′c (r1 , r1 ) = hc (r1 , r1 ), (2.20)
where h′c (r1 , r1 ) = (∂ h̃c (r1 , r12 )/∂r12 )r12 =0
R is the first-order derivative of the spherically averaged
2
correlation hole h̃c (r1 , r12 ) = (1/(4πr12 )) S(0,r12 ) hc (r1 , r1 +r12 )dr12 and S(0, r12 ) designates the
sphere centered at 0 and of radius r12 = |r12 |. The potential contribution to the correlation
energy can be written in terms of the correlation hole:
1 ρ(r1 )hc (r1 , r2 )

Z
Uc [ρ] = dr1 dr2 . (2.21)
2 R3 ×R3 |r1 − r2 |
In order to express the total correlation energy Ec [ρ] = Tc [ρ]+Uc [ρ] in a form similar to Eq. (2.21),
we need to introduce the adiabatic-connection formalism.
3
Therefore, the correlation hole is really a “hole” only in some region of space, and a “bump” in other regions.
11
2.2 The adiabatic connection
The idea of the adiabatic connection [40–42] (see, also, Ref. [43]) is to have a continuous
path between the non-interacting KS system and the physical system while keeping the ground-
state density constant. This allows one to obtain a convenient expression for the correlation
functional Ec [ρ] as an integral over this path. An infinity of such paths are possible, but the one
most often considered consists in switching on the electron-electron interaction linearly with a
coupling constant λ. The Hamiltonian along this adiabatic connection is
Ĥ λ = T̂ + λŴee + V̂ λ , (2.22)
where V̂ λ is the external local potential operator imposing that the ground-state density is the
same as the ground-state density of the physical system for all λ ∈ R. Of course, Eq. (2.22)
relies on a v-representability assumption, i.e. the external potential is assumed to exist for all
λ. The Hamiltonian (2.22) reduces to the KS non-interacting Hamiltonian for λ = 0 and to the
physical Hamiltonian for λ = 1.
Just as for the physical system, it is possible to define a universal functional associated with
the system of Eq. (2.22) for each value of the parameter λ,
F λ [ρ] = min hΨ|T̂ + λŴee |Ψi = hΨλ [ρ]|T̂ + λŴee |Ψλ [ρ]i, (2.23)
Ψ∈WρN
where Ψλ [ρ] denotes a minimizing wave function. This functional can be decomposed as
F λ [ρ] = Ts [ρ] + EHxc

λ
[ρ], (2.24)
λ [ρ] is the Hartree-exchange-correlation functional associated with the interaction
where EHxc
λ [ρ] = E λ [ρ] + E λ [ρ] + E λ [ρ], where the Hartree and
λŴee . One can write this functional as EHxc H x c
exchange contributions are simply linear in λ,
1 λ
Z
λ
EH [ρ] = ρ(r1 )ρ(r2 ) dr1 dr2 = λEH [ρ], (2.25)
2 R3 ×R3 |r1 − r2 |
and
Exλ [ρ] = hΦ[ρ]|λŴee |Φ[ρ]i − EH
λ
[ρ] = λEx [ρ]. (2.26)
The correlation contribution is nonlinear in λ,
Ecλ [ρ] = hΨλ [ρ]|T̂ + λŴee |Ψλ [ρ]i − hΦ[ρ]|T̂ + λŴee |Φ[ρ]i. (2.27)
We will assume that F λ [ρ] is of class C 1 as a function of λ for λ ∈ [0, 1] and that F λ=0 [ρ] =
Ts [ρ], the latter condition being guaranteed for nondegenerate KS systems [see footnote on the
definition of Ts [ρ] just before Eq. (1.14)]. Taking the derivative of Eq. (2.27) with respect to λ
and using the Hellmann-Feynman theorem for the wave function Ψλ [ρ]4 , we obtain
∂Ecλ [ρ]
= hΨλ [ρ]|Ŵee |Ψλ [ρ]i − hΦ[ρ]|Ŵee |Φ[ρ]i. (2.28)
∂λ
∂F λ [ρ] λ
4
In this context, the Hellmann-Feynman theorem states that in the derivative ∂λ
= h ∂Ψ∂λ[ρ] |T̂ +
λ
λŴee |Ψλ [ρ]i + hΨλ [ρ]|Ŵee |Ψλ [ρ]i + hΨλ [ρ]|T̂ + λŴee | ∂Ψ∂λ[ρ] i the first and third terms involving the derivative
of Ψλ [ρ] vanish. This is due to the fact that Ψλ [ρ] is obtained via the minimization of Eq. (2.23) and thus any
variation of Ψλ [ρ] which keeps the density constant (which is the case for a variation with respect to λ) gives a
vanishing variation of F λ [ρ].
12
Integrating over λ from 0 to 1, and using Ecλ=1 [ρ] = Ec [ρ] and Ecλ=0 [ρ] = 0, we arrive at the
adiabatic-connection formula for the correlation energy functional of the physical system
Z 1
Ec [ρ] = dλ hΨλ [ρ]|Ŵee |Ψλ [ρ]i − hΦ[ρ]|Ŵee |Φ[ρ]i. (2.29)
0
By introducing the correlation hole hλc (r1 , r2 ) associated to the wave function Ψλ [ρ], the adiabatic-
connection formula for the correlation energy can also be written as
1 1 ρ(r1 )hλc (r1 , r2 )
Z Z
Ec [ρ] = dλ dr1 dr2 , (2.30)
2 0 R3 ×R3 |r1 − r2 |
or, noting that hλc (r1 , r2 ) is the only quantity that depends on λ in Eq. (2.30), in a more compact
way,
1 ρ(r1 )h̄c (r1 , r2 )
Z
Ec [ρ] = dr1 dr2 , (2.31)
2 R3 ×R3 |r1 − r2 |
R1
where h̄c (r1 , r2 ) = 0 dλ hλc (r1 , r2 ) is the coupling-constant-integrated correlation hole. This
leads to the interpretation of Ec as the electrostatic interaction energy of an electron with its
coupling-constant-integrated correlation hole. As for the exchange energy, the correlation energy
functional can be written as Z
Ec [ρ] = ρ(r1 )εc [ρ](r1 )dr1 , (2.32)
R3
where εc [ρ](r1 ) is the correlation energy density per particle
1 h̄c (r1 , r2 )
Z
εc [ρ](r1 ) = dr2 , (2.33)
2 R3 |r1 − r2 |
which is a functional of the density. We can also define the correlation energy density ec [ρ](r) =
ρ(r)εc [ρ](r).
Finally, note that the sum-rule and cusp conditions of Eqs. (2.19) and (2.20) apply to the
λ-dependent correlation hole in the form
Z
3
∀r1 ∈ R , hλc (r1 , r2 )dr2 = 0, (2.34)
R3
and
∀r1 ∈ R3 , hλc ′ (r1 , r1 ) = λ hλc (r1 , r1 ). (2.35)
2.3 One-orbital and one-electron spatial regions

For systems composed of only one spin-↑ (or, symmetrically, one spin-↓) electron (e.g., the
hydrogen atom) with ground-state density ρ1e (r) = |ϕ1↑ (r)|2 where ϕ1↑ (r) is the unique occupied
KS orbital, the exchange hole in Eq. (2.12) simplifies to hx (r1 , r2 ) = −ρ(r2 ), and consequently
the exchange energy cancels out the Hartree energy:
Ex [ρ1e ] = −EH [ρ1e ]. (2.36)
Furthermore, the correlation energy vanishes:
Ec [ρ1e ] = 0. (2.37)
This must of course also be true for the spin-dependent version of the functionals introduced in
Section 1.3.3, i.e.
Ex [ρ1e , 0] = −EH [ρ1e ] (2.38)
13
and
Ec [ρ1e , 0] = 0. (2.39)
For systems composed of two opposite-spin electrons (e.g., the helium atom or the dihydrogen
molecule) in a unique doubly occupied KS orbital ϕ1 (r) = ϕ1↑ (r) = ϕ1↓ (r) with ground-state
density ρ↑↓ 2
2e (r) = 2|ϕ1 (r)| , the exchange hole simplifies to hx (r1 , r2 ) = −ρ(r2 )/2, and conse-
quently the exchange energy is equal to half the opposite of the Hartree energy:
1
Ex [ρ↑↓ ↑↓
2e ] = − EH [ρ2e ]. (2.40)
2
These are constraints for the exchange and correlation density functionals in the special cases
N = 1 and N = 2.
These special cases can be extended to more general systems. For systems with N ≥ 1
electrons containing a spatial region Ω↑1o in which, among the occupied KS orbitals, only one
spin-↑ (or, symmetrically, one spin-↓) orbital is not zero (or, more generally, takes non-negligible
values), we have again in this region
∀r1 , r2 ∈ Ω↑1o , hx (r1 , r2 ) = −ρ(r2 ), (2.41)
and therefore the contribution to the exchange energy density per particle coming from this
region must locally cancel out the contribution to the Hartree energy density per particle coming
from the same region,
Ω↑ Ω↑
∀r1 ∈ Ω↑1o , εx 1o (r1 ) = −εH1o (r1 ), (2.42)
where εΩ Ω
R R
H (r1 ) = (1/2) Ω ρ(r2 )/|r1 − r2 |dr2 and εx (r1 ) = (1/2) Ω hx (r1 , r2 )/|r1 − r2 |dr2 . Sim-
ilarly, for systems with N ≥ 2 electrons containing a spatial region Ω↑↓ 1o in which, among the
occupied KS orbitals, only one doubly occupied orbital is not zero, we have in this region
1
∀r1 , r2 ∈ Ω↑↓
1o , hx (r1 , r2 ) = − ρ(r2 ), (2.43)
2
and therefore the contribution to the exchange energy density per particle coming from this
region must locally be equal to half the opposite of the contribution to the Hartree energy
density per particle coming from the same region,
Ω↑↓ 1 Ω↑↓
∀r1 ∈ Ω↑↓ 1o 1o
1o , εx (r1 ) = − εH (r1 ). (2.44)
2
Thus, we see, particularly clearly for these Ω↑1o or Ω↑↓1o regions, that the Hartree functional
introduces a spurious self-interaction contribution which must be eliminated by the exchange
functional. Even though the concepts of Ω↑1o and Ω↑↓ 1o regions are formal, in practice they can
be approximately realized in chemical systems. For example, the unpaired electron in a radical
approximately corresponds to a Ω↑1o , and an electron pair in a single covalent bond, in a lone
pair, or in a core orbital approximately corresponds to a Ω↑↓
1o region.
We can also consider one-electron regions Ω1e that we define as5
∀r1 , r2 ∈ Ω1e , ∀λ ∈ (0, 1], ρλ2 (r1 , r2 ) = 0, (2.45)
where ρλ2 (r1 , r2 ) is the pair density associated to the wave function Ψλ [ρ] along the adiabatic
connection. This implies
∀r1 , r2 ∈ Ω1e , h̄xc (r1 , r2 ) = −ρ(r2 ), (2.46)
5
In the definition of Eq. (2.45) we exclude the point λ = 0 in order to allow for the possibility of a discontinuity
in λ there due to a degeneracy.
14
where h̄xc (r1 , r2 ) = hx (r1 , r2 ) + h̄c (r1 , r2 ) and, consequently, the contribution to the exchange-
correlation energy density per particle coming from this region must locally cancel out the
contribution to the Hartree energy density per particle coming from the same region,
Ω1e
∀r1 ∈ Ω1e , εΩ
xc (r1 ) = −εH (r1 ),
1e
(2.47)
where εΩ
R
xc (r1 ) = (1/2) Ω h̄xc (r1 , r2 )/|r1 −r2 |dr2 . For regions that are simultaneously one-electron
and one-orbital regions, this simply implies that the contribution to the correlation energy must
vanish,
Ω ∩Ω↑1o
∀r1 ∈ Ω1e ∩ Ω↑1o , εc 1e (r1 ) = 0, (2.48)
where εΩ
R
c (r1 ) = (1/2) Ω h̄c (r1 , r2 )/|r1 − r2 |dr2 , and we say that the correlation functional must
not introduce a self-interaction error. However, the definition of Ω1e regions also includes the
case of an electron entangled in several orbitals, such as the region around one hydrogen atom
in the dissociated dihydrogen molecule. In this latter case, the Hartree functional introduces
an additional spurious contribution (beyond the spurious self-interaction) which must be com-
pensated by a static correlation (or strong correlation) contribution in the exchange-correlation
functional.
2.4 Coordinate scaling

2.4.1 Uniform coordinate scaling
We consider a norm-preserving uniform scaling of the spatial coordinates in the N -electron

wave function along the adiabatic connection Ψλ [ρ] [introduced in Eq. (2.23)] while leaving
untouched the spin coordinates [44–46],
Ψλγ [ρ](r1 , σ1 , ..., rN , σN ) = γ 3N/2 Ψλ [ρ](γr1 , σ1 , , ..., γrN , σN ), (2.49)
where γ ∈ (0, +∞) is a scaling factor. The scaled wave function Ψλγ [ρ] yields the scaled density
ργ (r) = γ 3 ρ(γr), (2.50)

R R
with R3 ργ (r)dr = R3 ρ(r)dr = N , and minimizes hΨ|T̂ + λγ Ŵee |Ψi since
hΨλγ [ρ]|T̂ + λγ Ŵee |Ψλγ [ρ]i = γ 2 hΨλ [ρ]|T̂ + λŴee |Ψλ [ρ]i. (2.51)
We thus conclude that the scaled wave function at the density ρ and coupling constant λ corre-
sponds to the wave function at the scaled density ργ and coupling constant λγ,
Ψλγ [ρ] = Ψλγ [ργ ], (2.52)
or, equivalently,
Ψλ/γ λ
γ [ρ] = Ψ [ργ ], (2.53)
and that the universal density functional satisfies the scaling relation
F λγ [ργ ] = γ 2 F λ [ρ], (2.54)
or, equivalently,
F λ [ργ ] = γ 2 F λ/γ [ρ]. (2.55)
15
At λ = 0, we find the scaling relation of the KS wave function Φ[ρ] introduced in Section 1.3.1:
Φ[ργ ] = Φγ [ρ]. (2.56)
This directly leads to the scaling relation for the non-interacting kinetic density functional [see
Eq. (1.14)],
Ts [ργ ] = γ 2 Ts [ρ], (2.57)
for the Hartree density functional [see Eq. (1.19)],
EH [ργ ] = γEH [ρ], (2.58)
and for the exchange density functional [see Eq. (1.21)],
Ex [ργ ] = γEx [ρ]. (2.59)
However, the correlation density functional Ec [ρ] has the more complicated scaling (as F [ρ]),
Ecλ [ργ ] = γ 2 Ecλ/γ [ρ], (2.60)
and, in particular for λ = 1,

Ec [ργ ] = γ 2 Ec1/γ [ρ]. (2.61)
These scaling relations allow one to find the behavior of the density functionals in the high-
and low-density limits. In the high-density limit (γ → ∞), it can be shown from Eq. (2.61) that,
for nondegenerate KS systems, the correlation functional Ec [ρ] goes to a constant,
lim Ec [ργ ] = EcGL2 [ρ], (2.62)

γ→∞
where EcGL2 [ρ] is the second-order Görling–Levy (GL2) correlation energy [47, 48] (see Sec-
tion 7.2). This is also called the weak-correlation limit since in this limit the correlation energy
is negligible with respect to the exchange energy which is itself negligible with respect to the
non-interacting kinetic energy: |Ec [ργ ]| = O(γ 0 ) ≪ |Ex [ργ ]| = O(γ) ≪ Ts [ργ ] = O(γ 2 ). Equa-
tion (2.62) is an important constraint since atomic and molecular correlation energies are often
close to the high-density limit. For example, for the ground-state density of the helium atom,
we have Ec [ρ] = −0.0421 hartree and limγ→∞ Ec [ργ ] = −0.0467 hartree [49].
In the low-density limit (γ → 0), it can be shown from Eq. (2.55) that the Hartree-exchange-
correlation energy EHxc [ρ] goes to zero linearly in γ,
SCE
EHxc [ργ ] ∼ γ Wee [ρ], (2.63)
γ→0
SCE
where Wee [ρ] = inf hΨ|Ŵee |Ψi is the strictly-correlated-electron (SCE) functional [50–53].
Ψ∈WρN
This is also called the strong-interaction limit since in this limit the Hartree-exchange-correlation
energy dominates over the non-interacting kinetic energy: EHxc [ργ ] = O(γ) ≫ Ts [ργ ] = O(γ 2 ).
In this limit, the electrons strictly localize relatively to each other. In particular, for the uniform-
electron gas, this corresponds to the Wigner crystallization. Thus, in this limit, each electron
is within a one-electron region Ω1e [as defined in Eq. (2.45)]. For more information on the SCE
functional, see the chapter by Friesecke, Gerolin, and Gori-Giorgi in this volume.
16
2.4.2 Non-uniform coordinate scaling
We can also consider non-uniform one-dimensional or two-dimensional coordinate scalings of

the density [54, 55],
ρ(1)
γ (x, y, z) = γρ(γx, y, z) (2.64)
and
ρ(2) 2
γ (x, y, z) = γ ρ(γx, γy, z), (2.65)
which also preserve the number of the electrons. These non-uniform density scalings provide
constraints for the exchange and correlation functionals. In particular, in the non-uniform
one-dimensional high-density limit, the exchange functional remains finite and the correlation
functional vanishes [45, 56]
lim Ex [ρ(1)
γ ] > −∞ (2.66)
γ→∞
and
lim Ec [ρ(1)
γ ] = 0. (2.67)
γ→∞
Also, in the non-uniform two-dimensional low-density limit, we have [45, 56]:

1
lim Ex [ρ(2)
γ ] > −∞ (2.68)
γ→0 γ
and
1
lim Ec [ρ(2)
γ ] = 0. (2.69)
γ→0 γ
The conditions of Eqs. (2.66)-(2.69) are particularly useful because they also correspond to the
limit of rapidly varying densities [57].
2.5 Atoms in the limit of large nuclear charge

A practical realization of the uniform high-density limit is provided by atomic ions in the limit
of large nuclear charge, Z → ∞, at fixed electron number N (see Refs. [58–61]). In this limit, the
exact ground-state atomic density ρN,Z (r) becomes the density of the isoelectronic hydrogenic
(i.e., without electron-electron interaction) atom ρH N,Z (r) which obeys a simple scaling with Z:
ρN,Z (r) ∼ ρH 3 H
N,Z (r) = Z ρN,Z=1 (Zr). (2.70)
Z→∞
One can thus apply Eqs. (2.59) and (2.62) with γ = Z, which reveals that in an isoelectronic
series the exchange functional scales linearly with Z,
Ex [ρN,Z ] ∼ Ex [ρH
N,Z=1 ]Z, (2.71)
Z→∞
and, for nondegenerate KS systems, the correlation functional saturates to a constant,
lim Ec [ρN,Z ] = EcGL2 [ρH

N,Z=1 ]. (2.72)
Z→∞
17
Equations (2.71) and (2.72) are constraints for the exchange and correlation functionals, partic-
ularly relevant for highly ionized atoms but also for the core-electron regions of heavy atoms in
neutral systems.
Another very interesting limit is the one of large nuclear charge of neutral atoms, N = Z → ∞
(see, e.g., Ref. [62]). In this semiclassical limit, the exact ground-state atomic density ρN,Z (r)
tends to the Thomas–Fermi (TF) density of a neutral atom ρTF Z (r) which has a known scaling
with Z [63, 64]:
ρZ,Z (r) ∼ ρTF 2 TF
Z (r) = Z ρZ=1 (Z
1/3
r). (2.73)
Z→∞
In this limit, it was suggested that the exact exchange and correlation energies have the approx-
imate large-Z asymptotic expansions [65–67]
Ex [ρZ,Z ] ∼ −Ax Z 5/3 + Bx Z + · · · (2.74)

Z→∞
and
Ec [ρZ,Z ] ∼ −Ac Z ln Z + Bc Z + · · · , (2.75)
Z→∞
with the coefficients Ax = 0.220827, Ac = 0.020727, Bx ≈ 0.224, Bc ≈ 0.0372. Recently, it was
argued that there is in fact a missing term in Z ln Z in the expansion of the exchange energy in
Eq. (2.74) [68, 69].
2.6 Lieb–Oxford lower bound

Lieb and Oxford derived a lower bound for the indirect Coulomb energy (i.e., the two-particle
Coulomb potential energy beyond the Hartree energy) [70], which, when expressed in terms of
the exchange or exchange-correlation functional, takes the form [71]
Z
Ex [ρ] ≥ Exc [ρ] ≥ −CLO ρ(r)4/3 dr, (2.76)
R3
where the optimal (i.e., smallest) constant CLO (independent of the electron number N ) was
originally shown to be in the range 1.23 ≤ CLO ≤ 1.68 [70]. The range was later successively
narrowed to 1.4442 ≤ CLO ≤ 1.5765 [71–74]. This bound is approached only in the low-density
limit where the correlation energy becomes comparable to the exchange energy. Numerical
results suggest that for densities of most physical systems the Lieb–Oxford lower bound on the
exchange-correlation energy is far from being reached [75].
For two-electron densities, there is a specific tighter bound,
Z
Ex [ρ2e ] ≥ Exc [ρ2e ] ≥ −C2 ρ2e (r)4/3 dr, (2.77)
R3
with the best known constant C2 = 1.234 [70]. For one-electron densities, an even tighter bound
is known for the exchange functional [70, 76],
Z
Ex [ρ1e ] ≥ −C1 ρ1e (r)4/3 dr, (2.78)
R3
with the optimal constant C1 = 1.092. For two-electron spin-unpolarized densities, we have
Ex [ρ↑↓ ↑↓
2e ] = 2Ex [ρ1e ] with ρ1e = ρ2e /2, and Eq. (2.78) implies [77]
C1
Z
↑↓
Ex [ρ2e ] ≥ − 1/3 ρ↑↓
2e (r)
4/3
dr, (2.79)
2 R 3
which is a much tighter bound than the bounds of Eqs. (2.76) and (2.77).
18
3 Semilocal approximations for the exchange-correlation energy
We review here the different classes of semilocal approximations for the exchange-correlation
energy.
3.1 The local-density approximation

In the local-density approximation (LDA), introduced by Kohn and Sham [19], the exchange-
correlation functional is approximated as
Z
LDA
Exc [ρ] = eUEG
xc (ρ(r))dr, (3.1)
R3
where eUEG
xc (ρ) is the exchange-correlation energy density of the infinite uniform electron gas
(UEG) with the density ρ. The UEG represents a family of systems of interacting electrons
with an arbitrary spatially constant density ρ ∈ [0, +∞) that acts as a parameter. Thus, in the
LDA, the exchange-correlation energy density of an inhomogeneous system at a spatial point of
density ρ(r) is approximated as the exchange-correlation energy density of the UEG of the same
density.
In the spin-dependent version of LDA, sometimes specifically referred to as the local-spin-
density approximation (LSDA), the exchange-correlation functional is approximated as [31]
Z
LSDA
Exc [ρ↑ , ρ↓ ] = eUEG
xc (ρ↑ (r), ρ↓ (r))dr, (3.2)
R3
where eUEG
xc (ρ↑ , ρ↓ ) is the exchange-correlation energy density of the UEG with spin densities
ρ↑ and ρ↓ . For spin-unpolarized systems, we recover the spin-independent LDA as Exc LDA [ρ] =
LSDA
Exc [ρ/2, ρ/2].
The function eUEG
xc is a sum of exchange and correlation contributions, eUEG
xc = eUEG
x + eUEG
c ,
UEG
and it is convenient to introduce exchange and correlation energies per particle, εx UEG
and εc ,
such that eUEG
x = ρ εUEG and eUEG = ρ εUEG . The expression of the exchange energy per
x c c
particle of the spin-unpolarized UEG is
εUEG
x (ρ) = Cx ρ1/3 , (3.3)
where Cx = −(3/4)(3/π)1/3 , and the spin-polarized version is simply obtained from the spin-
scaling relation [Eq. (1.45)], leading to
εUEG
x (ρ↑ , ρ↓ ) = εUEG
x (ρ)φ4 (ζ), (3.4)
where ζ = (ρ↑ − ρ↓ )/ρ is the spin polarization and φ4 (ζ) is defined by the general spin-scaling
function
(1 + ζ)n/3 + (1 − ζ)n/3
φn (ζ) = . (3.5)
2
The LDA exchange functional is associated with the names of Dirac [78] and Slater [79]. For a
rigorous mathematical derivation of Eq. (3.3), see Ref. [80].
The correlation energy per particle εUEG
c (ρ↑ , ρ↓ ) of the UEG cannot be calculated analytically.
This quantity has been obtained numerically for a sample of densities and fitted to a parametrized
function satisfying the known high- and low-density expansions. Expressed in terms of the
Wigner–Seitz radius rs = (3/(4πρ))1/3 , the first terms of the high-density expansion (rs → 0)
have the form
εUEG
c (ρ↑ , ρ↓ ) = A(ζ) ln rs + B(ζ) + C(ζ)rs ln rs + O(rs ), (3.6)
19
with spin-unpolarized coefficients A(0) = (1 − ln 2)/π 2 , B(0) = −0.046921, C(0) = 0.009229,
and fully spin-polarized coefficients A(1) = A(0)/2, B(1) = −0.025738, C(1) = 0.004792. The
first terms of the low-density expansion (rs → +∞) have the form
!
a b c 1
εUEG
c (ρ↑ , ρ↓ ) = + + +O , (3.7)
rs rs3/2 rs2 5/2
rs
where the coefficients a = −0.895930, b = 1.325, and c = −0.365 are assumed to be independent
of ζ. The low-density limit of the UEG corresponds to the Wigner crystallization. For a recent
review of results on the UEG, see Ref. [81].
The two most used parametrizations are the one of Vosko, Wilk, and Nusair (VWN) [82]
and the more recent one of Perdew and Wang (PW92) [83] which we give here. In this
parametrization, the UEG correlation energy per particle is estimated using the approximate
spin-interpolation formula
f (ζ)
εPW92
c (ρ↑ , ρ↓ ) = εc (rs , 0) + αc (rs ) (1 − ζ 4 ) + [εc (rs , 1) − εc (rs , 0)]f (ζ)ζ 4 , (3.8)
f ′′ (0)
where εc (rs , ζ) is the UEG correlation energy per particle as a function of rs and ζ, f (ζ) =
[(1 + ζ)4/3 + (1 − ζ)4/3 − 2]/(24/3 − 2) is a spin-scaling function borrowed from the exchange
energy, and αc (rs ) = (∂ 2 εc (rs , ζ)/∂ζ 2 )ζ=0 is the spin stiffness. This spin-interpolation formula
was first proposed in the VWN parametrization based on a study of the ζ dependence of the
UEG correlation energy per particle at the random-phase approximation (RPA) level. A unique
parametrization function
 
1
G(rs , A, α1 , β1 , β2 , β3 , β4 ) = −2(1 + α1 rs )A ln 1 +  , (3.9)
1/2 3/2
2A β1 rs + β2 rs + β3 rs + β4 rs2
is then used for approximating εc (rs , 0), εc (rs , 1), and −αc (rs ), where
εc (rs , 0) = G(rs , A0 , α1,0 , β1,0 , β2,0 , β3,0 , β4,0 ), (3.10a)
εc (rs , 1) = G(rs , A1 , α1,1 , β1,1 , β2,1 , β3,1 , β4,1 ), (3.10b)
−αc (rs ) = G(rs , A2 , α1,2 , β1,2 , β2,2 , β3,2 , β4,2 ). (3.10c)
The form of G was chosen to reproduce the form of the high- and low-density expansions. The
parameters Ai , β1,i , and β2,i (with i ∈ {0, 1, 2}) are fixed by the first two terms of the high-
density expansion, while the parameters α1,i , β3,i , and β4,i are fitted to quantum Monte Carlo
(QMC) data [84] for εc (rs , 0) and εc (rs , 1), and to an estimation of −αc (rs ) extrapolated from
RPA data. The parameters are given in Table I of Ref. [83].
We now discuss the merits and deficiencies of the LDA. By construction, the LDA is of course
exact in the limit of uniform densities. More relevant to atomic and molecular systems is that
the LDA exchange and correlation energies are asymptotically exact in the limit of large nuclear
charge of neutral atoms N = Z → ∞. Indeed, in this semiclassical Thomas–Fermi limit, the
LDA gives the exact coefficients Ax and Ac of the leading terms in the asymptotic expansions
of Eqs. (2.74) and (2.74) [85]. However, the coefficients of the next terms are very different:
BxLDA ≈ 0 instead of Bx ≈ 0.224 and BcLDA ≈ −0.00451 instead of Bc ≈ 0.0372 [66].
Due to the scaling of the UEG exchange energy per particle,
εUEG
x (γ 3 ρ↑ , γ 3 ρ↓ ) = γεUEG
x (ρ↑ , ρ↓ ), (3.11)
20
the LDA exchange functional correctly scales linearly under uniform coordinate scaling of the
density [Eq. (2.59)]. Similarly, due the scaling of the UEG correlation energy per particle in the
low-density limit [Eq. (3.7)],
a
εUEG
c (γ 3 ρ↑ , γ 3 ρ↓ ) ∼ γ , (3.12)
γ→0 rs
the LDA correlation functional correctly scales linearly under uniform coordinate scaling to the
low-density limit [Eq. (2.63)]. However, from the behavior of εUEG c in the high-density limit
[Eq. (3.6)],
εUEG
c (γ 3 ρ↑ , γ 3 ρ↓ ) ∼ −A(ζ) ln γ, (3.13)
γ→∞
we see that the LDA correlation functional diverges logarithmically under uniform coordinate
scaling to the high-density limit whereas the exact correlation functional goes to a constant
for nondegenerate KS systems [Eq. (2.62)]. Consequently, in the limit of large nuclear charge,
Z → ∞, at fixed electron number N , the LDA exchange energy correctly scales linearly with
Z [Eq. (2.71)], albeit with an incorrect coefficient, and the LDA correlation energy does not
reproduce the exact saturation behavior [Eq. (2.72)] for a nondegenerate isoelectronic series but
incorrectly diverges [86]. Also, the LDA exchange and correlation functionals do not satisfy the
non-uniform scaling conditions of Eqs. (2.66)-(2.69), but instead both diverge in these limits.
The LDA can also be thought of as approximating the exchange and the (coupling-constant-
integrated) correlation holes of an inhomogeneous system in Eqs. (2.16) and (2.33) by the corre-
sponding exchange and correlation holes of the UEG. Namely, considering the spin-independent
version for simplicity, the LDA exchange hole is
hLDA
x (r1 , r2 ) = hUEG
x (ρ(r1 ), r12 ), (3.14)
with 2
9 j1 (kF r12 )
hUEG
x (ρ, r12 ) = −ρ , (3.15)
2 kF r12
where r12 = |r2 − r1 | is the interelectronic distance, kF = (3π 2 ρ)1/3 is the Fermi wave vector,
and j1 is the spherical Bessel function of the first kind. Similarly, the LDA correlation hole is
Z 1
h̄LDA
c (r ,
1 2r ) = h̄UEG
c (ρ(r 1 ), r12 ) = dλ hλ,UEG
c (ρ(r1 ), r12 ). (3.16)
0
Since the UEG is a physical system, the LDA exchange hole correctly fulfills the negativity
and sum-rule condition [Eqs. (2.11) and (2.13)] and the LDA correlation hole correctly fulfills
the sum-rule and electron-electron cusp condition [Eqs. (2.34) and (2.35)]. This constitutes
a significant merit of the LDA. However, because the LDA exchange hole hLDA x (r1 , r2 ) only
depends on ρ(r1 ) and not on ρ(r2 ), the LDA exchange functional does not entirely eliminate
the self-interaction contribution of the Hartree functional, in particular in one and two-electron
systems [Eqs. (2.36) or (2.38), and (2.40)], or in one-orbital spatial regions of many-electron
systems [Eqs. (2.42) and (2.44)]. Similarly, the LDA correlation functional does not vanish in
one-electron systems [Eq. (2.37) or (2.39)], or more generally in one-orbital one-electron regions
[Eq. (2.48)]. Thus, the LDA introduces a self-interaction error. Moreover, the LDA exchange-
correlation functional does not entirely cancel out the Hartree energy in entangled one-electron
spatial regions [Eq. (2.47], i.e. it introduces a static-correlation error.
Another deficiency of the LDA is that the (spin-independent) LDA exchange potential
δExLDA [ρ] 4
vxLDA (r) = = Cx ρ(r)1/3 , (3.17)
δρ(r) 3
21
decays exponentially at infinity for finite molecular systems (since the density ρ(r) decays ex-
ponentially), i.e. much too fast in comparison to the −1/|r| asymptotic behavior of the exact
exchange potential [Eq. (1.33)]. Since asymptotic spatial regions are dominated by the highest
occupied molecular orbital (HOMO) and are thus one-orbital regions (assuming the HOMO is
not degenerate), this is another signature of the incorrectness of the LDA exchange functional
in these one-orbital regions.
For a review of mathematical results on the LDA, see the chapter by Lewin, Lieb, and
Seiringer in this volume.
3.2 The gradient-expansion approximation

The next logical step beyond the LDA is the gradient-expansion approximation (GEA) [19],
in which the exchange-correlation functional is systematically expanded in the gradient and
higher-order derivatives of the density. One way of deriving the GEA is to start from the
UEG, introduce a weak and slowly-varying external potential δv(r), and expand the exchange-
correlation energy in terms of the gradients of the density (see, e.g., Refs. [87–90]). Alternatively,
one can perform a semiclassical expansion (i.e., an expansion in powers of the reduced Planck
constant ~) of the exact Exc [ρ] in terms of the gradients of the external potential and use the
mapping between the potential and the density to express it in terms of the gradients of the
density (see, e.g., Ref. [2]).
The spin-independent gradient expansion of the exchange functional is known up to fourth
order (GEA4) [90],
|∇ρ(r)|2
Z
GEA4 LDA (2)
Ex [ρ] = Ex [ρ] + Cx 4/3
dr
R3 ρ(r)
|∇2 ρ(r)|2 |∇ρ(r)|2 ∇2 ρ(r)
Z Z
(4) (4)
+Cx,1 dr + C x,2 dr, (3.18)
R3 ρ(r)2 R3 ρ(r)3
involving the density gradient ∇ρ(r) and Laplacian ∇2 ρ(r). Sham [91] obtained the second-
(2)
order coefficient Cx,S = −7/(432π(3π 2 )1/3 ) ≈ −0.001667. The calculation was done by starting
with the screened Yukawa interaction e−κr12 /r12 and taking the limit κ → 0 at the end of the
calculation. It was later shown that this calculation contains an order-of-limit problem and that
(2)
the correct Coulombic second-order coefficient is Cx = −5/(216π(3π 2 )1/3 ) ≈ −0.002382 [88,89].
(4) (4)
The fourth-order coefficients are Cx,1 = −73/(64800π 3 ) ≈ −0.000036, and Cx,2 ≈ 0.00009, where
the last one has been numerically estimated [90]. Note that each term in Eq. (3.18) correctly
fulfills the scaling relation of Eq. (2.59). The spin-dependent gradient exchange expansion is
simply obtained from the spin-scaling relation [Eq. (1.45)].
Similarly, Ma and Brueckner [87] obtained the spin-independent second-order gradient ex-
pansion (GEA2) of the correlation functional,
|∇ρ(r)|2
Z
GEA2 LDA
Ec [ρ] = Ec [ρ] + Cc(2) (rs (r)) dr, (3.19)
R3 ρ(r)4/3
(2)
with a second-order coefficient in the high-density limit Cc,MB (rs → 0) = 0.004235. It is be-
lieved [92] that this calculation contains a similar order-of-limit problem as in Sham’s coef-
(2)
ficient Cx,S , in such a way that these two coefficients must be combined to obtain the cor-
(2)
rect second-order exchange-correlation coefficient in the high-density limit Cxc (rs → 0) =
(2) (2)
Cx,S + Cc,MB (rs → 0). The correct second-order correlation coefficient in the high-density limit
(2) (2) (2)
is then Cc (rs → 0) = Cxc (rs → 0) − Cx = 0.004950. Similarly, the second-order correlation
22
(2) (2) (2) (2)
coefficient as a function of rs can be obtained by Cc (rs ) = Cxc (rs ) − Cx where Cxc (rs ) has
been parametrized in Ref. [93]. The spin-dependent generalization has the form [94, 95]
X Z ′ ∇ρσ (r) ∇ρσ′ (r)
EcGEA2 [ρ↑ , ρ↓ ] = EcLSDA [ρ↑ , ρ↓ ] + Ccσ,σ ,(2) (rs (r), ζ(r)) · dr, (3.20)
′ R3 ρσ (r)2/3 ρσ′ (r)2/3
σ,σ ∈{↑,↓}
σ,σ′ ,(2)
where the functions Cc (rs , ζ) have been numerically calculated in the high-density limit [94,
96].
The GEA should improve over the LDA for sufficiently slowly varying densities. Since the
spin-independent GEA2 exchange energy per particle has the form
εGEA2
x (ρ, ∇ρ) = ρ1/3 (Cx + Cx(2) x2 ), (3.21)
where x = |∇ρ|/ρ4/3 is a dimensionless reduced density gradient, the precise condition for
exchange is x ≪ 1. Unfortunately, for real systems like atoms and molecules, the reduced
density gradient x can be large in some regions of space. In particular, in the exponential
density tail, ρ(r) ∝ e−α|r| , the reduced density gradient diverges x(r) −−−−→ ∞. But this
|r|→∞ |r|→∞
is not as bad as it seems since ρ εGEA2
x goes to zero anyway in this limit. The situation is
more catastrophic for correlation. Indeed, in the high-density limit, the spin-independent GEA2
correlation energy per particle behaves as
εGEA2
c (γ 3 ρ, γ 4 ∇ρ) ∼ −A(0) ln γ + γ 1/2 Cc(2) (rs → 0)y 2 , (3.22)
γ→∞
where y = |∇ρ|/ρ7/6 is another reduced density gradient adapted to correlation. Therefore, in

this limit, the GEA2 correlation correction diverges to +∞ even faster than the LDA diverges
to −∞.
Another aspect of the deficiency of the GEA is that the corresponding GEA exchange and
correlation holes have unphysical long-range parts which break the negativity [Eq. (2.13)] and
sum-rule conditions [Eqs. (2.11) and (2.19)].
In practice, the GEA tends to deteriorate the results obtained at the LDA level. Truncated
gradient expansions should not be directly used but need to be resummed.
3.3 Generalized-gradient approximations

The failures of the GEA led to the development, which really started in the 1980s, of
generalized-gradient approximations (GGAs) with the generic form
Z
GGA
Exc [ρ↑ , ρ↓ ] = eGGA
xc (ρ↑ (r), ρ↓ (r), ∇ρ↑ (r), ∇ρ↓ (r))dr, (3.23)
R3
where eGGA
xc is some function. The GGAs are often called semilocal approximations in the
sense that eGGA
xc does not only use the local value of the spin densities ρ↑ (r) and ρ↓ (r) but also
“semilocal information” through its gradients6 ∇ρ↑ (r) and ∇ρ↓ (r).
Many GGA functionals have been proposed. They generally provide a big improvement
over LDA for molecular systems. However, their accuracy is still limited, in particular by self-
interaction and static-correlation errors. We review here some of the most widely used GGA
functionals.
6
For generality and simplicity, we consider here that the GGAs depend on the spin density gradients ∇ρ↑ and
∇ρ↓ , but due to rotational invariance GGAs actually depend only on the scalar quantities (∇ρ↑ )2 , (∇ρ↓ )2 , and
∇ρ↑ · ∇ρ↓ .
23
B88 exchange functional
In the Becke 88 (B88 or B) exchange functional [37], the exchange energy density is written
as
X βx2σ
eB88 UEG
x (ρ↑ , ρ↓ , ∇ρ↑ , ∇ρ↓ ) = ex (ρ↑ , ρ↓ ) − ρ4/3
σ , (3.24)
1 + 6βxσ sinh−1 (xσ )
σ∈{↑,↓}
4/3 4/3
where xσ = |∇ρσ |/ρσ . The fact that eB88x depends linearly on ρσ and nonlinearly only on
the dimensionless reduced density gradient xσ (r) guarantees the scaling relation of Eq. (2.59).
Using the exponential decay of the ground-state spin densities of Coulombic systems, ρσ (r) ∝
|r|→∞
e−ασ |r| , it can be verified that the chosen form for eB88
x satisfies the exact asymptotic behavior
of the exchange energy density per particle [Eq. (2.17)], although the corresponding exchange
potential does not satisfy the exact −1/r asymptotic behavior [Eq. (1.33)] [97]. For small xσ ,
eB88
x is correctly quadratic in xσ . The parameter β = 0.0042 was found by fitting to HF exchange
energies of rare-gas atoms. A very similar value of β can also be found by imposing the coefficient
Bx of the approximate large-Z asymptotic expansion of the exchange energy of neutral atoms
[Eq. (2.74)] [65]. It turns out that imposing the coefficient of the second-order gradient expansion
[Eq. (3.18)] would lead to a value of β about two times smaller and would greatly deteriorate
the accuracy of the functional for atoms and molecules.
LYP correlation functional

The Lee–Yang–Parr (LYP) [98] correlation functional is one of the rare functionals which have
not been constructed starting from LDA. It originates from the Colle–Salvetti [99] correlation-
energy approximation depending on the curvature of the HF hole. By using a gradient-expansion
approximation of the curvature of the HF hole, LYP turned the Colle-Salvetti expression into
a density functional depending on the density, the density gradient, and the Laplacian of the
density. The dependence on the Laplacian of the density can be exactly eliminated by an
integration by parts [100], giving the following correlation energy density
(
LYP 4 ρ↑ ρ↓
ec (ρ↑ , ρ↓ , ∇ρ↑ , ∇ρ↓ ) = −a − a b ω(ρ) ρ↑ ρ↓
1 + dρ−1/3 ρ
" #
X
5 δ(ρ)

δ(ρ) − 11 ρσ

47 7δ(ρ)

11/3 8/3 2 2 2
× 2 C F ρσ − − |∇ρσ | − |∇ρσ | + − |∇ρ|
2 18 9 ρ 18 18
σ∈{↑,↓}
)
2 2 2 2 2 2 2 2 2 2 2
− ρ |∇ρ| + ρ − ρ↑ |∇ρ↓ | + ρ − ρ↓ |∇ρ↑ | , (3.25)
3 3 3
where ω(ρ) = ρ−11/3 exp(−cρ−1/3 )/(1 + dρ−1/3 ), δ(ρ) = cρ−1/3 + dρ−1/3 /(1 + dρ−1/3 ), and
CF = (3/10)(3π 2 )2/3 . The parameters a = 0.04918, b = 0.132, c = 0.2533, and d = 0.349 were
obtained in the original Colle-Salvetti expression by a fit to Helium data. Note that the LYP
correlation energy vanishes for fully spin-polarized densities (ρ↑ = 0 or ρ↓ = 0) and therefore
correctly vanishes for one-electron systems [Eq. (2.39)].
PW91 exchange-correlation functional

The Perdew–Wang 91 (PW91) (see Refs. [71, 101, 102]) exchange-correlation functional is
based on a model of the exchange hole hx (r1 , r2 ) in Eq. (2.16) and of the coupling-constant-
integrated correlation hole h̄c (r1 , r2 ) in Eq. (2.33). The idea is to start from the GEA model
of these holes given as gradient expansions and remove the unrealistic long-range parts of these
24
holes to restore important constraints satisfied by the LDA. Specifically, the spurious positive
parts of the GEA exchange hole are removed to enforce the negativity condition of Eq. (2.13)
and a cutoff in |r1 − r2 | is applied to enforce the normalization condition of Eq. (2.11). Similarly,
a cutoff is applied on the GEA correlation hole to enforce the condition that the hole integrates
to zero [Eq. (2.34)]. The exchange and correlation energies per particle calculated from these
numerical holes are then fitted to functions of the density and density gradient chosen to satisfy
a number of exact constraints.
The spin-independent PW91 exchange energy density is written as
ePW91
x (ρ, ∇ρ) = eUEG
x (ρ)FxPW91 (s), (3.26)
where the so-called enhancement factor is

1 + 0.19645s sinh−1 (7.7956s) + [0.2743 − 0.1508 exp(−100s2 )]s2
FxPW91 (s) = , (3.27)
1 + 0.19645s sinh−1 (7.7956s) + 0.004s4
with the reduced density gradient s = |∇ρ|/(2kF ρ) = x/[2(3π 2 )1/3 ] where kF = (3π 2 ρ)1/3
is the Fermi wave vector. The spin-dependent PW91 exchange energy density is simply ob-
tained from the spin-scaling relation [Eq. (1.45)]: ePW91
x (ρ↑ , ρ↓ , ∇ρ↑ , ∇ρ↓ ) = [ePW91
x (2ρ↑ , 2∇ρ↑ )+
PW91
ex (2ρ↓ , 2∇ρ↓ )]/2. The enhancement factor Fx PW91 (s) satisfies the second-order gradient ex-
PW91 (2)
pansion [Eq. (3.18)], Fx (s) = 1+µ s +O(s ) with µ = −16π(π/3)2/3 Cx = 10/81 ≈ 0.1235,
2 4
the local Lieb–Oxford bound, FxPW91 (s) ≤ −CLO /(Cx 21/3 ) ≈ 1.804, which is a sufficient and nec-
essary condition for a spin-dependent GGA exchange functional to satisfy the Lieb–Oxford lower
bound [Eq. (2.76)] for all densities [77] (note however that 1.804 is not an optimal bound), and
the condition lims→∞ s1/2 FxPW91 (s) < ∞ which guarantees the non-uniform scaling finiteness
conditions of Eqs. (2.66) and (2.68) [57, 77].
The PW91 correlation energy density is written as
ePW91 (ρ↑ , ρ↓ , ∇ρ↑ , ∇ρ↓ ) = ρ εUEG (ρ↑ , ρ↓ ) + H PW91 (ρ↑ , ρ↓ , t) ,

c c (3.28)
where the gradient correction H PW91 (ρ↑ , ρ↓ , t) = H0 (ρ↑ , ρ↓ , t) + H1 (ρ↑ , ρ↓ , t) depends on another
reduced density 1/6
p gradient (adapted to correlation) t = |∇ρ|/(2φ2 (ζ)ks ρ) = y/[4φ2 (ζ)(3/π) ]
where ks = 4kF /π is the Thomas–Fermi screening wave vector and the spin-scaling function
φ2 (ζ) is defined by Eq. (3.5), with
2 1 + At2

3β 2α 2
H0 (ρ↑ , ρ↓ , t) = φ2 (ζ) ln 1 + t , (3.29)
2α β 1 + At2 + A2 t4
2α −1
A= exp(−2αεUEG
c (ρ↑ , ρ↓ )/(φ2 (ζ)3 β 2 )) − 1 , (3.30)
β
1/3
3 (2) 4 2 2 2
H1 (ρ↑ , ρ↓ , t) = 16 [Cxc (rs ) − Cc,MB (rs → 0) − Cx(2) ]φ2 (ζ)3 t2 e−100φ2 (ζ) ks t /kF , (3.31)
π
and Cxc (rs ) is taken from Ref. [93]. The function H0 (ρ↑ , ρ↓ , t) was chosen so that it fulfills the
second-order gradient expansion [Eq. (3.20)], H0 (ρ↑ , ρ↓ , t) = βφ2 (ζ)3 t2 + O(t4 ), using an approx-
imate ζ dependence [95] and the Ma-Brueckner high-density-limit second-order coefficient [87]
(2)
β = 16(3/π)1/3 Cc,MB (rs → 0) ≈ 0.06673, and so that it cancels the LDA correlation in the
large-t limit, limt→∞ H0 (ρ↑ , ρ↓ , t) = −εUEG
c (ρ↑ , ρ↓ ). The only fitted parameter is α = 0.09.
The function H1 (ρ↑ , ρ↓ , t) only serves to restore the correct second-order gradient expansion,
25
such that H PW91 (ρ↑ , ρ↓ , t) = 16(3/π)1/3 Cxc (rs )φ2 (ζ)3 t2 + O(t4 ), while keeping the large-t limit
unchanged.
PBE exchange-correlation functional

The Perdew–Burke–Ernzerhof (PBE) [103] exchange-correlation functional is a simplification
of the PW91 functional with no fitted parameters which gives almost the same energies. The
spin-independent PBE exchange energy density is written as
ePBE
x (ρ, ∇ρ) = eUEG
x (ρ)FxPBE (s), (3.32)
where the enhancement factor is

κ
FxPBE (s) = 1 + κ − . (3.33)
1 + µs2 /κ
The function FxPBE (s) has the second-order gradient expansion FxPBE (s) = 1 + µ s2 + O(s4 ),
(2)
and the parameter is chosen as µ = 16π(π/3)2/3 Cc,MB (rs → 0) ≈ 0.21951 so as to cancel the
correlation second-order gradient expansion. The second parameter κ is chosen so as to saturate
the local Lieb–Oxford bound, i.e. lims→∞ FxPBE (s) = 1 + κ = −CLO /(Cx 21/3 ) ≈ 1.804, leading
to κ = 0.804. The same exchange functional form was in fact proposed earlier in the Becke 86
(B86) functional [104] with empirical parameters (µ = 0.235, κ = 0.967).
A revised version of the PBE exchange functional, called revPBE, was proposed where the
local Lieb–Oxford bound constraint is relaxed and the parameter κ = 1.245 is found instead by
fitting to exchange-only total atomic energies for He and Ar, resulting in more accurate atomic
total energies and molecular atomization energies [105]. Another revised version of the PBE
exchange functional, called RPBE, was also proposed to achieve a similar improvement, while
still enforcing the local Lieb–Oxford bound, by changing the form of the enhancement factor to
FxRPBE (s) = 1 + κ(1 − exp(−µs2 /κ)) with the same parameters as in the original PBE [106].
The PBE correlation energy density is written as
ePBE (ρ↑ , ρ↓ , ∇ρ↑ , ∇ρ↓ ) = ρ εUEG (ρ↑ , ρ↓ ) + H PBE (ρ↑ , ρ↓ , t) ,

c c (3.34)
with the gradient correction
1 + At2

PBE 3 β 2
H (ρ↑ , ρ↓ , t) = A(0)φ2 (ζ) ln 1 + t (3.35)
A(0) 1 + At2 + A2 t4
and
β −1
A= exp(−εUEG
c (ρ↑ , ρ↓ )/(A(0)φ2 (ζ)3 )) − 1 . (3.36)
A(0)
As in the PW91 correlation functional, the function H PBE (ρ↑ , ρ↓ , t) has the second-order gradient
(2)
expansion H PBE (ρ↑ , ρ↓ , t) = βφ2 (ζ)3 t2 + O(t4 ) where β = 16(3/π)1/3 Cc,MB (rs → 0) ≈ 0.06673,
and it cancels the LDA correlation in the large-t limit, limt→∞ H PBE (ρ↑ , ρ↓ , t) = −εUEG
c (ρ↑ , ρ↓ ).
In contrast with the PW91 correlation functional, under uniform coordinate scaling to the high-
density limit, the PBE correlation functional correctly cancels out the logarithm divergence of
the LDA correlation functional [Eq. (3.13)], i.e. H PBE (γ 3 ρ↑ , γ 3 ρ↓ , γ 1/2 t) ∼ A(0)φ2 (ζ)3 ln γ
γ→∞
where A(0)φ2 (ζ)3 is a good approximation to the coefficient A(ζ) [95].
A variant of the PBE exchange-correlation functional, called PBEsol [107], targeted for
solid-state systems, was proposed where the correct second-order exchange gradient-expansion
(2)
coefficient is restored, i.e. µPBEsol = 16π(π/3)2/3 Cx = 10/81 ≈ 0.1235, and the second-order
26
correlation gradient-expansion coefficient βPBEsol = 0.046 is found by fitting to jellium surface
exchange-correlation energies.
B97-GGA exchange-correlation functional

The Becke 97 GGA (B97-GGA) exchange-correlation functional is the GGA part of the B97
hybrid functional [108] (see Section 4.1). The B97-GGA exchange energy density is
X
eB97-GGA
x (ρ↑ , ρ↓ , ∇ρ↑ , ∇ρ↓ ) = eUEG
x,σ (ρσ ) gx (xσ ), (3.37)
σ∈{↑,↓}
where eUEG UEG (ρ , 0) is the spin-σ contribution to the UEG exchange energy density
x,σ (ρσ ) = ex σ
4/3
and the gradient correction gx (xσ ) is a function of xσ = |∇ρσ |/ρσ ,
m
X
gx (xσ ) = cx,i ux (xσ )i , (3.38)
i=0
with ux (xσ ) = γx x2σ /(1 + γx x2σ ). The B97-GGA correlation energy density is written as the sum
of opposite- and same-spin contributions
X
eB97-GGA
c (ρ↑ , ρ↓ , ∇ρ↑ , ∇ρ↓ ) = eB97-GGA
c,↑↓ (ρ↑ , ρ↓ , ∇ρ↑ , ∇ρ↓ ) + eB97-GGA
c,σσ (ρσ , ∇ρσ ), (3.39)
σ∈{↑,↓}
where
eB97-GGA
c,↑↓ (ρ↑ , ρ↓ , ∇ρ↑ , ∇ρ↓ ) = eUEG
c,↑↓ (ρ↑ , ρ↓ ) gc,↑↓ (x↑↓ ) (3.40)
and
eB97-GGA
c,σσ (ρσ , ∇ρσ ) = eUEG
c,σσ (ρσ ) gc,σσ (xσ ). (3.41)
In these expressions, eUEG UEG (ρ , ρ ) − eUEG (ρ , 0) − eUEG (ρ , 0) and eUEG (ρ ) =

c,↑↓ (ρ↑ , ρ↓ ) = ec ↑ ↓ c ↑ c ↓ c,σσ α
eUEG
c (ρσ , 0) are estimations of the opposite- and same-spin contributions to the UEG correlation
energy
q density [109, 110]. The opposite-spin gradient correction is taken as a function of x↑↓ =
(x2↑ + x2↓ )/2,
m
c↑↓
X
↑↓ i
gc,↑↓ (x↑↓ ) = c,i uc (x↑↓ ) , (3.42)
i=0
with u↑↓ ↑↓ 2 ↑↓ 2
c (x↑↓ ) = γc x↑↓ /(1 + γc x↑↓ ), and the same-spin gradient correction is
m
X
gc,σσ (xσ ) = cσσ σσ i
c,i uc (xσ ) , (3.43)
i=0
↑↓
with uσσ σσ 2 σσ 2 σσ
c (xσ ) = γc xσ /(1 + γc xσ ). The parameters γx = 0.004, γc = 0.006, and γc = 0.2,
were roughly optimized on atomic exchange and correlation energies. The other parameters cx,i ,
c↑↓ σσ
c,i , cc,i for a polynomial degree m = 2 in Eqs. (3.38), (3.42), and (3.43) were optimized in the
B97 hybrid functional in the presence of a fraction of HF exchange energy (see Section 4.1).
The Hamprecht–Cohen–Tozer–Handy (HCTC) [111] exchange-correlation functional uses the
same form as the B97-GGA exchange-correlation functional but with a polynomial degree m =
4 and the parameters cx,i , c↑↓ σσ
c,i , cc,i were optimized without HF exchange a set of energetic
properties (atomic total energies, ionization energies, atomization energies), nuclear gradients,
and accurate exchange-correlation potentials.
27
3.4 Meta-generalized-gradient approximations
The meta-generalized-gradient approximations (meta-GGAs or mGGAs) are of the generic
form, in their spin-independent version,
Z
mGGA
Exc [ρ, τ ] = emGGA
xc (ρ(r), ∇ρ(r), ∇2 ρ(r), τ (r))dr, (3.44)
R3
i.e., they use more ingredients than the GGAs, namely the Laplacian of the density ∇2 ρ(r)
and/or the non-interacting positive kinetic energy density τ (r) associated with a single-determinant
wave function Φ,
N
Z
τ (r) = τΦ (r) = |∇r Φ(x, x2 , ..., xN )|2 dσdx2 ...dxN
2 {↑,↓}×(R3 ×{↑,↓})N−1
N
1X
= |∇ϕi (r)|2 , (3.45)
2
i=1
where {ϕi }i=1,...,N are the orbitals occupied in Φ. The meta-GGAs are considered as part of the
family of semilocal approximations, in the sense that τ (r) contains semilocal information with
respect to the orbitals.
Meta-GGAs can be viewed as implicit functionals of the density only, i.e. Exc mGGA [ρ, τ
Φ[ρ] ],
since τ (r) can be considered itself as an implicit functional of the density via the KS single-
determinant wave function Φ[ρ]. This view in which Exc mGGA [ρ] is a proper approximation
to the exchange-correlation density functional Exc [ρ] of the KS scheme is normally adopted
when constructing meta-GGAs approximations. However, the calculation of the functional
derivative of ExcmGGA [ρ] with respect to the density then requires the use of the complicated
optimized-effective-potential method (see Section 7). Therefore, in practical calculations, meta-

GGAs are usually reinterpreted as explicit functionals of a single-determinant wave function Φ,
mGGA [ρ , τ ], [112–119] or, in other words, approximations to an exact GKS exchange-
i.e. Exc Φ Φ
correlation functional (see Section 1.4).
In the latter approach, which we will here refer to as the meta-Kohn–Sham (mKS) scheme,
we introduce a functional Exc mKS [ρ, τ ] (to which meta-GGAs are approximations) defined for ρ
and τ simultaneously representable by a single-determinant wave function Φ ∈ S N and which

defines the GKS functional Exc S [Φ] = E mKS [ρ , τ ] [see Eq. (1.47)] giving the exact ground-state
xc Φ Φ
energy via Eq. (1.48),
n o
mKS
E0 = inf hΦ|T̂ + V̂ne |Φi + EH [ρΦ ] + Exc [ρΦ , τΦ ] , (3.46)
Φ∈S N
which, by taking variations with respect to the orbitals, gives the mKS equations:

1 2 mKS
− ∇ + vne (r) + vH (r) + vxc (r) ϕi (r) = εi ϕi (r). (3.47)
2
mKS (r) = v mKS (r) + v mKS (r) contains a usual local potential
Here, vxc xc,1 xc,2
mKS [ρ, τ ]
δExc
mKS
vxc,1 (r) = , (3.48)
δρ(r)
and a non-multiplicative operator [114–118]

mKS [ρ, τ ]

mKS 1 δExc
vxc,2 (r) =− ∇· ∇ , (3.49)
2 δτ (r)
28
evaluated with ρ(r) = N
P 2
PN 2
i=1 |ϕi (r)| and τ (r) = (1/2) i=1 |∇ϕi (r)| . Interestingly, the mKS
equations can be rewritten as a Schrödinger-like equation with a position-dependent mass
m(r) [120],
1 1 mKS
− ∇· ∇ + vne (r) + vH (r) + vxc,1 (r) ϕi (r) = εi ϕi (r), (3.50)
2 m(r)
mKS [ρ, τ ]/δτ (r) −1 . As in the KS scheme, the functional E mKS [ρ, τ ] is

where m(r) = 1 + δExc xc
decomposed into exchange and correlation contributions: Exc mKS [ρ, τ ] = E mKS [ρ, τ ] + E mKS [ρ, τ ].
x c
In the spin-dependent version of the mKS scheme, we consider a similar functional of the spin-
resolved densities and non-interacting positive kinetic energy densities Exc mKS [ρ , ρ , τ , τ ] and
↑ ↓ ↑ ↓
the spin-scaling relation of Eq. (1.45) is generalized to
1
ExmKS [ρ↑ , ρ↓ , τ↑ , τ↓ ] = ExmKS [2ρ↑ , 2τ↑ ] + ExmKS [2ρ↓ , 2τ↓ ] .

(3.51)
2
Correspondingly, the spin-dependent versions of the meta-GGAs are formulated in terms of the
spin-resolved quantities ρ↑ , ρ↓ , ∇ρ↑ , ∇ρ↓ , ∇2 ρ↑ , ∇2 ρ↓ , τ↑ , and τ↓ .
One motivation for the introduction of the variable τ (r) is that it appears in the expansion of
the spherically averaged exchange hole [entering in Eq. (2.16)] for small interelectronic distances
r12 [121], which for the case of a closed-shell system is
1 ρ(r1 )
Z
2 hx (r1 , r1 + r12 )dr12 = −
4πr12 S(0,r12 ) 2
|∇ρ(r1 )|2

1 1 2 2 4
− ∇ ρ(r1 ) − τ (r1 ) + r12 + O(r12 ), (3.52)
3 4 8ρ(r1 )
where S(0, r12 ) designates the sphere centered at 0 and of radius r12 = |r12 |. Thus τ (r) is needed
to describe the curvature of the exchange hole.
Another important motivation is that τ (r) is useful for identifying different types of spatial
regions of electronic systems [122]. This is done by comparing τ (r) with the von Weizsäcker
kinetic energy density,
|∇ρ(r)|2
τ W (r) = , (3.53)
8ρ(r)
which is the exact non-interacting kinetic energy density for one-electron systems and two-
electron spin-unpolarized systems, and, more generally, for one-orbital regions as introduced in
Section 2.3. For example, the indicator
τ W (r)
z(r) = , (3.54)
τ (r)
which takes its values in the range [0, 1] [123], identifies one-orbital regions (z = 1). A better
indicator is
τ (r) − τ W (r)
α(r) = , (3.55)
τ UEG (r)
where τ UEG (r) = (3/10)(3π 2 )2/3 ρ(r)5/3 is the non-interacting kinetic energy density of the UEG.
This indicator α(r) distinguishes one-orbital regions (α = 0), slowly varying density regions
(α ≈ 1), and regions of density overlap between closed shells that characterize noncovalent
bonds (α ≫ 1).
Nowadays, ∇2 ρ(r) is rarely used to construct meta-GGAs because it contains similar infor-
mation to τ (r), which can be seen by the second-order gradient expansion of τ (r) [124]:
1 |∇ρ(r)|2 1 2
τ GEA2 (r) = τ UEG (r) + + ∇ ρ(r). (3.56)
72 ρ(r) 6
29
In comparison to GGAs, meta-GGAs are more versatile and generally constitute an improve-
ment. Significantly, thanks to the use of τ , self-interaction errors in the correlation functional can
be essentially eliminated with meta-GGAs. They still suffer however from self-interaction errors
in the exchange functional. We now describe some of the most used meta-GGA functionals.
TPSS exchange-correlation functional

In the Tao–Perdew–Staroverov–Scuseria (TPSS) [125, 126] functional, the exchange energy
density is written as
eTPSS
x (ρ, ∇ρ, τ ) = eUEG
x (ρ)FxTPSS (s, z), (3.57)
where the enhancement factor is a function of s = |∇ρ|/(2kF ρ) and z = τ W /τ ,

κ
FxTPSS (s, z) = 1 + κ − , (3.58)
1+ xTPSS (s, z)/κ
with κ = 0.804 so as to saturate the local Lieb–Oxford bound (just like in the PBE exchange
functional) and
" s
2 1 3 2 1 4

10 z 146 73
xTPSS (s, z) = +c s 2
+ q̃ 2
b − q̃ b z + s
81 (1 + z 2 )2 2025 405 2 5 2
2 2 #
1 10 √ 10 3 √ 2
+ s4 + 2 e z + eµs6 / 1 + es2 , (3.59)
κ 81 81 5
and q̃b = (9/20)(α − 1)/[1 + bα(α − 1)]1/2 + 2s2 /3 (where α = (τ − τ W )/τ UEG = (5s2 /3)(z −1 − 1))
is a quantity that tends to the reduced density Laplacian q = ∇2 ρ/(4kF2 ρ) in the slowly varying
density limit [using Eq. (3.56)]. The function xTPSS (s, z) is chosen so as to satisfy the fourth-
order gradient expansion [Eq. (3.18)] which can be written in the form of the enhancement factor
FxGEA4 (s, z) = 1+(10/81)s2 +(146/2025)q 2 −(73/405)s2 q. The constant µ = 0.21951 is chosen to
retain the same large-s behavior of the PBE exchange functional, i.e. FxTPSS (s, z) ∼ FxPBE (s).
s→∞
The constants c = 1.59096 and e = 1.537 are chosen so as to eliminate the divergence of the
potential at the nucleus for a two-electron exponential density and to yield the correct exchange
energy (−0.3125 hartree) for the exact ground-state density of the hydrogen atom. Finally, the
constant b = 0.40 is chosen, quite arbitrarily, as the smallest value that makes FxTPSS (s, z) a
monotonically increasing function of s.
The TPSS correlation functional is constructed by making minor refinements to the previ-
ously developed Perdew–Kurth–Zupan–Blaha (PKZB) [127] meta-GGA correlation functional,
eTPSS
c (ρ↑ , ρ↓ , ∇ρ↑ , ∇ρ↓ , τ↑ , τ↓ ) = ρ εrevPKZB (ρ↑ , ρ↓ , ∇ρ↑ , ∇ρ↓ , τ↑ , τ↓ )
c revPKZB
(ρ↑ , ρ↓ , ∇ρ↑ , ∇ρ↓ , τ↑ , τ↓ )z 3 , (3.60)

× 1 + d εc
where the revised PKZB correlation energy per particle is
εrevPKZB (ρ↑ , ρ↓ , ∇ρ↑ , ∇ρ↓ , τ↑ , τ↓ ) = εPBE (ρ↑ , ρ↓ , ∇ρ↑ , ∇ρ↓ ) 1 + C(ζ, ξ)z 2

c c
X ρσ
− [1 + C(ζ, ξ)] z 2 ε̃PBE (ρ↑ , ρ↓ , ∇ρ↑ , ∇ρ↓ ), (3.61)
ρ c,σ
σ∈{↑,↓}
with ε̃PBE PBE (ρ , 0, ∇ρ , 0), εPBE (ρ , ρ , ∇ρ , ∇ρ )]

c,σ (ρ↑ , ρ↓ , ∇ρ↑ , ∇ρ↓ ) = max[εc σ σ c ↑ ↓ ↑ ↓
PBE
where εc (ρ↑ , ρ↓ , ∇ρ↑ , ∇ρ↓ ) is the PBE correlation energy per particle. Equation (3.61) con-
stitutes a one-electron self-interaction correction on the PBE correlation functional. Indeed, for
30
one-electron densities we have z = 1 and ζ = ±1, and the TPSS correlation energy correctly
vanishes [Eqs. (2.39)]. The TPSS correlation functional preserves many properties of the PBE
correlation functional: it has correct uniform coordinate scaling in the high- and low-density lim-
its, vanishing correlation energy in the large density-gradient limit, and the same second-order
gradient expansion (since the additional terms beyond PBE are at least in z 2 and thus only
change the fourth-order terms of the gradient expansion). The parameters d = 2.8 hartree−1
and C(0, 0) = 0.53 are chosen so as to recover the PBE surface correlation energy of jellium [128]
over the range of valence-electron bulk densities. The rest of the function is taken as
0.53 + 0.87ζ 2 + 0.50ζ 4 + 2.26ζ 6

C(ζ, ξ) = , (3.62)
(1 + ξ 2 [(1 + ζ)−4/3 + (1 − ζ)−4/3 ]/2)4
where ξ = |∇ζ|/(2kF ) is a reduced spin-polarization gradient. The function C(ζ, ξ) is chosen so

as to make the exchange-correlation energy independent of the spin polarization ζ in the low-
density limit [Eq. (2.63)] and to avoid that the self-interaction correction introduces additional
correlation energy density in the core-valence overlap region of monovalent atoms such as Li.
M06-L exchange-correlation functional

In the Minnesota 06 local (M06-L) exchange-correlation functional [129], the exchange energy
density is written as
X
eM06-L
x (ρ↑ , ρ↓ , ∇ρ↑ , ∇ρ↓ , τ↑ , τ↓ ) = ePBE UEG
x,σ (ρσ , ∇ρσ )f (wσ ) + ex,σ (ρσ )hx (xσ , Zσ ). (3.63)
σ∈{↑,↓}
The first term in Eq. (3.63), which has the same form as in the previously developed M05
exchange functional [130], contains the spin-σ PBE exchange energy density ePBE
x,σ (ρσ , ∇ρσ ) =
PBE
ex (ρσ , 0, ∇ρσ , 0) and the kinetic-energy density correction factor
11
X
f (wσ ) = ai wσi , (3.64)
i=0
5/3
where wσ = (τσUEG /τσ − 1)/(τσUEG /τσ + 1) with τσUEG = (3/10)(6π 2 )2/3 ρσ is an indicator of
the delocalization of the exchange hole [131]. The second term in Eq. (3.63), which has the
same form as in the VS98 exchange functional [132], contains the spin-σ UEG exchange energy
density eUEG UEG (ρ , 0) and the correction factor
x,σ (ρσ ) = ex σ
hx (xσ , Zσ ) = h(xσ , Zσ , dx,0 , dx,1 , dx,2 , dx,3 , dx,4 , αx ), (3.65)

4/3 5/3
where xσ = |∇ρσ |/ρσ and Zσ = 2(τσ − τσUEG )/ρσ and h is the parametrized function
d0 d1 x2 + d2 Z d3 x4 + d4 x2 Z
h(x, Z, d0 , d1 , d2 , d3 , d4 , α) = + + , (3.66)
γ(x, Z, α) γ(x, Z, α)2 γ(x, Z, α)3
with γ(x, Z, α) = 1 + α(x2 + Z).

The M06-L correlation energy is written as the sum of opposite- and same-spin contributions,
similarly to the B97-GGA correlation functional [Eq. (3.39)],
eM06-L
c (ρ↑ , ρ↓ , ∇ρ↑ , ∇ρ↓ , τ↑ , τ↓ ) = eM06-L
c,↑↓ (ρ↑ , ρ↓ , ∇ρ↑ , ∇ρ↓ , τ↑ , τ↓ )
X
+ eM06-L
c,σσ (ρσ , ∇ρσ , τσ ), (3.67)
σ∈{↑,↓}
31
where
eM06 UEG
c,↑↓ (ρ↑ , ρ↓ , ∇ρ↑ , ∇ρ↓ , τ↑ , τ↓ ) = ec,↑↓ (ρ↑ , ρ↓ ) [gc,↑↓ (x↑↓ ) + hc,↑↓ (x↑↓ , Z↑↓ )] , (3.68)
and
eM06-L UEG
c,σσ (ρσ , ∇ρσ ) = ec,σσ (ρσ ) [gc,σσ (xσ ) + hc,σσ (xσ , Zσ )] Dσ (zσ ), (3.69)
where the spin-decomposed UEG correlation energies eUEG UEG

c,↑↓ (ρ↑ , ρ↓ ) and ec,σσ (ρσ ) were already
defined after Eq. (3.41), and the gradient corrections gc,↑↓ (x↑↓ ) and gc,σσ (xσ ) are given in
Eqs. (3.42) and (3.43). The additional correction factors are
↑↓ ↑↓ ↑↓ ↑↓ ↑↓
h↑↓ ↑↓
c (x↑↓ , Z↑↓ ) = h(x↑↓ , Z↑↓ , dc,0 , dc,1 , dc,2 , dc,3 , dc,4 , αc ), (3.70)
q
where x↑↓ = (x2↑ + x2↓ )/2, Z↑↓ = Z↑ + Z↓ , and
hσσ σσ σσ σσ σσ σσ σσ
c (xσ , Zσ ) = h(xσ , Zσ , dc,0 , dc,1 , dc,2 , dc,3 , dc,4 , αc ). (3.71)
The factor Dσ (zσ ) = 1 − zσ , where zσ = τσW /τσ and τσW = |∇ρσ |2 /(8ρσ ), ensures that the
correlation energy correctly vanishes for one-electron systems [133].
The parameters γc↑↓ = 0.0031, and γcσσ = 0.06, were optimized on the correlation energies of
He and Ne. The parameters αx = 0.001867, α↑↓ σσ
c = 0.003050, and αc = 0.005151 were taken
↑↓ ↑↓
from Ref. [132]. The constraints a0 + dx,0 = 1, cc,0 + dc,0 = 1, and cc,0 + dσσ
σσ
c,0 = 1 are enforced to
obtain the correct UEG limit. The remaining 34 free parameters ai , c↑↓ σσ
c,i , cc,i for a polynomial de-
gree m = 4 in Eqs. (3.64), (3.42), and (3.43), and dx,i , d↑↓ σσ
c,i , dc,i in Eqs. (3.65), (3.70), and (3.71)
were optimized on a large set of diverse physicochemical properties concerning main-group ther-
mochemistry, reaction barrier heights, noncovalent interactions, electronic spectroscopy, and
transition metal bonding.
SCAN exchange-correlation functional

In the SCAN (strongly constrained and appropriately normed) [134] exchange-correlation
functional, the exchange energy density is written as
eSCAN
x (ρ, ∇ρ, τ ) = eUEG
x (ρ)FxSCAN (s, α), (3.72)
where the enhancement factor is a function of s = |∇ρ|/(2kF ρ) and α = (τ − τ W )/τ UEG ,
FxSCAN (s, α) = [h1x (s, α) + fx (α)(h0x − h1x (s, α))]gx (s), (3.73)
which interpolates between α = 0 and α ≈ 1, and extrapolates to α → ∞ using the function
fx (α) = exp[−c1x α/(1 − α)]θ(1 − α) − dx exp[c2x /(1 − α)]θ(α − 1), (3.74)
where θ is the Heaviside step function. The function gx (s) = 1−exp(−a1 s−1/2 ) is chosen to make
FxSCAN (s, α) vanish like s−1/2 as s → ∞, which guarantees the non-uniform scaling finiteness
conditions [Eqs. (2.66) and (2.68)] [57,77], and a1 = 4.9479 is taken to recover the exact exchange
energy of the hydrogen atom. For α ≈ 1 (slowly varying density regions), FxSCAN (s, α) ≈
h1x (s, α)gx (s) where h1x (s, α) is a PBE-like resummation of the fourth-order gradient expansion
[Eq. (3.18)],
k1
h1x (s, α) = 1 + k1 − , (3.75)
1 + xSCAN (s, α)/k1
32
where
2 /µ 2
xSCAN (s, α) = µs2 [1 + (b4 s2 /µ)e−|b4 |s ] + [b1 s2 + b2 (1 − α)e−b3 (1−α) ]2 , (3.76)
with µ = 10/81, b2 = (5913/405000)1/2 , b1 = (511/13500)/(2b2 ), b3 = 0.5, and b4 = µ2 /k1 −

1606/18225 − b21 . For α = 0 (one-orbital regions), FxSCAN (s, α = 0) = h0x gx (s) where h0x = 1.174
is chosen to saturate the local two-electron tight bound FxSCAN (s, α = 0) ≤ 1.174, which is a
sufficient and necessary condition for a meta-GGA exchange functional to satisfy the global tight
bound of Eq. (2.79) for all two-electron spin-unpolarized densities [77].
The SCAN correlation energy density is written as
eSCAN
c (ρ↑ , ρ↓ , ∇ρ↑ , ∇ρ↓ , τ↑ , τ↓ ) = ρ [ε1c (ρ↑ , ρ↓ , t) + fc (α)(ε0c (ρ↑ , ρ↓ , s) − ε1c (ρ↑ , ρ↓ , t))], (3.77)
which is again an interpolation between α = 0 and α = 1, and an extrapolation to α → ∞ using

the function
fc (α) = exp[−c1c α/(1 − α)]θ(1 − α) − dc exp[c2c /(1 − α)]θ(α − 1). (3.78)
For α = 1, the correlation energy par particle is taken as a revised version of the PBE correlation
energy per particle,
ε1c (ρ↑ , ρ↓ , t) = εUEG

c (ρ↑ , ρ↓ ) + H1SCAN (ρ↑ , ρ↓ , t) (3.79)
where
H1SCAN (ρ↑ , ρ↓ , t) = A(0)φ2 (ζ)3 ln 1 + w1 (1 − g(At2 )) ,

(3.80)
with t = |∇ρ|/(2φ2 (ζ)ks ρ), w1 = exp[−εUEGc (ρ↑ , ρ↓ )/(A(0)φ2 (ζ)3 )]−1, A = β(rs )/(A(0)w1 ), and
g(At ) = 1/(1+4At ) . The function has a second-order gradient expansion H1SCAN (ρ↑ , ρ↓ , t) =
2 2 1/4
β(rs )φ2 (ζ)t2 +O(t4 ) where the coefficient β(rs ) = 0.066725(1+0.1rs )/(1+0.1778rs ) is a rough fit
of the density dependence of the second-order gradient expansion correlation coefficient beyond
the Ma-Brueckner high-density-limit value and designed so that for rs → ∞ the second-order
gradient expansion terms for exchange and correlation cancel each other [135]. For α = 0, the
correlation energy par particle is constructed to be accurate for one- and two-electron systems
and is written as
ε0c (ρ↑ , ρ↓ , s) = [εLDA0

c (ρ) + H0SCAN (ρ, s)]Gc (ζ). (3.81)
The spin function Gc (ζ) = [1 − 2.3631(φ4 (ζ) − 1)](1 − ζ 12 ) is designed to make the correlation
energy vanish for one-electron densities (α = 0 and ζ = ±1) and to make the exchange-correlation
energy independent of ζ in the low-density limit [Eq. (2.63)]. Equation (3.81) includes a LDA-
type term [136]
b1c
εLDA0
c (ρ) = − 1/2
, (3.82)
1 + b2c rs + b3c rs
and a gradient correction
H0SCAN (ρ, s) = b1c ln [1 + w0 (1 − g∞ (ζ = 0, s))] , (3.83)
with w0 = exp(−εLDA0
c (ρ)/b1c )−1 and g∞ (ζ = 0, s) = limζ→0 limrs →∞ g(At) = 1/(1+0.512104s2 )1/4 .
The parameter b1c = 0.0285764 is determined so that the high-density limit of ε0c (ρ↑ , ρ↓ , s) repro-
duces the exact correlation energy of the Helium isoelectronic series in the large-nuclear charge
limit, i.e. limZ→∞ Ec [ρN =2,Z ] = EcGL2 [ρH
N =2,Z=1 ] = −0.0467 hartree [Eq. (2.72)]. The parameter
33
b3c = 0.125541 is determined to saturate the lower bound on the exchange-correlation energies
of two-electron densities [Eq. (2.77)]. The parameter b2c = 0.0889 is determined to reproduce
the exact exchange-correlation energy of the He atom.
The remaining seven parameters (k1 = 0.065, c1x = 0.667, c2x = 0.8, dx = 1.24, c1c = 0.64,
c2c = 1.5, and dc = 0.7) are determined by fitting to the approximate asymptotic expansions of
the exchange and correlation energies of neutral atoms in large nuclear charge limit [Eqs. (2.74)
and (2.75)], the binding energy curve of compressed Ar2 , and jellium surface exchange-correlation
energies.
4 Single-determinant hybrid approximations

4.1 Hybrid approximations
Based on arguments relying on the adiabatic-connection formalism, in 1993 Becke [137]
proposed to mix a fraction of the exact or Hartree-Fock (HF) exchange energy ExHF with GGA
functionals. In particular, he proposed a three-parameter hybrid (3H) approximation [138] of the
form, written here in its spin-independent version,
3H
Exc [Φ] = a ExHF [Φ] + b ExGGA [ρΦ ] + (1 − a − b) ExLDA [ρΦ ] + c EcGGA [ρΦ ] + (1 − c) EcLDA [ρΦ ], (4.1)
3H [Φ] is thought of as a functional of a
with empirical parameters a, b, and c. The functional Exc
N HF
single-determinant wave function Φ ∈ S since Ex [Φ] is itself a functional of Φ,
ExHF [Φ] = hΦ|Ŵee |Φi − EH [ρΦ ]

Nσ Z
Nσ X
1 X X ϕ∗iσ (r1 )ϕjσ (r1 )ϕ∗jσ (r2 )ϕiσ (r2 )
= − dr1 dr2 , (4.2)
2 R3 ×R3 |r1 − r2 |
σ∈{↑,↓} i=1 j=1
where {ϕiσ }i=1,...,Nσ are the orbitals occupied in Φ. In 1996, Becke proposed a simpler one-
parameter hybrid (1H) approximation [139],
1H
Exc [Φ] = a ExHF [Φ] + (1 − a) ExGGA [ρΦ ] + EcGGA [ρΦ ], (4.3)
where the fraction a of HF exchange has to be determined. For simplicity, we considered GGA
functionals ExGGA [ρΦ ] and EcGGA [ρΦ ] in Eq. (4.3) but we can more generally use meta-GGA
functionals ExmGGA [ρΦ , τΦ ] and EcmGGA [ρΦ , τΦ ].
These hybrid approximations should be considered as approximations of the GKS exchange-
S [Φ] in Eq. (1.47) with S[Φ] = a E HF [Φ]. The corresponding GKS
correlation functional Exc x
equations [Eq. (1.4)] then include the term
δS[Φ]
Z
HF
=a vx,σ (r, r′ )ϕiσ (r′ )dr′ , (4.4)
δϕ∗iσ (r) R3
HF (r, r′ ) is the nonlocal HF exchange potential

where vx,σ 7
Nσ
HF
X ϕjσ (r)ϕ∗jσ (r′ )
vx,σ (r, r′ ) =− . (4.5)
|r − r′ |
j=1
The main benefit of adding a fraction of HF exchange is to decrease the self-interaction

error (see Section 2.3) introduced by semilocal exchange functionals which tends to favor too
7
The possibility of combining a nonlocal HF potential with a local correlation potential was mentioned already
in 1965 in the paper by Kohn and Sham [19].
34
much delocalized electron densities over localized electron densities. The fraction of HF exchange
should however be small enough to keep the compensation of errors usually occurring between the
approximate semilocal exchange and correlation functionals. First, Becke used the value a = 0.5
in the so-called Becke Half-and-Half functional [137], but then fits to various experimental data
often repeatedly gave an optimal parameter a around 0.20-0.25. A rationale has been proposed
in favor of the value 0.25 [140]. By decreasing self-interaction errors in the exchange energy,
hybrid approximations are often a big improvement over semilocal approximations for molecular
systems with sufficiently large electronic gaps. However, for systems with small HOMO-LUMO
gaps, such as systems with stretched chemical bonds or with transition metal elements, they
tend to increase static-correlation errors.
An interesting extension of the hybrid approximations are the so-called local hybrids, which
use a position-dependent fraction a(r) of a (non-uniquely defined) HF exchange energy density
eHF
x (r) [141] (see, Ref. [142] for a recent review), and which belong to the wider family of hyper-
GGA functionals in which the correlation energy can also be expressed as a function of eHF
x (r) [6].
The local-hybrid approximations are much more flexible than the global hybrid approach exposed
in this Section but require more complicated and computationally expensive implementations.
For this reason, they have not been often used and we will not consider them any further here.
We now describe some of the most used hybrid approximations.
B3LYP exchange-correlation functional

The B3LYP exchange-correlation functional [143] is the most famous and widely used three-
parameter hybrid approximation [Eq. (4.1)]. It uses the B88 exchange functional and the LYP
correlation functional,
B3LYP
Exc [Φ] = a ExHF [Φ] + b ExB88 [ρ↑,Φ , ρ↓,Φ ] + (1 − a − b) ExLSDA [ρ↑,Φ , ρ↓,Φ ]
+c EcLYP [ρ↑,Φ , ρ↓,Φ ] + (1 − c) EcLSDA [ρ↑,Φ , ρ↓,Φ ], (4.6)
and the parameters a = 0.20, b = 0.72, and c = 0.81 were found by optimizing on a set of
atomization energies, ionization energies, proton affinities of small molecules and first-row total
atomic energies [138]. A caveat is that the VWN parametrization of the RPA correlation energy
(sometimes referred to as VWN3) of the UEG was actually used for EcLSDA [ρ↑ , ρ↓ ] instead of the
VWN parametrization of the accurate correlation energy (sometimes referred to as VWN5) of
the UEG [82].
B97 exchange-correlation functional

The Becke 97 (B97) exchange-correlation functional [108] is a GGA hybrid of the form
B97
Exc [Φ] = a ExHF [Φ] + (1 − a) ExB97-GGA [ρ↑,Φ , ρ↓,Φ ] + EcB97-GGA [ρ↑,Φ , ρ↓,Φ ], (4.7)
where the form of the B97-GGA exchange and correlation functionals were given in Eqs. (3.37)
and (3.39). The fraction of HF exchange a = 0.1943 and the remaining parameters cx,0 =
1.00459, cx,1 = 0.629639, cx,2 = 0.928509, c↑↓ ↑↓ ↑↓ σσ
c,0 = 0.9454, cc,1 = 0.7471, cc,2 = −4.5961, cc,0 =
σσ σσ
0.1737, cc,1 = 2.3487, and cc,2 = −2.4868 for a polynomial degree m = 2 in Eqs. (3.38), (3.42),
and (3.43) were optimized on a set of total energies, atomization energies, ionization energies, and
proton affinities. Note that, for xσ = 0, the UEG limit is not imposed, which would require the
parameters cx,0 , c↑↓ σσ
c,0 , and cc,0 to be all strictly equal to 1. With the above optimized parameters,
we see that it is nearly satisfied for the exchange energy and the opposite-spin correlation energy,
but very far from it for the same-spin correlation energy which is drastically reduced compared
to the LDA.
PBE0 exchange-correlation functional
35
The PBE0 exchange-correlation functional [144, 145] is a GGA hybrid using the PBE ex-
change and correlation functionals,
PBE0
Exc [Φ] = a ExHF [Φ] + (1 − a) ExPBE [ρ↑,Φ , ρ↓,Φ ] + EcPBE [ρ↑,Φ , ρ↓,Φ ], (4.8)
and the fraction of the HF exchange is fixed at a = 0.25 according to the rationale of Ref. [140].
This functional is also known under the name PBE1PBE. The “1” in the latter name emphasizes
that there is one parameter, a, while the “0” in the more common name PBE0 emphasizes that
this parameter is not found by fitting.
TPSSh exchange-correlation functional

The TPSSh exchange-correlation functional [146] is a meta-GGA hybrid using the TPSS
exchange and correlation functionals,
TPSSh
Exc [Φ] = a ExHF [Φ] + (1 − a) ExTPSS [ρ↑,Φ , ρ↓,Φ , τ↑,Φ , τ↓,Φ ] + EcTPSS [ρ↑,Φ , ρ↓,Φ , τ↑,Φ , τ↓,Φ ], (4.9)
and the fraction of the HF exchange a = 0.10 was determined by optimizing on a large set of
atomization energies.
M06 and M06-2X exchange-correlation functionals

The M06 exchange-correlation functional [147] is a meta-GGA hybrid using the M06-L ex-
change and correlation functionals,
M06
Exc [Φ] = a ExHF [Φ] + (1 − a) ExM06-L [ρ↑,Φ , ρ↓,Φ , τ↑,Φ , τ↓,Φ ] + EcM06-L [ρ↑,Φ , ρ↓,Φ , τ↑,Φ , τ↓,Φ ],(4.10)
and the parameters in the M06-L exchange and correlation functionals were reoptimized together
with the fraction of HF exchange a = 0.27 on the same large set of diverse physicochemical
properties used for the M06-L functional. In the M06-2X exchange-correlation functional the
fraction of HF exchange is doubled, i.e. a = 0.54, and the parameters were reoptimized with
the function hx (xσ , Zσ ) in Eq. (3.65) set to zero and excluding transition metal properties in
the training set. With this large fraction of HF exchange, the M06-2X functional is designed for
systems without transition metal elements.
4.2 Range-separated hybrid approximations

Based on earlier ideas of Savin [148] (exposed in details in Section 5.2), in 2001, Iikura,
Tsuneda, Yanai, and Hirao [149] proposed a long-range correction (LC) scheme in which the
exchange-correlation energy is written as, in its spin-independent version,
LC
Exc [Φ] = Exlr,µ,HF [Φ] + Exsr,µ,GGA [ρΦ ] + EcGGA [ρΦ ]. (4.11)
This scheme has also been referred to as the range-separated hybrid exchange (RSHX) scheme [150].
In Eq. (4.11), Exlr,µ,HF [Φ] is the HF exchange energy for a long-range electron-electron interac-
lr,µ
tion wee (r12 ) = erf(µr12 )/r12 (where erf is the error function and the parameter µ ∈ [0, +∞)
controls the range of the interaction),
Nσ Z
Nσ X
lr,µ,HF 1 X X
Ex [Φ] = − ϕ∗iσ (r1 )ϕjσ (r1 )ϕ∗jσ (r2 )ϕiσ (r2 )wee
lr,µ
(r12 )dr1 dr2 , (4.12)
2 R3 ×R3
σ∈{↑,↓} i=1 j=1
and Exsr,µ,GGA [ρ] is a GGA exchange energy functional for the complementary short-range inter-
sr,µ lr (r ). This latter functional can be thought of as an approximation
action wee (r12 ) = 1/r12 −wee 12
to the short-range exchange functional
1
Z
Exsr,µ [ρ] = sr,µ
ρ(r1 )hx (r1 , r2 )wee (r12 )dr1 dr2 , (4.13)
2 R3 ×R3
36
where hx (r1 , r2 ) is the KS exchange hole of Section 2.1. For µ = 0, the long-range HF exchange
energy vanishes, i.e. Exlr,µ=0,HF [Φ] = 0, and the short-range exchange functional reduces to the
standard KS exchange functional, i.e. Exsr,µ=0 [ρ] = Ex [ρ]. Reversely, for µ → ∞, the long-range
HF exchange energy reduces to the full-range HF exchange energy, i.e. Exlr,µ→∞,HF[Φ] = ExHF [Φ],
and the short-range exchange functional vanishes, i.e. Exsr,µ→∞ [ρ] = 0. Significantly, for large µ,
the short-range exchange functional becomes a local functional of the density [151, 152]:
π
Z
sr,µ
Ex [ρ] ∼ − 2 ρ(r)2 dr. (4.14)
µ→∞ 4µ R3
Like the hybrid approximations of Section 4.1, Eq. (4.11) should be considered as an approxi-
mation of the GKS exchange-correlation functional Exc S [Φ] in Eq. (1.47) with S[Φ] = E lr,µ,HF [Φ],
x
and the corresponding GKS equations [Eq. (1.4)] then includes a long-range nonlocal HF ex-
lr,µ,HF lr,µ
(r1 , r2 ) = − N ∗
P σ
change potential vx,σ j=1 ϕjσ (r1 )ϕjσ (r2 )wee (r12 ). Similarly to the hybrid
approximations, the introduction of a fraction of long-range HF exchange reduces the self-
interaction error (see, e.g., Ref. [153]). In addition, the short-range exchange part is easier
to approximate with semilocal density-functional approximations, as Eq. (4.14) strongly sug-
gests. In particular, the −1/r asymptotic behavior of the exchange potential [Eq. (1.33)], which
is difficult to satisfy with semilocal approximations, does not apply anymore to the short-range
exchange potential.
In 2004, Yanai, Tew, and Handy [154], introduced a more flexible scheme called the Coulomb-
attenuating method (CAM) [154] in which fractions of HF exchange are added at both short
range and long range,
CAM
Exc [Φ] = a Exsr,µ,HF [Φ]+b Exlr,µ,HF [Φ]+(1−a) Exsr,µ,GGA [ρΦ ]+(1−b) Exlr,µ,GGA [ρΦ ]+EcGGA [ρΦ ],
(4.15)
sr,µ,HF HF lr,µ,HF lr,µ,GGA
where Ex [Φ] = Ex [Φ]−Ex [Φ] is the short-range HF exchange energy and Ex =
ExGGA − Exsr,µ,GGA is a long-range GGA exchange energy. The reintroduction of HF exchange
at short range further reduces the self-interaction error and improves thermodynamic properties
such as atomization energies. Again, Eq. (4.15) should be considered as an approximation
of the GKS exchange-correlation functional Exc S [Φ] in Eq. (1.47) with S[Φ] = a E sr,µ,HF [Φ] +
x
b Exlr,µ,HF [Φ]. Other forms of modified electron-electron interactions are also possible (see, e.g.,
Refs. [152, 155, 156]).
The approximations in Eqs. (4.11) and (4.15) are usually collectively referred to as range-
separated hybrid approximations. Range-separated hybrids in the form of Eq. (4.15) are more
flexible than the hybrid approximations of Section 4.1, and consequently are potentially more
accurate, in particular for long-range electronic excitations. However, like the hybrid approxi-
mations, the presence of HF exchange tends to induce static-correlation errors for systems with
small HOMO-LUMO gaps.
The range-separation parameter µ (also sometimes denoted as ω) is generally chosen empir-
ically, e.g. by fitting to experimental data. In practice, a value around µ ≈ 0.3 − 0.5 bohr−1 ,
fixed for all systems, is often found to be optimal. It has also been proposed to adjust the
value of µ in each system, e.g. by requiring that the opposite of the HOMO energy be equal to
the ionization energy calculated by total energy differences [157–159]. These so-called optimally
tuned range-separated hybrids are well suited for the calculation of charge-transfer electronic
excitations but have the disadvantage of not being size consistent [160].
A natural idea is to use a position-dependent range-separation parameter µ(r) which allows
the range of the modified interaction to adapt to the local average electron-electron distance
in the diverse spatial regions of the system. These locally range-separated hybrids [161–163]
are promising but they induced computational complications and are still in the early stages of
development. We will thus not consider them any further here.
37
We now describe some of the most used approximations in the context of the range-separated
hybrids.
Short-range LDA exchange functional

The short-range LDA exchange functional [148, 151] can be obtained by using in Eq. (4.13)
the LDA exchange hole [Eq. (3.14)], which leads to
Z
sr,µ,LDA
Ex [ρ] = esr,µ,UEG
x (ρ(r))dr, (4.16)
R3
with the short-range UEG exchange energy density

8µ̃ √

sr,µ,UEG UEG 1 3 −1/(4µ̃2 ) 3
ex (ρ) = ex (ρ) 1 − π erf + (2µ̃ − 4µ̃ )e − 3µ̃ + 4µ̃ , (4.17)
3 2µ̃
where µ̃ = µ/(2kF ) is a dimensionless range-separation parameter. The spin-dependent version
is obtained from the same spin-scaling relation as in the standard case [Eq. (1.45)]. The short-
range LDA exchange functional becomes exact for large µ [Eq. (4.14)] and is the first building
block for constructing short-range exchange GGA functionals.
CAM-B3LYP exchange-correlation functional

The CAM-B3LYP exchange-correlation functional [154] uses Eq. (4.15) with short- and long-
range versions of the B88 exchange functional and the same correlation functional used in B3LYP
(i.e., 0.81 EcLYP + 0.19 EcLSDA ),
CAM-B3LYP
Exc [Φ] = a Exsr,µ,HF [Φ] + b Exlr,µ,HF [Φ] + (1 − a) Exsr,µ,B88 [ρ↑,Φ , ρ↓,Φ ]
+ (1 − b) Exlr,µ,B88 [ρ↑,Φ , ρ↓,Φ ] + 0.81EcLYP [ρ↑,Φ , ρ↓,Φ ] + 0.19EcLSDA [ρ↑,Φ , ρ↓,Φ ], (4.18)
where the parameters a = 0.19 and b = 0.65 were optimized on atomization energies and the
range-separation parameter µ = 0.33 bohr−1 was taken from Ref. [164] where it was optimized on
equilibrium distances of diatomics molecules. In this expression, the short-range B88 exchange
functional Exsr,µ,B88 is defined by using in Eq. (4.13) the following generic GGA model for the
exchange hole [149] (given here in its spin-independent version)
9 j1 (kGGA r12 ) 2

GGA
hx (ρ, ∇ρ, r12 ) = −ρ , (4.19)
2 kGGA r12
p
with kGGA = kF / eGGA x (ρ, ∇ρ)/eUEG
x (ρ). The exchange-hole model of Eq. (4.19) properly yields
the GGA exchange energy density eGGA x (ρ, ∇ρ) for µ = 0 and thus allows one to extend any
standard GGA exchange functional to a short-range GGA exchange functional. Note however
that it does not fulfill the sum rule [Eq. (2.11)]. The long-range B88 exchange functional is then
simply Exlr,µ,B88 = ExB88 − Exsr,µ,B88 .
LC-ωPBE exchange-correlation functional

The LC-ωPBE exchange-correlation functional [165, 166] uses a short-range version of the
PBE exchange functional as well as the standard PBE correlation functional,
LC-ωPBE
Exc [Φ] = Exlr,µ,HF [Φ] + Exsr,µ,PBE [ρ↑,Φ , ρ↓,Φ ] + EcPBE [ρ↑,Φ , ρ↓,Φ ]. (4.20)
The short-range PBE exchange functional is obtained by using in Eq. (4.13) the following GGA
exchange hole model constructed to yield the PBE exchange energy [167],
hPBE
x (ρ, ∇ρ, r12 ) = ρJ PBE (s, kF r12 ), (4.21)
38
where s = |∇ρ|/(2kF ρ) and
"
PBE A 1
J (s, u) = − 2 +
u 1 + (4/9)Au2
#
A 2 2 2 4 −Du2 2 2
+ B + C[1 + s F(s)]u + E[1 + s G(s)]u e e−s H(s)u . (4.22)
u2
Here, A, B, C, D, and E are constants chosen to obtain an oscillation-averaged UEG exchange
hole for s = 0, and F(s), G(s) and H(s) are functions determined so that the hole yields the PBE
exchange density for µ = 0, and satisfies the sum rule [Eq. (2.11)] and the small-r12 expansion
[Eq. (3.52)] using the gradient expansion of τ of Eq. (3.56). The range-separation parameter
is fixed at µ = ω = 0.4 bohr−1 which has been found to be close to optimal for atomization
energies, reaction barrier heights, and ionization energies [165].
ωB97X exchange-correlation functional

The ωB97X exchange-correlation functional [168] has the form of Eq. (4.15) with b = 1:
ωB97X
Exc [Φ] = a Exsr,µ,HF [Φ] + Exlr,µ,HF [Φ] + (1 − a) Exsr,µ,B97-GGA [ρ↑,Φ , ρ↓,Φ ]
+EcB97-GGA [ρ↑,Φ , ρ↓,Φ ]. (4.23)
The short-range B97-GGA exchange density is defined as
X
esr,µ,B97-GGA
x (ρ↑ , ρ↓ , ∇ρ↑ , ∇ρ↓ ) = esr,µ,UEG
x,σ (ρσ ) gx (xσ ), (4.24)
σ∈{↑,↓}
sr,µ,UEG
where ex,σ (ρσ ) = esr,µ,UEG
x (ρσ , 0) is the spin-σ contribution to the short-range UEG ex-
4/3
change energy density [Eq. (4.17)] and the gradient correction gx (xσ ) where xσ = |∇ρσ |/ρσ
has the same form as in Eq. (3.38) with polynomial degree m = 4. In Eq. (4.23), the correlation
functional has the same form as the B97-GGA correlation functional but again with polynomial
degree m = 4 in Eqs. (3.42) and (3.43). The fraction of short-range HF exchange a ≈ 0.16, the
range-separation parameter µ = ω = 0.3 bohr−1 , and the linear coefficients in Eqs. (3.38), (3.42),
and (3.43) were optimized on sets of atomic energies, atomization energies, ionization energies,
electron and proton affinities, reaction barrier heights, and noncovalent interactions, with the
constraints a + cx,0 = 1, c↑↓ σσ
c,0 = 1, and cc,0 = 1 to enforce the correct UEG limit.
HSE exchange-correlation functional

The Heyd–Scuseria–Ernzerhof (HSE) exchange-correlation functional [169] is of the form of
Eq. (4.15) with b = 0 (i.e., no long-range HF exchange),
HSE
Exc [Φ] = aExsr,µ,HF [Φ] + (1 − a)Exsr,µ,PBE [ρ↑,Φ , ρ↓,Φ ] + Exlr,µ,PBE [ρ↑,Φ , ρ↓,Φ ]
+EcPBE [ρ↑,Φ , ρ↓,Φ ], (4.25)
and involves the long-range PBE exchange functional Exlr,µ,PBE = ExPBE − Exsr,µ,PBE complemen-
tary to the short-range PBE exchange functional constructed from the PBE exchange hole model
[Eqs. (4.21) and (4.22)]. In order to reproduce reliable values for the band gap in semiconducting
solids, the range-separation parameter is fixed at µ = 0.15 bohr−1 , which is a very small value
compared to the other range-separated hybrids. It means that the range of electron-electron
distances covered by HF exchange is large, and the HSE functional could be thought of as a
regular hybrid approximation but with the very long-range contribution of the HF exchange
removed. This is particularly appropriate for solids since in these systems the very long-range
HF exchange is effectively balanced by the correlation effects (a phenomenon known as screen-
ing). The fraction of (short-range) HF exchange is fixed at a = 0.25 like in the PBE0 hybrid
functional.
39
5 Multideterminant hybrid approximations
5.1 Double-hybrid approximations
In 2006, Grimme [170] introduced a two-parameter double-hybrid (2DH) approximation,
written here in its spin-independent version,
2DH
Exc = ax ExHF [Φ] + (1 − ax ) ExGGA [ρΦ ] + (1 − ac )EcGGA [ρΦ ] + ac EcMP2 , (5.1)
mixing a fraction ax of the HF exchange energy with a GGA exchange functional, and a fraction
ac of the second-order Møller–Plesset (MP2) correlation energy EcMP2 with a GGA correlation
functional. In Eq. (5.1), the first three terms are first calculated in a self-consistent manner, and
then the last term EcMP2 is added perturbatively using the orbitals determined in the first step.
The expression of EcMP2 is [171]
N N
1 XX X X |hφi φj ||φa φb i|2
EcMP2 =− , (5.2)
4 εa + εb − εi − εj
i=1 j=1 a≥N +1 b≥N +1
where i, j and a, b run over occupied and virtual spin orbitals, respectively, εk are spin orbital en-
ergies, and hφi φj ||φa φb i = hφi φj |φa φb i − hφi φj |φb φa i are antisymmetrized two-electron integrals
with (in physicists’ notation)
φ∗p (x1 )φ∗q (x2 )φr (x1 )φs (x2 )
Z
hφp φq |φr φs i = dx1 dx2 . (5.3)
(R3 ×{↑,↓})2 |r1 − r2 |
Note that the notation in Eq. (5.2) assumes that the one-electron wave-function space is spanned
by a discrete set of spin orbitals. In the exact theory, the continuum limit of the set of virtual
spin orbitals is implied.
The rigorous framework underlying these double-hybrid approximations was established by
Sharkas, Toulouse, and Savin [172]. The idea is to decompose the universal density functional
of Eq. (1.8) as
λ
F [ρ] = min hΨ|T̂ + λŴee |Ψi + ĒHxc [ρ], (5.4)
Ψ∈WρN
λ [ρ] is a complementary density functional defined

where λ ∈ [0, 1] is a coupling constant and ĒHxc
λ [ρ] = E
to make Eq. (5.4) exact. From Eqs. (1.13) and (2.24), we see that ĒHxc λ
Hxc [ρ] − EHxc [ρ]
where EHxc [ρ] is the standard Hartree-exchange-correlation functional of the KS scheme and
λ [ρ] is the Hartree-exchange-correlation functional along the adiabatic connection. The
EHxc
Hartree and exchange contributions are simply linear in λ,
λ
ĒH [ρ] = (1 − λ)EH [ρ], (5.5)
Ēxλ [ρ] = (1 − λ)Ex [ρ], (5.6)

where EH [ρ] and Ex [ρ] are the standard Hartree and exchange functionals of the KS scheme.
Moreover, from the uniform coordinate scaling relation of Eq. (2.61), we have
Ēcλ [ρ] = Ec [ρ] − λ2 Ec [ρ1/λ ], (5.7)
where Ec [ρ] is the standard correlation functional of the KS scheme and ρ1/λ (r) = (1/λ)3 ρ(r/λ)
is the scaled density. The decomposition in Eq. (5.4) leads to the following expression of the
exact ground-state energy
n o
λ
E0 = inf hΨ|T̂ + V̂ne + λŴee |Ψi + ĒHxc [ρΨ ] , (5.8)
Ψ∈W N
40
where the infimum is over general multideterminant wave functions Ψ ∈ W N . This constitutes
a multideterminant extension of the KS scheme. Note that this multideterminant KS scheme
can trivially be extended to spin-dependent density functionals and functionals depending on
the kinetic-energy density [119].
The double-hybrid ansatz can be seen as a particular approximation within this multideter-
minant KS scheme [172]. To see this, we define a density-scaled one-parameter hybrid (DS1H)
approximation by restricting the minimization in Eq. (5.8) to single-determinant wave functions
Φ ∈ SN ,
n o
E0DS1H,λ = inf hΦ|T̂ + V̂ne + λŴee |Φi + ĒHxc λ
[ρΦ ] , (5.9)
Φ∈S N
obtaining an energy which necessarily depends on λ. A minimizing single-determinant wave

function Φλ must satisfy the self-consistent eigenvalue equation

T̂ + V̂ne + λV̂Hx [Φ ] + V̄ˆHxc
HF λ λ
[ρΦλ ] |Φλ i = E0λ |Φλ i, (5.10)
HF [Φλ ] is the nonlocal HF potential operator evaluated with the DS1H wave function
where V̂Hx
Φλ and V̄ˆHxc
λ [ρ ] is the local Hartree-exchange-correlation potential operator generated by the
Φλ
λ [ρ] and evaluated at the DS1H density ρ . If written explicitly in terms
energy functional ĒHxc Φλ
of spin orbitals, Eq. (5.10) would have the form of the GKS equations [Eq. (1.49)]. The DS1H
ground-state energy can be finally written as
E0DS1H,λ = hΦλ |T̂ + V̂ne |Φλ i + EH [ρΦλ ] + λExHF [Φλ ] + (1 − λ)Ex [ρΦλ ] + Ēcλ [ρΦλ ], (5.11)
where the full Hartree energy EH [ρ] has been recomposed. The exchange-correlation energy in
Eq. (5.11) is of the form of a hybrid approximation [Eq. (4.3)].
All that is missing in Eq. (5.11) is the correlation energy associated with the scaled interaction
λŴee . It can be calculated by a nonlinear Rayleigh-Schrödinger perturbation theory [173–175]
starting from the DS1H reference. Consider the following energy expression with the perturba-
tion parameter α ∈ [0, 1],
n o
E0λ,α = inf hΨ|T̂ + V̂ne + λV̂Hx HF λ λ
[Φ ] + αλŴ |Ψi + ĒHxc [ρΨ ] , (5.12)
Ψ∈W N

where λŴ = λ Ŵee − V̂Hx HF [Φλ ] is the scaled Møller–Plesset perturbation operator. For α = 0,
the stationary equation associated with Eq. (5.12) reduces to the DS1H eigenvalue equation
[Eq. (5.10)]. For α = 1, Eq. (5.12) reduces to Eq. (5.8), so E0λ,α=1 is the exact energy, indepen-
dently of λ. The sum of the zeroth-order energy and first-order energy correction gives simply
λ,(0) λ,(1)
the DS1H energy, E0DS1H,λ = E0 + E0 . Thanks to the existence of a Brillouin theorem
just like in standard Møller–Plesset perturbation theory (see Refs. [173–175]), only double exci-
tations contribute to the first-order wave-function correction Ψλ,(1) and the second-order energy
correction has a standard MP2 form
λ,(2)
E0 = λ2 hΦλ |Ŵ |Ψλ,(1) i = λ2 EcMP2 , (5.13)
where EcMP2 has the expression in Eq. (5.2) with DS1H spin orbitals and associated orbital
eigenvalues (which implicitly depend on λ). This second-order perturbation theory defines a
density-scaled one-parameter double-hybrid (DS1DH) approximation
λ,(2)
E0DS1DH,λ = E0DS1H,λ + E0 , (5.14)
41
which contains the exchange-correlation energy contribution
DS1DH,λ
Exc = λExHF [Φλ ] + (1 − λ)Ex [ρΦλ ] + Ēcλ [ρΦλ ] + λ2 EcMP2 . (5.15)
To make connection with the double-hybrid ansatz of Eq. (5.1), we can also define a one-
parameter double-hybrid (1DH) approximation, obtained by neglecting the density scaling in
the correlation functional, i.e. Ec [ρ1/λ ] ≈ Ec [ρ] in Eq. (5.7),
1DH,λ
Exc = λExHF [Φλ ] + (1 − λ)Ex [ρΦλ ] + (1 − λ2 )Ec [ρΦλ ] + λ2 EcMP2 , (5.16)
which, after using semilocal approximations for Ex [ρ] and Ec [ρ], has the form of Eq. (5.1) with
parameters ax = λ and ac = λ2 . In this rigorous formulation of the double-hybrid approxima-
tions, the fraction of HF exchange is thus connected to the fraction of MP2 correlation. Taking
into account approximately the scaling of the density in Ec [ρ1/λ ], it has also been proposed to use
ac = λ3 [176]. Fromager [177] also proposed an extension of this rigorous formulation in order
to justify the use of double-hybrid approximations with two parameters such that ac ≤ a2x = λ2 .
An essential advantage of double-hybrid approximations is that the presence of nonlocal MP2
correlation allows one to use a larger fraction of nonlocal HF exchange, which helps decreasing
the self-interaction error. This usually provides an improvement over hybrid approximations
for molecular systems with sufficiently large electronic gaps. However, a large fraction of HF
exchange and a fraction of MP2 correlation also generally means large static-correlation errors
in systems with small HOMO-LUMO gaps.
The first and still best known double-hybrid approximation is B2PLYP [170] which is based
on the B88 exchange functional and the LYP correlation functional,
B2PLYP
Exc = ax ExHF [Φ] + (1 − ax ) ExB88 [ρ↑,Φ , ρ↓,Φ ] + (1 − ac )EcLYP [ρ↑,Φ , ρ↓,Φ ] + ac EcMP2 , (5.17)
and the parameters ax = 0.53 and ac = 0.27 have been optimized on a set of atomization energies.
Interestingly, even though the two parameters have been optimized without any constraint, we
have ac ≈ a2x = 0.28 as predicted by Eq. (5.16).
It has also been proposed to use the spin-component-scaled (SCS) version of MP2 [178] to
construct spin-component-scaled double-hybrid approximations of the form [179, 180]
SCS-DH
Exc = ax ExHF [Φ] + (1 − ax ) ExGGA [ρΦ ] + (1 − ac )EcGGA [ρΦ ] + cOS Ec,OS
MP2 MP2
+ cSS Ec,SS ,
(5.18)
MP2 and
which contains four empirical parameters ax , ac , cOS , and cSS . In this expression, Ec,OS
MP2 are the opposite-spin (OS) and same-spin (SS) contributions to the MP2 correlation energy
Ec,SS
obtained by restricting the sums over i and j in Eq. (5.2) to spin orbitals of opposite and same
spins, respectively. Since in MP2 the same-spin component is usually overestimated relative
to the opposite-spin component, this SCS variant is a simple way to achieve higher accuracy
without increasing computational cost.
For reviews on different flavors of double hybrids and their assessments, the reader may
consult Refs. [181–184]. It has also been proposed to construct double-hybrid approxima-
tions where the MP2 correlation term is extended to a higher-order correlation method such
as RPA [185–189] or coupled-cluster [190, 191]. More generally, the multideterminant extension
of the KS scheme of Eq. (5.8) allows one to define hybrids combining any wave-function method
with density functionals. For example, a multiconfiguration hybrid approximation based on
Eq. (5.8) which combines a multiconfiguration self-consistent-field (MCSCF) wave function with
density functionals has been proposed in the goal of tackling strongly correlated systems [192].
This approach has also been used to combine valence-bond (VB) theory [193] or variational
two-electron reduced-density-matrix theory [194] with DFT.
42
5.2 Range-separated double-hybrid approximations
5.2.1 Range-separated one-parameter double-hybrid approximations
In 2005, Ángyán, Gerber, Savin, and Toulouse [173] introduced what could be called the first
range-separated one-parameter double-hybrid approximation, i.e. combining HF exchange and
MP2 correlation with density functionals using a one-parameter decomposition of the electron-
electron interaction. This is based on the range-separated multideterminant extension of the
KS scheme introduced earlier by Savin [148] (see, also, Refs. [152, 155, 195]) and which actually
predates and inspired the multideterminant extension of the KS scheme of Eq. (5.8).
The idea is to decompose the universal density functional of Eq. (1.8) as
lr,µ sr,µ
F [ρ] = min hΨ|T̂ + Ŵee |Ψi + ĒHxc [ρ], (5.19)
Ψ∈WρN
lr,µ
where Ŵee is the long-range electron-electron interaction operator (associated with the pair po-
lr,µ
tential wee (r12 ) = erf(µr12 )/r12 as already used in the range-separated hybrids of Section 4.2)
sr,µ
and ĒHxc [ρ] is the complementary short-range density functional defined to make Eq. (5.19)
exact. As before, the parameter µ ∈ [0, +∞) controls the range of the separation. The comple-
sr,µ lr,µ
mentary short-range functional can be written as ĒHxc [ρ] = EHxc [ρ] − EHxc [ρ], where EHxc [ρ]
lr,µ
is the standard Hartree-exchange-correlation functional of the KS scheme and EHxc [ρ] is the
lr,µ
Hartree-exchange-correlation functional associated with the long-range interaction wee (r12 ).
It is often convenient to decompose the short-range functional as (see Refs. [196–198] for an
alternative decomposition)
sr,µ sr,µ
ĒHxc [ρ] = EH [ρ] + Exsr,µ [ρ] + Ēcsr,µ [ρ], (5.20)
sr,µ
where EH [ρ] is the short-range Hartree functional,
1
Z
sr,µ sr,µ
EH [ρ] = ρ(r1 )ρ(r2 )wee (r12 )dr1 dr2 , (5.21)
2 R3 ×R3
sr,µ lr,µ
with the short-range interaction wee (r12 ) = 1/r12 − wee (r12 ), Exsr,µ [ρ] is the short-range ex-
change functional [Eq. (4.13)] which can also be written as
sr,µ
Exsr [ρ] = hΦ[ρ]|Ŵee
sr,µ
|Φ[ρ]i − EH [ρ], (5.22)
with the KS single-determinant wave function Φ[ρ], and Ēcsr,µ [ρ] is the complementary short-
range correlation functional. Just like for Eq. (5.8), the decomposition in Eq. (5.19) leads to the
following expression of the exact ground-state energy
n o
lr,µ sr,µ
E0 = inf hΨ|T̂ + V̂ne + Ŵee |Ψi + ĒHxc [ρΨ ] , (5.23)
Ψ∈W N
where the infimum is over general multideterminant wave functions Ψ ∈ W N .

To obtain a MP2/DFT hybrid scheme, we proceed analogously to Section 5.1. First, we define
the following range-separated hybrid (RSH) approximation by restricting the minimization in
Eq. (5.23) to single-determinant wave functions Φ ∈ S N ,
n o
E0RSH,µ = inf hΦ|T̂ + V̂ne + Ŵee lr,µ sr,µ
|Φi + ĒHxc [ρΦ ] , (5.24)
Φ∈S N
obtaining an energy which necessarily depends on µ. A minimizing single-determinant wave

function Φµ must satisfy the self-consistent eigenvalue equation

T̂ + V̂ne + V̂Hx [Φ ] + V̄ˆHxc
lr,µ,HF µ sr,µ
[ρΦµ ] |Φµ i = E0µ |Φµ i, (5.25)
43
lr,µ,HF µ
where V̂Hx [Φ ] is the nonlocal long-range HF potential operator evaluated with the RSH
wave function Φµ and V̄ˆHxc
sr,µ
[ρΦµ ] is the local short-range Hartree-exchange-correlation potential
sr,µ
operator generated by the energy functional ĒHxc [ρ] and evaluated at the RSH density ρΦµ . The
RSH ground-state energy can be finally written as
E0RSH,µ = hΦµ |T̂ + V̂ne |Φµ i + EH [ρΦµ ] + Exlr,µ,HF [Φµ ] + Exsr,µ [ρΦµ ] + Ēcsr,µ [ρΦµ ], (5.26)
where the full Hartree energy EH [ρ] has been recomposed. The exchange-correlation energy in
Eq. (5.26) has a similar form as in the LC scheme of Eq. (4.11).
To calculate the missing long-range correlation energy in Eq. (5.26), we can define a nonlinear
Rayleigh-Schrödinger perturbation theory [173–175] starting from the RSH reference. We start
from the following energy expression with the perturbation parameter α ∈ [0, 1],
n o
E0µ,α = lr,µ,HF µ
inf hΨ|T̂ + V̂ne + V̂Hx lr,µ
[Φ ] + αŴ lr,µ |Ψi + ĒHxc [ρΨ ] , (5.27)
Ψ∈W N

lr,µ lr,µ,HF µ
where Ŵ lr,µ = Ŵee − V̂Hx [Φ ] is the long-range Møller–Plesset perturbation operator.
For α = 0, the stationary equation associated with Eq. (5.27) reduces to the RSH eigenvalue
equation [Eq. (5.25)]. For α = 1, Eq. (5.27) reduces to Eq. (5.23), so E0µ,α=1 is the exact energy,
independently of µ. The sum of the zeroth-order energy and first-order energy correction gives
µ,(0) µ,(1)
simply the RSH energy, E0RSH,µ = E0 + E0 . As in Section 5.1, only double excitations con-
tribute to the first-order wave-function correction Ψµ,(1) and the second-order energy correction
has a standard MP2 form
µ,(2)
E0 = hΦµ |Ŵ lr,µ |Ψµ,(1) i = Eclr,µ,MP2 , (5.28)
where Eclr,µ,MP2 has the same expression as in Eq. (5.2) with RSH spin orbitals and associated
orbital eigenvalues (which implicitly depend on µ) but using the long-range two-electron integrals
Z
lr,µ
hφp φq |φr φs i = φ∗p (x1 )φ∗q (x2 )φr (x1 )φs (x2 )wee
lr,µ
(r12 )dx1 dx2 , (5.29)
(R3 ×{↑,↓})2
instead of the standard two-electron integrals of Eq. (5.3). This second-order perturbation theory
defines a RSH+MP2 approximation,
E0RSH+MP2,µ = E0RSH,µ + Eclr,µ,MP2 , (5.30)
which contains the exchange-correlation energy contribution

RSH+MP2,µ
Exc = Exlr,µ,HF [Φµ ] + Exsr,µ [ρΦµ ] + Ēcsr,µ [ρΦµ ] + Eclr,µ,MP2 . (5.31)
When using semilocal density-functional approximations for the short-range functionals Exsr,µ [ρ]
and Ēcsr,µ [ρ], the RSH+MP2 exchange-correlation energy expression of Eq. (5.31) thus consti-
tutes range-separated double-hybrid approximations similar to the double hybrids of Section 5.1.
The optimal value for the range-separation parameter is often around µ ≈ 0.5 bohr−1 [150,199].
This scheme has the advantage of dropping the long-range part of both the exchange and cor-
relation density functionals which are usually not well described by semilocal density-functional
approximations. Moreover, using a long-range MP2 correlation energy has the advantage of lead-
ing to a correct qualitative description of London dispersion interaction energies [173, 200–202],
while displaying a fast convergence with the one-electron basis size [203]. Similar to the SCS
double hybrids [Eq. (5.18)], a SCS variant of the RSH+MP2 scheme has also been proposed [204].
The range-separated multideterminant extension of the KS scheme of Eq. (5.23) allows one
to define various hybrid schemes combining any wave-function method with density functionals.
44
For example, one can go beyond second order by using long-range coupled-cluster [205–208]
or random-phase approximations [207, 209–212]. To describe strongly correlated systems, one
can also use for the long-range part wave-function methods such as configuration interaction
(CI) [213–216], MCSCF [217–219], density-matrix renormalization group (DMRG) [220], or mul-
tireference perturbation theory [221]. Density-matrix functional theory (DMFT) [222–224] and
Green-function methods [225, 226] have also been used for the long-range part.
We now consider the approximations used for Exsr,µ [ρ] and Ēcsr,µ [ρ]. In Section 4.2, we have
already described the short-range exchange LDA [Eq. (4.16)] and some short-range exchange
GGAs for Exsr,µ [ρ]. Here, we describe the short-range LDA correlation functional and another
short-range GGA exchange-correlation functional.
Short-range LDA correlation functional

The complementary short-range LDA (or LSDA) correlation functional is
Z
sr,µ,LSDA
Ēc [ρ↑ , ρ↓ ] = ēcsr,µ,UEG (ρ↑ (r), ρ↓ (r))dr, (5.32)
R3
where ēsr,µ,UEG
c (ρ↑ , ρ↓ ) = ρ ε̄sr,µ,UEG
c (ρ↑ , ρ↓ ) is the complementary short-range UEG correlation
sr,µ,UEG
energy density. In this expression, ε̄c (ρ↑ , ρ↓ ) is defined by
ε̄csr,µ,UEG (ρ↑ , ρ↓ ) = εUEG

c (ρ↑ , ρ↓ ) − εclr,µ,UEG (ρ↑ , ρ↓ ), (5.33)
where εUEG
c (ρ↑ , ρ↓ ) and εclr,µ,UEG (ρ↑ , ρ↓ ) are the correlation energies per particle of the UEG
with the standard Coulomb and long-range electron-electron interactions, respectively. A simple
spin-independent parametrization of ε̄sr,µ,UEGc was given in Ref. [227]. A better spin-dependent
parametrization was constructed in Ref. [228] which uses the PW92 parametrization for εUEG c (ρ↑ , ρ↓ )
[Eq. (3.8)] and the following parametrization for εclr,µ,UEG (ρ↑ , ρ↓ ) in terms of rs = (3/(4πρ))1/3
and ζ = (ρ↑ − ρ↓ )/ρ:
εclr,µ,UEG (ρ↑ , ρ↓ ) =
√
µ r
φ2 (ζ)3 Q φ2 (ζ)s + a1 (rs , ζ)µ3 + a2 (rs , ζ)µ4 + a3 (rs , ζ)µ5 + a4 (rs , ζ)µ6 + a5 (rs , ζ)µ8
. (5.34)
(1 + b0 (rs )2 µ2 )4
In this expression, φ2 (ζ) is a spin-scaling function defined by Eq. (3.5), Q(x) is a function
determined from the small-µ and/or small-rs limit,
1 + ax + bx2 + cx3

2 ln(2) − 2
Q(x) = ln , (5.35)
π2 1 + ax + dx2
with a = 5.84605, c = 3.91744, d = 3.44851, b = d − 3πα/[4 ln(2) − 4], α = 4/(9π)1/3 , and the
functions ai (rs , ζ) are
a1 (rs , ζ) = 4b0 (rs )6 C3 (rs , ζ) + b0 (rs )8 C5 (rs , ζ), (5.36a)
a2 (rs , ζ) = 4b0 (rs )6 C2 (rs , ζ) + b0 (rs )8 C4 (rs , ζ) + 6b0 (rs )4 εPW92

c (rs , ζ), (5.36b)
a3 (rs , ζ) = b0 (rs )8 C3 (rs , ζ), (5.36c)
a4 (rs , ζ) = b0 (rs )8 C2 (rs , ζ) + 4b0 (rs )6 εPW92
c (rs , ζ), (5.36d)
a5 (rs , ζ) = b0 (rs )8 εPW92
c (rs , ζ), (5.36e)
45
where εPW92
c (rs , ζ) is the PW92 parametrization of the UEG correlation energy per particle.
The functions Ci (rs , ζ) are determined from the large-µ limit,
3gc (0, rs , ζ)
C2 (rs , ζ) = − , (5.37a)
8rs3
g(0, rs , ζ)
C3 (rs , ζ) = − √ , (5.37b)
2πrs3
9[gc′′ (0, rs , ζ) + (1 − ζ 2 )D2 (rs )]
C4 (rs , ζ) = − , (5.37c)
64rs3
9[g′′ (0, rs , ζ) + (1 − ζ 2 )D3 (rs )]
C5 (rs , ζ) = − √ , (5.37d)
40 2πrs3
where g(0, rs , ζ) is the on-top pair-distribution function8 of the Coulombic UEG and g ′′ (0, rs , ζ)
is its second-order derivative with respect to r12 at r12 = 0, and similarly for their correlation
parts gc (0, rs , ζ) = g(0, rs , ζ)−(1−ζ 2 )/2 and gc′′ (0, rs , ζ) = g ′′ (0, rs , ζ)−φ8 (ζ)/(5α2 rs2 ) with φ8 (ζ)
defined by Eq. (3.5). The ζ-dependence of the latter quantities is assumed to be exchange-like, i.e.
2 g ′′ (0, r /ζ 1/3 , ζ = 1) + ζ 2 g ′′ (0, r /ζ 1/3 , ζ =
g(0, rs , ζ) ≈ (1 − ζ 2 )g(0, rs , ζ = 0) and g′′ (0, rs , ζ) ≈ ζ+ s + − s −
1) where ζ± = (1 ± ζ)/2. The on-top pair-distribution function has been parametrized in
Ref. [229] as
g(0, rs , ζ = 0) = (1 − Brs + Crs2 + Drs3 + Ers4 )e−F rs , (5.38)
with B = 0.7317 − d, C = 0.08193, D = −0.01277, E = 0.001859, and F = 0.7524. The
remaining functions were determined by fitting to QMC data:
b0 (rs ) = 0.784949rs , (5.39a)
25/3 1 − 0.02267rs
g′′ (0, rs , ζ = 1) = , (5.39b)
5α2 rs2 1 + 0.4319rs + 0.04rs2
e−0.547rs
D2 (rs ) = (−0.388rs + 0.676rs2 ), (5.39c)
rs2
e−0.31rs
D3 (rs ) = (−4.95rs + rs2 ). (5.39d)
rs3
Short-range PBE(GWS) exchange-correlation functional

The Goll-Werner-Stoll (GWS) variant of the short-range PBE exchange-correlation func-
tional [205,206] is a slight modification of the short-range PBE functional developed in Ref. [230].
The exchange energy density is
exsr,µ,PBE(GWS) (ρ, ∇ρ) = esr,µ,UEG

x (ρ)Fx (s, µ̃), (5.40)
with an enhancement factor of the same form as in the standard PBE exchange functional,
κ
Fx (s, µ̃) = 1 + κ − , (5.41)
1 + b(µ̃)s2 /κ
with s = |∇ρ|/(2kF ρ) and µ̃ = µ/(2kF ). In this expression, κ = 0.840, as in the stan-

dard PBE exchange functional, to saturate the local Lieb–Oxford bound (for µ = 0) and
8
For a general system, the pair-distribution function g(r1 , r2 ) is defined from the pair density ρ2 (r1 , r2 )
[Eq. (2.1)] as ρ2 (r1 , r2 ) = ρ(r1 )ρ(r2 )g(r1 , r2 ). The on-top value is the value at electron coalescence, i.e. for
r1 = r2 .
46
2
b(µ̃) = bPBE [bT (µ̃)/bT (0)]e−αx µ̃ where bPBE = 0.21951 is the second-order gradient-expansion
coefficient of the standard PBE exchange functional, and bT (µ̃) is a function coming from the
second-order GEA of the short-range exchange energy [196, 230],
2
T −c1 (µ̃) + c2 (µ̃)e1/(4µ̃ )
b (µ̃) = , (5.42)
c3 (µ̃) + 54c4 (µ̃)e1/(4µ̃2 )
with c1 (µ̃) = 1 + 22µ̃2 + 144µ̃4 , c2 (µ̃) = 2µ̃2 (−7 + 72µ̃2 ), c3 (µ̃) = −864µ̃4 (−1 + 2µ̃2 ), and
√
c4 (µ̃) = µ̃2 [−3 − 24µ̃2 + 32µ̃4 + 8µ̃ π erf(1/(2µ̃))]. Finally, αx = 19.0 is a damping parameter
optimized for the He atom.
Similarly, the correlation energy density has the same form as the standard PBE correlation
functional,
ēcsr,µ,PBE(GWS) (ρ↑ , ρ↓ , ∇ρ↑ , ∇ρ↓ ) = ρ ε̄sr,µ,UEG (ρ↑ , ρ↓ ) + H µ (ρ↑ , ρ↓ , t) ,

c (5.43)
with t = |∇ρ|/(2φ2 (ζ)ks ρ) and the gradient correction

1 + A(µ)t2

µ 3 β(µ) 2
H (ρ↑ , ρ↓ , t) = A(0)φ2 (ζ) ln 1 + t , (5.44)
A(0) 1 + A(µ)t2 + A(µ)2 t4
where
β(µ) −1
A(µ) = exp(−ε̄csr,µ,UEG (ρ↑ , ρ↓ )/(A(0)φ2 (ζ)3 )) − 1 , (5.45)
A(0)
and
!αc
ε̄csr,µ,UEG (ρ↑ , ρ↓ )
β(µ) = β PBE , (5.46)
ε̄csr,µ=0,UEG (ρ↑ , ρ↓ )
and the value of A(0) is given after Eq. (3.6). In Eq. (5.46), β = 0.066725 is the second-order
gradient coefficient of the standard PBE correlation functional and αc = 2.78 is a damping
parameter optimized for the He atom.
For µ = 0, this short-range PBE exchange-correlation functional reduces to the standard
PBE exchange-correlation functional and for large µ it reduces to the short-range LDA exchange-
correlation functional.
5.2.2 Range-separated two-parameter double-hybrid approximations
In 2018, Kalai and Toulouse [231] introduced what we will call range-separated two-parameter
double-hybrid approximations, combining HF exchange and MP2 correlation with density func-
tionals using a two-parameter decomposition of the electron-electron in a way reminiscent of the
CAM decomposition [Eq. (4.15)] (see, also, Refs. [208,232]). This is based on a multideterminant
extension of the KS scheme which generalizes the schemes of Sections 5.1 and 5.2.1.
We first decompose the universal density functional of Eq. (1.8) as
lr,µ sr,µ sr,µ,λ
F [ρ] = min hΨ|T̂ + Ŵee + λŴee |Ψi + ĒHxc [ρ], (5.47)
Ψ∈WρN
where the parameter µ ∈ [0, +∞) controls the range of the separation as always, the parameter
λ ∈ [0, 1] corresponds to the fraction of the short-range electron-electron interaction in the
sr,µ,λ
wave-function part, and ĒHxc [ρ] is the complementary short-range density functional defined
to make this decomposition exact. As before, the latter functional can be decomposed as
sr,µ,λ sr,µ,λ
ĒHxc [ρ] = EH [ρ] + Exsr,µ,λ [ρ] + Ēcsr,µ,λ [ρ]. (5.48)
47
The Hartree and exchange contributions are linear in λ,
sr,µ,λ sr,µ
EH [ρ] = (1 − λ)EH [ρ], (5.49)
Exsr,µ,λ [ρ] = (1 − λ)Exsr,µ [ρ], (5.50)

sr,µ
where EH [ρ] and Exsr,µ [ρ] are the short-range Hartree and exchange functionals introduced in
Section 5.2.1, and the correlation contribution can be written as
Ēcsr,µ,λ [ρ] = Ec [ρ] − Ecµ,λ [ρ], (5.51)
where Ec [ρ] is the standard KS correlation functional and Ecµ,λ [ρ] is the correlation functional
lr,µ sr,µ
associated with the interaction wee (r12 ) + λwee (r12 ). The exact ground-state energy can then
be expressed as
n o
lr,µ sr,µ sr,µ,λ
E0 = inf hΨ|T̂ + V̂ne + Ŵee + λŴee |Ψi + ĒHxc [ρΨ ] , (5.52)
Ψ∈W N
which constitutes a generalization of Eqs. (5.8) and (5.23).

To obtain a MP2/DFT hybrid scheme, we proceed in full analogy to Sections 5.1 and 5.2.1.
First, we define the following single-determinant range-separated two-parameter hybrid (RS2H)
approximation,
n o
E0RS2H,µ,λ = inf hΦ|T̂ + V̂ne + Ŵee lr,µ sr,µ
+ λŴee sr,µ,λ
|Φi + ĒHxc [ρΦ ] , (5.53)
Φ∈S N
and use it as a reference for defining a perturbation theory similarly to Eqs. (5.12) and (5.27).
At second order, we obtain
E0RS2H+MP2,µ,λ = E0RS2H,µ,λ + Ecµ,λ,MP2 , (5.54)
where Ecµ,λ,MP2 is the MP2 correlation energy expression evaluated with RS2H spin orbitals and
lr,µ
orbital eigenvalues, and the two-electron integrals associated with the interaction wee (r12 ) +
sr,µ
λwee (r12 ). This RS2H+MP2 scheme thus contains the exchange-correlation energy contribu-
tion
RS2H+MP2,µ,λ
Exc = Exlr,µ,HF [Φµ,λ ] + λExsr,µ,HF [Φµ,λ ] + (1 − λ)Exsr,µ [ρΦµ,λ ]
+Ēcsr,µ,λ [ρΦµ,λ ] + Ecµ,λ,MP2 , (5.55)
where Φµ,λ is a minimizing single-determinant wave function in Eq. (5.53).

A good approximation for the λ-dependence of the complementary correlation functional
Ēcsr,µ,λ [ρ]
is [231]
√
Ēcsr,µ,λ [ρ] ≈ Ēcsr,µ [ρ] − λ2 Ēcsr,µ λ
[ρ], (5.56)
where Ēcsr,µ [ρ] is the short-range correlation functional introduced in Section 5.2.1. In particular,
the λ-dependence in Eq. (5.56) is correct both in the high-density limit, for a non-degenerate KS
system, and in the low-density limit. Thanks to Eqs. (5.50) and (5.56), the semilocal density-
functional approximations for Exsr,µ [ρ] and Ēcsr,µ [ρ] of Section 5.2.1 can be reused here without
developing new ones. In Ref. [231], the short-range PBE(GWS) exchange and correlation func-
tionals were used, and the optimal parameters µ = 0.46 bohr−1 and λ = 0.58 were found on
small sets of atomization energies and reaction barrier heights, i.e. values similar to the ones
usually used separately in range-separated hybrids and double hybrids.
48
The RS2H+MP2 scheme improves a bit over the RSH+MP2 scheme of Section 5.2.1, partic-
ularly for interaction energies of hydrogen-bonded systems. Even if the presence of short-range
MP2 correlation deteriorates in principle the convergence rate with the one-electron basis size,
in practice the fraction of pure short-range MP2 correlation (λ2 ≈ 0.34) is small enough to keep
a fast basis convergence. Accuracy can be improved, particularly for dispersion interactions, by
supplanting the MP2 term by coupled-cluster or random-phase approximations [233]. Like for
the approach of Section 5.2.1, many wave-function methods could be used in the general scheme
of Eq. (5.52).
6 Semiempirical dispersion corrections and nonlocal van der Waals den-

sity functionals
Among the previously considered exchange-correlation approximations, only the range-separated
double hybrids of Section 5.2, thanks to their long-range nonlocal correlation component, are
capable of fully describing London dispersion interactions, crucial for describing weakly bonded
systems. To improve the other approximations (semilocal functionals, single-determinant hy-
brids, double hybrids without range separation) for weakly bonded systems, it has been pro-
posed to add to them a semiempirical dispersion correction or a nonlocal van der Waals density
functional. We now describe these approaches.
6.1 Semiempirical dispersion corrections

To explicitly account for London dispersion interactions, it has been proposed in the 2000s
to add to the standard approximate functionals a semiempirical dispersion correction of the
form [234–237]
Nn X
Nn
X C6αβ
Edisp = −s f (Rαβ ) 6 , (6.1)
α=1 β=1
Rαβ
β>α
where Rαβ is the distance between each pair of atoms and C6αβ is the London dispersion coefficient
between these atoms. Here, f (Rαβ ) is a damping function which tends to 1 at large Rαβ and
tends to zero at small Rαβ , e.g.
1
f (Rαβ ) = , (6.2)
−d(Rαβ /RvdW
αβ −1)
1+e
with the sum of tabulated atomic van der Waals radii Rαβ vdW = RvdW + RvdW and a constant d,
α β
and s is a scaling parameter that can be adjusted for each approximate functional. The disper-
sion coefficient C6αβ for any pair of atoms is empirically calculated from tabulated same-atom
dispersion coefficients C6αα and/or atomic polarizabilities. This approach was named “DFT-D”
by Grimme [236].
The last version of DFT-D (referred to as DFT-D3) also includes C8αβ two-body terms and
C9αβγ three-body terms [238]. There have also been various proposals to make the determination
of dispersion coefficients less empirical, such as the scheme of Becke and Johnson [239] based
on the exchange-hole dipole moment, the scheme of Tkatchenko and Scheffler [240] based on a
Hirshfeld atomic partitioning, or the scheme of Sato and Nakai [241] based on the local-response
approximation [242].
The “DFT-D” approach provides a big and inexpensive improvement for the description of
weakly bonded systems. One limitation is that the semiempirical dispersion correction, being
49
just a force field in its simplest variant, affects only the molecular geometry of the system but
not directly its electronic structure. Some of the most used DFT-D functionals are:
• The PBE-D exchange-correlation functional [237], based on the PBE functional with a
scaling parameter s = 0.75;
• The B97-D exchange-correlation functional [237], based on the B97-GGA functional with
a scaling parameter s = 1.25 and reoptimized linear coefficients in Eqs. (3.38), (3.42),
and (3.43) in the presence of the semiempirical dispersion correction;
• The B3LYP-D exchange-correlation functional [237], based on the B3LYP hybrid func-
tional with a scaling parameter s = 1.05;
• The ωB97X-D exchange-correlation functional [243], based on the ωB97X range-separated

hybrid functional with a scaling parameter s = 1, a modified damping function, and reop-
timized parameters in ωB97X in the presence of the semiempirical dispersion correction.
The semiempirical dispersion correction can also be added to double-hybrid approximations.

For example, B2PLYP-D [244] is based on the B2PLYP double hybrid with a scaling parameter
s = 0.55. The scaling parameter is small since the fraction of MP2 correlation in B2PLYP
already partially takes into account dispersion interactions. It has also been proposed to add
a semiempirical dispersion correction to the SCS version of the double hybrids [Eq. (5.18)], re-
sulting in a family of dispersion-corrected spin-component-scaled double-hybrid (DSD) approx-
imations [179, 180, 245]. An example of double hybrid is this latter family is DSD-BLYP [179]
which uses the B88 exchange functional and the LYP correlation functional.
6.2 Nonlocal van der Waals density functionals

Another approach to describe dispersion interactions is to add to the standard approximate
functionals a so-called nonlocal van der Waals density functional of the form [246–250]
1
Z
nl
Ec [ρ] = ρ(r1 )ρ(r2 )φ(r1 , r2 )dr1 dr2 , (6.3)
2 R3 ×R3
where φ(r1 , r2 ) is a correlation kernel. Two main families of such nonlocal correlation functionals
exist: the “van der Waals density functionals” (vdW-DF) [246,249] and the Vydrov-Van Voorhis
(VV) functionals [247, 248, 250].
We will only describe the last version of the VV functionals, i.e. the VV10 nonlocal correla-
tion functional [250]. In this functional, the correlation kernel is taken as
3
φVV10 (r1 , r2 ) = − + βδ(r1 − r2 ), (6.4)
2g(r1 , r12 )g(r2 , r12 )(g(r1 , r12 ) + g(r2 , r12 ))
where r12 = |r2 − r1 | is the interelectronic distance, β is a constant determining the local (delta-
distribution) part of the kernel, and the function g is defined as
2
g(r, r12 ) = ω0 (r)r12 + κ(r). (6.5)
q
ω (r)2
In Eq. (6.5), ω0 (r) = ωg (r)2 + p 3 involves the square of the local plasma frequency ωp (r)2 =
4πρ(r) and the square of the local band gap ωg (r)2 = C|∇ρ(r)|4 /ρ(r)4 where C is an adjustable
parameter controlling the large-r12 asymptotic dispersion coefficients, and κ(r) = b kF (r)2 /ωp (r)
where kF (r) = (3π 2 ρ(r))1/3 is the local Fermi wave vector and b is an adjustable parameter
50
controlling the short-range damping of the large-r12 asymptote. As expected for dispersion
interactions, in the large-r12 limit, φVV10 (r1 , r2 ) behaves as 1/r12
6 :
3
φVV10 (r1 , r2 ) ∼ − 6 . (6.6)
r12 →∞ 2ω0 (r1 )ω0 (r2 )(ω0 (r1 ) + ω0 (r2 ))r12
The constant β = (3/b2 )3/4 /16 is chosen to make Ecnl [ρ] vanish in the uniform density limit,
leaving thus this limit unchanged when Ecnl [ρ] is added to another density functional. The ad-
justable parameters C ≈ 0.009 and b ≈ 6 are found by optimization of C6 dispersion coefficients
and of weak intermolecular interaction energies, respectively, the precise values depending on
which exchange-correlation functional the VV10 correction is added to.
Nonlocal van der Waals density functionals are necessarily more computationally expensive
than semiempirical dispersion corrections. However, they have the advantage of being less
empirical and, since they are functionals of the density, of impacting directly on the electronic
structure of the system. The VV10 nonlocal functional has been incorporated in a number of
recently developed exchange-correlation functionals, for example:
• The ωB97X-V exchange-correlation functional [251], based on the ωB97X range-separated

hybrid [Eq. (4.23)] with reoptimized linear coefficients in Eq. (3.38) with polynomial degree
m = 2 and in Eqs. (3.42) and (3.43) with polynomial degree m = 1, as well as reoptimized
VV10 parameters C = 0.01 and b = 6.0;
• The ωB97M-V exchange-correlation functional [252], based on the ωB97X range-separated
hybrid [Eq. (4.23)] but with more general and combinatorially optimized meta-GGA ex-
change and correlation enhancement factors and the same VV10 parameters C = 0.01 and
b = 6.0 as in ωB97X-V.
7 Orbital-dependent exchange-correlation density functionals

We discuss here some exchange-correlation density functionals explicitly depending on the
KS orbitals (for a review, see Ref. [253]). Since the KS orbitals are themselves functionals of
the density, these exchange-correlation expressions are thus implicit functionals of the density
(for notational simplicity, this dependence on the density of the orbitals and other intermediate
quantities will not be explicitly indicated). In fact, the single-determinant and multideterminant
hybrid approximations of Sections 4 and 5 already belong to this family, with the caveat that the
orbitals are obtained with a nonlocal potential. In this section, we are concerned with orbital-
dependent exchange-correlation energy functionals with orbitals obtained with a local potential,
i.e. staying within the KS scheme9 . These approximations tend to be more computationally
involved than the approximations previously seen and have thus been much less used so far.
7.1 Exact exchange

The exact exchange (EXX) energy functional [Eq. (1.21)] can be expressed in terms of the
KS orbitals,
Nσ Z
Nσ X
1 X X ϕ∗iσ (r1 )ϕjσ (r1 )ϕ∗jσ (r2 )ϕiσ (r2 )
Ex [ρ] = − dr1 dr2 , (7.1)
2 R3 ×R3 |r1 − r2 |
σ∈{↑,↓} i=1 j=1
9
The boundary between the various single-determinant and multideterminant hybrids of Sections 4 and 5 and
the orbital-dependent functionals of the present Section is however thin. For example, it is possible to optimize
the orbitals using a local potential in hybrids or range-separated hybrids [254–256], and in double hybrids or
range-separated double hybrids [257, 258].
51
and has exactly the same form as the HF exchange [Eq. (4.2)], but the orbitals used in both
expressions are in general different.
Since the exact exchange energy in Eq. (7.1) is not an explicit functional of the density, the
corresponding exchange potential vx (r) = δEx [ρ]/δρ(r) cannot be calculated directly. We can
however find an workable equation for vx (r) by first considering the functional derivative of Ex [ρ]
with respect to the KS potential vs (r) and then applying the chain rule:
δEx [ρ] δEx [ρ] δρ(r′ ) ′
Z
= ′
dr . (7.2)
δvs (r) R3 δρ(r ) δvs (r)
Introducing the non-interacting KS static linear-response function χ0 (r′ , r) = δρ(r′ )/δvs (r), we
can rewrite Eq. (7.2) as
δEx [ρ]
Z
vx (r′ )χ0 (r′ , r)dr′ = , (7.3)
R 3 δvs (r)
which is known as the optimized-effective-potential (OEP) equation for the exact-exchange po-
tential [48, 259, 260].
Using first-order perturbation theory on the KS system, explicit expressions in terms of the
orbitals can be derived for χ0 (r′ , r) and δEx [ρ]/δvs (r). The expression of χ0 (r′ , r) is
Nσ
′
X X X ϕ∗iσ (r′ )ϕ∗aσ (r)ϕiσ (r)ϕaσ (r′ )
χ0 (r , r) = − + c.c. , (7.4)
εaσ − εiσ
σ∈{↑,↓} i=1 a≥Nσ +1
where i and a run over occupied and virtual spatial orbitals, respectively, and c.c. stands for
the complex conjugate. The expression of δEx [ρ]/δvs (r) is
Nσ
Nσ X
δEx [ρ] X X X ϕaσ (r)ϕ∗iσ (r)
= hϕaσ ϕjσ |ϕjσ ϕiσ i + c.c. , (7.5)
δvs (r) εaσ − εiσ
σ∈{↑,↓} i=1 j=1 a≥Nσ +1
where hϕaσ ϕjσ |ϕjσ ϕiσ i are two-electron integrals over KS spatial orbitals:
ϕ∗aσ (r1 )ϕ∗jσ (r2 )ϕjσ (r1 )ϕiσ (r2 )
Z
hϕaσ ϕjσ |ϕjσ ϕiσ i = dr1 dr2 . (7.6)
R3 ×R3 |r1 − r2 |
Applying this OEP method with the EXX energy (and no correlation energy functional)
is an old idea [259, 261], but reasonably efficient calculations for molecules have been possible
only relatively recently [262, 263]. The EXX occupied orbitals turn out to be very similar to
the HF occupied orbitals, and thus the EXX ground-state properties are also similar to the HF
ones. However, the EXX virtual orbitals (which see a −1/r asymptotic potential for a neutral
system) tend to be much less diffuse than the HF virtual orbitals (which see an exponentially
decaying potential for a neutral system), and may be more adapted for calculating excited-state
properties.
7.2 Second-order Görling–Levy perturbation theory

In 1993, Görling and Levy [47, 48] developed a perturbation theory in terms of the coupling
constant λ of the adiabatic connection (Section 2.2) which provides an explicit orbital-dependent
second-order approximation for the correlation energy functional. The Hamiltonian along the
adiabatic connection [Eq. (2.22)] can be written as
Ĥ λ = T̂ + λŴee + V̂ λ
= Ĥs + λ(Ŵee − V̂Hx ) − V̂cλ , (7.7)
52
where Ĥs = Ĥ λ=0 = T̂ + V̂s is the KS non-interacting reference Hamiltonian (which will be
assumed to have a nondegenerate ground state). Equation (7.7) was obtained by decomposing
the potential operator keeping the density constant as V̂ λ = V̂s − λV̂Hx − V̂cλ where V̂s = V̂ λ=0
is the KS potential operator, λV̂Hx is the Hartree-exchange potential operator which is linear in
(2)
λ, and V̂cλ is the correlation potential which starts at second order in λ, i.e. V̂cλ = λ2 V̂c + · · · .
Using a complete set of orthonormal eigenfunctions Φn and eigenvalues En of the KS Hamiltonian,
Ĥs |Φn i = En |Φn i, the normalized ground-state wave function of the Hamiltonian Ĥ λ can be
expanded as Ψλ = Φ + λΨ(1) + · · · where Φ = Φ0 is the ground-state KS single-determinant
wave function and Ψ(1) is its first-order correction given by
X hΦn |Ŵee − V̂Hx |Φi
|Ψ(1) i = − |Φn i. (7.8)
En − E0
n6=0
Using the expression in Eq. (2.27), the correlation energy functional can also be expanded
in powers of λ:
Ecλ [ρ] = hΨλ |T̂ + λŴee |Ψλ i − hΦ|T̂ + λŴee |Φi.
= Ec(0) + λEc(1) + λ2 Ec(2) + · · · . (7.9)
(0)
Since Ψλ=0 = Φ, the zeroth-order term vanishes: Ec
= 0. Using the expression of the first-
order derivative of Ecλ with respect to λ in Eq. (2.28), i.e. ∂Ecλ /∂λ = hΨλ |Ŵee |Ψλ i − hΦ|Ŵee |Φi,
(1)
we find that the first-order term vanishes as well: Ec = 0. The second-order term corresponds
to the second-order Görling–Levy (GL2) correlation energy and is given by
EcGL2 [ρ] = Ec(2) = hΦ|Ŵee |Ψ(1) i = hΦ|Ŵee − VHx |Ψ(1) i, (7.10)
where the second equality comes from the factR that hΦ|V̂Hx |Ψ(1) i = 0 since it is the derivative
with respect to λ at λ = 0 of hΨλ |V̂Hx |Ψλ i = R3 vHx (r)ρ(r)dr which does not depend on λ by
virtue of the fact that the density ρ(r) is constant along the adiabatic connection. Using the
last expression in Eq. (7.10) allows one to express the GL2 correlation energy as
X |hΦ|Ŵee − V̂Hx |Φn i|2
EcGL2 [ρ] = − . (7.11)
En − E0
n6=0
It is instructive to decompose the GL2 correlation energy as

EcGL2 [ρ] = EcMP2 + EcS , (7.12)
where EcMP2 is a MP2-like correlation energy evaluated with KS spin orbitals,
N N
1 XX X X |hφi φj ||φa φb i|2
EcMP2 = − , (7.13)
i=1 j=1 a≥N +1 b≥N +1
and EcSis the contribution coming from the single excitations (which does not vanish here,
contrary to HF-based MP2 perturbation theory),
N
X X |hφi |V̂ HF − V̂x |φa i|2
x
EcS =− , (7.14)
εa − εi
i=1 a≥N +1
involving the difference between the integrals over the nonlocal HF exchange HF
PN R potential∗hφi |V̂x |φa i =
− j=1 hφi φj |φj φa i and over the local KS exchange potential hφi |V̂x |φa i = R3 ×{↑,↓} φi (x)vx (r)φa (x)dx.
Calculations of the GL2 correlation energy using either a non-self-consistent post-EXX im-
plementation or a more complicated OEP self-consistent procedure have been tested (see, e.g.,
Refs. [264–266]) but the results are often disappointing. It is preferable to go beyond second
order with random-phase approximations in the adiabatic-connection fluctuation-dissipation ap-
proach.
53
7.3 Adiabatic-connection fluctuation-dissipation approach
7.3.1 Exact adiabatic-connection fluctuation-dissipation expression
Using the adiabatic-connection formula of Eq. (2.29), the correlation energy functional can
be written as
ρλ2,c (r1 , r2 )
Z 1
1 1
Z Z
λ λ
Ec [ρ] = dλ hΨ |Ŵee |Ψ i − hΦ|Ŵee |Φi = dλ dr1 dr2 , (7.15)
0 2 0 R3 ×R3 |r1 − r2 |
where ρλ2,c (r1 , r2 ) = ρλ2 (r1 , r2 ) − ρ2,KS (r1 , r2 ) is the correlation part of the pair density along
the adiabatic connection. The pair density ρλ2 (r1 , r2 ) can be expressed with the pair-density
operator ρ̂2 (r1 , r2 ) = ρ̂(r1 )ρ̂(r2 ) − δ(r1 − r2 )ρ̂(r1 ) where ρ̂(r) is the density operator,
ρλ2 (r1 , r2 ) = hΨλ |ρ̂2 (r1 , r2 )|Ψλ i = hΨλ |ρ̂(r1 )ρ̂(r2 )|Ψλ i − δ(r1 − r2 )hΨλ |ρ̂(r1 )|Ψλ i, (7.16)
and the KS pair density ρ2,KS (r1 , r2 ) simply corresponds to the case λ = 0,
ρ2,KS (r1 , r2 ) = ρλ=0

2 (r1 , r2 ) = hΦ|ρ̂(r1 )ρ̂(r2 )|Φi − δ(r1 − r2 )hΦ|ρ̂(r1 )|Φi. (7.17)
Since the density does not change with λ, i.e. hΨλ |ρ̂(r)|Ψλ i = hΦ|ρ̂(r)|Φi = ρ(r), the correlation
pair density needed in Eq. (7.15) can thus be expressed as
ρλ2,c (r1 , r2 ) = hΨλ |ρ̂(r1 )ρ̂(r2 )|Ψλ i − hΦ|ρ̂(r1 )ρ̂(r2 )|Φi. (7.18)
We would like to calculate ρλ2,c (r1 , r2 ) without having to calculate the complicated many-
body wave function Ψλ . For this, we consider the retarded dynamic linear-response function
along the adiabatic connection in frequency space (or the so-called Lehmann representation)
X hΨλ |ρ̂(r1 )|Ψλ ihΨλ |ρ̂(r2 )|Ψλ i hΨλ |ρ̂(r2 )|Ψλn ihΨλn |ρ̂(r1 )|Ψλ i
n n
χλ (r1 , r2 ; ω) = − , (7.19)
ω − ωnλ + i0+ ω + ωnλ + i0+
n6=0
where the sums are over all eigenstates Ψλn of the Hamiltonian Ĥ λ , i.e. Ĥ λ |Ψλn i = Enλ |Ψλn i,
except the ground state Ψλ = Ψλ0 , and ωnλ = Enλ − E0λ are the corresponding excitation energies.
By contour integrating χλ (r1 , r2 , ω) around the right half ω-complex plane, we arrive at the
(zero-temperature) fluctuation-dissipation theorem,
Z +∞
dω
nλ2,c (r1 , r2 ) = − [χλ (r1 , r2 , iω) − χ0 (r1 , r2 , iω)], (7.20)
−∞ 2π
which relates ground-state correlations in the time-independent system ρλ2,c (r1 , r2 ) to the linear
response of the system due to a time-dependent external perturbation χλ (r1 , r2 , ω).
Combining Eqs. (7.15) and (7.20), we finally obtain the exact adiabatic-connection fluctuation-
dissipation (ACFD) formula for the correlation energy [40, 42] (see, also, Ref. [43]):
1 1
Z +∞
dω χλ (r1 , r2 ; iω) − χ0 (r1 , r2 ; iω)
Z Z
Ec [ρ] = − dλ dr1 dr2 . (7.21)
2 0 −∞ 2π R3 ×R3 |r1 − r2 |
The usefulness of the ACFD formula is due to the fact that there are practical ways of directly
calculating χλ (r1 , r2 ; ω) without having to calculate the many-body wave function Ψλ . In linear-
response time-dependent density-functional theory (TDDFT), one can find a Dyson-like equation
for χλ (r1 , r2 ; ω),
Z
λ
χλ (r1 , r2 ; ω) = χ0 (r1 , r2 ; ω) + χ0 (r1 , r3 ; ω)fHxc (r3 , r4 ; ω)χλ (r4 , r2 ; ω)dr3 dr4 , (7.22)
R3 ×R3
54
λ (r , r ; ω) is the Hartree-exchange-correlation kernel associated to the Hamiltonian
where fHxc 3 4
H λ . Here, Eq. (7.22) will be considered as the definition for fHxcλ . In principle, the exact
correlation energy can be obtained with Eqs. (7.21) and (7.22). In practice, however, we need
λ .
to use an approximation for fHxc
7.3.2 Random-phase approximations
In the direct random-phase approximation (dRPA, also just referred to as RPA, or some-
times as time-dependent Hartree), only the Hartree part of the kernel, which is linear in λ and
independent from ω, is retained [40, 267],
dRPA,λ
fHxc (r1 , r2 ; ω) = fHλ (r1 , r2 ) = λwee (r1 , r2 ), (7.23)
where wee (r1 , r2 ) = 1/|r1 − r2 | is the Coulomb interaction, and the corresponding dRPA linear-
response function then satisfies the equation
Z
χdRPA
λ (r ,
1 2r ; ω) = χ (r ,
0 1 2 r ; ω) + λ χ0 (r1 , r3 ; ω)wee (r3 , r4 )χdRPA
λ (r4 , r2 ; ω)dr3 dr4 .(7.24)
R3 ×R3
The physical contents of this approximation can be seen by iterating Eq. (7.24) which generates
an infinite series,
Z
dRPA
χλ (r1 , r2 ; ω) = χ0 (r1 , r2 ; ω) + λ χ0 (r1 , r3 ; ω)wee (r3 , r4 )χ0 (r4 , r2 ; ω)dr3 dr4
R3 ×R3
Z
+λ2 χ0 (r1 , r3 ; ω)wee (r3 , r4 )χ0 (r4 , r5 ; ω)wee (r5 , r6 )χ0 (r6 , r2 ; ω)dr3 dr4 dr5 dr6 + · · · ,
R3 ×R3 ×R3 ×R3
(7.25)
which, after plugging it into Eq. (7.21), leads to the dRPA correlation energy as the following
perturbation expansion10
" Z
1 1
Z +∞
dω χ0 (r1 , r3 ; iω)χ0 (r4 , r2 ; iω)
Z
dRPA
Ec [ρ] = − dλ λ dr1 dr2 dr3 dr4
2 0 −∞ 2π 3 3 3
R ×R ×R ×R 3 |r1 − r2 | |r3 − r4 |
#
χ0 (r1 , r3 ; iω)χ0 (r4 , r5 ; iω)χ0 (r6 , r2 ; iω)
Z
2
+λ dr1 dr2 dr3 dr4 dr5 dr6 + · · · . (7.26)
R3 ×R3 ×R3 ×R3 ×R3 ×R3 |r1 − r2 | |r3 − r4 | |r5 − r6 |
Using now the Lehmann representation [Eq. (7.19)] of the KS dynamic linear-response function
in terms of the KS orbitals and their energies,
Nσ
"
X X X ϕ∗iσ (r1 )ϕaσ (r1 )ϕ∗aσ (r2 )ϕiσ (r2 )
χ0 (r1 , r2 ; ω) =
ω − (εaσ − εiσ ) + i0+
σ∈{↑,↓} i=1 a≥Nσ +1
#
ϕ∗iσ (r2 )ϕaσ (r2 )ϕ∗aσ (r1 )ϕiσ (r1 )
− , (7.27)
ω + (εaσ − εiσ ) + i0+
10 dRPA
Using
R +∞the operator viewpoint, the series in Eq. (7.26) can be formally summed in the form Ec [ρ] =
1/(4π) −∞ dω Tr[ln(1 − χ0 (iω)wee ) + χ0 (iω)wee ] (see, e.g., Ref. [268]).
55
one can obtain, after quite some work,
N N
1 XX X X |hφi φj |φa φb i|2
EcdRPA [ρ] = −
i=1 j=1 a≥N +1 b≥N +1
N X N
N X
X X X X hφi φj |φa φb ihφj φk |φb φc ihφk φi |φc φa i
+ + ··· .
(εa + εb − εi − εj )(εa + εc − εi − εk )
i=1 j=1 k=1 a≥N +1 b≥N +1 c≥N +1
(7.28)
The dRPA correlation energy is the sum of all the direct terms (i.e., no exchange terms) of
the perturbation expansion up to infinite order. In the language of diagrammatic perturbation
theory, we say that the dRPA correlation energy is the sum of all direct ring diagrams. Of course,
Eq. (7.28) is not the way to calculate the dRPA correlation energy in practice. This is done by
solving the Dyson equation [Eq. (7.24)] without explicitly expanding in powers of λ, e.g. using
matrix equations from linear-response TDDFT [211, 269] or coupled-cluster theory [207, 270].
Most dRPA correlation energy (combined with the EXX energy) calculations are done
in a non-self-consistent way, but self-consistent OEP dRPA calculations have also been per-
formed [271, 272]. One of the main advantage of dRPA is that it accounts for long-range dis-
persion interactions [273–275]. However, it shows large self-interaction errors. To overcome the
latter drawback and improve the general accuracy, one can add exchange and beyond terms in
various ways (see, e.g., Refs. [209, 211, 268, 276–288]). This remains an active area of research.
For reviews on random-phase approximations, the reader may consult Refs. [289–291].
Acknowledgements
I thank Éric Cancès and Gero Friesecke for discussions and comments on the manuscript.
This review chapter grew out of my lecture notes for DFT courses given in several summer
schools (ISTPC summer schools in June 2015 and 2017 in France, ICS summer school in July
2015 in France, and ESQC summer school in September 2019 in Italy). This work has received
funding from the European Research Council (ERC) under the European Union’s Horizon 2020
research and innovation programme (grant agreement EMC2 No 810367).
References
[1] R. G. Parr and W. Yang, Density-Functional Theory of Atoms and Molecules (Oxford
University Press, New York, 1989).
[2] R. M. Dreizler and E. K. U. Gross, Density Functional Theory (Springer-Verlag, Berlin,
1990).
[3] C. Fiolhais, F. Nogueira and M. A. L. Marques, eds., A Primer in Density Functional
Theory, Vol. 620 of Lecture Notes in Physics (Springer, Berlin, 2003).
[4] E. Engel and R. M. Dreizler, Density Functional Theory: An Advanced Course, Theoretical
and Mathematical Physics (Springer-Verlag, Berlin Heidelberg, 2011).
[5] T. Tsuneda, Density Functional Theory in Quantum Chemistry (Springer, Tokyo, 2014).
[6] J. P. Perdew and K. Schmidt, AIP Conf. Proc. 577, 1 (2001).
[7] G. E. Scuseria and V. N. Staroverov, in Theory and Applications of Computational Chem-
istry: The First 40 years, edited by C. E. Dykstra, G. Frenking, K. S. Kim and G. E.
Scuseria (Elsevier, 2005), pp. 669–724.
56
[8] K. Capelle, Braz. J. Phys. 36, 1318 (2006).
[9] A. J. Cohen, P. Mori-Sánchez and W. Yang, Chem. Rev. 112, 289 (2012).
[10] K. Burke, J. Chem. Phys. 136, 150901 (2012).
[11] V. N. Staroverov, in A Matter of Density: Exploring the Electron Density Concept in the
Chemical, Biological, and Materials Sciences, edited by N. Sukumar (John Wiley & Sons,
Hoboken, NJ, 2013), pp. 125–156.
[12] A. D. Becke, J. Chem. Phys. 140, 18A301 (2014).
[13] R. O. Jones, Rev. Mod. Phys. 87, 897 (2015).
[14] H. S. Yu, S. L. Li and D. G. Truhlar, J. Chem. Phys. 145, 130901 (2016).
[15] N. Mardirossian and M. Head-Gordon, Mol. Phys. 115, 2315 (2017).
[16] P. Hohenberg and W. Kohn, Phys. Rev. 136, B 864 (1964).
[17] M. Levy, Proc. Natl. Acad. Sci. U.S.A. 76, 6062 (1979).
[18] E. H. Lieb, Int. J. Quantum Chem. 24, 243 (1983).
[19] W. Kohn and L. J. Sham, Phys. Rev. 140, A1133 (1965).
[20] T. L. Gilbert, Phys. Rev. B 12, 2111 (1975).
[21] J. E. Harriman, Phys. Rev. A 24, 680 (1981).
[22] H. Englisch and R. Englisch, Phys. Stat. Sol. 123, 711 (1984).
[23] H. Englisch and R. Englisch, Phys. Stat. Sol. 124, 373 (1984).
[24] I. Lindgren and S. Salomonson, Adv. Quantum Chem. 43, 95 (2003).
[25] P. E. Lammert, Int. J. Quantum Chem. 107, 1943 (2007).
[26] S. Kvaal, U. Ekström, A. M. Teale and T. Helgaker, J. Chem. Phys. 140, 18A518 (2014).
[27] T. Helgaker, Principles of Density-Functional Theory, Lecture at the GdR CORREL Mini-
School on Mathematics in Electronic Structure Theory, Université Pierre et Marie Curie,
Paris, 2017, https://trygvehelgaker.no/Presentations/Paris_2017.pdf.
[28] A. Laestadius, M. Penz, E. I. Tellgren, M. Ruggenthaler, S. Kvaal and T. Helgaker, J.
Chem. Phys. 149, 164103 (2018).
[29] T. Grabo, T. Kreibich, S. Kurth and E. K. U. Gross, in Strong Coulomb Correlation in
Electronic Structure: Beyond the Local Density Approximation, edited by V. Anisimov
(Gordon & Breach, Tokyo, 2000).
[30] C.-O. Almbladh and U. von Barth, Phys. Rev. B 31, 3231 (1985).
[31] U. von Barth and L. Hedin, J. Phys. C 5, 1629 (1972).
[32] A. K. Rajagopal and J. Callaway, Phys. Rev. B 7, 1912 (1973).
[33] J. P. Perdew and A. Zunger, Phys. Rev. B 23, 5048 (1981).
[34] G. L. Oliver and J. P. Perdew, Phys. Rev. A 20, 397 (1979).
[35] A. Seidl, A. Görling, P. Vogl, J. A. Majewski and M. Levy, Phys. Rev. B 53, 3764 (1996).
[36] N. H. March, Phys. Rev. A 36, 5077 (1987).
[37] A. D. Becke, Phys. Rev. A 38, 3098 (1988).
[38] J. C. Kimball, Phys. Rev. A 7, 1648 (1973).
[39] D. P. Tew, J. Chem. Phys. 129, 014104 (2008).
[40] D. C. Langreth and J. P. Perdew, Solid State Commun. 17, 1425 (1975).
[41] O. Gunnarsson and B. I. Lundqvist, Phys. Rev. B 13, 4274 (1976).
[42] D. C. Langreth and J. P. Perdew, Phys. Rev. B 15, 2884 (1977).
[43] J. Harris and R. O. Jones, J. Phys. F 4, 1170 (1974).
57
[44] M. Levy and J. P. Perdew, Phys. Rev. A 32, 2010 (1985).
[45] M. Levy, Phys. Rev. A 43, 4637 (1991).
[46] M. Levy, in Density Functional Theory, edited by E. Gross and R. Dreizler (Plenum Press,
New York, 1995).
[47] A. Görling and M. Levy, Phys. Rev. B 47, 13105 (1993).
[48] A. Görling and M. Levy, Phys. Rev. A 50, 196 (1994).
[49] C.-J. Huang and C. J. Umrigar, Phys. Rev. A 56, 290 (1997).
[50] M. Seidl, J. P. Perdew and M. Levy, Phys. Rev. A 59, 51 (1999).
[51] M. Seidl, Phys. Rev. A 60, 4387 (1999).
[52] M. Seidl, P. Gori-Giorgi and A. Savin, Phys. Rev. A 75, 042511 (2007).
[53] P. Gori-Giorgi and M. Seidl, Phys. Chem. Chem. Phys. 12, 14405 (2010).
[54] H. Ou-Yang and M. Levy, Phys. Rev. A 42, 155 (1990).
[55] M. Levy and H. Ou-Yang, Phys. Rev. A 42, 651 (1990).
[56] A. Görling and M. Levy, Phys. Rev. A 45, 1509 (1992).
[57] M. Levy and J. P. Perdew, Phys. Rev. B 48, 11638 (1993).
[58] S. Ivanov and M. Levy, J. Phys. Chem. A 102, 3151 (1998).
[59] V. N. Staroverov, G. E. Scuseria, J. P. Perdew, J. Tao and E. R. Davidson, Phys. Rev. A
70, 012502 (2004).
[60] G. Friesecke and B. D. Goddard, SIAM J. Math. Analysis 41, 631 (2009).
[61] G. Friesecke and B. D. Goddard, Phys. Rev. A 81, 032516 (2010).
[62] A. D. Kaplan, B. Santra, P. Bhattarai, K. Wagle, S. T. u. R. Chowdhury, P. Bhetwal,
J. Yu, H. Tang, K. Burke, M. Levy and J. P. Perdew, J. Chem. Phys. 153, 074114 (2020).
[63] E. H. Lieb and B. Simon, Phys. Rev. Lett. 31, 681 (1973).
[64] E. H. Lieb and B. Simon, Adv. Math. 23, 22 (1977).
[65] P. Elliott and K. Burke, Can. J. Chem. 87, 1485 (2009).
[66] K. Burke, A. Cancio, T. Gould and S. Pittalis, J. Chem. Phys. 145, 054112 (2016).
[67] A. Cancio, G. P. Chen, B. T. Krull and K. Burke, J. Chem. Phys. 149, 084116 (2018).
[68] T. J. Daas, D. P. Kooi, A. J. A. F. Grooteman, M. Seidl and P. Gori-Giorgi, J. Chem.
Theory Comput. 18, 1584 (2022).
[69] N. Argaman, J. Redd, A. C. Cancio and K. Burke, arxiv.org/abs/2204.01030 (2022).
[70] E. H. Lieb and S. Oxford, Int. J. Quantum. Chem. 19, 427 (1981).
[71] J. P. Perdew, in Electronic Structure of Solids ’91 , edited by P. Ziesche and H. Eschrig
(Akademie Verlag, berlin, 1991).
[72] G. K.-L. Chan and N. C. Handy, Phys. Rev. A 59, 3075 (1999).
[73] C. Cotar and M. Petrache, arxiv.org/abs/1707.07664 (2019).
[74] M. Lewin, E. H. Lieb and R. Seiringer, arxiv.org/abs/2203.12473 (2022).
[75] M. M. Odashima and K. Capelle, J. Chem. Phys. 127, 054106 (2007).
[76] S. R. Gadre, L. J. Bartolotti and N. C. Handy, J. Chem. Phys. 72, 1034 (1980).
[77] J. P. Perdew, A. Ruzsinszky, J. Sun and K. Burke, J. Chem. Phys. 140, 18A533 (2014).
[78] P. A. M. Dirac, Proc. Cambridge Phil. Soc. 26, 376 (1930).
[79] J. C. Slater, Phys. Rev. 81, 385 (1951).
[80] G. Friesecke, Commun. Math. Phys. 184, 143 (1997).
[81] P.-F. Loos and P. M. W. Gill, WIREs Comput. Mol. Sci. 6, 410 (2016), doi:
10.1002/wcms.1257.
58
[82] S. J. Vosko, L. Wilk and M. Nusair, Can. J. Phys. 58, 1200 (1980).
[83] J. P. Perdew and Y. Wang, Phys. Rev. B 45, 13244 (1992).
[84] D. M. Ceperley and B. J. Alder, Phys. Rev. Lett. 45, 566 (1980).
[85] J. P. Perdew, L. A. Constantin, E. Sagvolden and K. Burke, Phys. Rev. Lett. 97, 223002
(2006).
[86] J. P. Perdew, E. R. McMullen and A. Zunger, Phys. Rev. A 23, 2785 (1981).
[87] S.-K. Ma and K. A. Brueckner, Phys. Rev. 165, 165 (1968).
[88] L. Kleinman and S. Lee, Phys. Rev. B 37, 4634 (1988).
[89] E. Engel and S. H. Vosko, Phys. Rev. B 42, 4940 (1990).
[90] P. S. Svendsen and U. von Barth, Phys. Rev. B 54, 17402 (1996).
[91] L. J. Sham, in Computational Methods in Band Theory, edited by P. Marcus, J. F. Janak
and A. R. Williams (Plenum, New York, 1971).
[92] D. C. Langreth and S. H. Vosko, Phys. Rev. Lett. 59, 497 (1987).
[93] M. Rasolt and D. J. W. Geldart, Phys. Rev. B 34, 1325 (1986).
[94] M. Rasolt, Phys. Rev. B 16, 3234 (1977).
[95] Y. Wang and J. Perdew, Phys. Rev. B 44, 13298 (1991).
[96] M. Rasolt and H. L. Davis, Phys. Lett. 86A, 45 (1981).
[97] E. Engel, J. A. Chevary, L. D. Macdonald and S. H. Vosko, Z. Phys. D 23, 7 (1992).
[98] C. Lee, W. Yang and R. G. Parr, Phys. Rev. B 37, 785 (1988).
[99] R. Colle and O. Salvetti, Theor. Chim. Acta 37, 329 (1975).
[100] B. Miehlich, A. Savin, H. Stoll and H. Preuss, Chem. Phys. Lett. 157, 200 (1989).
[101] J. P. Perdew, J. A. Chevary, S. H. Vosko, K. A. Jackson, M. R. Pederson, D. J. Singh and
C. Fiolhais, Phys. Rev. B 46, 6671 (1992).
[102] K. Burke, J. P. Perdew and Y. Wang, Derivation of a generalized gradient approximation:
The PW91 density functional (Plenum, NY, 1997), p. 81.
[103] J. P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett. 77, 3865 (1996).
[104] A. D. Becke, J. Chem. Phys. 84, 4524 (1986).
[105] Y. Zhang and W. Yang, Phys. Rev. Lett. 80, 890 (1998).
[106] B. Hammer, L. B. Hansen and J. K. Nørskov, Phys. Rev. B 59, 7413 (1999).
[107] J. P. Perdew, A. Ruzsinszky, G. I. Csonka, O. A. Vydrov, G. E. Scuseria, L. A. Constantin,
X. Zhou and K. Burke, Phys. Rev. Lett. 100, 136406 (2008).
[109] H. Stoll, C. Pavlidou and H. Preuss, Theor. Chim. Acta 49, 143 (1978).
[110] H. Stoll, E. Golka and H. Preuss, Theor. Chim. Acta 55, 29 (1980).
[111] F. A. Hamprecht, A. J. Cohen, D. J. Tozer and N. C. Handy, J. Chem. Phys. 109, 6264
(1998).
[112] R. Neumann, R. H. Nobes and N. C. Handy, Mol. Phys. 87, 1 (1996).
[113] R. Neumann and N. C. Handy, Chem. Phys. Lett. 252, 19 (1996).
[114] C. Adamo, M. Ernzerhof and G. E. Scuseria, J. Chem. Phys. 112, 2643 (2000).
[115] A. V. Arbuznikov, M. Kaupp, V. G. Malkin, R. Reviakine and O. L. Malkina, Phys. Chem.
Chem. Phys. 4, 5467 (2002).
[116] F. Furche and J. P. Perdew, J. Chem. Phys. 124, 044103 (2006).
[117] J. Sun, M. Marsman, G. I. Csonka, A. Ruzsinszky, P. Hao, Y.-S. Kim and G. K. J. P.
Perdew, Phys. Rev. B 84, 035117 (2011).
59
[118] F. Zahariev, S. S. Leang and M. S. Gordon, J. Chem. Phys. 138, 244108 (2013).
[119] S. M. O. Souvi, K. Sharkas and J. Toulouse, J. Chem. Phys. 140, 084107 (2014).
[120] F. G. Eich and M. Hellgren, J. Chem. Phys. 141, 224107 (2014).
[121] A. D. Becke, Int. J. Quantum. Chem. 23, 1915 (1983).
[122] J. Sun, B. Xiao, Y. Fang, R. Haunschild, P. Hao, A. Ruzsinszky, G. I. Csonka, G. E.
Scuseria and J. P. Perdew, Phys. Rev. Lett. 111, 106401 (2013).
[123] S. Kurth, J. P. Perdew and P. Blaha, Int. J. Quantum Chem. 75, 889 (1999).
[124] M. Brack, B. K. Jennings and Y. H. Chu, Phys. Lett. 65A, 1 (1976).
[125] J. Tao, J. P. Perdew, V. N. Staroverov and G. E. Scuseria, Phys. Rev. Lett. 91, 146401
(2003).
[126] J. P. Perdew, J. Tao, V. N. Staroverov and G. E. Scuseria, J. Chem. Phys. 120, 6898
(2004).
[127] J. P. Perdew, S. Kurth, A. Zupan and P. Blaha, Phys. Rev. Lett. 82, 2544 (1999).
[128] N. D. Lang and W. Kohn, Phys. Rev. B 1, 4555 (1970).
[129] Y. Zhao and D. G. Truhlar, J. Chem. Phys. 125, 194101 (2006).
[130] Y. Zhao, N. E. Schultz and D. G. Truhlar, J. Chem. Phys. 123, 161103 (2005).
[132] T. Van Voorhis and G. E. Scuseria, J. Chem. Phys. 109, 400 (1998).
[134] J. Sun, A. Ruzsinszky and J. P. Perdew, Phys. Rev. Lett. 115, 036402 (2015).
[135] J. P. Perdew, A. Ruzsinszky, G. I. Csonka, L. A. Constantin and J. Sun, Phys. Rev. Lett.
103, 026403 (2009).
[136] J. Sun, J. P. Perdew, Z. Yang and H. Peng, J. Chem. Phys. 144, 191101 (2016).
[140] J. P. Perdew, M. Ernzerhof and K. Burke, J. Chem. Phys 105, 9982 (1996).
[141] J. Jaramillo, G. E. Scuseria and M. Ernzerhof, J. Chem. Phys. 118, 1068 (2003).
[142] T. M. Maier, A. V. Arbuznikov and M. Kaupp, WIREs Comput. Mol. Sci. 9, e1378 (2019).
[143] P. J. Stephens, F. J. Devlin, C. F. Chabalowski and M. J. Frisch, J. Phys. Chem. 98, 11623
(1994).
[144] C. Adamo and V. Barone, J. Chem. Phys. 110, 6158 (1999).
[145] M. Ernzerhof and G. E. Scuseria, J. Chem. Phys. 110, 5029 (1999).
[146] V. N. Staroverov, G. E. Scuseria, J. Tao and J. P. Perdew, J. Chem. Phys. 119, 12129
(2003).
[147] Y. Zhao and D. G. Truhlar, Theor. Chem. Acc. 120, 215 (2008).
[148] A. Savin, in Recent Developments of Modern Density Functional Theory, edited by J. M.
Seminario (Elsevier, Amsterdam, 1996), pp. 327–357.
[149] H. Iikura, T. Tsuneda, T. Yanai and K. Hirao, J. Chem. Phys. 115, 3540 (2001).
[150] I. C. Gerber and J. G. Ángyán, Chem. Phys. Lett. 415, 100 (2005).
[151] P. M. W. Gill, R. D. Adamson and J. A. Pople, Mol. Phys. 88, 1005 (1996).
[152] J. Toulouse, F. Colonna and A. Savin, Phys. Rev. A 70, 062505 (2004).
[153] B. Mussard and J. Toulouse, Mol. Phys. 115, 161 (2017).
60
[154] T. Yanai, D. P. Tew and N. C. Handy, Chem. Phys. Lett. 393, 51 (2004).
[155] A. Savin and H.-J. Flad, Int. J. Quantum. Chem. 56, 327 (1995).
[156] T. M. Henderson, A. F. Izmaylov, G. E. Scuseria and A. Savin, J. Chem. Phys. 127, 221103
(2007).
[157] T. Stein, L. Kronik and R. Baer, J. Am. Chem. Soc. 131, 2818 (2009).
[158] T. Stein, L. Kronik and R. Baer, J. Chem. Phys. 131, 244119 (2009).
[159] R. Baer, E. Livshits and U. Salzner, Annu. Rev. Phys. Chem. 61, 85 (2010).
[160] A. Karolewski, L. Kronik and S. Kümmel, J. Chem. Phys. 138, 204115 (2013).
[161] A. V. Krukau, G. E. Scuseria, J. P. Perdew and A. Savin, J. Chem. Phys. 129, 124103
(2008).
[162] T. Aschebrock and S. Kümmel, J. Chem. Phys. 151, 154108 (2019).
[163] S. Klawohn and H. Bahmann, J. Chem. Theory Comput. 16, 953 (2020).
[164] Y. Tawada, T. Tsuneda, S. Yanagisawa, T. Yanai and K. Hirao, J. Chem. Phys. 120, 8425
(2004).
[165] O. A. Vydrov, J. Heyd, A. V. Krukau and G. E. Scuseria, J. Chem. Phys. 125, 074106
(2006).
[166] O. A. Vydrov and G. E. Scuseria, J. Chem. Phys. 125, 234109 (2006).
[167] M. Ernzerhof and J. P. Perdew, J. Chem. Phys. 109, 3313 (1998).
[168] J.-D. Chai and M. Head-Gordon, J. Chem. Phys. 128, 084106 (2008).
[169] J. Heyd, G. E. Scuseria and M. Ernzerhof, J. Chem. Phys. 118, 8207 (2003).
[170] S. Grimme, J. Chem. Phys. 124, 034108 (2006).
[171] A. Szabo and N. S. Ostlund, Modern Quantum Chemistry: Introduction to Advanced
Electronic Structure Theory (Dover, New York, 1996).
[172] K. Sharkas, J. Toulouse and A. Savin, J. Chem. Phys. 134, 064113 (2011).
[173] J. G. Ángyán, I. C. Gerber, A. Savin and J. Toulouse, Phys. Rev. A 72, 012510 (2005).
[174] E. Fromager and H. J. A. Jensen, Phys. Rev. A 78, 022504 (2008).
[175] J. G. Ángyán, Phys. Rev. A 78, 022510 (2008).
[176] J. Toulouse, K. Sharkas, E. Brémond and C. Adamo, J. Chem. Phys. 135, 101102 (2011).
[177] E. Fromager, J. Chem. Phys. 135, 244106 (2011).
[178] S. Grimme, J. Chem. Phys. 118, 9095 (2003).
[179] S. Kozuch, D. Gruzman and J. M. L. Martin, J. Phys. Chem. C 114, 20801 (2010).
[180] S. Kozuch and J. M. L. Martin, Phys. Chem. Chem. Phys. 13, 20104 (2011).
[181] J. C. Sancho-Garcı́a and C. Adamo, Phys. Chem. Chem. Phys. 15, 14581 (2013).
[182] L. Goerigk and S. Grimme, WIREs Comput. Mol. Sci. 4, 576 (2014).
[183] N. Q. Su and X. Xu, WIREs Comput. Mol. Sci. 6, 721 (2016).
[184] N. Mehta, M. Casanova-Páez and L. Goerigk, Phys. Chem. Chem. Phys. 20, 23175 (2018).
[185] A. Ruzsinszky, J. P. Perdew and G. I. Csonka, J. Chem. Theory Comput. 6, 127 (2010).
[186] S. Ahnen, A.-S. Hehn, K. D. Vogiatzis, M. A. Trachsel, S. Leutwyler and W. Klopper,
Chem. Phys. 441, 17 (2014).
[187] P. D. Mezei, G. I. Csonka, A. Ruzsinszky and M. Kállay, J. Chem. Theory Comput. 11,
4615 (2015).
[188] S. Grimme and M. Steinmetz, Phys. Chem. Chem. Phys. 18, 20926 (2016).
[189] P. D. Mezei, G. I. Csonka, A. Ruzsinszky and M. Kállay, J. Chem. Theory Comput. 13,
796 (2017).
61
[190] A. J. Garza, I. W. Bulik, T. M. Henderson and G. E. Scuseria, J. Chem. Phys. 142, 044109
(2015).
[191] B. Chan, L. Goerigk and L. Radom, J. Comput. Chem. 37, 183 (2016).
[192] K. Sharkas, A. Savin, H. J. A. Jensen and J. Toulouse, J. Chem. Phys. 137, 044104 (2012).
[193] F. Ying, C. Zhou, P. Zheng, J. Luan, P. Su and W. Wu, Front. Chem. 7, 225 (2019).
[194] M. Mostafanejad, M. D. Liebenthal, and A. E. DePrince III, J. Chem. Theory Comput.
16, 2274 (2020).
[195] A. Savin, in Recent Advances in Density Functional Theory, edited by D. P. Chong (World
Scientific, 1996), pp. 129–148.
[196] J. Toulouse, Ph.D. thesis, Université Pierre et Marie Curie (Paris 6) (2005), tel.archives-
ouvertes.fr/tel-00550772.
[197] J. Toulouse, P. Gori-Giorgi and A. Savin, Theor. Chem. Acc. 114, 305 (2005).
[198] A. Stoyanova, A. M. Teale, J. Toulouse, T. Helgaker and E. Fromager, J. Chem. Phys.
139, 134113 (2013).
[199] B. Mussard, P. Reinhardt, J. G. Ángyán and J. Toulouse, J. Chem. Phys. 142, 154123
(2015), Erratum: J. Chem. Phys. 142, 219901 (2015).
[200] I. C. Gerber and J. G. Ángyán, Chem. Phys. Lett. 416, 370 (2005).
[201] I. C. Gerber and J. G. Ángyán, J. Chem. Phys. 126, 044103 (2007).
[202] D. E. Taylor, J. G. Ángyán, G. Galli, C. Zhang, F. Gygi, K. Hirao, J. W. Song, K. Rahul,
O. A. von Lilienfeld, R. Podeszwa, I. W. Bulik, T. M. Henderson, G. E. Scuseria,
J. Toulouse, R. Peverati, D. G. Truhlar and K. Szalewicz, J. Chem. Phys. 145, 124105
(2016).
[203] O. Franck, B. Mussard, E. Luppi and J. Toulouse, J. Chem. Phys. 142, 074107 (2015).
[204] G. Sansone, B. Civalleri, D. Usvyat, J. Toulouse, K. Sharkas and L. Maschio, J. Chem.
Phys. 143, 102811 (2015).
[205] E. Goll, H.-J. Werner and H. Stoll, Phys. Chem. Chem. Phys. 7, 3917 (2005).
[206] E. Goll, H.-J. Werner, H. Stoll, T. Leininger, P. Gori-Giorgi and A. Savin, Chem. Phys.
329, 276 (2006).
[207] J. Toulouse, W. Zhu, A. Savin, G. Jansen and J. G. Ángyán, J. Chem. Phys. 135, 084119
(2011).
[208] A. J. Garza, I. W. Bulik, T. M. Henderson and G. E. Scuseria, Phys. Chem. Chem. Phys.
17, 22412 (2015).
[209] J. Toulouse, I. C. Gerber, G. Jansen, A. Savin and J. G. Ángyán, Phys. Rev. Lett. 102,
096404 (2009).
[210] B. G. Janesko, T. M. Henderson and G. E. Scuseria, J. Chem. Phys. 130, 081105 (2009).
[211] J. Toulouse, W. Zhu, J. G. Ángyán and A. Savin, Phys. Rev. A 82, 032502 (2010).
[212] J. Paier, B. G. Janesko, T. M. Henderson, G. E. Scuseria, A. Grüneis and G. Kresse, J.
Chem. Phys. 132, 094103 (2010).
[213] T. Leininger, H. Stoll, H.-J. Werner and A. Savin, Chem. Phys. Lett. 275, 151 (1997).
[214] R. Pollet, A. Savin, T. Leininger and H. Stoll, J. Chem. Phys. 116, 1250 (2002).
[215] D. Casanova, J. Chem. Phys. 148, 124118 (2018).
[216] A. Ferté, E. Giner and J. Toulouse, J. Chem. Phys. 150, 084103 (2019).
[217] E. Fromager, J. Toulouse and H. J. A. Jensen, J. Chem. Phys. 126, 074111 (2007).
[218] E. Fromager, F. Réal, P. Wåhlin, U. Wahlgren and H. J. A. Jensen, J. Chem. Phys. 131,
054107 (2009).
62
[219] E. D. Hedegård, J. Toulouse and H. J. A. Jensen, J. Chem. Phys. 148, 214103 (2018).
[220] E. D. Hedegård, S. Knecht, J. S. Kielberg, H. J. A. Jensen and M. Reiher, J. Chem. Phys.
142, 224108 (2015).
[221] E. Fromager, R. Cimiraglia and H. J. A. Jensen, Phys. Rev. A 81, 024502 (2010).
[222] K. Pernal, Phys. Rev. A 81, 052511 (2010).
[223] D. R. Rohr, J. Toulouse and K. Pernal, Phys. Rev. A 82, 052502 (2010).
[224] D. R. Rohr and K. Pernal, J. Chem. Phys. 135, 074104 (2011).
[225] E. Rebolini and J. Toulouse, J. Chem. Phys. 144, 094107 (2016).
[226] A. A. Kananenka and D. Zgid, J. Chem. Theory Comput. 13, 5317 (2017).
[227] J. Toulouse, A. Savin and H.-J. Flad, Int. J. Quantum Chem. 100, 1047 (2004).
[228] S. Paziani, S. Moroni, P. Gori-Giorgi and G. B. Bachelet, Phys. Rev. B 73, 155111 (2006).
[229] P. Gori-Giorgi and J. P. Perdew, Phys. Rev. B 64, 155102 (2001).
[230] J. Toulouse, F. Colonna and A. Savin, J. Chem. Phys. 122, 014110 (2005).
[231] C. Kalai and J. Toulouse, J. Chem. Phys. 148, 164105 (2018).
[232] Y. Cornaton and E. Fromager, Int. J. Quantum Chem. 114, 1199 (2014).
[233] C. Kalai, B. Mussard and J. Toulouse, J. Chem. Phys. 151, 074102 (2019).
[234] M. Elstner, P. Hobza, T. Frauenheim, S. Suhai and E. Kaxiras, J. Chem. Phys. 114, 5149
(2001).
[235] Q. Wu and W. Yang, J. Chem. Phys. 116, 515 (2002).
[236] S. Grimme, J. Comput. Chem. 25, 1463 (2004).
[237] S. Grimme, J. Comput. Chem. 27, 1787 (2006).
[238] S. Grimme, J. Antony, S. Ehrlich and H. Krieg, J. Chem. Phys. 132, 154104 (2010).
[239] A. D. Becke and E. R. Johnson, J. Chem. Phys. 127, 154108 (2007).
[240] A. Tkatchenko and M. Scheffler, Phys. Rev. Lett. 102, 073005 (2009).
[241] T. Sato and H. Nakai, J. Chem. Phys. 133, 194101 (2010).
[242] J. F. Dobson and B. P. Dinte, Phys. Rev. Lett. 76, 1780 (1996).
[243] J.-D. Chai and M. Head-Gordon, Phys. Chem. Chem. Phys. 10, 6615 (2008).
[244] T. Schwabe and S. Grimme, Phys. Chem. Chem. Phys. 9, 3397 (2007).
[245] S. Kozuch and J. M. L. Martin, J. Comput. Chem. 34, 2327 (2013).
[246] M. Dion, H. Rydberg, E. Schröder, D. C. Langreth and B. I. Lundqvist, Phys. Rev. Lett.
92, 246401 (2004).
[247] O. A. Vydrov and T. Van Voorhis, J. Chem. Phys. 130, 104105 (2009).
[248] O. A. Vydrov and T. Van Voorhis, Phys. Rev. Lett. 103, 063004 (2009).
[249] K. Lee, E. D. Murray, L. Kong, B. I. Lundqvist and D. C. Langreth, Phys. Rev. B 82,
081101 (2010).
[250] O. A. Vydrov and T. Van Voorhis, J. Chem. Phys. 133, 244103 (2010).
[251] N. Mardirossian and M. Head-Gordon, Phys. Chem. Chem. Phys. 16, 9904 (2014).
[252] N. Mardirossian and M. Head-Gordon, J. Chem. Phys. 144, 214110 (2016).
[253] S. Kümmel and L. Kronik, Rev. Mod. Phys. 80, 3 (2008).
[254] A. V. Arbuznikov, J. Struct. Chem. 48, S1 (2007).
[255] J. Kim, K. Hong, S.-Y. Hwang, S. Ryu, S. Choi and W. Y. Kim, Phys. Chem. Chem.
Phys. 19, 10177 (2017).
[256] S. Śmiga and L. A. Constantin, J. Phys. Chem. A 124, 5606 (2020).
63
[257] S. Śmiga, O. Franck, B. Mussard, A. Buksztel, I. Grabowski, E. Luppi and J. Toulouse,
J. Chem. Phys. 145, 144102 (2016).
[258] S. Śmiga, I. Grabowski, M. Witkowski, B. Mussard and J. Toulouse, J. Chem. Theory
Comput. 16, 211 (2020).
[259] J. D. Talman and W. F. Shadwick, Phys. Rev. A 14, 36 (1976).
[260] A. Görling and M. Levy, Int. J. Quantum Chem. Symp. 56, 93 (1995).
[261] R. T. Sharp and G. K. Horton, Phys. Rev. 90, 317 (1953).
[262] S. Ivanov, S. Hirata and R. J. Bartlett, Phys. Rev. Lett. 83, 5455 (1999).
[263] A. Görling, Phys. Rev. Lett. 83, 5459 (1999).
[264] I. Grabowski, S. Hirata, S. Ivanov and R. J. Bartlett, J. Chem. Phys. 116, 4415 (2002).
[265] E. Engel, in A Primer in Density Functional Theory, edited by C. Fiolhais, F. Nogueira
and M. A. L. Marques (Springer, Berlin, 2003), Vol. 620 of Lecture Notes in Physics, pp.
56–122.
[266] P. Mori-Sánchez, Q. Wu and W. Yang, J. Chem. Phys. 123, 062204 (2005).
[267] D. C. Langreth and J. P. Perdew, Phys. Rev. B 21, 5469 (1980).
[268] B. Mussard, D. Rocca, G. Jansen and J. G. Ángyán, J. Chem. Theory Comput. 12, 2191
(2016).
[269] F. Furche, Phys. Rev. B 64, 195120 (2001).
[270] G. E. Scuseria, T. M. Henderson and D. C. Sorensen, J. Chem. Phys. 129, 231101 (2008).
[271] M. Hellgren, D. R. Rohr and E. K. U. Gross, J. Chem. Phys. 136, 034106 (2012).
[272] P. Bleiziffer, A. Heßelmann and A. Görling, J. Chem. Phys. 139, 084113 (2013).
[273] J. F. Dobson, K. McLennan, A. Rubio, J. Wang, T. Gould, H. M. Le and B. P. Dinte,
Aust. J. Chem. 54, 513 (2001).
[274] J. F. Dobson, J. Wang, B. P. Dinte, K. McLennan and H. M. Le, Int. J. Quantum. Chem.
101, 579 (2005).
[275] J. F. Dobson and T. Gould, J. Phys.: Condens. Matter 24, 073201 (2012).
[276] A. Grüneis, M. Marsman, J. Harl, L. Schimka and G. Kresse, J. Chem. Phys. 131, 154115
(2009).
[277] M. Hellgren and U. von Barth, J. Chem. Phys. 132, 044101 (2010).
[278] A. Heßelmann and A. Görling, Mol. Phys. 108, 359 (2010).
[279] J. G. Ángyán, R.-F. Liu, J. Toulouse and G. Jansen, J. Chem. Theory Comput. 7, 3116
(2011).
[280] A. Heßelmann and A. Görling, Phys. Rev. Lett. 106, 093001 (2011).
[281] A. Heßelmann, Phys. Rev. A 012517, 85 (2012).
[282] J. E. Bates and F. Furche, J. Chem. Phys. 139, 171103 (2013).
[283] N. Colonna, M. Hellgren, and S. de Gironcoli, Phys. Rev. B 90, 125150 (2014).
[284] J. Erhard, P. Bleiziffer and A. Görling, Phys. Rev. Lett. 117, 143002 (2016).
[285] M. Hellgren, N. Colonna, and S. de Gironcoli, Phys. Rev. B 98, 045117 (2018).
[286] J. E. Bates, N. Sengupta, J. Sensenig and A. Ruzsinszky, J. Chem. Theory Comput. 14,
2979 (2018).
[287] F. Hummel, A. Grüneis, G. Kresse and P. Ziesche, J. Chem. Theory Comput. 15, 3223
(2019).
[288] A. Görling, Phys. Rev. B 99, 235120 (2019).
[289] H. Eshuis, J. Bates and F. Furche, Theor. Chem. Acc. 131, 1084 (2012).
64
[290] X. Ren, P. Rinke, C. Joas and M. Scheffler, J. Mater. Sci. 47, 7447 (2012).
[291] G. P. Chen, V. K. Voora, M. M. Agee, S. G. Balasubramani and F. Furche, Annu. Rev.
Phys. Chem. 68, 421 (2017).
65

Review of Approximations For The Exchange-Correlation Energy in Density-Functional Theory

Uploaded by

Copyright:

Available Formats

Review of Approximations For The Exchange-Correlation Energy in Density-Functional Theory

Uploaded by

Document Information

Original Description:

Original Title

Copyright

Available Formats

Share this document

Share or Embed Document

Sharing Options

Did you find this document useful?

Is this content inappropriate?

Copyright:

Available Formats

Review of Approximations For The Exchange-Correlation Energy in Density-Functional Theory

Uploaded by

Copyright:

Available Formats

Chapter in the book Density Functional Theory edited by Éric Cancès and Gero Friesecke.

Review of approximations for the exchange-correlation energy in

Institut Universitaire de France, F-75005 Paris, France

In this chapter, we provide a review of ground-state Kohn–Sham density-functional theory

4 Single-determinant hybrid approximations 34

5 Multideterminant hybrid approximations 40

7 Orbital-dependent exchange-correlation density functionals 51

1 Basics of density-functional theory

HΨ(x1 , x2 , ..., xN ) = EΨ(x1 , x2 , ..., xN ), (1.2)

Ĥ = T̂ + Ŵee + V̂ne , (1.4)

1.2 The universal density functional

F [ρ] = T [ρ] + Wee [ρ], (1.12)

1.3 The Kohn–Sham scheme

F [ρ] = Ts [ρ] + EHxc [ρ], (1.13)

EHxc [ρ] = EH [ρ] + Exc [ρ], (1.18)

where EH [ρ] is the Hartree energy functional,

Ex [ρ] = hΦ[ρ]|Ŵee |Φ[ρ]i − EH [ρ], (1.21)

and Ec [ρ] is the correlation energy functional,

where γσ , for σ ∈ {↑, ↓}, is the spin-dependent one-particle KS density matrix,

where the density is expressed in terms of the orbitals as

vs (r) = vne (r) + vHxc (r) (1.30)

vHxc (r) = vH (r) + vxc (r), (1.31)

δExc [ρ]/δρ(r) is the exchange-correlation potential. Likewise, following the decomposition of

vxc (r) = vx (r) + vc (r), (1.32)

whereas the correlation potential decays faster [30].

1.3.3 Extension to spin density-functional theory

ρ↓ ∈ D N↓ , the universal density functional is now defined as [33],

WρN↑ ,ρ↓ = {Ψ ∈ W N , ρ↑,Ψ = ρ↑ , ρ↓,Ψ = ρ↓ }. (1.36)

In Eq. (1.35), Ψ[ρ↑ , ρ↓ ] designates a minimizing wave function.

F [ρ↑ , ρ↓ ] = Ts [ρ↑ , ρ↓ ] + EH [ρ] + Exc [ρ↑ , ρ↓ ], (1.37)

where Ts [ρ↑ , ρ↓ ] is defined as

Ts [ρ↑ , ρ↓ ] = min hΦ|T̂ |Φi = hΦ[ρ↑ , ρ↓ ]|T̂ |Φ[ρ↑ , ρ↓ ]i, (1.38)

SρN↑ ,ρ↓ = {Φ ∈ S N , ρ↑,Φ = ρ↑ , ρ↓,Φ = ρ↓ }. (1.39)

and the spin density,

Exc [ρ↑ , ρ↓ ] = Ex [ρ↑ , ρ↓ ] + Ec [ρ↑ , ρ↓ ], (1.44)

where S[Φ] is any functional of a single-determinant wave function Φ ∈ S N leading to a mini-

we can express the exact ground-state energy as

2 Exact expressions and constraints for the Kohn–Sham exchange and

Mirroring the decomposition of the Hartree-exchange-correlation energy performed in the KS

ρ2 (r1 , r2 ) = ρ(r1 )ρ(r2 ) + ρ2,xc (r1 , r2 ). (2.3)

The product of the densities ρ(r1 )ρ(r

ρ2,xc (r1 , r2 ) = ρ(r1 )hxc (r1 , r2 ), (2.4)

ρcond (r1 , r2 ) = ρ(r2 ) + hxc (r1 , r2 ). (2.5)

The positivity of ρ2 (r1 , r2 ) implies that

hxc (r1 , r2 ) ≥ −ρ(r2 ). (2.6)

Moreover, from Eq. (2.5), we have the following sum rule:

ρ2,KS (r1 , r2 ) = ρ(r1 )ρ(r2 ) + ρ2,x (r1 , r2 ), (2.8)

ρ2,x (r1 , r2 ) = ρ(r1 )hx (r1 , r2 ), (2.9)

1 ρ(r1 )hx (r1 , r2 )

where εx [ρ](r1 ) is the exchange energy density per particle

The correlation hole is defined as the difference

hc (r1 , r2 ) = hxc (r1 , r2 ) − hx (r1 , r2 ), (2.18)

1 ρ(r1 )hc (r1 , r2 )

F λ [ρ] = Ts [ρ] + EHxc

2.3 One-orbital and one-electron spatial regions

Ex [ρ1e ] = −EH [ρ1e ]. (2.36)

Furthermore, the correlation energy vanishes:

∀r1 , r2 ∈ Ω↑1o , hx (r1 , r2 ) = −ρ(r2 ), (2.41)