# Review of approximations for the exchange-correlation energy in density-
functional theory

Julien Toulouse
Laboratoire de Chimie Théorique (LCT), Sorbonne Université and CNRS, F-75005
Paris, France
Institut Universitaire de France, F-75005 Paris, France

(24 February, 2021)

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 approximations for the
exchange-correlation energy: semilocal approximations, single-determinant
hybrid approximations, multideterminant hybrid approximations, dispersion-
corrected approximations, 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, 2, 3, 4, 5] and the review articles
of Refs. [6, 7, 8, 9, 10, 11, 12, 13, 14, 15].

###### Contents

  1. 1 Basics of density-functional theory
    1. 1.1 The many-body problem
    2. 1.2 The universal density functional
    3. 1.3 The Kohn-Sham scheme
      1. 1.3.1 Decomposition of the universal functional
      2. 1.3.2 The Kohn-Sham equations
      3. 1.3.3 Extension to spin density-functional theory
    4. 1.4 The generalized Kohn-Sham scheme
  2. 2 Exact expressions and constraints for the Kohn-Sham exchange and correlation functionals
    1. 2.1 The exchange and correlation holes
    2. 2.2 The adiabatic connection
    3. 2.3 One-orbital and one-electron spatial regions
    4. 2.4 Coordinate scaling
      1. 2.4.1 Uniform coordinate scaling
      2. 2.4.2 Non-uniform coordinate scaling
    5. 2.5 Atoms in the limit of large nuclear charge
    6. 2.6 Lieb-Oxford lower bound
  3. 3 Semilocal approximations for the exchange-correlation energy
    1. 3.1 The local-density approximation
    2. 3.2 The gradient-expansion approximation
    3. 3.3 Generalized-gradient approximations
    4. 3.4 Meta-generalized-gradient approximations
  4. 4 Single-determinant hybrid approximations
    1. 4.1 Hybrid approximations
    2. 4.2 Range-separated hybrid approximations
  5. 5 Multideterminant hybrid approximations
    1. 5.1 Double-hybrid approximations
    2. 5.2 Range-separated double-hybrid approximations
      1. 5.2.1 Range-separated one-parameter double-hybrid approximations
      2. 5.2.2 Range-separated two-parameter double-hybrid approximations
  6. 6 Semiempirical dispersion corrections and nonlocal van der Waals density functionals
    1. 6.1 Semiempirical dispersion corrections
    2. 6.2 Nonlocal van der Waals density functionals
  7. 7 Orbital-dependent exchange-correlation density functionals
    1. 7.1 Exact exchange
    2. 7.2 Second-order Görling-Levy perturbation theory
    3. 7.3 Adiabatic-connection fluctuation-dissipation approach
      1. 7.3.1 Exact adiabatic-connection fluctuation-dissipation expression
      2. 7.3.2 Random-phase approximations
  8. Acknowledgements

##  1 Basics of density-functional theory

###  1.1 The many-body problem

We consider a $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,

$H=-\frac{1}{2}\sum_{i=1}^{N}\nabla^{2}_{\mathbf{r}_{i}}+\frac{1}{2}\sum_{i=1}^{N}\sum_{\begin{subarray}{c}j=1\\\
i\neq
j\end{subarray}}^{N}\frac{1}{|\mathbf{r}_{i}-\mathbf{r}_{j}|}+\sum_{i=1}^{N}v_{\text{ne}}(\mathbf{r}_{i}),$
(1.1)

where $\nabla^{2}_{\mathbf{r}_{i}}=\Delta_{\mathbf{r}_{i}}$ is the Laplacian
with respect to the electron coordinate $\mathbf{r}_{i}$ and
$v_{\text{ne}}(\mathbf{r}_{i})=-\sum_{\alpha=1}^{N_{\text{n}}}Z_{\alpha}/|\mathbf{r}_{i}-\mathbf{R}_{\alpha}|$
is the nuclei-electron interaction depending of the positions
$\\{\mathbf{R}_{\alpha}\\}$ and charges $\\{Z_{\alpha}\\}$ of the
$N_{\text{n}}$ nuclei. The stationary electronic states are determined by the
time-independent Schrödinger equation

$H\Psi(\mathbf{x}_{1},\mathbf{x}_{2},...,\mathbf{x}_{N})=E\Psi(\mathbf{x}_{1},\mathbf{x}_{2},...,\mathbf{x}_{N}),$
(1.2)

where $\Psi(\mathbf{x}_{1},\mathbf{x}_{2},...,\mathbf{x}_{N})$ is a wave
function written with space-spin coordinates
$\mathbf{x}_{i}=(\mathbf{r}_{i},\sigma_{i})\in\mathbb{R}^{3}\times\\{\uparrow,\downarrow\\}$
(with $\\{\uparrow,\downarrow\\}\cong\mathbb{Z}_{2}$ 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 Schrödinger equation (1.2) can be rewritten in a
representation-independent formalism

$\hat{H}|\Psi\rangle=E|\Psi\rangle,$ (1.3)

where the Hamiltonian is formally written as

$\hat{H}=\hat{T}+\hat{W}_{\text{ee}}+\hat{V}_{\text{ne}},$ (1.4)

with the kinetic-energy operator $\hat{T}$, the electron-electron interaction
operator $\hat{W}_{\text{ee}}$, and the nuclei-electron interaction operator
$\hat{V}_{\text{ne}}$.

The quantity of primary interest is the ground-state energy $E_{0}$. The
variational theorem establishes that $E_{0}$ can be expressed as an infimum

$E_{0}=\inf_{\Psi\in{\cal W}^{N}}\langle\Psi|\hat{H}|\Psi\rangle,$ (1.5)

where the search is over the set of $N$-electron antisymmetric normalized wave
functions $\Psi$ having a finite kinetic energy

${\cal
W}^{N}=\Big{\\{}\Psi\in\bigwedge^{N}L^{2}(\mathbb{R}^{3}\times\\{\uparrow,\downarrow\\};\mathbb{C}),\Psi\in
H^{1}((\mathbb{R}^{3}\times\\{\uparrow,\downarrow\\})^{N};\mathbb{C}),\langle\Psi|\Psi\rangle=1\Big{\\}},$
(1.6)

where $\bigwedge^{N}$ is the $N$-fold antisymmetrized tensor product, $L^{2}$
and $H^{1}$ are the standard Lebesgue 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 its first-order derivative), and
$\langle\Psi|\Psi\rangle$ designates the $L^{2}$ inner product. Density-
functional theory (DFT) is based on a reformulation of the variational theorem
in terms of the one-electron density defined as111In Eq. (1.7), an integration
over a spin coordinate $\sigma$ just means a sum over the two values
$\sigma\in\\{\uparrow,\downarrow\\}$.

$\rho_{\Psi}(\mathbf{r})=N\int_{\\{\uparrow,\downarrow\\}\times(\mathbb{R}^{3}\times\\{\uparrow,\downarrow\\})^{N-1}}\left|\Psi(\mathbf{x},\mathbf{x}_{2},...,\mathbf{x}_{N})\right|^{2}\text{d}\sigma\text{d}\mathbf{x}_{2}...\text{d}\mathbf{x}_{N},$
(1.7)

which is normalized to the electron number,
$\int_{\mathbb{R}^{3}}\rho_{\Psi}(\mathbf{r})\text{d}\mathbf{r}=N$.

###  1.2 The universal density functional

Building on the work on Hohenberg and Kohn [16], Levy [17] and Lieb [18]
proposed to define the following universal density functional $F[\rho]$ using
a constrained-search approach

$F[\rho]=\min_{\Psi\in{\cal
W}^{N}_{\rho}}\langle\Psi|\hat{T}+\hat{W}_{\text{ee}}|\Psi\rangle=\langle\Psi[\rho]|\hat{T}+\hat{W}_{\text{ee}}|\Psi[\rho]\rangle,$
(1.8)

where the minimization is done over the set of $N$-electron wave functions
$\Psi$ yielding the fixed density $\rho$ [via Eq. (1.7)]

${\cal W}^{N}_{\rho}=\\{\Psi\in{\cal W}^{N},\rho_{\Psi}=\rho\\}.$ (1.9)

In Eq. (1.8), for a given density $\rho$, $\Psi[\rho]$ denotes a minimizing
wave function, which is known to exists [18] but is possibly not unique. This
so-called Levy-Lieb functional $F[\rho]$ is defined on the set of
$N$-representable densities [18]

$\displaystyle{\cal D}^{N}$ $\displaystyle=$
$\displaystyle\\{\rho\;|\;\exists\Psi\in{\cal
W}^{N}\;\text{s.t.}\;\rho_{\Psi}=\rho\\}$ (1.10) $\displaystyle=$
$\displaystyle\\{\rho\in L^{1}(\mathbb{R}^{3})\;|\;\rho\geq
0,\int_{\mathbb{R}^{3}}\rho(\mathbf{r})\text{d}\mathbf{r}=N,\sqrt{\rho}\in
H^{1}(\mathbb{R}^{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

$\displaystyle E_{0}$ $\displaystyle=$ $\displaystyle\inf_{\rho\in{\cal
D}^{N}}\left\\{F[\rho]+\int_{\mathbb{R}^{3}}v_{\text{ne}}(\mathbf{r})\rho(\mathbf{r})\text{d}\mathbf{r}\right\\},$
(1.11)

and if a minimizer exists then it is a ground-state density
$\rho_{0}(\mathbf{r})$ for the potential $v_{\text{ne}}(\mathbf{r})$. Hence,
the ground-state energy can be in principle obtained by minimizing over the
density $\rho$, 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 $\Psi$. However, the explicit expression of $F[\rho]$
in terms of the density is not known, and the direct approximations for
$F[\rho]$ that have been tried so far turn out not to be accurate enough.

If there is a unique wave function $\Psi[\rho]$ (up to a phase factor) in Eq.
(1.8), we can define kinetic and potential contributions to $F[\rho]$

$F[\rho]=T[\rho]+W_{\text{ee}}[\rho],$ (1.12)

where $T[\rho]=\langle\Psi[\rho]|\hat{T}|\Psi[\rho]\rangle$ and
$W_{\text{ee}}[\rho]=\langle\Psi[\rho]|\hat{W}_{\text{ee}}|\Psi[\rho]\rangle$.
The kinetic-energy functional $T[\rho]$ 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
directly $F[\rho]$ can be circumvented by decomposing $F[\rho]$ as

$F[\rho]=T_{\text{s}}[\rho]+E_{\text{Hxc}}[\rho],$ (1.13)

where $T_{\text{s}}[\rho]$ is the non-interacting kinetic-energy functional
which can be defined with a constrained search222It is also possible to define
the non-interacting kinetic-energy functional analogously to the Levy-Lieb
functional in Eq. (1.8) by minimizing over wave functions $\Psi\in{\cal
W}^{N}_{\rho}$, i.e. $T_{\text{s},\text{LL}}[\rho]=\min_{\Psi\in{\cal
W}^{N}_{\rho}}\langle\Psi|\hat{T}|\Psi\rangle$ [18]. In this case, the
corresponding minimizing KS wave function can generally be a linear
combination of Slater determinants. However, we often have
$T_{\text{s},\text{LL}}[\rho]=T_{\text{s}}[\rho]$, in particular for densities
$\rho$ that come from a non-interacting ground-state wave function which is
not degenerate. In this chapter, we will normally assume this nondegeneracy
condition.

$T_{\text{s}}[\rho]=\min_{\Phi\in{\cal
S}^{N}_{\rho}}\langle\Phi|\hat{T}|\Phi\rangle=\langle\Phi[\rho]|\hat{T}|\Phi[\rho]\rangle,$
(1.14)

where the minimization is over the set of $N$-electron single-determinant wave
functions $\Phi$ yielding the fixed density $\rho$

${\cal S}^{N}_{\rho}=\\{\Phi\in{\cal S}^{N},\rho_{\Phi}=\rho\\},$ (1.15)

where ${\cal S}^{N}$ is the set of $N$-electron single-determinant wave
functions build from orthonormal spin orbitals

${\cal
S}^{N}=\\{\Phi=\phi_{1}\wedge\phi_{2}\wedge...\wedge\phi_{N}\;|\;\forall
i\;\phi_{i}\in
H^{1}(\mathbb{R}^{3}\times\\{\uparrow,\downarrow\\};\mathbb{C}),\;\forall
i,j\;\langle\phi_{i}|\phi_{j}\rangle=\delta_{i,j}\\},$ (1.16)

where $\phi_{1}\wedge\phi_{2}\wedge...\wedge\phi_{N}$ designates the
normalized $N$-fold antisymmetrized tensor product of $N$ spin orbitals. The
functional $T_{\text{s}}[\rho]$ is defined over the entire set of
$N$-representable densities ${\cal D}^{N}$ since any $N$-representable density
can be obtained from a single-determinant wave function [20, 21, 18]. In Eq.
(1.14), for a given density $\rho$, $\Phi[\rho]$ 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
$E_{\text{Hxc}}[\rho]$ 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 $T_{\text{s}}[\rho]$ by reformulating the minimization over
densities in Eq. (1.11) as a minimization over single-determinant wave
functions $\Phi$

$E_{0}=\inf_{\Phi\in{\cal
S}^{N}}\left\\{\langle\Phi|\hat{T}+\hat{V}_{\text{ne}}|\Phi\rangle+E_{\text{Hxc}}[\rho_{\Phi}]\right\\},$
(1.17)

and if a minimum exists then any minimizing single-determinant wave function
in Eq. (1.17) gives a ground-state density $\rho_{0}(\mathbf{r})$. Thus, the
exact ground-state energy can in principle be obtained 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 $\Phi$ and therefore it still represents a tremendous simplification
over the usual variational theorem involving a multideterminant wave function
$\Psi$. The advantage of Eq. (1.17) over Eq. (1.11) is that a major part of
the kinetic energy can be treated explicitly with the single-determinant wave
function $\Phi$, and only $E_{\text{Hxc}}[\rho]$ needs to be approximated as
an explicit functional of the density.

In practice, $E_{\text{Hxc}}[\rho]$ is decomposed as

$E_{\text{Hxc}}[\rho]=E_{\text{H}}[\rho]+E_{\text{xc}}[\rho],$ (1.18)

where $E_{\text{H}}[\rho]$ is the Hartree energy functional

$E_{\text{H}}[\rho]=\frac{1}{2}\int_{\mathbb{R}^{3}\times\mathbb{R}^{3}}\frac{\rho(\mathbf{r}_{1})\rho(\mathbf{r}_{2})}{|\mathbf{r}_{1}-\mathbf{r}_{2}|}\text{d}\mathbf{r}_{1}\text{d}\mathbf{r}_{2},$
(1.19)

representing the classical electrostatic repulsion energy for the charge
distribution $\rho(\mathbf{r})$ and which is calculated exactly, and
$E_{\text{xc}}[\rho]$ is the exchange-correlation energy functional that
remains to approximate. If there is a unique KS wave function $\Phi[\rho]$ (up
to a phase factor), we can further decompose $E_{\text{xc}}[\rho]$ as

$E_{\text{xc}}[\rho]=E_{\text{x}}[\rho]+E_{\text{c}}[\rho],$ (1.20)

where $E_{\text{x}}[\rho]$ is the exchange energy functional

$\displaystyle E_{\text{x}}[\rho]$ $\displaystyle=$
$\displaystyle\langle\Phi[\rho]|\hat{W}_{\text{ee}}|\Phi[\rho]\rangle-
E_{\text{H}}[\rho],$ (1.21)

and $E_{\text{c}}[\rho]$ is the correlation energy functional

$E_{\text{c}}[\rho]=\langle\Psi[\rho]|\hat{T}+\hat{W}_{\text{ee}}|\Psi[\rho]\rangle-\langle\Phi[\rho]|\hat{T}+\hat{W}_{\text{ee}}|\Phi[\rho]\rangle=T_{\text{c}}[\rho]+U_{\text{c}}[\rho],$
(1.22)

which contains a kinetic contribution
$T_{\text{c}}[\rho]=\langle\Psi[\rho]|\hat{T}|\Psi[\rho]\rangle-\langle\Phi[\rho]|\hat{T}|\Phi[\rho]\rangle$
and a potential contribution
$U_{\text{c}}[\rho]=\langle\Psi[\rho]|\hat{W}_{\text{ee}}|\Psi[\rho]\rangle-\langle\Phi[\rho]|\hat{W}_{\text{ee}}|\Phi[\rho]\rangle$.
Using the fact that $\Phi[\rho]$ is a single-determinant wave function, it can
be shown that the exchange functional can be expressed as

$E_{\text{x}}[\rho]=-\frac{1}{2}\sum_{\sigma\in\\{\uparrow,\downarrow\\}}\int_{\mathbb{R}^{3}\times\mathbb{R}^{3}}\frac{|\gamma_{\sigma}(\mathbf{r}_{1},\mathbf{r}_{2})|^{2}}{|\mathbf{r}_{1}-\mathbf{r}_{2}|}\text{d}\mathbf{r}_{1}\text{d}\mathbf{r}_{2},$
(1.23)

where $\gamma_{\sigma}$, for $\sigma\in\\{\uparrow,\downarrow\\}$, is the
spin-dependent one-particle KS density matrix

$\gamma_{\sigma}(\mathbf{r},\mathbf{r}^{\prime})=N\int_{(\mathbb{R}^{3}\times\\{\uparrow,\downarrow\\})^{N-1}}\Phi[\rho](\mathbf{r}^{\prime},\sigma,\mathbf{x}_{2},...,\mathbf{x}_{N})^{*}\;\Phi[\rho](\mathbf{r},\sigma,\mathbf{x}_{2},...,\mathbf{x}_{N})\text{d}\mathbf{x}_{2}...\text{d}\mathbf{x}_{N},$
(1.24)

which shows that $E_{\text{x}}[\rho]\leq 0$. Moreover, from the variational
definition of $F[\rho]$, we see that $E_{\text{c}}[\rho]\leq 0$.

####  1.3.2 The Kohn-Sham equations

The single-determinant wave function $\Phi$ in Eq. (1.17) is constructed from
a set of $N$ orthonormal occupied spin-orbitals $\\{\phi_{i}\\}_{i=1,...,N}$.
To enforce ${S}_{z}$ spin symmetry, each spin-orbital is factorized as
$\phi_{i}(\mathbf{x})=\varphi_{i}(\mathbf{r})\chi_{\sigma_{i}}(\sigma)$ where
$\varphi_{i}\in H^{1}(\mathbb{R}^{3})$ is a spatial orbital and
$\chi_{\sigma_{i}}$ is a spin function from $\\{\uparrow,\downarrow\\}$ to
$\\{0,1\\}$ such that
$\forall\sigma_{i},\sigma\in\\{\uparrow,\downarrow\\},\,\chi_{\sigma_{i}}(\sigma)=\delta_{\sigma_{i},\sigma}$
($\sigma_{i}$ is the spin of the spin-orbital $i$). Alternatively, when this
is convenient, we will sometimes reindex the spatial orbitals,
$\\{\varphi_{i}\\}\longrightarrow\\{\varphi_{i\sigma}\\}$, including
explicitly the spin $\sigma$ 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

$E[\\{\varphi_{i}\\}]=\frac{1}{2}\sum_{i=1}^{N}\int_{\mathbb{R}^{3}}|\nabla\varphi_{i}(\mathbf{r})|^{2}\text{d}\mathbf{r}+\int_{\mathbb{R}^{3}}v_{\text{ne}}(\mathbf{r})\rho(\mathbf{r})\text{d}\mathbf{r}+E_{\text{Hxc}}[\rho],$
(1.25)

where the density is expressed in terms of the orbitals as

$\rho(\mathbf{r})=\sum_{i=1}^{N}\left|\varphi_{i}(\mathbf{r})\right|^{2}.$
(1.26)

The minimization over $\Phi$ can then be recast into a minimization of
$E[\\{\varphi_{i}\\}]$ with respect to the spatial orbitals
$\\{\varphi_{i}\\}$ with the constraint of keeping the orbitals
orthonormalized. The stationary condition with respect to variations of
$\varphi_{i}(\mathbf{r})$ leads to the KS equations [19]

$\left(-\frac{1}{2}\nabla^{2}+v_{\text{ne}}(\mathbf{r})+v_{\text{Hxc}}(\mathbf{r})\right)\varphi_{i}(\mathbf{r})=\varepsilon_{i}\varphi_{i}(\mathbf{r}),$
(1.27)

where $\varepsilon_{i}$ is the Lagrange multiplier associated to the
normalization condition of $\varphi_{i}$ and $v_{\text{Hxc}}(\mathbf{r})$ is
the Hartree-exchange-correlation potential defined as the functional
derivative of $E_{\text{Hxc}}[\rho]$ with respect to $\rho(\mathbf{r})$

$v_{\text{Hxc}}(\mathbf{r})=\frac{\delta
E_{\text{Hxc}}[\rho]}{\delta\rho(\mathbf{r})},$ (1.28)

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

$\displaystyle
h_{\text{s}}(\mathbf{r})=-\frac{1}{2}\nabla^{2}+v_{\text{s}}(\mathbf{r}),$
(1.29)

where

$\displaystyle
v_{\text{s}}(\mathbf{r})=v_{\text{ne}}(\mathbf{r})+v_{\text{Hxc}}(\mathbf{r}),$
(1.30)

is the KS potential, and $\varepsilon_{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
$\\{\varphi_{i}\\}_{i=1,...,N}$ through the density [Eq. (1.26)]. The operator
$h_{\text{s}}(\mathbf{r})$ defines the KS system which is a system of $N$ non-
interacting electrons in an effective external potential
$v_{\text{s}}(\mathbf{r})$ ensuring that the density $\rho(\mathbf{r})$ in Eq.
(1.26) is the same as the exact ground-state density $\rho_{0}(\mathbf{r})$ of
the physical system of $N$ interacting electrons. The exact ground-state
energy $E_{0}$ is then obtained by injecting the KS orbitals in Eq. (1.25).
The other (unoccupied) eigenfunctions in Eq. (1.27) define virtual KS orbitals
$\\{\varphi_{a}\\}_{a\geq N+1}$.

Note that to define the potential $v_{\text{Hxc}}(\mathbf{r})$ in Eq. (1.28) a
form of differentiability of the functional $E_{\text{Hxc}}[\rho]$, also
referred to as $v$-representability of the density, has been assumed.
Justifying this is in fact subtle and has been debated [22, 23, 24, 25, 26,
27, 28] (see also the chapter by Kvaal in the volume). Here, we will simply
assume that a form of differentiability of $E_{\text{Hxc}}[\rho]$ holds on at
least a restricted set of densities that allows one to define the potential
$v_{\text{Hxc}}(\mathbf{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 $v_{\text{s}}(\mathbf{r})$ tends to a constant as
$|\mathbf{r}|\to\infty$ and we choose this constant to be zero. Note also that
the assumption of the existence of the KS potential $v_{\text{s}}(\mathbf{r})$
in Eq. (1.29), which does not depend on spin coordinates, implies that each
spin-orbital must indeed has a definite ${S}_{z}$ spin value.

Following the decomposition of $E_{\text{Hxc}}[\rho]$ in Eq. (1.18), the
potential $v_{\text{Hxc}}(\mathbf{r})$ is decomposed as

$\displaystyle
v_{\text{Hxc}}(\mathbf{r})=v_{\text{H}}(\mathbf{r})+v_{\text{xc}}(\mathbf{r}),$
(1.31)

where $v_{\text{H}}(\mathbf{r})=\delta
E_{\text{H}}[\rho]/\delta\rho(\mathbf{r})=\int_{\mathbb{R}^{3}}\rho(\mathbf{r}^{\prime})/|\mathbf{r}-\mathbf{r}^{\prime}|\text{d}\mathbf{r}^{\prime}$
is the Hartree potential and $v_{\text{xc}}(\mathbf{r})=\delta
E_{\text{xc}}[\rho]/\delta\rho(\mathbf{r})$ is the exchange-correlation
potential. Likewise, following the decomposition of $E_{\text{xc}}[\rho]$ in
Eq. (1.20), and assuming that both $E_{\text{x}}[\rho]$ and
$E_{\text{c}}[\rho]$ are differentiable with respect to $\rho$, the potential
$v_{\text{xc}}(\mathbf{r})$ can be further decomposed as

$\displaystyle
v_{\text{xc}}(\mathbf{r})=v_{\text{x}}(\mathbf{r})+v_{\text{c}}(\mathbf{r}),$
(1.32)

where $v_{\text{x}}(\mathbf{r})=\delta
E_{\text{x}}[\rho]/\delta\rho(\mathbf{r})$ is the exchange potential and
$v_{\text{c}}(\mathbf{r})=\delta E_{\text{c}}[\rho]/\delta\rho(\mathbf{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 $v_{\text{x}}(\mathbf{r})$ instead of the nonlocal HF exchange one,
and an additional correlation potential. One important feature of the exchange
potential is its long-range asymptotic behavior (see, e.g., Ref. [29])

$\displaystyle
v_{\text{x}}(\mathbf{r})\underset{|\mathbf{r}|\to\infty}{\sim}-\frac{1}{|\mathbf{r}|},$
(1.33)

whereas the correlation potential decays faster.

####  1.3.3 Extension to spin density-functional theory

For dealing with an external magnetic field, DFT has been extended from the
total density to spin-resolved densities [30, 31]. 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 $\rho_{\sigma,\Psi}$ with
$\sigma\in\\{\uparrow,\downarrow\\}$ associated to a wave function $\Psi$ is
defined as

$\rho_{\sigma,\Psi}(\mathbf{r})=N\int_{(\mathbb{R}^{3}\times\\{\uparrow,\downarrow\\})^{N-1}}\left|\Psi(\mathbf{x},\mathbf{x}_{2},...,\mathbf{x}_{N})\right|^{2}\text{d}\mathbf{x}_{2}...\text{d}\mathbf{x}_{N},$
(1.34)

and integrates to the number of $\sigma$-spin electrons $N_{\sigma}$, i.e.
$\int_{\mathbb{R}^{3}}\rho_{\sigma,\Psi}(\mathbf{r})\text{d}\mathbf{r}=N_{\sigma}$.
For $\rho_{\uparrow}\in{\cal D}^{N_{\uparrow}}$ and $\rho_{\downarrow}\in{\cal
D}^{N_{\downarrow}}$, the universal density functional is now defined as [32]

$F[\rho_{\uparrow},\rho_{\downarrow}]=\min_{\Psi\in{\cal
W}^{N}_{\rho_{\uparrow},\rho_{\downarrow}}}\langle\Psi|\hat{T}+\hat{W}_{\text{ee}}|\Psi\rangle=\langle\Psi[\rho_{\uparrow},\rho_{\downarrow}]|\hat{T}+\hat{W}_{\text{ee}}|\Psi[\rho_{\uparrow},\rho_{\downarrow}]\rangle,$
(1.35)

where the search is over the set of normalized antisymmetric wave functions
$\Psi$ with $N=N_{\uparrow}+N_{\downarrow}$ electrons and yielding the fixed
spin densities $\rho_{\uparrow}$ and $\rho_{\downarrow}$

${\cal W}^{N}_{\rho_{\uparrow},\rho_{\downarrow}}=\\{\Psi\in{\cal
W}^{N},\;\rho_{\uparrow,\Psi}=\rho_{\uparrow},\;\rho_{\downarrow,\Psi}=\rho_{\downarrow}\\}.$
(1.36)

In Eq. (1.35), $\Psi[\rho_{\uparrow},\rho_{\downarrow}]$ designates a
minimizing wave function.

A spin-dependent KS scheme is obtained by decomposing
$F[\rho_{\uparrow},\rho_{\downarrow}]$ as

$F[\rho_{\uparrow},\rho_{\downarrow}]=T_{\text{s}}[\rho_{\uparrow},\rho_{\downarrow}]+E_{\text{H}}[\rho]+E_{\text{xc}}[\rho_{\uparrow},\rho_{\downarrow}],$
(1.37)

where $T_{\text{s}}[\rho_{\uparrow},\rho_{\downarrow}]$ is defined as

$T_{\text{s}}[\rho_{\uparrow},\rho_{\downarrow}]=\min_{\Phi\in{\cal
S}^{N}_{\rho_{\uparrow},\rho_{\downarrow}}}\langle\Phi|\hat{T}|\Phi\rangle=\langle\Phi[\rho_{\uparrow},\rho_{\downarrow}]|\hat{T}|\Phi[\rho_{\uparrow},\rho_{\downarrow}]\rangle,$
(1.38)

with a constrained search over the set of single-determinant wave functions
$\Phi$ yielding the fixed spin densities $\rho_{\uparrow}$ and
$\rho_{\downarrow}$

${\cal S}^{N}_{\rho_{\uparrow},\rho_{\downarrow}}=\\{\Phi\in{\cal
S}^{N},\;\rho_{\uparrow,\Phi}=\rho_{\uparrow},\;\rho_{\downarrow,\Phi}=\rho_{\downarrow}\\},$
(1.39)

$\Phi[\rho_{\uparrow},\rho_{\downarrow}]$ denotes a minimizing KS single-
determinant wave function, $E_{\text{H}}[\rho]$ is the Hartree energy which is
a functional of the total density $\rho=\rho_{\uparrow}+\rho_{\downarrow}$
only [Eq. (1.19)], and $E_{\text{xc}}[\rho_{\uparrow},\rho_{\downarrow}]$ is
the spin-dependent exchange-correlation energy functional. The ground-state
energy is then obtained as

$E_{0}=\inf_{\Phi\in{\cal
S}^{N}}\left\\{\langle\Phi|\hat{T}+\hat{V}_{\text{ne}}|\Phi\rangle+E_{\text{H}}[\rho_{\Phi}]+E_{\text{xc}}[\rho_{\uparrow,\Phi},\rho_{\downarrow,\Phi}]\right\\}.$
(1.40)

Writing the spatial orbitals of the spin-unrestricted determinant as
$\\{\varphi_{i\sigma}\\}_{i=1,...,N}$ (with the index explicitly including the
spin $\sigma$ now for clarity), we have now the spin-dependent KS equations

$\left(-\frac{1}{2}\nabla^{2}+v_{\text{ne}}(\mathbf{r})+v_{\text{H}}(\mathbf{r})+v_{\text{xc},\sigma}(\mathbf{r})\right)\varphi_{i\sigma}(\mathbf{r})=\varepsilon_{i\sigma}\varphi_{i\sigma}(\mathbf{r}),$
(1.41)

with the spin-dependent exchange-correlation potential

$\displaystyle v_{\text{xc},\sigma}(\mathbf{r})=\frac{\delta
E_{\text{xc}}[\rho_{\uparrow},\rho_{\downarrow}]}{\delta\rho_{\sigma}(\mathbf{r})},$
(1.42)

and the spin density

$\displaystyle\rho_{\sigma}(\mathbf{r})=\sum_{i=1}^{N_{\sigma}}\left|\varphi_{i\sigma}(\mathbf{r})\right|^{2}.$
(1.43)

As before, if there is a unique KS wave function
$\Phi[\rho_{\uparrow},\rho_{\downarrow}]$ (up to a phase factor), we can
decompose $E_{\text{xc}}[\rho_{\uparrow},\rho_{\downarrow}]$ into exchange and
correlation contributions

$E_{\text{xc}}[\rho_{\uparrow},\rho_{\downarrow}]=E_{\text{x}}[\rho_{\uparrow},\rho_{\downarrow}]+E_{\text{c}}[\rho_{\uparrow},\rho_{\downarrow}],$
(1.44)

with
$E_{\text{x}}[\rho_{\uparrow},\rho_{\downarrow}]=\langle\Phi[\rho_{\uparrow},\rho_{\downarrow}]|\hat{W}_{\text{ee}}|\Phi[\rho_{\uparrow},\rho_{\downarrow}]\rangle-
E_{\text{H}}[\rho]$. It turns out that the spin-dependent exchange functional
$E_{\text{x}}[\rho_{\uparrow},\rho_{\downarrow}]$ can be exactly expressed in
terms of the spin-independent exchange functional $E_{\text{x}}[\rho]$ [33]

$E_{\text{x}}[\rho_{\uparrow},\rho_{\downarrow}]=\frac{1}{2}\left(E_{\text{x}}[2\rho_{\uparrow}]+E_{\text{x}}[2\rho_{\downarrow}]\right),$
(1.45)

which is known as the spin-scaling relation and stems directly from the fact
the $\uparrow$\- and $\downarrow$-spin electrons are uncoupled in the exchange
energy [see Eq. (1.23)]. Therefore, any approximation for $E_{\text{x}}[\rho]$
can be easily extended to an approximation for
$E_{\text{x}}[\rho_{\uparrow},\rho_{\downarrow}]$. Unfortunately, there is no
such relation for the spin-dependent correlation functional
$E_{\text{c}}[\rho_{\uparrow},\rho_{\downarrow}]$.

Obviously, in the spin-unpolarized case, i.e.
$\rho_{\uparrow}=\rho_{\downarrow}=\rho/2$, this spin-dependent formalism
reduces to the spin-independent one.

###  1.4 The generalized Kohn-Sham scheme

An important extension of the KS scheme is the so-called generalized Kohn-Sham
(GKS) scheme [34], which recognizes that the universal density functional
$F[\rho]$ of Eq. (1.8) can be decomposed in other ways than in the KS
decomposition of Eq. (1.13). In particular, we can decompose $F[\rho]$ as

$\displaystyle F[\rho]=\min_{\Phi\in{\cal
S}^{N}_{\rho}}\left\\{\langle\Phi|\hat{T}|\Phi\rangle+E_{\text{H}}[\rho_{\Phi}]+S[\Phi]\right\\}+\bar{S}[\rho],$
(1.46)

where $S[\Phi]$ is any functional of a single-determinant wave function
$\Phi\in{\cal S}^{N}$ leading to a minimum in Eq. (1.46), and $\bar{S}[\rho]$
is the corresponding complementary density functional that makes Eq. (1.46)
exact. Defining the $S$-dependent GKS exchange-correlation functional as

$\displaystyle E_{\text{xc}}^{S}[\Phi]=S[\Phi]+\bar{S}[\rho_{\Phi}],$ (1.47)

we can express the exact ground-state energy as

$\displaystyle E_{0}=\inf_{\Phi\in{\cal
S}^{N}}\left\\{\langle\Phi|\hat{T}+\hat{V}_{\text{ne}}|\Phi\rangle+E_{\text{H}}[\rho_{\Phi}]+E_{\text{xc}}^{S}[\Phi]\right\\},$
(1.48)

and if a minimum exists then any minimizing single-determinant wave function
in Eq. (1.48) gives a ground-state density $\rho_{0}(\mathbf{r})$. Similarly
to the KS equations [Eq. (1.27)], Eq. (1.48) leads to the one-electron GKS
equations

$\left(-\frac{1}{2}\nabla^{2}+v_{\text{ne}}(\mathbf{r})+v_{\text{H}}(\mathbf{r})+v_{\bar{S}}(\mathbf{r})\right)\varphi_{i\sigma}(\mathbf{r})+\frac{\delta
S[\Phi]}{\delta\varphi_{i\sigma}^{*}(\mathbf{r})}=\varepsilon_{i\sigma}\varphi_{i\sigma}(\mathbf{r}),$
(1.49)

where $v_{\bar{S}}(\mathbf{r})=\delta\bar{S}[\rho]/\delta\rho(\mathbf{r})$ is
a local potential and $\delta S[\Phi]/\delta\varphi_{i\sigma}^{*}(\mathbf{r})$
generates a one-electron (possibly nonlocal) operator.

In the special case $S[\Phi]=0$, we recover the KS exchange-correlation
density functional

$\displaystyle E_{\text{xc}}^{S=0}[\Phi]=E_{\text{xc}}[\rho_{\Phi}].$ (1.50)

Due to the freedom in the choice of $S[\Phi]$, there is an infinity of GKS
exchange-correlation functionals $E_{\text{xc}}^{S}[\Phi]$ giving the exact
ground-state energy via Eq. (1.48). This freedom and the fact that $\Phi$
carries more information than $\rho_{\Phi}$ gives the possibility to design
more accurate approximations for the exchange-correlation energy.

Of course, by starting from the density functional
$F[\rho_{\uparrow},\rho_{\downarrow}]$ in Eq. (1.35), this GKS scheme can be
extended to the spin-dependent case, leading to GKS exchange-correlation
functionals of the form
$E_{\text{xc}}^{S}[\Phi]=S[\Phi]+\bar{S}[\rho_{\uparrow,\Phi},\rho_{\downarrow,\Phi}]$.

##  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
$\Psi[\rho]$ defined in Eq. (1.8)

$\rho_{2}(\mathbf{r}_{1},\mathbf{r}_{2})=N(N-1)\int_{\\{\uparrow,\downarrow\\}^{2}\times(\mathbb{R}^{3}\times\\{\uparrow,\downarrow\\})^{N-2}}\left|\Psi[\rho](\mathbf{x}_{1},\mathbf{x}_{2},...,\mathbf{x}_{N})\right|^{2}\text{d}\sigma_{1}\text{d}\sigma_{2}\text{d}\mathbf{x}_{3}...\text{d}\mathbf{x}_{N},$
(2.1)

which is a functional of the density, and is normalized to the number of
electron pairs,
$\int_{\mathbb{R}^{3}\times\mathbb{R}^{3}}\rho_{2}(\mathbf{r}_{1},\mathbf{r}_{2})\text{d}\mathbf{r}_{1}\text{d}\mathbf{r}_{2}=N(N-1)$.
It is proportional to the probability density of finding two electrons at
positions $(\mathbf{r}_{1},\mathbf{r}_{2})$ with all the other electrons being
anywhere. The pair density is useful to express the expectation value of the
electron-electron interaction operator

$\langle\Psi[\rho]|\hat{W}_{\text{ee}}|\Psi[\rho]\rangle=\frac{1}{2}\int_{\mathbb{R}^{3}\times\mathbb{R}^{3}}\frac{\rho_{2}(\mathbf{r}_{1},\mathbf{r}_{2})}{|\mathbf{r}_{1}-\mathbf{r}_{2}|}\text{d}\mathbf{r}_{1}\text{d}\mathbf{r}_{2}.$
(2.2)

Mirroring the decomposition of the Hartree-exchange-correlation energy
performed in the KS scheme [Eq. (1.18)], the pair density can be decomposed as

$\rho_{2}(\mathbf{r}_{1},\mathbf{r}_{2})=\rho(\mathbf{r}_{1})\rho(\mathbf{r}_{2})+\rho_{2,\text{xc}}(\mathbf{r}_{1},\mathbf{r}_{2}).$
(2.3)

The product of the densities $\rho(\mathbf{r}_{1})\rho(\mathbf{r}_{2})$
corresponds to the case of independent electrons [up to a change of
normalization, i.e.
$\int_{\mathbb{R}^{3}\times\mathbb{R}^{3}}\rho(\mathbf{r}_{1})\rho(\mathbf{r}_{2})\text{d}\mathbf{r}_{1}\text{d}\mathbf{r}_{2}=N^{2}$
instead of $N(N-1)$] and the exchange-correlation pair density
$\rho_{2,\text{xc}}(\mathbf{r}_{1},\mathbf{r}_{2})$ represents the
modification of the pair density due to exchange and correlation effects
between the electrons. It can be further written as

$\rho_{2,\text{xc}}(\mathbf{r}_{1},\mathbf{r}_{2})=\rho(\mathbf{r}_{1})h_{\text{xc}}(\mathbf{r}_{1},\mathbf{r}_{2}),$
(2.4)

where $h_{\text{xc}}(\mathbf{r}_{1},\mathbf{r}_{2})$ is the exchange-
correlation hole. Introducing the conditional density
$\rho^{\text{cond}}(\mathbf{r}_{1},\mathbf{r}_{2})=\rho_{2}(\mathbf{r}_{1},\mathbf{r}_{2})/\rho(\mathbf{r}_{1})$
of the remaining $N-1$ electrons at $\mathbf{r}_{2}$ given that one electron
has been found at $\mathbf{r}_{1}$, the exchange-correlation hole can be
interpreted as the modification of
$\rho^{\text{cond}}(\mathbf{r}_{1},\mathbf{r}_{2})$ due to exchange and
correlation effects

$\rho^{\text{cond}}(\mathbf{r}_{1},\mathbf{r}_{2})=\rho(\mathbf{r}_{2})+h_{\text{xc}}(\mathbf{r}_{1},\mathbf{r}_{2}).$
(2.5)

The positivity of $\rho_{2}(\mathbf{r}_{1},\mathbf{r}_{2})$ implies that

$h_{\text{xc}}(\mathbf{r}_{1},\mathbf{r}_{2})\geq-\rho(\mathbf{r}_{2}).$ (2.6)

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

$\forall\mathbf{r}_{1}\in\mathbb{R}^{3},\;\int_{\mathbb{R}^{3}}h_{\text{xc}}(\mathbf{r}_{1},\mathbf{r}_{2})\text{d}\mathbf{r}_{2}=-1.$
(2.7)

We can separate the exchange and correlation contributions in the exchange-
correlation hole. For this, consider the pair density
$\rho_{2,\text{KS}}(\mathbf{r}_{1},\mathbf{r}_{2})$ associated with the KS
single-determinant wave function $\Phi[\rho]$ defined in Eq. (1.14). It can be
decomposed as

$\rho_{2,\text{KS}}(\mathbf{r}_{1},\mathbf{r}_{2})=\rho(\mathbf{r}_{1})\rho(\mathbf{r}_{2})+\rho_{2,\text{x}}(\mathbf{r}_{1},\mathbf{r}_{2}),$
(2.8)

where $\rho_{2,\text{x}}(\mathbf{r}_{1},\mathbf{r}_{2})$ is the exchange pair
density, which is further written as

$\rho_{2,\text{x}}(\mathbf{r}_{1},\mathbf{r}_{2})=\rho(\mathbf{r}_{1})h_{\text{x}}(\mathbf{r}_{1},\mathbf{r}_{2}),$
(2.9)

where $h_{\text{x}}(\mathbf{r}_{1},\mathbf{r}_{2})$ is the exchange hole. Just
like the exchange-correlation hole, the exchange hole satisfies the conditions

$h_{\text{x}}(\mathbf{r}_{1},\mathbf{r}_{2})\geq-\rho(\mathbf{r}_{2}),$ (2.10)

and

$\forall\mathbf{r}_{1}\in\mathbb{R}^{3},\;\int_{\mathbb{R}^{3}}h_{\text{x}}(\mathbf{r}_{1},\mathbf{r}_{2})\text{d}\mathbf{r}_{2}=-1.$
(2.11)

Moreover, since the exchange hole can be written as [compare with Eq. (1.23)]

$h_{\text{x}}(\mathbf{r}_{1},\mathbf{r}_{2})=-\frac{1}{\rho(\mathbf{r}_{1})}\sum_{\sigma\in\\{\uparrow,\downarrow\\}}|\gamma_{\sigma}(\mathbf{r}_{1},\mathbf{r}_{2})|^{2},$
(2.12)

where
$\gamma_{\sigma}(\mathbf{r}_{1},\mathbf{r}_{2})=\sum_{i=1}^{N_{\sigma}}\varphi_{i\sigma}^{*}(\mathbf{r}_{2})\varphi_{i\sigma}(\mathbf{r}_{1})$
is the spin-dependent one-particle KS density matrix, it thus appears that it
is always non-positive

$h_{\text{x}}(\mathbf{r}_{1},\mathbf{r}_{2})\leq 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

$E_{\text{x}}[\rho]=\frac{1}{2}\int_{\mathbb{R}^{3}\times\mathbb{R}^{3}}\frac{\rho(\mathbf{r}_{1})h_{\text{x}}(\mathbf{r}_{1},\mathbf{r}_{2})}{|\mathbf{r}_{1}-\mathbf{r}_{2}|}\text{d}\mathbf{r}_{1}\text{d}\mathbf{r}_{2},$
(2.14)

leading to the interpretation of $E_{\text{x}}$ as the electrostatic
interaction energy of an electron and its exchange hole. It is useful to write
the exchange energy functional as

$\displaystyle
E_{\text{x}}[\rho]=\int_{\mathbb{R}^{3}}\rho(\mathbf{r}_{1})\varepsilon_{\text{x}}[\rho](\mathbf{r}_{1})\text{d}\mathbf{r}_{1},$
(2.15)

where $\varepsilon_{\text{x}}[\rho](\mathbf{r}_{1})$ is the exchange energy
density per particle

$\displaystyle\varepsilon_{\text{x}}[\rho](\mathbf{r}_{1})=\frac{1}{2}\int_{\mathbb{R}^{3}}\frac{h_{\text{x}}(\mathbf{r}_{1},\mathbf{r}_{2})}{|\mathbf{r}_{1}-\mathbf{r}_{2}|}\text{d}\mathbf{r}_{2},$
(2.16)

which is itself a functional of the density. It is also convenient to define
the exchange energy density
$e_{\text{x}}[\rho](\mathbf{r})=\rho(\mathbf{r})\varepsilon_{\text{x}}[\rho](\mathbf{r})$.
For finite systems, we have the exact asymptotic behavior [35, 36]

$\displaystyle\varepsilon_{\text{x}}[\rho](\mathbf{r})\underset{|\mathbf{r}|\to+\infty}{\sim}-\frac{1}{2|\mathbf{r}|}.$
(2.17)

The correlation hole is defined as the difference

$h_{\text{c}}(\mathbf{r}_{1},\mathbf{r}_{2})=h_{\text{xc}}(\mathbf{r}_{1},\mathbf{r}_{2})-h_{\text{x}}(\mathbf{r}_{1},\mathbf{r}_{2}),$
(2.18)

and, from Eqs. (2.7) and (2.11), satisfies the sum rule

$\forall\mathbf{r}_{1}\in\mathbb{R}^{3},\;\int_{\mathbb{R}^{3}}h_{\text{c}}(\mathbf{r}_{1},\mathbf{r}_{2})\text{d}\mathbf{r}_{2}=0,$
(2.19)

which implies that the correlation hole has negative and positive
contributions333Therefore, the correlation hole is really a “hole” only in
some region of space, and a “bump” in other regions.. In contrast with the
exchange hole which is a smooth function of the interelectronic coordinate
$\mathbf{r}_{12}=\mathbf{r}_{2}-\mathbf{r}_{1}$, the correlation hole
satisfies the electron-electron cusp condition (i.e., it has a derivative
discontinuity in $\mathbf{r}_{12}$) [37]

$\forall\mathbf{r}_{1}\in\mathbb{R}^{3},\;h_{\text{c}}^{\prime}(\mathbf{r}_{1},\mathbf{r}_{1})=h_{\text{c}}(\mathbf{r}_{1},\mathbf{r}_{1}),$
(2.20)

where
$h_{\text{c}}^{\prime}(\mathbf{r}_{1},\mathbf{r}_{1})=(\partial\tilde{h}_{\text{c}}(\mathbf{r}_{1},r_{12})/\partial
r_{12})_{r_{12}=0}$ is the first-order derivative of the spherically averaged
correlation hole $\tilde{h}_{\text{c}}(\mathbf{r}_{1},r_{12})=(1/(4\pi
r_{12}^{2}))\int_{S(\mathbf{0},r_{12})}h_{\text{c}}(\mathbf{r}_{1},\mathbf{r}_{1}+\mathbf{r}_{12})\text{d}\mathbf{r}_{12}$
and $S(\mathbf{0},r_{12})$ designates the sphere centered at $\mathbf{0}$ and
of radius $r_{12}=|\mathbf{r}_{12}|$. The potential contribution to the
correlation energy can be written in terms of the correlation hole

$U_{\text{c}}[\rho]=\frac{1}{2}\int_{\mathbb{R}^{3}\times\mathbb{R}^{3}}\frac{\rho(\mathbf{r}_{1})h_{\text{c}}(\mathbf{r}_{1},\mathbf{r}_{2})}{|\mathbf{r}_{1}-\mathbf{r}_{2}|}\text{d}\mathbf{r}_{1}\text{d}\mathbf{r}_{2}.$
(2.21)

In order to express the total correlation energy
$E_{\text{c}}[\rho]=T_{\text{c}}[\rho]+U_{\text{c}}[\rho]$ in a form similar
to Eq. (2.21), we need to introduce the adiabatic-connection formalism.

###  2.2 The adiabatic connection

The idea of the adiabatic connection [38, 39, 40] (see, also, Ref. [41]) 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
$E_{\text{c}}[\rho]$ 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 $\lambda$. The
Hamiltonian along this adiabatic connection is

$\hat{H}^{\lambda}=\hat{T}+\lambda\hat{W}_{\text{ee}}+\hat{V}^{\lambda},$
(2.22)

where $\hat{V}^{\lambda}$ 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 $\lambda\in\mathbb{R}$. Of course, Eq. (2.22) relies
on a $v$-representability assumption, i.e. the external potential is assumed
to exist for all $\lambda$. The Hamiltonian (2.22) reduces to the KS non-
interacting Hamiltonian for $\lambda=0$ and to the physical Hamiltonian for
$\lambda=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 $\lambda$

$\displaystyle F^{\lambda}[\rho]$ $\displaystyle=$
$\displaystyle\min_{\Psi\in{\cal
W}^{N}_{\rho}}\langle\Psi|\hat{T}+\lambda\hat{W}_{\text{ee}}|\Psi\rangle=\langle\Psi^{\lambda}[\rho]|\hat{T}+\lambda\hat{W}_{\text{ee}}|\Psi^{\lambda}[\rho]\rangle,$
(2.23)

where $\Psi^{\lambda}[\rho]$ denotes a minimizing wave function. This
functional can be decomposed as

$\displaystyle
F^{\lambda}[\rho]=T_{\text{s}}[\rho]+E_{\text{Hxc}}^{\lambda}[\rho],$ (2.24)

where $E_{\text{Hxc}}^{\lambda}[\rho]$ is the Hartree-exchange-correlation
functional associated with the interaction $\lambda\hat{W}_{\text{ee}}$. It
can be decomposed as
$E_{\text{Hxc}}^{\lambda}[\rho]=E_{\text{H}}^{\lambda}[\rho]+E_{\text{x}}^{\lambda}[\rho]+E_{\text{c}}^{\lambda}[\rho]$,
where the Hartree and exchange contributions are simply linear in $\lambda$

$E_{\text{H}}^{\lambda}[\rho]=\frac{1}{2}\int_{\mathbb{R}^{3}\times\mathbb{R}^{3}}\rho(\mathbf{r}_{1})\rho(\mathbf{r}_{2})\frac{\lambda}{|\mathbf{r}_{1}-\mathbf{r}_{2}|}\text{d}\mathbf{r}_{1}\text{d}\mathbf{r}_{2}=\lambda
E_{\text{H}}[\rho],$ (2.25)

and

$E_{\text{x}}^{\lambda}[\rho]=\langle\Phi[\rho]|\lambda\hat{W}_{\text{ee}}|\Phi[\rho]\rangle-
E_{\text{H}}^{\lambda}[\rho]=\lambda E_{\text{x}}[\rho],$ (2.26)

and a correlation contribution, which is nonlinear in $\lambda$,

$E_{\text{c}}^{\lambda}[\rho]=\langle\Psi^{\lambda}[\rho]|\hat{T}+\lambda\hat{W}_{\text{ee}}|\Psi^{\lambda}[\rho]\rangle-\langle\Phi[\rho]|\hat{T}+\lambda\hat{W}_{\text{ee}}|\Phi[\rho]\rangle.$
(2.27)

We will assume that $F^{\lambda}[\rho]$ is of class $C^{1}$ as a function of
$\lambda$ for $\lambda\in[0,1]$ and that
$F^{\lambda=0}[\rho]=T_{\text{s}}[\rho]$, the latter condition being
guaranteed for nondegenerate KS systems [see footnote on the definition of
$T_{\text{s}}[\rho]$ just before Eq. (1.14)]. Taking the derivative of Eq.
(2.27) with respect to $\lambda$ and using the Hellmann-Feynman theorem for
the wave function $\Psi^{\lambda}[\rho]$444In this context, the Hellmann-
Feynman theorem states that in the derivative $\frac{\partial
F^{\lambda}[\rho]}{\partial\lambda}=\langle\frac{\partial\Psi^{\lambda}[\rho]}{\partial\lambda}|\hat{T}+\lambda\hat{W}_{\text{ee}}|\Psi^{\lambda}[\rho]\rangle+\langle\Psi^{\lambda}[\rho]|\hat{W}_{\text{ee}}|\Psi^{\lambda}[\rho]\rangle+\langle\Psi^{\lambda}[\rho]|\hat{T}+\lambda\hat{W}_{\text{ee}}|\frac{\partial\Psi^{\lambda}[\rho]}{\partial\lambda}\rangle$
the first and third terms involving the derivative of $\Psi^{\lambda}[\rho]$
vanish. This is due to the fact that $\Psi^{\lambda}[\rho]$ is obtained via
the minimization of Eq. (2.23) and thus any variation of
$\Psi^{\lambda}[\rho]$ which keeps the density constant (which is the case for
a variation with respect to $\lambda$) gives a vanishing variation of
$F^{\lambda}[\rho]$., we obtain

$\displaystyle\frac{\partial E_{\text{c}}^{\lambda}[\rho]}{\partial\lambda}$
$\displaystyle=$
$\displaystyle\langle\Psi^{\lambda}[\rho]|\hat{W}_{\text{ee}}|\Psi^{\lambda}[\rho]\rangle-\langle\Phi[\rho]|\hat{W}_{\text{ee}}|\Phi[\rho]\rangle,$
(2.28)

Integrating over $\lambda$ from $0$ to $1$, and using
$E_{\text{c}}^{\lambda=1}[\rho]=E_{\text{c}}[\rho]$ and
$E_{\text{c}}^{\lambda=0}[\rho]=0$, we arrive at the adiabatic-connection
formula for the correlation energy functional of the physical system

$E_{\text{c}}[\rho]=\int_{0}^{1}\text{d}\lambda\;\langle\Psi^{\lambda}[\rho]|\hat{W}_{\text{ee}}|\Psi^{\lambda}[\rho]\rangle-\langle\Phi[\rho]|\hat{W}_{\text{ee}}|\Phi[\rho]\rangle.$
(2.29)

By introducing the correlation hole
$h_{\text{c}}^{\lambda}(\mathbf{r}_{1},\mathbf{r}_{2})$ associated to the wave
function $\Psi^{\lambda}[\rho]$, the adiabatic-connection formula for the
correlation energy can also be written as

$E_{\text{c}}[\rho]=\frac{1}{2}\int_{0}^{1}\text{d}\lambda\int_{\mathbb{R}^{3}\times\mathbb{R}^{3}}\frac{\rho(\mathbf{r}_{1})h_{\text{c}}^{\lambda}(\mathbf{r}_{1},\mathbf{r}_{2})}{|\mathbf{r}_{1}-\mathbf{r}_{2}|}\text{d}\mathbf{r}_{1}\text{d}\mathbf{r}_{2},$
(2.30)

or, noting that $h_{\text{c}}^{\lambda}(\mathbf{r}_{1},\mathbf{r}_{2})$ is the
only quantity that depends on $\lambda$ in Eq. (2.30), in a more compact way,

$E_{\text{c}}[\rho]=\frac{1}{2}\int_{\mathbb{R}^{3}\times\mathbb{R}^{3}}\frac{\rho(\mathbf{r}_{1})\bar{h}_{\text{c}}(\mathbf{r}_{1},\mathbf{r}_{2})}{|\mathbf{r}_{1}-\mathbf{r}_{2}|}\text{d}\mathbf{r}_{1}\text{d}\mathbf{r}_{2},$
(2.31)

where
$\bar{h}_{\text{c}}(\mathbf{r}_{1},\mathbf{r}_{2})=\int_{0}^{1}\text{d}\lambda\;h_{\text{c}}^{\lambda}(\mathbf{r}_{1},\mathbf{r}_{2})$
is the coupling-constant-integrated correlation hole. It leads to the
interpretation of $E_{\text{c}}$ 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

$E_{\text{c}}[\rho]=\int_{\mathbb{R}^{3}}\rho(\mathbf{r}_{1})\varepsilon_{\text{c}}[\rho](\mathbf{r}_{1})\text{d}\mathbf{r}_{1},$
(2.32)

where $\varepsilon_{\text{c}}[\rho](\mathbf{r}_{1})$ is the correlation energy
density per particle

$\varepsilon_{\text{c}}[\rho](\mathbf{r}_{1})=\frac{1}{2}\int_{\mathbb{R}^{3}}\frac{\bar{h}_{\text{c}}(\mathbf{r}_{1},\mathbf{r}_{2})}{|\mathbf{r}_{1}-\mathbf{r}_{2}|}\text{d}\mathbf{r}_{2},$
(2.33)

which is a functional of the density. We can also define the correlation
energy density
$e_{\text{c}}[\rho](\mathbf{r})=\rho(\mathbf{r})\varepsilon_{\text{c}}[\rho](\mathbf{r})$.

Finally, note that the sum-rule and cusp conditions of Eqs. (2.19) and (2.20)
apply to the $\lambda$-dependent correlation hole in the form

$\forall\mathbf{r}_{1}\in\mathbb{R}^{3},\;\int_{\mathbb{R}^{3}}h_{\text{c}}^{\lambda}(\mathbf{r}_{1},\mathbf{r}_{2})\text{d}\mathbf{r}_{2}=0,$
(2.34)

and

$\forall\mathbf{r}_{1}\in\mathbb{R}^{3},\;h_{\text{c}}^{\lambda\;\prime}(\mathbf{r}_{1},\mathbf{r}_{1})=\lambda\;h_{\text{c}}^{\lambda}(\mathbf{r}_{1},\mathbf{r}_{1}).$
(2.35)

###  2.3 One-orbital and one-electron spatial regions

For systems composed of only one spin-$\uparrow$ (or, symmetrically, one
spin-$\downarrow$) electron (e.g., the hydrogen atom) with ground-state
density $\rho_{\text{1e}}(\mathbf{r})=|\varphi_{1\uparrow}(\mathbf{r})|^{2}$
where $\varphi_{1\uparrow}(\mathbf{r})$ is the unique occupied KS orbital, the
exchange hole in Eq. (2.12) simplifies to
$h_{\text{x}}(\mathbf{r}_{1},\mathbf{r}_{2})=-\rho(\mathbf{r}_{2})$, and
consequently the exchange energy cancels out the Hartree energy

$E_{\text{x}}[\rho_{\text{1e}}]=-E_{\text{H}}[\rho_{\text{1e}}].$ (2.36)

Besides, the correlation energy vanishes

$E_{\text{c}}[\rho_{\text{1e}}]=0.$ (2.37)

This is of course also true for the spin-dependent version of the functionals
introduced in Section 1.3.3, i.e.

$E_{\text{x}}[\rho_{\text{1e}},0]=-E_{\text{H}}[\rho_{\text{1e}}].$ (2.38)

and

$E_{\text{c}}[\rho_{\text{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
$\varphi_{1}(\mathbf{r})=\varphi_{1\uparrow}(\mathbf{r})=\varphi_{1\downarrow}(\mathbf{r})$
with ground-state density
$\rho_{\text{2e}}^{\uparrow\downarrow}(\mathbf{r})=2|\varphi_{1}(\mathbf{r})|^{2}$,
the exchange hole simplifies to
$h_{\text{x}}(\mathbf{r}_{1},\mathbf{r}_{2})=-\rho(\mathbf{r}_{2})/2$, and
consequently the exchange energy is equal to half the opposite of the Hartree
energy

$E_{\text{x}}[\rho_{\text{2e}}^{\uparrow\downarrow}]=-\frac{1}{2}E_{\text{H}}[\rho_{\text{2e}}^{\uparrow\downarrow}].$
(2.40)

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\geq 1$ electrons containing a spatial region
$\Omega_{\text{1o}}^{\uparrow}$ in which, among the occupied KS orbitals, only
one spin-$\uparrow$ (or, symmetrically, one spin-$\downarrow$) orbital is not
zero (or, more generally, takes non-negligible values), we have again in this
region

$\forall\mathbf{r}_{1},\mathbf{r}_{2}\in\Omega_{\text{1o}}^{\uparrow},\;h_{\text{x}}(\mathbf{r}_{1},\mathbf{r}_{2})=-\rho(\mathbf{r}_{2}),$
(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

$\displaystyle\forall\mathbf{r}_{1}\in\Omega_{\text{1o}}^{\uparrow},\;\varepsilon_{\text{x}}^{\Omega_{\text{1o}}^{\uparrow}}(\mathbf{r}_{1})=-\varepsilon_{\text{H}}^{\Omega_{\text{1o}}^{\uparrow}}(\mathbf{r}_{1}),$
(2.42)

where
$\varepsilon_{\text{H}}^{\Omega}(\mathbf{r}_{1})=(1/2)\int_{\Omega}\rho(\mathbf{r}_{2})/|\mathbf{r}_{1}-\mathbf{r}_{2}|\text{d}\mathbf{r}_{2}$
and
$\varepsilon_{\text{x}}^{\Omega}(\mathbf{r}_{1})=(1/2)\int_{\Omega}h_{\text{x}}(\mathbf{r}_{1},\mathbf{r}_{2})/|\mathbf{r}_{1}-\mathbf{r}_{2}|\text{d}\mathbf{r}_{2}$.
Similarly, for systems with $N\geq 2$ electrons containing a spatial region
$\Omega_{\text{1o}}^{\uparrow\downarrow}$ in which, among the occupied KS
orbitals, only one doubly occupied orbital is not zero, we have in this region

$\forall\mathbf{r}_{1},\mathbf{r}_{2}\in\Omega_{\text{1o}}^{\uparrow\downarrow},\;h_{\text{x}}(\mathbf{r}_{1},\mathbf{r}_{2})=-\frac{1}{2}\rho(\mathbf{r}_{2}),$
(2.43)

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

$\displaystyle\forall\mathbf{r}_{1}\in\Omega_{\text{1o}}^{\uparrow\downarrow},\;\varepsilon_{\text{x}}^{\Omega_{\text{1o}}^{\uparrow\downarrow}}(\mathbf{r}_{1})=-\frac{1}{2}\varepsilon_{\text{H}}^{\Omega_{\text{1o}}^{\uparrow\downarrow}}(\mathbf{r}_{1}).$
(2.44)

Thus, we see that, in $\Omega_{\text{1o}}^{\uparrow}$ or
$\Omega_{\text{1o}}^{\uparrow\downarrow}$ regions, the Hartree functional
introduces a spurious self-interaction contribution which must be eliminated
by the exchange functional. Even though the concepts of
$\Omega_{\text{1o}}^{\uparrow}$ and $\Omega_{\text{1o}}^{\uparrow\downarrow}$
regions are formal, in practice they can be approximately realized in chemical
systems. For example, the unpair electron in a radical approximately
corresponds to a $\Omega_{\text{1o}}^{\uparrow}$, and an electron pair in a
single covalent bond, in a lone pair, or in a core orbital approximately
corresponds to a $\Omega_{\text{1o}}^{\uparrow\downarrow}$ region.

We can also consider one-electron regions $\Omega_{\text{1e}}$ that we define
as555For generally, in the definition of Eq. (2.45) we exclude the point
$\lambda=0$ in order to allow for the possibility of a discontinuity in
$\lambda$ there due to a degeneracy.

$\displaystyle\forall\mathbf{r}_{1},\mathbf{r}_{2}\in\Omega_{\text{1e}},\;\forall\lambda\in(0,1],\;\rho_{2}^{\lambda}(\mathbf{r}_{1},\mathbf{r}_{2})=0,$
(2.45)

where $\rho_{2}^{\lambda}(\mathbf{r}_{1},\mathbf{r}_{2})$ is the pair density
associated to the wave function $\Psi^{\lambda}[\rho]$ along the adiabatic
connection. This implies

$\forall\mathbf{r}_{1},\mathbf{r}_{2}\in\Omega_{\text{1e}},\;\bar{h}_{\text{xc}}(\mathbf{r}_{1},\mathbf{r}_{2})=-\rho(\mathbf{r}_{2}),$
(2.46)

where
$\bar{h}_{\text{xc}}(\mathbf{r}_{1},\mathbf{r}_{2})=h_{\text{x}}(\mathbf{r}_{1},\mathbf{r}_{2})+\bar{h}_{\text{c}}(\mathbf{r}_{1},\mathbf{r}_{2})$
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

$\displaystyle\forall\mathbf{r}_{1}\in\Omega_{\text{1e}},\;\varepsilon_{\text{xc}}^{\Omega_{\text{1e}}}(\mathbf{r}_{1})=-\varepsilon_{\text{H}}^{\Omega_{\text{1e}}}(\mathbf{r}_{1}),$
(2.47)

where
$\varepsilon_{\text{xc}}^{\Omega}(\mathbf{r}_{1})=(1/2)\int_{\Omega}\bar{h}_{\text{xc}}(\mathbf{r}_{1},\mathbf{r}_{2})/|\mathbf{r}_{1}-\mathbf{r}_{2}|\text{d}\mathbf{r}_{2}$.
For regions that are simultaneously one-electron and one-orbital regions, this
simply implies that the contribution to the correlation energy must vanish

$\displaystyle\forall\mathbf{r}_{1}\in\Omega_{\text{1e}}\cap\Omega_{\text{1o}}^{\uparrow},\;\varepsilon_{\text{c}}^{\Omega_{\text{1e}}\cap\Omega_{\text{1o}}^{\uparrow}}(\mathbf{r}_{1})=0,$
(2.48)

where
$\varepsilon_{\text{c}}^{\Omega}(\mathbf{r}_{1})=(1/2)\int_{\Omega}\bar{h}_{\text{c}}(\mathbf{r}_{1},\mathbf{r}_{2})/|\mathbf{r}_{1}-\mathbf{r}_{2}|\text{d}\mathbf{r}_{2}$,
and we say that the correlation functional must not introduce a self-
interaction error. However, the definition of $\Omega_{\text{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 compensated
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
$\Psi^{\lambda}[\rho]$ [introduced in Eq. (2.23)] while leaving untouched the
spin coordinates [42, 43, 44]

$\Psi_{\gamma}^{\lambda}[\rho](\mathbf{r}_{1},\sigma_{1},...,\mathbf{r}_{N},\sigma_{N})=\gamma^{3N/2}\Psi^{\lambda}[\rho](\gamma\mathbf{r}_{1},\sigma_{1},,...,\gamma\mathbf{r}_{N},\sigma_{N}),$
(2.49)

where $\gamma\in(0,+\infty)$ is a scaling factor. The scaled wave function
$\Psi_{\gamma}^{\lambda}[\rho]$ yields the scaled density

$\rho_{\gamma}(\mathbf{r})=\gamma^{3}\rho(\gamma\mathbf{r}),$ (2.50)

with
$\int_{\mathbb{R}^{3}}\\!\rho_{\gamma}(\mathbf{r})\text{d}\mathbf{r}=\int_{\mathbb{R}^{3}}\\!\rho(\mathbf{r})\text{d}\mathbf{r}=N$
and minimizes
$\langle\Psi|\hat{T}+\lambda\gamma\hat{W}_{\text{ee}}|\Psi\rangle$ since

$\langle\Psi^{\lambda}_{\gamma}[\rho]|\hat{T}+\lambda\gamma\hat{W}_{\text{ee}}|\Psi^{\lambda}_{\gamma}[\rho]\rangle=\gamma^{2}\langle\Psi^{\lambda}[\rho]|\hat{T}+\lambda\hat{W}_{\text{ee}}|\Psi^{\lambda}[\rho]\rangle.$
(2.51)

We thus conclude that the scaled wave function at the density $\rho$ and
coupling constant $\lambda$ corresponds to the wave function at the scaled
density $\rho_{\gamma}$ and coupling constant $\lambda\gamma$

$\Psi_{\gamma}^{\lambda}[\rho]=\Psi^{\lambda\gamma}[\rho_{\gamma}],$ (2.52)

or, equivalently,

$\Psi_{\gamma}^{\lambda/\gamma}[\rho]=\Psi^{\lambda}[\rho_{\gamma}],$ (2.53)

and that the universal density functional satisfies the scaling relation

$F^{\lambda\gamma}[\rho_{\gamma}]=\gamma^{2}F^{\lambda}[\rho],$ (2.54)

or, equivalently,

$F^{\lambda}[\rho_{\gamma}]=\gamma^{2}F^{\lambda/\gamma}[\rho].$ (2.55)

At $\lambda=0$, we find the scaling relation of the KS wave function
$\Phi[\rho]$ introduced in Section 1.3.1

$\Phi[\rho_{\gamma}]=\Phi_{\gamma}[\rho].$ (2.56)

This directly leads to the scaling relation for the non-interacting kinetic
density functional [see Eq. (1.14)]

$T_{\text{s}}[\rho_{\gamma}]=\gamma^{2}T_{\text{s}}[\rho],$ (2.57)

for the Hartree density functional [see Eq. (1.19)]

$E_{\text{H}}[\rho_{\gamma}]=\gamma E_{\text{H}}[\rho],$ (2.58)

and for the exchange density functional [see Eq. (1.21)]

$E_{\text{x}}[\rho_{\gamma}]=\gamma E_{\text{x}}[\rho].$ (2.59)

However, the correlation density functional $E_{\text{c}}[\rho]$ has the more
complicated scaling (as $F[\rho]$)

$E_{\text{c}}^{\lambda}[\rho_{\gamma}]=\gamma^{2}E_{\text{c}}^{\lambda/\gamma}[\rho],$
(2.60)

and, in particular for $\lambda=1$,

$E_{\text{c}}[\rho_{\gamma}]=\gamma^{2}E_{\text{c}}^{1/\gamma}[\rho].$ (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
($\gamma\to\infty$), it can be shown from Eq. (2.61) that, for nondegenerate
KS systems, the correlation functional $E_{\text{c}}[\rho]$ goes to a constant

$\lim_{\gamma\to\infty}E_{\text{c}}[\rho_{\gamma}]=E_{\text{c}}^{\text{GL2}}[\rho],$
(2.62)

where $E_{\text{c}}^{\text{GL2}}[\rho]$ is the second-order Görling-Levy (GL2)
correlation energy [45, 46] (see Section 7.2). This is also called the weak-
correlation limit since in this limit the correlation energy is negligible
with respect to exchange energy which is itself negligible with respect to the
non-interacting kinetic energy:
$|E_{\text{c}}[\rho_{\gamma}]|=O(\gamma^{0})\ll|E_{\text{x}}[\rho_{\gamma}]|=O(\gamma)\ll
T_{\text{s}}[\rho_{\gamma}]=O(\gamma^{2})$. Equation (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 $E_{\text{c}}[\rho]=-0.0421$ hartree and
$\lim_{\gamma\to\infty}E_{\text{c}}[\rho_{\gamma}]=-0.0467$ hartree [47].

In the low-density limit ($\gamma\to 0$), it can be shown from Eq. (2.55) that
the Hartree-exchange-correlation $E_{\text{Hxc}}[\rho]$ goes to zero linearly
in $\gamma$

$E_{\text{Hxc}}[\rho_{\gamma}]\underset{\gamma\to
0}{\sim}\gamma\;W_{\text{ee}}^{\text{SCE}}[\rho],$ (2.63)

where $\displaystyle W_{\text{ee}}^{\text{SCE}}[\rho]=\inf_{\Psi\in{\cal
W}^{N}_{\rho}}\langle\Psi|\hat{W}_{\text{ee}}|\Psi\rangle$ is the strictly-
correlated-electron (SCE) functional [48, 49, 50, 51]. This is also called the
strong-interaction limit since in this limit the Hartree-exchange-correlation
energy dominates over the non-interacting kinetic energy:
$E_{\text{Hxc}}[\rho_{\gamma}]=O(\gamma)\gg
T_{\text{s}}[\rho_{\gamma}]=O(\gamma^{2})$. In this limit, the electrons
strictly localize relatively to each other. In particular, for the uniform-
electron gas, it corresponds to the Wigner crystallization. Thus, in this
limit, each electron is within a one-electron region $\Omega_{\text{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.

####  2.4.2 Non-uniform coordinate scaling

We can also consider non-uniform one-dimensional or two-dimensional coordinate
scalings of the density [52, 53]

$\displaystyle\rho_{\gamma}^{(1)}(x,y,z)=\gamma\rho(\gamma x,y,z),$ (2.64)

and

$\displaystyle\rho_{\gamma}^{(2)}(x,y,z)=\gamma^{2}\rho(\gamma x,\gamma y,z),$
(2.65)

which also preserve the number of the electrons. These non-uniform density
scaling 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
[43, 54]

$\displaystyle\lim_{\gamma\to\infty}E_{\text{x}}[\rho_{\gamma}^{(1)}]>-\infty,$
(2.66)
$\displaystyle\lim_{\gamma\to\infty}E_{\text{c}}[\rho_{\gamma}^{(1)}]=0.$
(2.67)

Also, in the non-uniform two-dimensional low-density limit, we have [43, 54]

$\displaystyle\lim_{\gamma\to
0}\frac{1}{\gamma}E_{\text{x}}[\rho_{\gamma}^{(2)}]>-\infty,$ (2.68)

and

$\displaystyle\lim_{\gamma\to
0}\frac{1}{\gamma}E_{\text{c}}[\rho_{\gamma}^{(2)}]=0.$ (2.69)

The conditions of Eqs. (2.66)-(2.69) are particularly useful because they also
correspond to the limit of rapidly varying densities [55].

###  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\to\infty$, at fixed
electron number $N$ (see Refs. [56, 57, 58, 59]). In this limit, the exact
ground-state atomic density $\rho_{N,Z}(\mathbf{r})$ becomes the density of
the isoelectronic hydrogenic (i.e., without electron-electron interaction)
atom $\rho_{N,Z}^{\text{H}}(\mathbf{r})$ which a simple scaling with $Z$

$\rho_{N,Z}(\mathbf{r})\underset{Z\to\infty}{\sim}\rho_{N,Z}^{\text{H}}(\mathbf{r})=Z^{3}\rho_{N,Z=1}^{\text{H}}(Z\mathbf{r}).$
(2.70)

One can thus apply Eqs. (2.59) and (2.62) with $\gamma=Z$, which reveals that
in an isoelectronic series the exchange functional scales linearly with $Z$

$E_{\text{x}}[\rho_{N,Z}]\underset{Z\to\infty}{\sim}E_{\text{x}}[\rho_{N,Z=1}^{\text{H}}]Z,$
(2.71)

and, for nondegenerate KS systems, the correlation functional saturates to a
constant

$\lim_{Z\to\infty}E_{\text{c}}[\rho_{N,Z}]=E_{\text{c}}^{\text{GL2}}[\rho_{N,Z=1}^{\text{H}}].$
(2.72)

Equations (2.71) and (2.72) are constraints for the exchange and correlation
functionals, particularly 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\to\infty$ (see, e.g., Ref. [60]). In this semiclassical limit, the
exact ground-state atomic density $\rho_{N,Z}(\mathbf{r})$ tends to the
Thomas-Fermi (TF) density of a neutral atom $\rho_{Z}^{\text{TF}}(\mathbf{r})$
which has a known scaling with $Z$ [61, 62]

$\rho_{Z,Z}(\mathbf{r})\underset{Z\to\infty}{\sim}\rho_{Z}^{\text{TF}}(\mathbf{r})=Z^{2}\rho_{Z=1}^{\text{TF}}(Z^{1/3}\mathbf{r}).$
(2.73)

In this limit, the exact exchange and correlation energies have the large-$Z$
asymptotic expansions [63, 64, 65]

$E_{\text{x}}[\rho_{Z,Z}]\underset{Z\to\infty}{\sim}-A_{\text{x}}Z^{5/3}+B_{\text{x}}Z+\cdots,$
(2.74) $E_{\text{c}}[\rho_{Z,Z}]\underset{Z\to\infty}{\sim}-A_{\text{c}}Z\ln
Z+B_{\text{c}}Z+\cdots,$ (2.75)

with the coefficients $A_{\text{x}}=0.220827$, $A_{\text{c}}=0.020727$,
$B_{\text{x}}\approx 0.224$, $B_{\text{c}}\approx 0.0372$.

###  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) [66],
which, when expressed in terms of the exchange or exchange-correlation
functional, takes the form [67]

$E_{\text{x}}[\rho]\geq E_{\text{xc}}[\rho]\geq-
C_{\text{LO}}\int_{\mathbb{R}^{3}}\rho(\mathbf{r})^{4/3}\text{d}\mathbf{r},$
(2.76)

where the optimal (i.e., smallest) constant $C_{\text{LO}}$ (independent of
the electron number $N$) was originally shown to be in the range $1.23\leq
C_{\text{LO}}\leq 1.68$ [66]. The range was later narrowed to $1.43\leq
C_{\text{LO}}\leq 1.6358$ [67, 68]. This bound is approached only in the low-
density limit where the correlation energy becomes comparable to the exchange
energy.

For two-electron densities, there is a specific tighter bound

$E_{\text{x}}[\rho_{\text{2e}}]\geq E_{\text{xc}}[\rho_{\text{2e}}]\geq-
C_{2}\int_{\mathbb{R}^{3}}\rho_{\text{2e}}(\mathbf{r})^{4/3}\text{d}\mathbf{r},$
(2.77)

with the best known constant $C_{2}=1.234$ [66]. For one-electron densities,
an even tighter bound is known for the exchange functional [69, 66]

$E_{\text{x}}[\rho_{\text{1e}}]\geq-
C_{1}\int_{\mathbb{R}^{3}}\rho_{\text{1e}}(\mathbf{r})^{4/3}\text{d}\mathbf{r},$
(2.78)

with the optimal constant $C_{1}=1.092$. For two-electron spin-unpolarized
densities, we have
$E_{\text{x}}[\rho_{\text{2e}}^{\uparrow\downarrow}]=2E_{\text{x}}[\rho_{\text{1e}}]$
with $\rho_{\text{1e}}=\rho_{\text{2e}}^{\uparrow\downarrow}/2$, and Eq.
(2.78) implies [70]

$E_{\text{x}}[\rho_{\text{2e}}^{\uparrow\downarrow}]\geq-\frac{C_{1}}{2^{1/3}}\int_{\mathbb{R}^{3}}\rho_{\text{2e}}^{\uparrow\downarrow}(\mathbf{r})^{4/3}\text{d}\mathbf{r},$
(2.79)

which is a much tighter bound than the bounds of Eqs. (2.76) and (2.77).

##  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 by

$E_{\text{xc}}^{\text{LDA}}[\rho]=\int_{\mathbb{R}^{3}}e_{\text{xc}}^{\text{UEG}}(\rho(\mathbf{r}))\text{d}\mathbf{r},$
(3.1)

where $e_{\text{xc}}^{\text{UEG}}(\rho)$ is the exchange-correlation energy
density of the infinite uniform electron gas (UEG) with the density $\rho$.
The UEG represents a family of systems of interacting electrons with an
arbitrary spatially constant density $\rho\in[0,+\infty)$ that acts as a
parameter. Thus, in the LDA, the exchange-correlation energy density of an
inhomogeneous system at a spatial point of density $\rho(\mathbf{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 [30]

$E_{\text{xc}}^{\text{LDA}}[\rho_{\uparrow},\rho_{\downarrow}]=\int_{\mathbb{R}^{3}}e_{\text{xc}}^{\text{UEG}}(\rho_{\uparrow}(\mathbf{r}),\rho_{\downarrow}(\mathbf{r}))\text{d}\mathbf{r},$
(3.2)

where $e_{\text{xc}}^{\text{UEG}}(\rho_{\uparrow},\rho_{\downarrow})$ is the
exchange-correlation energy density of the UEG with spin densities
$\rho_{\uparrow}$ and $\rho_{\downarrow}$. For spin-unpolarized systems, we
recover the spin-independent LDA as
$E_{\text{xc}}^{\text{LDA}}[\rho]=E_{\text{xc}}^{\text{LDA}}[\rho/2,\rho/2]$.

The function $e_{\text{xc}}^{\text{UEG}}$ is a sum of exchange and correlation
contributions,
$e_{\text{xc}}^{\text{UEG}}=e_{\text{x}}^{\text{UEG}}+e_{\text{c}}^{\text{UEG}}$,
and it is convenient to introduce exchange and correlation energies per
particle, $\varepsilon_{\text{x}}^{\text{UEG}}$ and
$\varepsilon_{\text{c}}^{\text{UEG}}$, such that
$e_{\text{x}}^{\text{UEG}}=\rho\;\varepsilon_{\text{x}}^{\text{UEG}}$ and
$e_{\text{c}}^{\text{UEG}}=\rho\;\varepsilon_{\text{c}}^{\text{UEG}}$. The
expression of the exchange energy per particle of the spin-unpolarized UEG is

$\varepsilon_{\text{x}}^{\text{UEG}}(\rho)=C_{\text{x}}\;\rho^{1/3},$ (3.3)

where $C_{\text{x}}=-(3/4)(3/\pi)^{1/3}$, and the spin-polarized version is
simply obtained from the spin-scaling relation [Eq. (1.45)], leading to

$\varepsilon_{\text{x}}^{\text{UEG}}(\rho_{\uparrow},\rho_{\downarrow})=\varepsilon_{\text{x}}^{\text{UEG}}(\rho)\phi_{4}(\zeta),$
(3.4)

where $\zeta=(\rho_{\uparrow}-\rho_{\downarrow})/\rho$ is the spin
polarization and $\phi_{4}(\zeta)$ is defined by the general spin-scaling
function

$\phi_{n}(\zeta)=\frac{(1+\zeta)^{n/3}+(1-\zeta)^{n/3}}{2}.$ (3.5)

The LDA exchange functional is associated with the names of Dirac [71] and
Slater [72]. For a rigorous mathematical derivation of Eq. (3.3), see Ref.
[73].

The correlation energy per particle
$\varepsilon_{\text{c}}^{\text{UEG}}(\rho_{\uparrow},\rho_{\downarrow})$ 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 $r_{\text{s}}=(3/(4\pi\rho))^{1/3}$, the first terms
of the high-density expansion ($r_{\text{s}}\to 0$) have the form

$\varepsilon_{\text{c}}^{\text{UEG}}(\rho_{\uparrow},\rho_{\downarrow})=A(\zeta)\ln
r_{\text{s}}+B(\zeta)+C(\zeta)r_{\text{s}}\ln r_{\text{s}}+O(r_{\text{s}}),$
(3.6)

with spin-unpolarized coefficients $A(0)=(1-\ln 2)/\pi^{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$, and the first terms of the low-density
expansion ($r_{\text{s}}\to+\infty$) have the form

$\varepsilon_{\text{c}}^{\text{UEG}}(\rho_{\uparrow},\rho_{\downarrow})=\frac{a}{r_{\text{s}}}+\frac{b}{r_{\text{s}}^{3/2}}+\frac{c}{r_{\text{s}}^{2}}+O\left(\frac{1}{r_{\text{s}}^{5/2}}\right),$
(3.7)

where the coefficients $a=-0.895930$, $b=1.325$, and $c=-0.365$ are assumed to
be independent of $\zeta$. The low-density limit of the UEG corresponds to the
Wigner crystallization. For a recent review on results on the UEG, see Ref.
[74].

The two most used parametrizations are the one of Vosko, Wilk, and Nusair
(VWN) [75] and the more recent one of Perdew and Wang (PW92) [76] which we
give here. In this parametrization, the UEG correlation energy per particle is
estimated using an approximate spin-interpolation formula

$\displaystyle\varepsilon_{\text{c}}^{\text{PW92}}(\rho_{\uparrow},\rho_{\downarrow})=\varepsilon_{\text{c}}(r_{\text{s}},0)+\alpha_{\text{c}}(r_{\text{s}})\frac{f(\zeta)}{f^{\prime\prime}(0)}(1-\zeta^{4})+[\varepsilon_{\text{c}}(r_{\text{s}},1)-\varepsilon_{\text{c}}(r_{\text{s}},0)]f(\zeta)\zeta^{4},$
(3.8)

where $\varepsilon_{\text{c}}(r_{\text{s}},\zeta)$ is the UEG correlation
energy per particle as a function of $r_{\text{s}}$ and $\zeta$,
$f(\zeta)=[(1+\zeta)^{4/3}+(1-\zeta)^{4/3}-2]/(2^{4/3}-2)$ is a spin-scaling
function borrowed from the exchange energy, and
$\alpha_{\text{c}}(r_{\text{s}})=(\partial^{2}\varepsilon_{\text{c}}(r_{\text{s}},\zeta)/\partial\zeta^{2})_{\zeta=0}$
is the spin stiffness. This spin-interpolation formula was first proposed in
the VWN parametrization based on a study of the $\zeta$ dependence of the UEG
correlation energy per particle at the random-phase approximation (RPA) level.
A unique parametrization function

$\displaystyle
G(r_{\text{s}},A,\alpha_{1},\beta_{1},\beta_{2},\beta_{3},\beta_{4})=-2(1+\alpha_{1}r_{\text{s}})A\ln\left[1+\frac{1}{2A\left(\beta_{1}r_{\text{s}}^{1/2}+\beta_{2}r_{\text{s}}+\beta_{3}r_{s}^{3/2}+\beta_{4}r_{\text{s}}^{2}\right)}\right],$
(3.9)

is then used for approximating $\varepsilon_{\text{c}}(r_{\text{s}},0)$,
$\varepsilon_{\text{c}}(r_{\text{s}},1)$, and
$-\alpha_{\text{c}}(r_{\text{s}})$

$\displaystyle\varepsilon_{\text{c}}(r_{\text{s}},0)=G(r_{\text{s}},A_{0},\alpha_{1,0},\beta_{1,0},\beta_{2,0},\beta_{3,0},\beta_{4,0}),$
(3.10a)
$\displaystyle\varepsilon_{\text{c}}(r_{\text{s}},1)=G(r_{\text{s}},A_{1},\alpha_{1,1},\beta_{1,1},\beta_{2,1},\beta_{3,1},\beta_{4,1}),$
(3.10b)
$\displaystyle-\alpha_{\text{c}}(r_{\text{s}})=G(r_{\text{s}},A_{2},\alpha_{1,2},\beta_{1,2},\beta_{2,2},\beta_{3,2},\beta_{4,2}).$
(3.10c)

The form of $G$ was chosen to reproduce the form of the high- and low-density
expansions. The parameters $A$’s, $\beta_{1}$’s, and $\beta_{2}$’s are fixed
by the first two terms of the high-density expansion, while the parameters
$\alpha_{1}$’s, $\beta_{3}$’s, and $\beta_{4}$’s are fitted to quantum Monte
Carlo (QMC) data [77] for $\varepsilon_{\text{c}}(r_{\text{s}},0)$ and
$\varepsilon_{\text{c}}(r_{\text{s}},1)$, and to an estimation of
$-\alpha_{\text{c}}(r_{\text{s}})$ extrapolated from RPA data. The parameters
are given in Table I of Ref. [76].

We know 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\to\infty$. Indeed, in this semiclassical Thomas-Fermi limit, the LDA
gives the exact coefficients $A_{\text{x}}$ and $A_{\text{c}}$ of the leading
terms in the asymptotic expansions of Eqs. (2.74) and (2.74) [78]. However,
the coefficients of the next terms are very incorrect:
$B_{\text{x}}^{\text{LDA}}\approx 0$ instead of $B_{\text{x}}\approx 0.224$
and $B_{\text{c}}^{\text{LDA}}\approx-0.00451$ instead of $B_{\text{c}}\approx
0.0372$ [64].

Due to the scaling of the UEG exchange energy per particle

$\varepsilon_{\text{x}}^{\text{UEG}}(\gamma^{3}\rho_{\uparrow},\gamma^{3}\rho_{\downarrow})=\gamma\varepsilon_{\text{x}}^{\text{UEG}}(\rho_{\uparrow},\rho_{\downarrow}),$
(3.11)

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)]

$\varepsilon_{\text{c}}^{\text{UEG}}(\gamma^{3}\rho_{\uparrow},\gamma^{3}\rho_{\downarrow})\underset{\gamma\to
0}{\sim}\gamma\;\frac{a}{r_{\text{s}}},$ (3.12)

the LDA correlation functional correctly scales linearly under uniform
coordinate scaling to the low-density limit [Eq. (2.63)]. However, from the
behavior of $\varepsilon_{\text{c}}^{\text{UEG}}$ in the high-density limit
[Eq. (3.6)],

$\varepsilon_{\text{c}}^{\text{UEG}}(\gamma^{3}\rho_{\uparrow},\gamma^{3}\rho_{\downarrow})\underset{\gamma\to\infty}{\sim}-A(\zeta)\;\ln\gamma,$
(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\to\infty$, 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
[79]. 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 corresponding exchange and correlation holes of
the UEG. Namely, considering the spin-independent version for simplicity, the
LDA exchange hole is

$h_{\text{x}}^{\text{LDA}}(\mathbf{r}_{1},\mathbf{r}_{2})=h_{\text{x}}^{\text{UEG}}(\rho(\mathbf{r}_{1}),r_{12}),$
(3.14)

with

$h_{\text{x}}^{\text{UEG}}(\rho,r_{12})=-\rho\;\frac{9}{2}\left(\frac{j_{1}(k_{\text{F}}r_{12})}{k_{\text{F}}r_{12}}\right)^{2},$
(3.15)

where $r_{12}=|\mathbf{r}_{2}-\mathbf{r}_{1}|$ is the interelectronic
distance, $k_{\text{F}}=(3\pi^{2}\rho)^{1/3}$ is the Fermi wave vector, and
$j_{1}$ is the spherical Bessel function of the first kind. Similarly, the LDA
correlation hole is

$\bar{h}_{\text{c}}^{\text{LDA}}(\mathbf{r}_{1},\mathbf{r}_{2})=\bar{h}_{\text{c}}^{\text{UEG}}(\rho(\mathbf{r}_{1}),r_{12})=\int_{0}^{1}\text{d}\lambda\;h_{\text{c}}^{\lambda,\text{UEG}}(\rho(\mathbf{r}_{1}),r_{12}).$
(3.16)

Since the UEG is a physical system, the LDA exchange hole correctly fulfils
the negativity and sum-rule condition [Eqs. (2.11) and (2.13)] and the LDA
correlation hole correctly fulfils 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
$h_{\text{x}}^{\text{LDA}}(\mathbf{r}_{1},\mathbf{r}_{2})$ only depends on
$\rho(\mathbf{r}_{1})$ and not on $\rho(\mathbf{r}_{2})$, the LDA exchange
functional does not entirely eliminate the self-interaction contribution of
the Hartree functional in one and two-electron systems [Eqs. (2.36) or (2.38),
and (2.40)], or more generally 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

$v_{\text{x}}^{\text{LDA}}(\mathbf{r})=\frac{\delta
E_{\text{x}}^{\text{LDA}}[\rho]}{\delta\rho(\mathbf{r})}=\frac{4}{3}C_{\text{x}}\;\rho(\mathbf{r})^{1/3},$
(3.17)

decays exponentially at infinity for finite molecular systems (since the
density $\rho(\mathbf{r})$ decays exponentially), i.e. much too fast in
comparison to the $-1/|\mathbf{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 $\delta v(\mathbf{r})$, and expand the exchange-
correlation energy in terms of the gradients of the density (see, e.g., Refs.
[80, 81, 82, 83]). Alternatively, one can perform a semiclassical expansion
(i.e., an expansion in powers of the reduced Planck constant $\hbar$) of the
exact $E_{\text{xc}}[\rho]$ 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) [83]

$\displaystyle E_{\text{x}}^{\text{GEA4}}[\rho]$ $\displaystyle=$
$\displaystyle
E_{\text{x}}^{\text{LDA}}[\rho]+C_{\text{x}}^{(2)}\int_{\mathbb{R}^{3}}\frac{|\nabla\rho(\mathbf{r})|^{2}}{\rho(\mathbf{r})^{4/3}}\text{d}\mathbf{r}$
(3.18)
$\displaystyle+C_{\text{x},1}^{(4)}\int_{\mathbb{R}^{3}}\frac{|\nabla^{2}\rho(\mathbf{r})|^{2}}{\rho(\mathbf{r})^{2}}\text{d}\mathbf{r}+C_{\text{x},2}^{(4)}\int_{\mathbb{R}^{3}}\frac{|\nabla\rho(\mathbf{r})|^{2}\,\nabla^{2}\rho(\mathbf{r})}{\rho(\mathbf{r})^{3}}\text{d}\mathbf{r},$

involving the density gradient $\nabla\rho(\mathbf{r})$ and Laplacian
$\nabla^{2}\rho(\mathbf{r})$. Sham [84] obtained the second-order coefficient
$C_{\text{x},\text{S}}^{(2)}=-7/(432\pi(3\pi^{2})^{1/3})\approx-0.001667$. The
calculation was done by starting with the screened Yukawa interaction
$e^{-\kappa r_{12}}/r_{12}$ and taking the limit $\kappa\to 0$ at the end of
the calculation. It was later shown that this calculation contains an order-
of-limit problem and that the correct Coulombic second-order coefficient is
$C_{\text{x}}^{(2)}=-5/(216\pi(3\pi^{2})^{1/3})\approx-0.002382$ [81, 82]. The
fourth-order coefficients are
$C_{\text{x},1}^{(4)}=-73/(64800\pi^{3})\approx-0.000036$, and
$C_{\text{x},2}^{(4)}\approx 0.00009$, where the last one has been numerically
estimated [83]. Note that each term in Eq. (3.18) correctly fulfils 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 [80] obtained the spin-independent second-order
gradient expansion (GEA2) of the correlation functional

$E_{\text{c}}^{\text{GEA2}}[\rho]=E_{\text{c}}^{\text{LDA}}[\rho]+\int_{\mathbb{R}^{3}}C_{\text{c}}^{(2)}(r_{\text{s}}(\mathbf{r}))\frac{|\nabla\rho(\mathbf{r})|^{2}}{\rho(\mathbf{r})^{4/3}}\text{d}\mathbf{r},$
(3.19)

with a second-order coefficient in the high-density limit
$C_{\text{c},\text{MB}}^{(2)}(r_{\text{s}}=0)=0.004235$. It is believed [85]
that this calculation contains a similar order-of-limit problem as in Sham’s
coefficient $C_{\text{x},\text{S}}^{(2)}$, in such a way that these two
coefficients must be combined to obtain the correct second-order exchange-
correlation coefficient in the high-density limit
$C_{\text{xc}}^{(2)}(r_{\text{s}}=0)=C_{\text{x},\text{S}}^{(2)}+C_{\text{c},\text{MB}}^{(2)}(r_{\text{s}}=0)$.
The correct second-order correlation coefficient in the high-density limit is
then
$C_{\text{c}}^{(2)}(r_{\text{s}}=0)=C_{\text{xc}}^{(2)}(r_{\text{s}}=0)-C_{\text{x}}^{(2)}=0.004950$.
Similarly, the second-order correlation coefficient as a function of
$r_{\text{s}}$ can be obtained by
$C_{\text{c}}^{(2)}(r_{\text{s}})=C_{\text{xc}}^{(2)}(r_{\text{s}})-C_{\text{x}}^{(2)}$
where $C_{\text{xc}}^{(2)}(r_{\text{s}})$ has been parametrized in Ref. [86].
The spin-dependent generalization has the form [87, 88]

$E_{\text{c}}^{\text{GEA2}}[\rho_{\uparrow},\rho_{\downarrow}]=E_{\text{c}}^{\text{LDA}}[\rho_{\uparrow},\rho_{\downarrow}]+\sum_{\sigma,\sigma^{\prime}\in\\{\uparrow,\downarrow\\}}\int_{\mathbb{R}^{3}}C_{\text{c}}^{\sigma,\sigma^{\prime},(2)}(r_{\text{s}}(\mathbf{r}),\zeta(\mathbf{r}))\frac{\nabla\rho_{\sigma}(\mathbf{r})}{\rho_{\sigma}(\mathbf{r})^{2/3}}\cdot\frac{\nabla\rho_{\sigma^{\prime}}(\mathbf{r})}{\rho_{\sigma^{\prime}}(\mathbf{r})^{2/3}}\text{d}\mathbf{r},$
(3.20)

where the functions
$C_{\text{c}}^{\sigma,\sigma^{\prime},(2)}(r_{\text{s}},\zeta)$ have been
numerically calculated in the high-density limit [87, 89].

The GEA should improve over the LDA for sufficiently slowly varying densities.
Since the spin-independent GEA2 exchange energy per particle has the form

$\varepsilon_{\text{x}}^{\text{GEA2}}(\rho,\nabla\rho)=\rho^{1/3}(C_{\text{x}}+C_{\text{x}}^{(2)}x^{2}),$
(3.21)

where $x=|\nabla\rho|/\rho^{4/3}$ is a dimensionless reduced density gradient,
the precise condition for exchange is $x\ll 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,
$\rho(\mathbf{r})\underset{|\mathbf{r}|\to\infty}{\propto}e^{-\alpha|\mathbf{r}|}$,
the reduced density gradient diverges
$x(\mathbf{r})\xrightarrow[|\mathbf{r}|\to\infty]{}\infty$. But this is not as
bad as it seems since $\rho\;\varepsilon_{\text{x}}^{\text{GEA2}}$ 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

$\varepsilon_{\text{c}}^{\text{GEA2}}(\gamma^{3}\rho,\gamma^{4}\nabla\rho)\underset{\gamma\to\infty}{\sim}-A(0)\ln\gamma+\gamma^{1/2}C_{\text{c}}^{(2)}(r_{\text{s}}=0)y^{2},$
(3.22)

where $y=|\nabla\rho|/\rho^{7/6}$ is another reduced density gradient adapted
to correlation. Therefore, in this limit, the GEA2 correlation correction
diverges to $+\infty$ even faster than the LDA diverges to $-\infty$.

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)].
Consequently, truncated gradient expansions should not be directly used but
need to be resummed.

###  3.3 Generalized-gradient approximations

The failures of the GEA lead to the development, which really started in the
1980s, of generalized-gradient approximations (GGAs) with the generic form

$E_{\text{xc}}^{\text{GGA}}[\rho_{\uparrow},\rho_{\downarrow}]=\int_{\mathbb{R}^{3}}e_{\text{xc}}^{\text{GGA}}(\rho_{\uparrow}(\mathbf{r}),\rho_{\downarrow}(\mathbf{r}),\nabla\rho_{\uparrow}(\mathbf{r}),\nabla\rho_{\downarrow}(\mathbf{r}))\text{d}\mathbf{r},$
(3.23)

where $e_{\text{xc}}^{\text{GGA}}$ is some function. The GGAs are often called
semilocal approximations in the sense that $e_{\text{xc}}^{\text{GGA}}$ does
not only use the local value of the spin densities
$\rho_{\uparrow}(\mathbf{r})$ and $\rho_{\downarrow}(\mathbf{r})$ but also
“semilocal information” through its gradients666For generality and simplicity,
we consider here that the GGAs depend on the spin density gradients
$\nabla\rho_{\uparrow}$ and $\nabla\rho_{\downarrow}$, but due to rotational
invariance GGAs actually depend only on the scalar quantities
$(\nabla\rho_{\uparrow})^{2}$, $(\nabla\rho_{\downarrow})^{2}$, and
$\nabla\rho_{\uparrow}\cdot\nabla\rho_{\downarrow}$.
$\nabla\rho_{\uparrow}(\mathbf{r})$ and $\nabla\rho_{\downarrow}(\mathbf{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.

B88 exchange functional

In the Becke 88 (B88 or B) exchange functional [36], the exchange energy
density is written as

$\displaystyle
e_{\text{x}}^{\text{B88}}(\rho_{\uparrow},\rho_{\downarrow},\nabla\rho_{\uparrow},\nabla\rho_{\downarrow})=e_{\text{x}}^{\text{UEG}}(\rho_{\uparrow},\rho_{\downarrow})-\sum_{\sigma\in\\{\uparrow,\downarrow\\}}\rho_{\sigma}^{4/3}\frac{\beta
x_{\sigma}^{2}}{1+6\beta x_{\sigma}\sinh^{-1}(x_{\sigma})},$ (3.24)

where $x_{\sigma}=|\nabla\rho_{\sigma}|/\rho_{\sigma}^{4/3}$. The fact that
$e_{\text{x}}^{\text{B88}}$ depends linearly on $\rho_{\sigma}^{4/3}$ and
nonlinearly only on the dimensionless reduced density gradient
$x_{\sigma}(\mathbf{r})$ imposes the scaling relation of Eq. (2.59). Using the
exponential decay of the ground-state spin densities of Coulombic systems,
$\rho_{\sigma}(\mathbf{r})\underset{|\mathbf{r}|\to\infty}{\propto}e^{-\alpha_{\sigma}|\mathbf{r}|}$,
it can be verified that the chosen form for $e_{\text{x}}^{\text{B88}}$
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)] [90]. For small
$x_{\sigma}$, $e_{\text{x}}^{\text{B88}}$ is correctly quadratic in
$x_{\sigma}$. The parameter $\beta=0.0042$ was found by fitting to HF exchange
energies of rare-gas atoms. A very similar value of $\beta$ can also be found
by imposing the correct coefficient $B_{\text{x}}$ of the large-$Z$ asymptotic
expansion of the exchange energy of neutral atoms [Eq. (2.74)] [63]. It turns
out that imposing the coefficient of the second-order gradient expansion [Eq.
(3.18)] would lead to a value of $\beta$ about twice smaller and would greatly
deteriorate the accuracy of the functional for atoms and molecules.

LYP correlation functional

The Lee-Yang-Parr (LYP) [91] correlation functional is one of the rare
functionals which have not been constructed starting from LDA. It originates
from the Colle-Salvetti [92] 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 [93], giving the following
correlation energy density

$\displaystyle
e_{\text{c}}^{\text{LYP}}(\rho_{\uparrow},\rho_{\downarrow},\nabla\rho_{\uparrow},\nabla\rho_{\downarrow})=-a\frac{4}{1+d\rho^{-1/3}}\frac{\rho_{\uparrow}\rho_{\downarrow}}{\rho}-a\,b\,\omega(\rho)\Biggl{\\{}\rho_{\uparrow}\rho_{\downarrow}\phantom{xxxxxxxxxxxxxxxxxxxxxx}$
$\displaystyle\times\Bigg{[}\sum_{\sigma\in\\{\uparrow,\downarrow\\}}\left(2^{11/3}C_{\text{F}}\rho_{\sigma}^{8/3}-\left(\frac{5}{2}-\frac{\delta(\rho)}{18}\right)|\nabla\rho_{\sigma}|^{2}-\frac{\delta(\rho)-11}{9}\frac{\rho_{\sigma}}{\rho}|\nabla\rho_{\sigma}|^{2}\right)+\left(\frac{47}{18}-\frac{7\delta(\rho)}{18}\right)|\nabla\rho|^{2}\Bigg{]}$
$\displaystyle-\frac{2}{3}\rho^{2}|\nabla\rho|^{2}+\left(\frac{2}{3}\rho^{2}-\rho_{\uparrow}^{2}\right)|\nabla\rho_{\downarrow}|^{2}+\left(\frac{2}{3}\rho^{2}-\rho_{\downarrow}^{2}\right)|\nabla\rho_{\uparrow}|^{2}\Biggl{\\}},\phantom{xxxxx}$
(3.25)

where $\omega(\rho)=\rho^{-11/3}\exp(-c\rho^{-1/3})/(1+d\rho^{-1/3})$,
$\delta(\rho)=c\rho^{-1/3}+d\rho^{-1/3}/(1+d\rho^{-1/3})$, and
$C_{\text{F}}=(3/10)(3\pi^{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 ($\rho_{\uparrow}=0$ or
$\rho_{\downarrow}=0$) and therefore correctly vanishes for one-electron
systems [Eq. (2.39)].

PW91 exchange-correlation functional

The Perdew-Wang 91 (PW91) (see Refs. [67, 94, 95]) exchange-correlation
functional is based on a model of the exchange hole
$h_{\text{x}}(\mathbf{r}_{1},\mathbf{r}_{2})$ in Eq. (2.16) and of the
coupling-constant-integrated correlation hole
$\bar{h}_{\text{c}}(\mathbf{r}_{1},\mathbf{r}_{2})$ 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 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 $|\mathbf{r}_{1}-\mathbf{r}_{2}|$ 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

$\displaystyle
e_{\text{x}}^{\text{PW91}}(\rho,\nabla\rho)=e_{\text{x}}^{\text{UEG}}(\rho)F_{\text{x}}^{\text{PW91}}(s),$
(3.26)

where the so-called enhancement factor is

$\displaystyle
F_{\text{x}}^{\text{PW91}}(s)=\frac{1+0.19645s\sinh^{-1}(7.7956s)+[0.2743-0.1508\exp(-100s^{2})]s^{2}}{1+0.19645s\sinh^{-1}(7.7956s)+0.004s^{4}},$
(3.27)

with the reduced density gradient
$s=|\nabla\rho|/(2k_{\text{F}}\rho)=x/[2(3\pi^{2})^{1/3}]$ where
$k_{\text{F}}=(3\pi^{2}\rho)^{1/3}$ is the Fermi wave vector. The spin-
dependent PW91 exchange energy density is simply obtained from the spin-
scaling relation [Eq. (1.45)]:
$e_{\text{x}}^{\text{PW91}}(\rho_{\uparrow},\rho_{\downarrow},\nabla\rho_{\uparrow},\nabla\rho_{\downarrow})=[e_{\text{x}}^{\text{PW91}}(2\rho_{\uparrow},2\nabla\rho_{\uparrow})+e_{\text{x}}^{\text{PW91}}(2\rho_{\downarrow},2\nabla\rho_{\downarrow})]/2$.
The enhancement factor $F_{\text{x}}^{\text{PW91}}(s)$ satisfies the second-
order gradient expansion [Eq. (3.18)],
$F_{\text{x}}^{\text{PW91}}(s)=1+\mu\;s^{2}+O(s^{4})$ with
$\mu=-16\pi(\pi/3)^{2/3}C_{\text{x}}^{(2)}=10/81\approx 0.1235$, the local
Lieb-Oxford bound, $F_{\text{x}}^{\text{PW91}}(s)\leq-
C_{\text{LO}}/(C_{\text{x}}2^{1/3})\approx 1.804$, which is a sufficient and
necessary condition for a spin-dependent GGA exchange functional to satisfy
the Lieb-Oxford lower bound [Eq. (2.76)] for all densities [70] (note however
that $1.804$ is not an optimal bound), and the condition
$\lim_{s\to\infty}s^{1/2}F_{\text{x}}^{\text{PW91}}(s)<\infty$ which
guarantees the non-uniform scaling finiteness conditions of Eqs. (2.66) and
(2.68) [55, 70].

The PW91 correlation energy density is written as

$\displaystyle
e_{\text{c}}^{\text{PW91}}(\rho_{\uparrow},\rho_{\downarrow},\nabla\rho_{\uparrow},\nabla\rho_{\downarrow})=\rho\left[\varepsilon_{\text{c}}^{\text{UEG}}(\rho_{\uparrow},\rho_{\downarrow})+H^{\text{PW91}}(\rho_{\uparrow},\rho_{\downarrow},t)\right],$
(3.28)

where the gradient correction
$H^{\text{PW91}}(\rho_{\uparrow},\rho_{\downarrow},t)=H_{0}(\rho_{\uparrow},\rho_{\downarrow},t)+H_{1}(\rho_{\uparrow},\rho_{\downarrow},t)$
depends on another reduced density gradient (adapted to correlation)
$t=|\nabla\rho|/(2\phi_{2}(\zeta)k_{\text{s}}\rho)=y/[4\phi_{2}(\zeta)(3/\pi)^{1/6}]$
where $k_{\text{s}}=\sqrt{4k_{\text{F}}/\pi}$ is the Thomas-Fermi screening
wave vector and the spin-scaling function $\phi_{2}(\zeta)$ is defined by Eq.
(3.5), with

$\displaystyle
H_{0}(\rho_{\uparrow},\rho_{\downarrow},t)=\phi_{2}(\zeta)^{3}\frac{\beta^{2}}{2\alpha}\ln\left[1+\frac{2\alpha}{\beta}t^{2}\frac{1+{\cal
A}t^{2}}{1+{\cal A}t^{2}+{\cal A}^{2}t^{4}}\right],$ (3.29)
$\displaystyle{\cal
A}=\frac{2\alpha}{\beta}\left[\exp(-2\alpha\varepsilon_{\text{c}}^{\text{UEG}}(\rho_{\uparrow},\rho_{\downarrow})/(\phi_{2}(\zeta)^{3}\beta^{2}))-1\right]^{-1},$
(3.30) $\displaystyle
H_{1}(\rho_{\uparrow},\rho_{\downarrow},t)=16\left(\frac{3}{\pi}\right)^{1/3}[C_{\text{xc}}(r_{\text{s}})-C_{\text{c},\text{MB}}^{(2)}(r_{\text{s}}=0)-C_{\text{x}}^{(2)}]\phi_{2}(\zeta)^{3}t^{2}e^{-100\phi_{2}(\zeta)^{4}k_{\text{s}}^{2}t^{2}/k_{\text{F}}^{2}},$
(3.31)

and $C_{\text{xc}}(r_{\text{s}})$ is taken from Ref. [86]. The function
$H_{0}(\rho_{\uparrow},\rho_{\downarrow},t)$ was chosen so that it fulfils the
second-order gradient expansion [Eq. (3.20)],
$H_{0}(\rho_{\uparrow},\rho_{\downarrow},t)=\beta\phi_{2}(\zeta)^{3}t^{2}+O(t^{4})$,
using an approximate $\zeta$ dependence [88] and the Ma-Brueckner high-
density-limit second-order coefficient [80]
$\beta=16(3/\pi)^{1/3}C_{\text{c},\text{MB}}^{(2)}(r_{\text{s}}=0)\approx
0.06673$, and so that it cancels the LDA correlation in the large-$t$ limit,
$\lim_{t\to\infty}H_{0}(\rho_{\uparrow},\rho_{\downarrow},t)=-\varepsilon_{\text{c}}^{\text{UEG}}(\rho_{\uparrow},\rho_{\downarrow})$.
The only fitted parameter is $\alpha=0.09$. The function
$H_{1}(\rho_{\uparrow},\rho_{\downarrow},t)$ only serves to restore the
correct second-order gradient expansion, such that
$H^{\text{PW91}}(\rho_{\uparrow},\rho_{\downarrow},t)=16(3/\pi)^{1/3}C_{\text{xc}}(r_{\text{s}})\phi_{2}(\zeta)^{3}t^{2}+O(t^{4})$,
while keeping the large-$t$ limit unchanged.

PBE exchange-correlation functional

The Perdew-Burke-Ernzerhof (PBE) [96] 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

$\displaystyle
e_{\text{x}}^{\text{PBE}}(\rho,\nabla\rho)=e_{\text{x}}^{\text{UEG}}(\rho)F_{\text{x}}^{\text{PBE}}(s),$
(3.32)

where the enhancement factor is

$\displaystyle F_{\text{x}}^{\text{PBE}}(s)=1+\kappa-\frac{\kappa}{1+\mu
s^{2}/\kappa}.$ (3.33)

The function $F_{\text{x}}^{\text{PBE}}(s)$ has the second-order gradient
expansion $F_{\text{x}}^{\text{PBE}}(s)=1+\mu\;s^{2}+O(s^{4})$, and the
parameter is chosen as
$\mu=16\pi(\pi/3)^{2/3}C_{\text{c},\text{MB}}^{(2)}(r_{\text{s}}=0)\approx
0.21951$ so as to cancel the correlation second-order gradient expansion. The
second parameter $\kappa$ is chosen so as to saturate the local Lieb-Oxford
bound, i.e.
$\lim_{s\to\infty}F_{\text{x}}^{\text{PBE}}(s)=1+\kappa=-C_{\text{LO}}/(C_{\text{x}}2^{1/3})\approx
1.804$, leading to $\kappa=0.804$. The same exchange functional form was in
fact proposed earlier in the Becke 86 (B86) functional [97] with empirical
parameters ($\mu=0.235$, $\kappa=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
$\kappa=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 [98]. 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
$F_{\text{x}}^{\text{RPBE}}(s)=1+\kappa(1-\exp(-\mu s^{2}/\kappa))$ with the
same parameters as in the original PBE [99].

The PBE correlation energy density is written as

$\displaystyle
e_{\text{c}}^{\text{PBE}}(\rho_{\uparrow},\rho_{\downarrow},\nabla\rho_{\uparrow},\nabla\rho_{\downarrow})=\rho\left[\varepsilon_{\text{c}}^{\text{UEG}}(\rho_{\uparrow},\rho_{\downarrow})+H^{\text{PBE}}(\rho_{\uparrow},\rho_{\downarrow},t)\right],$
(3.34)

with the gradient correction

$\displaystyle
H^{\text{PBE}}(\rho_{\uparrow},\rho_{\downarrow},t)=A(0)\phi_{2}(\zeta)^{3}\ln\left[1+\frac{\beta}{A(0)}t^{2}\frac{1+{\cal
A}t^{2}}{1+{\cal A}t^{2}+{\cal A}^{2}t^{4}}\right],$ (3.35)

and

$\displaystyle{\cal
A}=\frac{\beta}{A(0)}\left[\exp(-\varepsilon_{\text{c}}^{\text{UEG}}(\rho_{\uparrow},\rho_{\downarrow})/(A(0)\phi_{2}(\zeta)^{3}))-1\right]^{-1},$
(3.36)

As in the PW91 correlation functional, the function
$H^{\text{PBE}}(\rho_{\uparrow},\rho_{\downarrow},t)$ has the second-order
gradient expansion
$H^{\text{PBE}}(\rho_{\uparrow},\rho_{\downarrow},t)=\beta\phi_{2}(\zeta)^{3}t^{2}+O(t^{4})$
where
$\beta=16(3/\pi)^{1/3}C_{\text{c},\text{MB}}^{(2)}(r_{\text{s}}=0)\approx
0.06673$, and it cancels the LDA correlation in the large-$t$ limit,
$\lim_{t\to\infty}H^{\text{PBE}}(\rho_{\uparrow},\rho_{\downarrow},t)=-\varepsilon_{\text{c}}^{\text{UEG}}(\rho_{\uparrow},\rho_{\downarrow})$.
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^{\text{PBE}}(\gamma^{3}\rho_{\uparrow},\gamma^{3}\rho_{\downarrow},\gamma^{1/2}t)\underset{\gamma\to\infty}{\sim}A(0)\phi_{2}(\zeta)^{3}\;\ln\gamma$
where $A(0)\phi_{2}(\zeta)^{3}$ is a good approximation to the coefficient
$A(\zeta)$ [88].

A variant of the PBE exchange-correlation functional, called PBEsol [100],
targeted for solid-state systems was proposed where the correct second-order
exchange gradient-expansion coefficient is restored, i.e.
$\mu_{\text{PBEsol}}=16\pi(\pi/3)^{2/3}C_{\text{x}}^{(2)}=10/81\approx
0.1235$, and the second-order correlation gradient-expansion coefficient
$\beta_{\text{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 [101] (see Section 4.1). The B97-GGA exchange energy
density is

$\displaystyle
e_{\text{x}}^{\text{B97-GGA}}(\rho_{\uparrow},\rho_{\downarrow},\nabla\rho_{\uparrow},\nabla\rho_{\downarrow})=\sum_{\sigma\in\\{\uparrow,\downarrow\\}}e_{\text{x},\sigma}^{\text{UEG}}(\rho_{\sigma})\;g_{\text{x}}(x_{\sigma}),$
(3.37)

where
$e_{\text{x},\sigma}^{\text{UEG}}(\rho_{\sigma})=e_{\text{x}}^{\text{UEG}}(\rho_{\sigma},0)$
is the spin-$\sigma$ contribution to the UEG exchange energy density and the
gradient correction $g_{\text{x}}(x_{\sigma})$ is a function of
$x_{\sigma}=|\nabla\rho_{\sigma}|/\rho_{\sigma}^{4/3}$

$\displaystyle
g_{\text{x}}(x_{\sigma})=\sum_{i=0}^{m}c_{\text{x},i}\;u_{\text{x}}(x_{\sigma})^{i},$
(3.38)

with
$u_{\text{x}}(x_{\sigma})=\gamma_{\text{x}}x_{\sigma}^{2}/(1+\gamma_{\text{x}}x_{\sigma}^{2})$.
The B97-GGA correlation energy density is written as the sum of opposite- and
same-spin contributions

$\displaystyle
e_{\text{c}}^{\text{B97-GGA}}(\rho_{\uparrow},\rho_{\downarrow},\nabla\rho_{\uparrow},\nabla\rho_{\downarrow})=e_{\text{c},\uparrow\downarrow}^{\text{B97-GGA}}(\rho_{\uparrow},\rho_{\downarrow},\nabla\rho_{\uparrow},\nabla\rho_{\downarrow})+\sum_{\sigma\in\\{\uparrow,\downarrow\\}}e_{\text{c},\sigma\sigma}^{\text{B97-GGA}}(\rho_{\sigma},\nabla\rho_{\sigma}),$
(3.39)

where

$\displaystyle
e_{\text{c},\uparrow\downarrow}^{\text{B97-GGA}}(\rho_{\uparrow},\rho_{\downarrow},\nabla\rho_{\uparrow},\nabla\rho_{\downarrow})=e_{\text{c},\uparrow\downarrow}^{\text{UEG}}(\rho_{\uparrow},\rho_{\downarrow})\;g_{\text{c},\uparrow\downarrow}(x_{\uparrow\downarrow}),$
(3.40)

and

$\displaystyle
e_{\text{c},\sigma\sigma}^{\text{B97-GGA}}(\rho_{\sigma},\nabla\rho_{\sigma})=e_{\text{c},\sigma\sigma}^{\text{UEG}}(\rho_{\sigma})\;g_{\text{c},\sigma\sigma}(x_{\sigma}).$
(3.41)

In these expressions,
$e_{\text{c},\uparrow\downarrow}^{\text{UEG}}(\rho_{\uparrow},\rho_{\downarrow})=e_{\text{c}}^{\text{UEG}}(\rho_{\uparrow},\rho_{\downarrow})-e_{\text{c}}^{\text{UEG}}(\rho_{\uparrow},0)-e_{\text{c}}^{\text{UEG}}(\rho_{\downarrow},0)$
and
$e_{\text{c},\sigma\sigma}^{\text{UEG}}(\rho_{\alpha})=e_{\text{c}}^{\text{UEG}}(\rho_{\sigma},0)$
are estimations of the opposite- and same-spin contributions to the UEG
correlation energy density [102, 103]. The opposite-spin gradient correction
is taken as a function of
$x_{\uparrow\downarrow}=\sqrt{(x_{\uparrow}^{2}+x_{\downarrow}^{2})/2}$

$\displaystyle
g_{\text{c},\uparrow\downarrow}(x_{\uparrow\downarrow})=\sum_{i=0}^{m}c_{\text{c},i}^{\uparrow\downarrow}\;u_{\text{c}}^{\uparrow\downarrow}(x_{\uparrow\downarrow})^{i},$
(3.42)

with
$u_{\text{c}}^{\uparrow\downarrow}(x_{\uparrow\downarrow})=\gamma_{\text{c}}^{\uparrow\downarrow}x_{\uparrow\downarrow}^{2}/(1+\gamma_{\text{c}}^{\uparrow\downarrow}x_{\uparrow\downarrow}^{2})$,
and the same-spin gradient correction is

$\displaystyle
g_{\text{c},\sigma\sigma}(x_{\sigma})=\sum_{i=0}^{m}c_{\text{c},i}^{\sigma\sigma}\;u_{\text{c}}^{\sigma\sigma}(x_{\sigma})^{i},$
(3.43)

with
$u_{\text{c}}^{\sigma\sigma}(x_{\sigma})=\gamma_{\text{c}}^{\sigma\sigma}x_{\sigma}^{2}/(1+\gamma_{\text{c}}^{\sigma\sigma}x_{\sigma}^{2})$.
The parameters $\gamma_{\text{x}}=0.004$,
$\gamma_{\text{c}}^{\uparrow\downarrow}=0.006$, and
$\gamma_{\text{c}}^{\sigma\sigma}=0.2$, were roughly optimized on atomic
exchange and correlation energies. The other parameters $c_{\text{x},i}$,
$c_{\text{c},i}^{\uparrow\downarrow}$, $c_{\text{c},i}^{\sigma\sigma}$ 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) [104] exchange-correlation functional
uses the same form as the B97-GGA exchange-correlation functional but with a
polynomial degree $m=4$ and the parameters $c_{\text{x},i}$,
$c_{\text{c},i}^{\uparrow\downarrow}$, $c_{\text{c},i}^{\sigma\sigma}$ were
optimized without HF exchange a set of energetic properties (atomic total
energies, ionization energies, atomization energies), nuclear gradients, and
accurate exchange-correlation potentials.

###  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,

$E_{\text{xc}}^{\text{mGGA}}[\rho,\tau]=\int_{\mathbb{R}^{3}}e_{\text{xc}}^{\text{mGGA}}(\rho(\mathbf{r}),\nabla\rho(\mathbf{r}),\nabla^{2}\rho(\mathbf{r}),\tau(\mathbf{r}))\text{d}\mathbf{r},$
(3.44)

i.e., they use more ingredients than the GGAs, namely the Laplacian of the
density $\nabla^{2}\rho(\mathbf{r})$ and/or the non-interacting positive
kinetic energy density $\tau(\mathbf{r})$ associated to a single-determinant
wave function $\Phi$

$\displaystyle\tau(\mathbf{r})=\tau_{\Phi}(\mathbf{r})$ $\displaystyle=$
$\displaystyle\frac{N}{2}\int_{\\{\uparrow,\downarrow\\}\times(\mathbb{R}^{3}\times\\{\uparrow,\downarrow\\})^{N-1}}\left|\nabla_{\mathbf{r}}\Phi(\mathbf{x},\mathbf{x}_{2},...,\mathbf{x}_{N})\right|^{2}\text{d}\sigma\text{d}\mathbf{x}_{2}...\text{d}\mathbf{x}_{N}$
(3.45) $\displaystyle=$
$\displaystyle\frac{1}{2}\sum_{i=1}^{N}\left|\nabla\varphi_{i}(\mathbf{r})\right|^{2},$

where $\\{\varphi_{i}\\}_{i=1,...,N}$ are the orbitals occupied in $\Phi$. The
meta-GGAs are considered as part of the family of semilocal approximations, in
the sense that $\tau(\mathbf{r})$ contains semilocal information with respect
to the orbitals.

Meta-GGAs can be viewed as implicit functionals of the density only, i.e.
$E_{\text{xc}}^{\text{mGGA}}[\rho,\tau_{\Phi[\rho]}]$, since
$\tau(\mathbf{r})$ can be considered itself as an implicit functional of the
density via the KS single-determinant wave function $\Phi[\rho]$. This view in
which $E_{\text{xc}}^{\text{mGGA}}[\rho]$ is a proper approximation to the
exchange-correlation density functional $E_{\text{xc}}[\rho]$ of the KS scheme
is normally adopted when constructing meta-GGAs approximations. However, the
calculation of the functional derivative of
$E_{\text{xc}}^{\text{mGGA}}[\rho]$ 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 $\Phi$, i.e.
$E_{\text{xc}}^{\text{mGGA}}[\rho_{\Phi},\tau_{\Phi}]$, [105, 106, 107, 108,
109, 110, 111, 112] or, in other words, approximations to an exact GKS
exchange-correlation functional (see Section 1.4).

In the latter approach, that we will here refer to as the meta-Kohn-Sham (mKS)
scheme, we introduce a functional $E_{\text{xc}}^{\text{mKS}}[\rho,\tau]$ (to
which meta-GGAs are approximations) defined for $\rho$ and $\tau$
simultaneously representable by a single-determinant wave function
$\Phi\in{\cal S}^{N}$ and which defines the GKS functional
$E_{\text{xc}}^{S}[\Phi]=E_{\text{xc}}^{\text{mKS}}[\rho_{\Phi},\tau_{\Phi}]$
[see Eq. (1.47)] giving the exact ground-state energy via Eq. (1.48)

$\displaystyle E_{0}=\inf_{\Phi\in{\cal
S}^{N}}\left\\{\langle\Phi|\hat{T}+\hat{V}_{\text{ne}}|\Phi\rangle+E_{\text{H}}[\rho_{\Phi}]+E_{\text{xc}}^{\text{mKS}}[\rho_{\Phi},\tau_{\Phi}]\right\\},$
(3.46)

which, by taking variations with respect to the orbitals, gives the mKS
equations

$\left(-\frac{1}{2}\nabla^{2}+v_{\text{ne}}(\mathbf{r})+v_{\text{H}}(\mathbf{r})+v_{\text{xc}}^{\text{mKS}}(\mathbf{r})\right)\varphi_{i}(\mathbf{r})=\varepsilon_{i}\varphi_{i}(\mathbf{r}),$
(3.47)

where
$v_{\text{xc}}^{\text{mKS}}(\mathbf{r})=v_{\text{xc},1}^{\text{mKS}}(\mathbf{r})+v_{\text{xc},2}^{\text{mKS}}(\mathbf{r})$
contains a usual local potential

$v_{\text{xc},1}^{\text{mKS}}(\mathbf{r})=\frac{\delta
E_{\text{xc}}^{\text{mKS}}[\rho,\tau]}{\delta\rho(\mathbf{r})},$ (3.48)

and a non-multiplicative operator [107, 108, 109, 110, 111]

$v_{\text{xc},2}^{\text{mKS}}(\mathbf{r})=-\frac{1}{2}\nabla\cdot\left(\frac{\delta
E_{\text{xc}}^{\text{mKS}}[\rho,\tau]}{\delta\tau(\mathbf{r})}\nabla\right),$
(3.49)

evaluated with $\rho(\mathbf{r})=\sum_{i=1}^{N}|\varphi_{i}(\mathbf{r})|^{2}$
and $\tau(\mathbf{r})=(1/2)\sum_{i=1}^{N}|\nabla\varphi_{i}(\mathbf{r})|^{2}$.
Interestingly, the mKS equations can be rewritten as a Schrödinger-like
equation with a position-dependent mass $m(\mathbf{r})$ [113]

$\left(-\frac{1}{2}\nabla\cdot\frac{1}{m(\mathbf{r})}\nabla+v_{\text{ne}}(\mathbf{r})+v_{\text{H}}(\mathbf{r})+v_{\text{xc},1}^{\text{mKS}}(\mathbf{r})\right)\varphi_{i}(\mathbf{r})=\varepsilon_{i}\varphi_{i}(\mathbf{r}),$
(3.50)

where $m(\mathbf{r})=\left(1+\delta
E_{\text{xc}}^{\text{mKS}}[\rho,\tau]/\delta\tau(\mathbf{r})\right)^{-1}$. As
in the KS scheme, the functional $E_{\text{xc}}^{\text{mKS}}[\rho,\tau]$ is
decomposed into exchange and correlation contributions:
$E_{\text{xc}}^{\text{mKS}}[\rho,\tau]=E_{\text{x}}^{\text{mKS}}[\rho,\tau]+E_{\text{c}}^{\text{mKS}}[\rho,\tau]$.
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
$E_{\text{xc}}^{\text{mKS}}[\rho_{\uparrow},\rho_{\downarrow},\tau_{\uparrow},\tau_{\downarrow}]$
and the spin-scaling relation of Eq. (1.45) is generalized to

$E_{\text{x}}^{\text{mKS}}[\rho_{\uparrow},\rho_{\downarrow},\tau_{\uparrow},\tau_{\downarrow}]=\frac{1}{2}\left(E_{\text{x}}^{\text{mKS}}[2\rho_{\uparrow},2\tau_{\uparrow}]+E_{\text{x}}^{\text{mKS}}[2\rho_{\downarrow},2\tau_{\downarrow}]\right).$
(3.51)

Correspondingly, the spin-dependent version of the meta-GGAs are formulated in
terms of the spin-resolved quantities $\rho_{\uparrow}$, $\rho_{\downarrow}$,
$\nabla\rho_{\uparrow}$, $\nabla\rho_{\downarrow}$,
$\nabla^{2}\rho_{\uparrow}$, $\nabla^{2}\rho_{\downarrow}$, $\tau_{\uparrow}$,
and $\tau_{\downarrow}$.

One motivation for the introduction of the variable $\tau(\mathbf{r})$ is that
it appears in the expansion of the spherically averaged exchange hole
[entering in Eq. (2.16)] for small interelectronic distances $r_{12}$ [114],
which for the case of a closed-shell system is

$\displaystyle\frac{1}{4\pi
r_{12}^{2}}\int_{S(\mathbf{0},r_{12})}h_{\text{x}}(\mathbf{r}_{1},\mathbf{r}_{1}+\mathbf{r}_{12})\text{d}\mathbf{r}_{12}=-\frac{\rho(\mathbf{r}_{1})}{2}\phantom{xxxxxxxxxxxxxxxxxxx}$
$\displaystyle-\frac{1}{3}\left(\frac{1}{4}\nabla^{2}\rho(\mathbf{r}_{1})-\tau(\mathbf{r}_{1})+\frac{|\nabla\rho(\mathbf{r}_{1})|^{2}}{8\rho(\mathbf{r}_{1})}\right)r_{12}^{2}+O(r_{12}^{4}),$
(3.52)

where $S(\mathbf{0},r_{12})$ designates the sphere centered at $\mathbf{0}$
and of radius $r_{12}=|\mathbf{r}_{12}|$. Thus $\tau(\mathbf{r})$ is needed to
describe the curvature of the exchange hole.

Another important motivation is that $\tau(\mathbf{r})$ is useful to identify
different types of spatial regions of electronic systems [115]. This is done
by comparing $\tau(\mathbf{r})$ with the von Weizsäcker kinetic energy density

$\tau^{\text{W}}(\mathbf{r})=\frac{|\nabla\rho(\mathbf{r})|^{2}}{8\rho(\mathbf{r})},$
(3.53)

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

$z(\mathbf{r})=\frac{\tau^{\text{W}}(\mathbf{r})}{\tau(\mathbf{r})},$ (3.54)

which takes its values in the range $[0,1]$ [116], identifies one-orbital
regions ($z=1$). A better indicator is

$\alpha(\mathbf{r})=\frac{\tau(\mathbf{r})-\tau^{\text{W}}(\mathbf{r})}{\tau^{\text{UEG}}(\mathbf{r})},$
(3.55)

where
$\tau^{\text{UEG}}(\mathbf{r})=(3/10)(3\pi^{2})^{2/3}\rho(\mathbf{r})^{5/3}$
is the non-interacting kinetic energy density of the UEG. This indicator
$\alpha(\mathbf{r})$ distinguishes one-orbital regions ($\alpha=0$), slowly
varying density regions ($\alpha\approx 1$), and regions of density overlap
between closed shells that characterize noncovalent bonds ($\alpha\gg 1$).

Nowadays, $\nabla^{2}\rho(\mathbf{r})$ is rarely used to construct meta-GGAs
because it contains similar information than $\tau(\mathbf{r})$, which can be
seen by the second-order gradient expansion of $\tau(\mathbf{r})$ [117]

$\tau^{\text{GEA2}}(\mathbf{r})=\tau^{\text{UEG}}(\mathbf{r})+\frac{1}{72}\frac{|\nabla\rho(\mathbf{r})|^{2}}{\rho(\mathbf{r})}+\frac{1}{6}\nabla^{2}\rho(\mathbf{r}).$
(3.56)

In comparison to GGAs, meta-GGAs are more versatile and generally constitute
an improvement. Significantly, thanks to the use of $\tau$, 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) [118, 119] functional, the
exchange energy density is written as

$\displaystyle
e_{\text{x}}^{\text{TPSS}}(\rho,\nabla\rho,\tau)=e_{\text{x}}^{\text{UEG}}(\rho)F_{\text{x}}^{\text{TPSS}}(s,z),$
(3.57)

where the enhancement factor is a function of
$s=|\nabla\rho|/(2k_{\text{F}}\rho)$ and $z=\tau^{\text{W}}/\tau$

$\displaystyle
F_{\text{x}}^{\text{TPSS}}(s,z)=1+\kappa-\frac{\kappa}{1+x^{\text{TPSS}}(s,z)/\kappa},$
(3.58)

with $\kappa=0.804$ so as to saturate the local Lieb-Oxford bound (just like
in the PBE exchange functional) and

$\displaystyle x^{\text{TPSS}}(s,z)$ $\displaystyle=$
$\displaystyle\Biggl{[}\left(\frac{10}{81}+c\frac{z^{2}}{(1+z^{2})^{2}}\right)s^{2}+\frac{146}{2025}\tilde{q}_{b}^{2}-\frac{73}{405}\tilde{q}_{b}\sqrt{\frac{1}{2}\left(\frac{3}{5}z\right)^{2}+\frac{1}{2}s^{4}}$
(3.59)
$\displaystyle+\frac{1}{\kappa}\left(\frac{10}{81}\right)^{2}s^{4}+2\sqrt{e}\frac{10}{81}\left(\frac{3}{5}z\right)^{2}+e\mu
s^{6}\Biggl{]}/\left(1+\sqrt{e}s^{2}\right)^{2},$

and $\tilde{q}_{b}=(9/20)(\alpha-1)/[1+b\alpha(\alpha-1)]^{1/2}+2s^{2}/3$
(where $\alpha=(\tau-\tau^{\text{W}})/\tau^{\text{UEG}}=(5s^{2}/3)(z^{-1}-1)$)
is a quantity that tends to the reduced density Laplacian
$q=\nabla^{2}\rho/(4k_{\text{F}}^{2}\rho)$ in the slowly varying density limit
[using Eq. (3.56)]. The function $x^{\text{TPSS}}(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
$F_{\text{x}}^{\text{GEA4}}(s,z)=1+(10/81)s^{2}+(146/2025)q^{2}-(73/405)s^{2}q$.
The constant $\mu=0.21951$ is chosen to retain the same large-$s$ behavior of
the PBE exchange functional, i.e.
$F_{\text{x}}^{\text{TPSS}}(s,z)\underset{s\to\infty}{\sim}F_{\text{x}}^{\text{PBE}}(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
$F_{\text{x}}^{\text{TPSS}}(s,z)$ a monotonically increasing function of $s$.

The TPSS correlation functional is constructed by making minor refinements to
the previously developed Perdew-Kurth-Zupan-Blaha (PKZB) [120] meta-GGA
correlation functional

$\displaystyle
e_{\text{c}}^{\text{TPSS}}(\rho_{\uparrow},\rho_{\downarrow},\nabla\rho_{\uparrow},\nabla\rho_{\downarrow},\tau_{\uparrow},\tau_{\downarrow})$
$\displaystyle=$
$\displaystyle\rho\;\varepsilon_{\text{c}}^{\text{revPKZB}}(\rho_{\uparrow},\rho_{\downarrow},\nabla\rho_{\uparrow},\nabla\rho_{\downarrow},\tau_{\uparrow},\tau_{\downarrow})$
(3.60)
$\displaystyle\times\left[1+d\;\varepsilon_{\text{c}}^{\text{revPKZB}}(\rho_{\uparrow},\rho_{\downarrow},\nabla\rho_{\uparrow},\nabla\rho_{\downarrow},\tau_{\uparrow},\tau_{\downarrow})z^{3}\right],$

where the revised PKZB correlation energy per particle is

$\displaystyle\varepsilon_{\text{c}}^{\text{revPKZB}}(\rho_{\uparrow},\rho_{\downarrow},\nabla\rho_{\uparrow},\nabla\rho_{\downarrow},\tau_{\uparrow},\tau_{\downarrow})=\varepsilon_{\text{c}}^{\text{PBE}}(\rho_{\uparrow},\rho_{\downarrow},\nabla\rho_{\uparrow},\nabla\rho_{\downarrow})\left[1+C(\zeta,\xi)z^{2}\right]$
$\displaystyle-\left[1+C(\zeta,\xi)\right]z^{2}\sum_{\sigma\in\\{\uparrow,\downarrow\\}}\frac{\rho_{\sigma}}{\rho}\tilde{\varepsilon}_{\text{c},\sigma}^{\text{PBE}}(\rho_{\uparrow},\rho_{\downarrow},\nabla\rho_{\uparrow},\nabla\rho_{\downarrow}),$
(3.61)

with
$\tilde{\varepsilon}_{\text{c},\sigma}^{\text{PBE}}(\rho_{\uparrow},\rho_{\downarrow},\nabla\rho_{\uparrow},\nabla\rho_{\downarrow})=\max[\varepsilon_{\text{c}}^{\text{PBE}}(\rho_{\sigma},0,\nabla\rho_{\sigma},0),\varepsilon_{\text{c}}^{\text{PBE}}(\rho_{\uparrow},\rho_{\downarrow},\nabla\rho_{\uparrow},\nabla\rho_{\downarrow})]$
where
$\varepsilon_{\text{c}}^{\text{PBE}}(\rho_{\uparrow},\rho_{\downarrow},\nabla\rho_{\uparrow},\nabla\rho_{\downarrow})$
is the PBE correlation energy per particle. Equation (3.61) constitutes a one-
electron self-interaction correction on the PBE correlation functional.
Indeed, for one-electron densities we have $z=1$ and $\zeta=\pm 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 limits,
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 [121] over the
range of valence-electron bulk densities. The rest of the function is taken as

$\displaystyle
C(\zeta,\xi)=\frac{0.53+0.87\zeta^{2}+0.50\zeta^{4}+2.26\zeta^{6}}{(1+\xi^{2}[(1+\zeta)^{-4/3}+(1-\zeta)^{-4/3}]/2)^{4}},$
(3.62)

where $\xi=|\nabla\zeta|/(2k_{\text{F}})$ is a reduced spin-polarization
gradient. The function $C(\zeta,\xi)$ is chosen so as to make the exchange-
correlation energy independent on the spin polarization $\zeta$ 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 [122], the
exchange energy density is written as

$\displaystyle
e_{\text{x}}^{\text{M06-L}}(\rho_{\uparrow},\rho_{\downarrow},\nabla\rho_{\uparrow},\nabla\rho_{\downarrow},\tau_{\uparrow},\tau_{\downarrow})=\sum_{\sigma\in\\{\uparrow,\downarrow\\}}e_{\text{x},\sigma}^{\text{PBE}}(\rho_{\sigma},\nabla\rho_{\sigma})f(w_{\sigma})+e_{\text{x},\sigma}^{\text{UEG}}(\rho_{\sigma})h_{\text{x}}(x_{\sigma},Z_{\sigma}).$
(3.63)

The first term in Eq. (3.63), which has the same form as in the previously
developed M05 exchange functional [123], contains the spin-$\sigma$ PBE
exchange energy density
$e_{\text{x},\sigma}^{\text{PBE}}(\rho_{\sigma},\nabla\rho_{\sigma})=e_{\text{x}}^{\text{PBE}}(\rho_{\sigma},0,\nabla\rho_{\sigma},0)$
and the kinetic-energy density correction factor

$\displaystyle f(w_{\sigma})=\sum_{i=0}^{11}a_{i}w_{\sigma}^{i},$ (3.64)

where
$w_{\sigma}=(\tau_{\sigma}^{\text{UEG}}/\tau_{\sigma}-1)/(\tau_{\sigma}^{\text{UEG}}/\tau_{\sigma}+1)$
with $\tau_{\sigma}^{\text{UEG}}=(3/10)(6\pi^{2})^{2/3}\rho_{\sigma}^{5/3}$ is
an indicator of the delocalization of the exchange hole [124]. The second term
in Eq. (3.63), which has the same form as in the VS98 exchange functional
[125], contains the spin-$\sigma$ UEG exchange energy density
$e_{\text{x},\sigma}^{\text{UEG}}(\rho_{\sigma})=e_{\text{x}}^{\text{UEG}}(\rho_{\sigma},0)$
and the correction factor

$\displaystyle
h_{\text{x}}(x_{\sigma},Z_{\sigma})=h(x_{\sigma},Z_{\sigma},d_{\text{x},0},d_{\text{x},1},d_{\text{x},2},d_{\text{x},3},d_{\text{x},4},\alpha_{\text{x}}),$
(3.65)

where $x_{\sigma}=|\nabla\rho_{\sigma}|/\rho_{\sigma}^{4/3}$ and
$Z_{\sigma}=2(\tau_{\sigma}-\tau_{\sigma}^{\text{UEG}})/\rho_{\sigma}^{5/3}$
and $h$ is the parametrized function

$\displaystyle
h(x,Z,d_{0},d_{1},d_{2},d_{3},d_{4},\alpha)=\frac{d_{0}}{\gamma(x,Z,\alpha)}+\frac{d_{1}x^{2}+d_{2}Z}{\gamma(x,Z,\alpha)^{2}}+\frac{d_{3}x^{4}+d_{4}x^{2}Z}{\gamma(x,Z,\alpha)^{3}},$
(3.66)

with $\gamma(x,Z,\alpha)=1+\alpha(x^{2}+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)],

$\displaystyle
e_{\text{c}}^{\text{M06-L}}(\rho_{\uparrow},\rho_{\downarrow},\nabla\rho_{\uparrow},\nabla\rho_{\downarrow},\tau_{\uparrow},\tau_{\downarrow})$
$\displaystyle=$ $\displaystyle
e_{\text{c},\uparrow\downarrow}^{\text{M06-L}}(\rho_{\uparrow},\rho_{\downarrow},\nabla\rho_{\uparrow},\nabla\rho_{\downarrow},\tau_{\uparrow},\tau_{\downarrow})$
(3.67)
$\displaystyle+\sum_{\sigma\in\\{\uparrow,\downarrow\\}}e_{\text{c},\sigma\sigma}^{\text{M06-L}}(\rho_{\sigma},\nabla\rho_{\sigma},\tau_{\sigma}),$

where

$\displaystyle
e_{\text{c},\uparrow\downarrow}^{\text{M06}}(\rho_{\uparrow},\rho_{\downarrow},\nabla\rho_{\uparrow},\nabla\rho_{\downarrow},\tau_{\uparrow},\tau_{\downarrow})=e_{\text{c},\uparrow\downarrow}^{\text{UEG}}(\rho_{\uparrow},\rho_{\downarrow})\;\left[g_{\text{c},\uparrow\downarrow}(x_{\uparrow\downarrow})+h_{\text{c},\uparrow\downarrow}(x_{\uparrow\downarrow},Z_{\uparrow\downarrow})\right],$
(3.68)

and

$\displaystyle
e_{\text{c},\sigma\sigma}^{\text{M06-L}}(\rho_{\sigma},\nabla\rho_{\sigma})=e_{\text{c},\sigma\sigma}^{\text{UEG}}(\rho_{\sigma})\;\left[g_{\text{c},\sigma\sigma}(x_{\sigma})+h_{\text{c},\sigma\sigma}(x_{\sigma},Z_{\sigma})\right]D_{\sigma}(z_{\sigma}),$
(3.69)

where the spin-decomposed UEG correlation energies
$e_{\text{c},\uparrow\downarrow}^{\text{UEG}}(\rho_{\uparrow},\rho_{\downarrow})$
and $e_{\text{c},\sigma\sigma}^{\text{UEG}}(\rho_{\sigma})$ were already
defined after Eq. (3.41), and the gradient corrections
$g_{\text{c},\uparrow\downarrow}(x_{\uparrow\downarrow})$ and
$g_{\text{c},\sigma\sigma}(x_{\sigma})$ are given in Eqs. (3.42) and (3.43).
The additional correction factors are

$\displaystyle
h_{\text{c}}^{\uparrow\downarrow}(x_{\uparrow\downarrow},Z_{\uparrow\downarrow})=h(x_{\uparrow\downarrow},Z_{\uparrow\downarrow},d_{\text{c},0}^{\uparrow\downarrow},d_{\text{c},1}^{\uparrow\downarrow},d_{\text{c},2}^{\uparrow\downarrow},d_{\text{c},3}^{\uparrow\downarrow},d_{\text{c},4}^{\uparrow\downarrow},\alpha_{\text{c}}^{\uparrow\downarrow}),$
(3.70)

where $x_{\uparrow\downarrow}=\sqrt{(x_{\uparrow}^{2}+x_{\downarrow}^{2})/2}$,
$Z_{\uparrow\downarrow}=Z_{\uparrow}+Z_{\downarrow}$, and

$\displaystyle
h_{\text{c}}^{\sigma\sigma}(x_{\sigma},Z_{\sigma})=h(x_{\sigma},Z_{\sigma},d_{\text{c},0}^{\sigma\sigma},d_{\text{c},1}^{\sigma\sigma},d_{\text{c},2}^{\sigma\sigma},d_{\text{c},3}^{\sigma\sigma},d_{\text{c},4}^{\sigma\sigma},\alpha_{\text{c}}^{\sigma\sigma}).$
(3.71)

The factor $D_{\sigma}(z_{\sigma})=1-z_{\sigma}$ where
$z_{\sigma}=\tau^{\text{W}}_{\sigma}/\tau_{\sigma}$ and
$\tau^{\text{W}}_{\sigma}=|\nabla\rho_{\sigma}|^{2}/(8\rho_{\sigma})$ ensures
that the correlation energy correctly vanishes for one-electron systems [126].

The parameters $\gamma_{\text{c}}^{\uparrow\downarrow}=0.0031$, and
$\gamma_{\text{c}}^{\sigma\sigma}=0.06$, were optimized on the correlation
energies of He and Ne. The parameters $\alpha_{\text{x}}=0.001867$,
$\alpha_{\text{c}}^{\uparrow\downarrow}=0.003050$, and
$\alpha_{\text{c}}^{\sigma\sigma}=0.005151$ were taken from Ref. [125]. The
constraints $a_{0}+d_{\text{x},0}=1$,
$c_{\text{c},0}^{\uparrow\downarrow}+d_{\text{c},0}^{\uparrow\downarrow}=1$,
and $c_{\text{c},0}^{\sigma\sigma}+d_{\text{c},0}^{\sigma\sigma}=1$ are
enforced to obtain the correct UEG limit. The remaining 34 free parameters
$a_{i}$, $c_{\text{c},i}^{\uparrow\downarrow}$,
$c_{\text{c},i}^{\sigma\sigma}$ for a polynomial degree $m=4$ in Eqs. (3.64),
(3.42), and (3.43), and $d_{\text{x},i}$,
$d_{\text{c},i}^{\uparrow\downarrow}$, $d_{\text{c},i}^{\sigma\sigma}$ in Eqs.
(3.65), (3.70), and (3.71) were optimized on a large set of diverse
physicochemical properties concerning main-group thermochemistry, reaction
barrier heights, noncovalent interactions, electronic spectroscopy, and
transition metal bonding.

SCAN exchange-correlation functional

In the SCAN (strongly constrained and appropriately normed) [127] exchange-
correlation functional, the exchange energy density is written as

$\displaystyle
e_{\text{x}}^{\text{SCAN}}(\rho,\nabla\rho,\tau)=e_{\text{x}}^{\text{UEG}}(\rho)F_{\text{x}}^{\text{SCAN}}(s,\alpha),$
(3.72)

where the enhancement factor is a function of
$s=|\nabla\rho|/(2k_{\text{F}}\rho)$ and
$\alpha=(\tau-\tau^{\text{W}})/\tau^{\text{UEG}}$

$\displaystyle
F_{\text{x}}^{\text{SCAN}}(s,\alpha)=[h_{\text{x}}^{1}(s,\alpha)+f_{\text{x}}(\alpha)(h_{\text{x}}^{0}-h_{\text{x}}^{1}(s,\alpha))]g_{\text{x}}(s).$
(3.73)

which interpolates between $\alpha=0$ and $\alpha\approx 1$, and extrapolate
to $\alpha\to\infty$ using the function

$\displaystyle
f_{\text{x}}(\alpha)=\exp[-c_{1\text{x}}\alpha/(1-\alpha)]\theta(1-\alpha)-d_{\text{x}}\exp[\text{c}_{2\text{x}}/(1-\alpha)]\theta(\alpha-1),$
(3.74)

where $\theta$ is the Heaviside step function. The function
$g_{\text{x}}(s)=1-\exp(-a_{1}s^{-1/2})$ is chosen to make
$F_{\text{x}}^{\text{SCAN}}(s,\alpha)$ vanish like $s^{-1/2}$ as $s\to\infty$,
which guarantees the non-uniform scaling finiteness conditions [Eqs. (2.66)
and (2.68)] [55, 70], and $a_{1}=4.9479$ is taken to recover the exact
exchange energy of the hydrogen atom. For $\alpha\approx 1$ (slowly varying
density regions), $F_{\text{x}}^{\text{SCAN}}(s,\alpha)\approx
h_{\text{x}}^{1}(s,\alpha)g_{\text{x}}(s)$ where $h_{\text{x}}^{1}(s,\alpha)$
is a PBE-like resummation of the fourth-order gradient expansion [Eq. (3.18)]

$\displaystyle
h_{\text{x}}^{1}(s,\alpha)=1+k_{1}-\frac{k_{1}}{1+x^{\text{SCAN}}(s,\alpha)/k_{1}},$
(3.75)

where

$\displaystyle x^{\text{SCAN}}(s,\alpha)=\mu
s^{2}[1+(b_{4}s^{2}/\mu)e^{-|b_{4}|s^{2}/\mu}]+[b_{1}s^{2}+b_{2}(1-\alpha)e^{-b_{3}(1-\alpha)^{2}}]^{2},$
(3.76)

with $\mu=10/81$, $b_{2}=(5913/405000)^{1/2}$, $b_{1}=(511/13500)/(2b_{2})$,
$b_{3}=0.5$, and $b_{4}=\mu^{2}/k_{1}-1606/18225-b_{1}^{2}$. For $\alpha=0$
(one-orbital regions),
$F_{\text{x}}^{\text{SCAN}}(s,\alpha=0)=h_{\text{x}}^{0}g_{\text{x}}(s)$ where
$h_{\text{x}}^{0}=1.174$ is chosen to saturate the local two-electron tight
bound $F_{\text{x}}^{\text{SCAN}}(s,\alpha=0)\leq 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 [70].

The SCAN correlation energy density is written as

$\displaystyle
e_{\text{c}}^{\text{SCAN}}(\rho_{\uparrow},\rho_{\downarrow},\nabla\rho_{\uparrow},\nabla\rho_{\downarrow},\tau_{\uparrow},\tau_{\downarrow})=\rho\;[\varepsilon_{\text{c}}^{1}(\rho_{\uparrow},\rho_{\downarrow},t)+f_{\text{c}}(\alpha)(\varepsilon_{\text{c}}^{0}(\rho_{\uparrow},\rho_{\downarrow},s)-\varepsilon_{\text{c}}^{1}(\rho_{\uparrow},\rho_{\downarrow},t))],$
(3.77)

which is again an interpolation between $\alpha=0$ and $\alpha=1$, and an
extrapolation to $\alpha\to\infty$ using the function

$\displaystyle
f_{\text{c}}(\alpha)=\exp[-c_{1\text{c}}\alpha/(1-\alpha)]\theta(1-\alpha)-d_{\text{c}}\exp[\text{c}_{2\text{c}}/(1-\alpha)]\theta(\alpha-1).$
(3.78)

For $\alpha=1$, the correlation energy par particle is taken as a revised
version of the PBE correlation energy per particle

$\displaystyle\varepsilon_{\text{c}}^{1}(\rho_{\uparrow},\rho_{\downarrow},t)=\varepsilon_{\text{c}}^{\text{UEG}}(\rho_{\uparrow},\rho_{\downarrow})+H_{1}^{\text{SCAN}}(\rho_{\uparrow},\rho_{\downarrow},t)$
(3.79)

where

$\displaystyle
H_{1}^{\text{SCAN}}(\rho_{\uparrow},\rho_{\downarrow},t)=A(0)\phi_{2}(\zeta)^{3}\ln\left[1+w_{1}(1-g({\cal
A}t^{2}))\right],$ (3.80)

with $t=|\nabla\rho|/(2\phi_{2}(\zeta)k_{\text{s}}\rho)$,
$w_{1}=\exp[-\varepsilon_{\text{c}}^{\text{UEG}}(\rho_{\uparrow},\rho_{\downarrow})/(A(0)\phi_{2}(\zeta)^{3})]-1$,
${\cal A}=\beta(r_{\text{s}})/(A(0)w_{1})$, and $g({\cal A}t^{2})=1/(1+4{\cal
A}t^{2})^{1/4}$. The function has a second-order gradient expansion
$H_{1}^{\text{SCAN}}(\rho_{\uparrow},\rho_{\downarrow},t)=\beta(r_{\text{s}})\phi_{2}(\zeta)t^{2}+O(t^{4})$
where the coefficient
$\beta(r_{\text{s}})=0.066725(1+0.1r_{\text{s}})/(1+0.1778r_{\text{s}})$ 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 $r_{\text{s}}\to\infty$ the second-order gradient
expansion terms for exchange and correlation cancel each other [128]. For
$\alpha=0$, the correlation energy par particle is constructed to be accurate
for one- and two-electron systems and is written as

$\displaystyle\varepsilon_{\text{c}}^{0}(\rho_{\uparrow},\rho_{\downarrow},s)=[\varepsilon_{\text{c}}^{\text{LDA0}}(\rho)+H_{0}^{\text{SCAN}}(\rho,s)]G_{\text{c}}(\zeta).$
(3.81)

The spin function
$G_{\text{c}}(\zeta)=[1-2.3631(\phi_{4}(\zeta)-1)](1-\zeta^{12})$ is designed
to make the correlation energy vanish for one-electron densities ($\alpha=0$
and $\zeta=\pm 1$) and to make the exchange-correlation energy independent of
$\zeta$ in the low-density limit [Eq. (2.63)]. Equation (3.81) includes a LDA-
type term [129]

$\displaystyle\varepsilon_{\text{c}}^{\text{LDA0}}(\rho)=-\frac{b_{1\text{c}}}{1+b_{2\text{c}}r_{\text{s}}^{1/2}+b_{3\text{c}}r_{\text{s}}},$
(3.82)

and a gradient correction

$\displaystyle
H_{0}^{\text{SCAN}}(\rho,s)=b_{1\text{c}}\ln\left[1+w_{0}(1-g_{\infty}(\zeta=0,s))\right],$
(3.83)

with $w_{0}=\exp(-\varepsilon_{\text{c}}^{\text{LDA0}}(\rho)/b_{1\text{c}})-1$
and $g_{\infty}(\zeta=0,s)=\lim_{\zeta\to
0}\lim_{r_{\text{s}}\to\infty}g({\cal A}t)=1/(1+0.512104s^{2})^{1/4}$. The
parameter $b_{1\text{c}}=0.0285764$ is determined so that the high-density