1. Introduction

One of the fundamental problems in condensed-matter physics and quantum chemistry is the theoretical study of electronic properties of systems ranging from atoms and molecules to complex materials. Since electrons are governed by the laws of quantum mechanics, all these systems are fully described by the Schrödinger equation. Analytic solutions of the Schrödinger equation are obtainable for very simple systems only. However, for systems with large number of atoms, the electrostatic repulsion between the electrons makes its numerical resolution very difficult. In that case, it is natural to consider the many-electron wave function $\psi (\mathbf{r} _{1}, \mathbf{r}_{2}, ..., \mathbf{r}_{i}, ...)$ with appropriate approximations. The earliest and widely used approximation was that of Hartree ^[1], which expresses the wave function of system as a product of one-electron wave functions, so that the problem reduces to a one-electron Schrödinger equation. Then, considerable improvement of the energy computation was made by incorporating the exchange effects with the so-called Hartree-Fock approximation ^[2], which replaces the product of one-electron wave functions by a linear combination of orbitals ^[3].

When the adiabatic approximation is made that permits to separate the movements of ions and electron, the electronic Schrödinger equation, describing the state of the electrons interacting in the presence of ions, reads:

$\left[\sum\limits_{i}\left(-\frac{\hbar ^{2}}{2m}\nabla _{i}^{2}\right) +\frac{1}{% 2}\sum\limits_{i}\sum\limits_{j\neq i}U_{ee}\left( \left\vert \mathbf{r}_{i}-\mathbf{r}% _{j}\right\vert \right) +\sum\limits_{i}\sum\limits_{\alpha }U_{ei}\left( \left\vert \mathbf{R}_{\alpha }-\mathbf{r}_{i}\right\vert \right) \right] \psi _{el} = E_{el}\psi _{el}\, ,$

(1)

or under the compact form:

$\left[T+U_{ee}+U_{ei}\right] \psi _{el} = E_{el}\psi _{el}\, ,$

(2)

where the two first terms in the Hamiltonian are the kinetic and electron-electron interaction operators, respectively, and the third term represents the potential energy operator of the electron $i$ in the field of all ions.

Within the Hartree approximation, the $N$-electron wave function $\psi _{el}$ is the product of one-electron wave functions:

$\psi _{el}(\mathbf{r}_{1}, \mathbf{r}_{2}, ..., \mathbf{r}_{N}) = \prod% \limits_{i}\psi _{i}(\mathbf{r}_{i}) = \psi _{1}(\mathbf{r}_{1})\psi _{2}(% \mathbf{r}_{2})...\psi _{N}(\mathbf{r}_{N})\, .$

(3)

Consequently, the electronic energy $E_{el},$ as a sum of one-electron energies, is given by:

$\begin{eqnarray} E_{el} & = &\sum\limits_{i}E_{i} = \left\langle \psi _{el}\left\vert H\right\vert \psi _{el}\right\rangle \notag \\ && \notag \\ & = &\int \!\cdots \!\int \psi _{1}^{\ast }(\mathbf{r}_{1})...\psi _{N}^{\ast }(\mathbf{r}_{N})\left[T+U_{ei}\right] \psi _{1}(\mathbf{r}_{1})...\psi _{N}(\mathbf{r}_{N})d\mathbf{r}_{1}...d\mathbf{r}_{N} \notag \\ &&+\frac{1}{2}\sum\limits_{i}\sum\limits_{j\neq i}\int \!\cdots \!\int \psi _{i}^{\ast }(% \mathbf{r}_{i})\psi _{j}^{\ast }(\mathbf{r}_{j})\left[\frac{e^{2}}{4\pi \varepsilon _{0}r_{ij}}\right] \psi _{i}(\mathbf{r}_{i})\psi _{j}(\mathbf{r}% _{j})d\mathbf{r}_{i}d\mathbf{r}_{j}\, . \label{EnHart} \end{eqnarray}$

(4)

The last term in the previous equation is the so-called Coulomb energy. The simplification achieved in this term comes from the fact that certain factors are equal to unity as a consequence of normalization of the wave function. Now, the implementation of the Hartree approximation amounts to define the Hartree potential $W^{H}\left(r_{i}\right)$ representing the potential energy of the electron $i$ in the field of all other electrons:

$W^{H}\left( r_{i}\right) = \frac{1}{2}\sum\limits_{j\neq i}\int \frac{e^{2}}{4\pi \varepsilon _{0}}\frac{\left\vert \psi _{j}(r_{j})\right\vert ^{2}}{r_{ij}}d% \mathbf{r}_{j}\, ,$

(5)

in order to reduce the Schrödinger equation for multiple electrons (Eq. (1)) to the sum of one-electron Schrödinger equations:

$\left[-\frac{\hbar ^{2}}{2m}\nabla _{i}^{2}+\sum\limits_{\alpha }U_{ei}\left( \left\vert \mathbf{R}_{\alpha }-\mathbf{r}_{i}\right\vert \right) +W^{H}\left( r_{i}\right) \right] \psi _{i} = E_{i}\psi _{i}\, ,$

(6)

whereby the wave function of the system $\psi _{el}$ is equal to the product of one-electron wave functions ($\psi _{el} = \prod\limits_{i}\psi _{i}$) and the electronic energy of the system is equal to the sum of energies of each electron ($E_{el} = \sum_{i}E_{i}$).

Within the Hartree-Fock approximation in place of the Hartree approximation, the $N$-electron wave function $\psi _{el}$ is constructed by taking a linear combination of products of one-electron wave functions written in the determinantal form ^[4]:

$\psi _{el}(\mathbf{r}_{1}, \mathbf{r}_{2}, ..., \mathbf{r}_{N}) = \frac{1}{\sqrt{% N!}}\sum\limits_{\mathbf{P}}(-1)^{p}\mathbf{P}\psi _{1}(\mathbf{r}_{1})\psi _{2}(% \mathbf{r}_{2})...\psi _{N}(\mathbf{r}_{N})\, ,$

(7)

where $\mathbf{P}$ is the permutation operator. As a result, an extra term known as the exchange energy }must be added onto equation (4), which comes from the exchange of electronic coordinates in the expanded determinantal wave function. The exchange energy reads:

$-\frac{1}{2}\sum\limits_{i}\sum\limits_{j\neq i}\int \int \psi _{i}^{\ast }(\mathbf{r}% _{j})\psi _{j}^{\ast }(\mathbf{r}_{i})\left[\frac{e^{2}}{4\pi \varepsilon _{0}r_{ij}}\right] \psi _{i}(\mathbf{r}_{j})\psi _{j}(\mathbf{r}_{i})d% \mathbf{r}_{i}d\mathbf{r}_{j}\, ,$

(8)

where the double sum is over the functions of same spin. Similarly to Eq. (5), its operator in the one-electron Schrödinger equation for electron $i$ in the field of all other electrons is:

$-\sum\limits_{j\neq i}\int \frac{e^{2}}{4\pi \varepsilon _{0}}\frac{\psi _{j}^{\ast }(\mathbf{r}_{j})\psi _{i}(\mathbf{r}_{j})\psi _{j}(\mathbf{r}_{i})}{% r_{ij}\psi _{i}(\mathbf{r}_{i})}d\mathbf{r}_{j}\, .$

(9)

Beside these considerations, a different description of many-electron systems was made by Thomas ^[5] and Fermi ^[6], who noted that the energy of electrons moving in the potential of nuclei can be directly calculated from the electron density $n(\mathbf{r})$ rather than from the wave function. If $N$ is the number of electrons in the system, the electron density is defined as:

$n(\mathbf{r}) = N\int \!\cdots \!\int \psi _{el}^{\ast }(\mathbf{r}, \mathbf{r}% _{2}, ..., \mathbf{r}_{N})\psi _{el}(\mathbf{r}, \mathbf{r}_{2}, ..., \mathbf{r}% _{N})d\mathbf{r}_{2}...d\mathbf{r}_{N}\, ,$

(10)

with the normalization condition $\int n(\mathbf{r})d\mathbf{r} = N$ and the limit $n(\mathbf{r\rightarrow \infty })\rightarrow 0.$ The term $n(\mathbf{r}%)d\mathbf{r}$ represents the elemental probability of finding any electron in the volume element $d\mathbf{r}$ at position $\mathbf{r.}$ As a result, the electronic energy $E_{el}$ may be calculated in terms of the electron density as:

$E_{el} = \left\langle \psi _{el}\left\vert H\right\vert \psi _{el}\right\rangle = \int \left[T+U_{ee}+U_{ei}\right] n(\mathbf{r})d\mathbf{% r}\, .$

(11)

Initially the authors ^[5,6] assumed that the motions of electrons were uncorrelated and that the corresponding kinetic energy could be described by a local approximation electron density. Shortly later, Dirac ^[7] proposed that the exchange effects be included by incorporating a term obtained from exchange energy density in the homogeneous system. It was demonstrated that the calculations based on the homogeneous system could be satisfactory even if the electron density is far from homogeneity. The Thomas-Fermi approach is a variational method providing an alternative to the Schrödinger equation for calculating the ground state electronic energy associated to the electron density $n_{0}(\mathbf{r}).$ For systems with many electrons, the advantage of the Thomas-Fermi (TF) method over the Hartree-Fock method is that, instead of solving the many-electron wave function $\psi (\mathbf{r}_{1}, \mathbf{r}_{2}, ..., \mathbf{r}_{N})$ in the Schrödinger equation (with $3N$ space coordinates and $N$ spin coordinates), the TF method deals with the electron density $n(\mathbf{r})$ depending on $3$ coordinates, whatever the number of electrons in the system, plus one spin coordinate possibly.

Even if it was rigorously demonstrated by Hohenberg and Kohn ^[8] and Kohn and Sham ^[9] that the ground state energy $E_{el}$ (and related properties) of a system can be linked to the electron density $n(% \mathbf{r})$, the Thomas-Fermi and Hartree-Fock methods can be regarded as forerunners of the density-functional theory (DFT). Hohenberg and Kohn proved that the ground-state energy is uniquely defined by a functional^* $E[n(\mathbf{r})]$ of the electron density profile, and Kohn and Sham postulated that many-electron system can be equivalently represented by a non-interacting reference system where each electron moves independently in the effective potential $v_{KS}(r)$ consisting of the Coulomb repulsion between the electrons and the potential arising from the exchange and correlation effects, plus the potential of nuclei whose the pseudopotential method is a representative application. The main problem in the initial density-functional theory was to determine the functionals of the kinetic energy and the exchange energy accounting for the electron-electron interaction. In the modern version of the DFT, self-consistent equations are solved for a set of orbitals whose electron density $n(\mathbf{r})$ is exactly that of the real system. In practical calculations, the exchange and correlation (xc) contributions are approximated. The simplest xc approximation is the LDA (local density approximation) extensively used in 1970s. The GGA (generalized gradient approximations) was developed later and, in the early 1990s, hybrid functionals were introduced by Becke ^[10], replacing a fraction of the conventional GGA approximation with HF (Hartree-Fock) exchange. Among the hybrid functionals, the most popular approximation in use in chemistry today, with about 80% of the occurrences in the literature, is the B3LYP (exchange of Becke ^[11] and correlation of Lee-Yang-Parr ^[12]).

^*While a function is a relation that assigns a number $f(x)$ to a number $x$, a functional is a relation associating a number $F[f(\mathbf{r})]$ with a function $f(\mathbf{r})$.

This article provides a tutorial overview of the DFT. It starts with a presentation of the simple atomic model of Thomas-Fermi—even if it has a limited success in describing the properties of real systems—because it is a good introduction to the DFT formalism. Then, it sketches the outline of fundamentals of the density-functional theory based on the works of Hohenberg, Kohn and Sham. As an application, the Thomas-Fermi approximation is rederived properly from the DFT and possible improvements are mentioned. A discussion of the exchange and correlation effects in the uniform electron gas is included giving rise to the local density approximation (LDA), and the generalized-gradient approximations (GGA), to treat the inhomogeneous electron gas, is also briefly presented. Lastly, mention is made of the main hybrid functionals and the software packages successfully applied to diverse materials of chemical, physical and biological interest.

---

2. Atomic Model of Thomas and Fermi

The Thomas-Fermi (TF) method has been developed independently by Thomas and Fermi according to a semi-classical approach, which uses conventional functions instead of quantum operators and borrows at the same time certain ideas from quantum mechanics like the electron density. The TF method, which circumvents the concept of wave function, constitutes a good introduction to the density-functional theory.

Contrary to the Bohr model, in the TF atomic model the electrons around the nucleus are assumed to constitute a free electron gas at $0$ K. It is well known that the average energy per electron in the free electron gas is very large compared to $k_{B}T,$ so that the quantum states associated with each energy level are so dense that it is possible to define the density of states:

$g(E) = \frac{V}{4\pi ^{2}}\left( \frac{2m}{\hbar ^{2}}\right) ^{\!\!3/2}\!\!E^{1/2}\, .$

(12)

As a consequence, the TF method is better suited for atoms with large atomic number $Z$ (the number of electrons in the system is $N = Z)$. Taken into account that each quantum state can be occupied by two electrons with opposite spins, the number of electrons in neutral atoms below the upper limit of energy $\mu$ is:

$Z = 2\int_{0}^{\mu }g(E)dE = 2\frac{V}{4\pi ^{2}}\left( \frac{2m}{\hbar ^{2}}% \right) ^{\!\!3/2}\frac{2}{3}\mu ^{3/2}\, .$

(13)

This relation allows to write the chemical potential $\mu$ in terms of the electron density $n = \frac{Z}{V}$:

$\mu = \left( 3\pi ^{2}n\right) ^{2/3}\left( \frac{\hbar ^{2}}{2m}\right) = % \frac{\hbar ^{2}k_{F}^{2}}{2m} = \frac{p_{F}^{2}}{2m}\, ,$

(14)

where $p_{F}( = \hbar k_{F})$ is the momentum and $k_{F}\left[ = \left(3\pi ^{2}n\right) ^{1/3}\right]$ the wavevector at the Fermi level.

In the semi-classical atomic TF model, it is common to represent the kinetic energy, per electron, by $\mu$ and to take in consideration its dependence on the position $r$ in the atom through the electron density $% n(r)$, according to the relation $\mu (r) = \left(3\pi ^{2}n(r)\right)^{2/3}\left(\frac{\hbar ^{2}}{2m}\right)$. Now we consider the attractive potential energy $v(r) = -e\varphi (r)$ of an electron at the position $r$ from the nucleus, resulting from the electrostatic potential $% \varphi (r)$. Thus, the total energy, $E =$ $\mu (r)+v(r),$ of an electron at the position $r$, is:

$E = \frac{\hbar^{2}}{2m}\left[3\pi^{2}n(r)\right]^{2/3}-e\varphi (r)\, .$

(15)

It should be stressed that the chemical potential must remain constant over the atom to avoid a concentration of the electrons in a particular region, contrary to what happens in the Bohr model. In fact, $E$ is a constant negative or equal to zero since the electrons are not permitted to escape from the atom. For a neutral atom, it is stated that $E = 0$ to guarantee that $n(r) = 0$ and $\varphi (r) = 0$ when $r\rightarrow \infty$. Eq. (15) is the central equation of the Thomas-Fermi model, which links the electron density $n(r)$ to the potential energy $v(r),$ i.e., to the electrostatic potential $\varphi (r),$ in place of the Schrödinger equation, namely:

$n(r) = \frac{1}{\left(3\pi ^{2}\right)}\! \left(\frac{2me}{\hbar^{2}}% \right)^{\!\!3/2} \!\!\varphi^{3/2}(r)\, .$

(16)

This equation is the integral form of the TF one; it is possible to transform it into a differential form by using the Poisson equation, which links the electrostatic potential $\varphi (\mathbf{r})$ to the density of charge $\rho (\mathbf{r}) = -en(\mathbf{r})$ as^*:

$\nabla ^{2}\varphi (\mathbf{r}) = -\frac{\rho (\mathbf{r})}{\varepsilon _{0}} = % \frac{en(\mathbf{r})}{\varepsilon _{0}}.$

(17)

^*Here, the Poisson equation is written in the rationalized M.K.S. unit system (SI units) as $\nabla ^{2}\varphi (\mathbf{r}) = -\frac{\rho (\mathbf{r})}{ \varepsilon _{0}}$, rather than in the unrationalized Gaussian unit system where $\nabla ^{2}\varphi (\mathbf{r}) = -4\pi \rho (\mathbf{r}).$

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

$\nabla ^{2}\varphi (\mathbf{r}) = \frac{en(\mathbf{r})}{\varepsilon _{0}} = % \frac{e}{3\pi ^{2}\varepsilon _{0}}\left( \frac{2me}{\hbar ^{2}}\right) ^{\!\!3/2}\varphi ^{3/2}(\mathbf{r}).$

(18)

By using the radial part of the Laplacian in spherical coordinates, and considering the change of variable $u(r) = \frac{r}{Ze}\varphi (r)$, the TF equation reads:

$\begin{eqnarray} \frac{1}{r^{2}}\frac{d}{dr}\left( r^{2}\frac{d\varphi (r)}{dr}\right) & = &% \frac{Ze}{r}\frac{d^{2}u(r)}{dr^{2}} = \frac{e}{3\pi ^{2}\varepsilon _{0}}% \left( \frac{2me}{\hbar ^{2}}\right) ^{\!\!3/2}\left[\frac{Zeu(r)}{r}\right] ^{3/2}, \notag \\ \frac{d^{2}u(r)}{dr^{2}} & = &\frac{e\left( Ze\right) ^{1/2}}{3\pi ^{2}\varepsilon _{0}}\left( \frac{2me}{\hbar ^{2}}\right) ^{\!\!3/2}\left[ \frac{u^{3}(r)}{r}\right] ^{1/2}. \label{9} \end{eqnarray}$

(19)

To solve this differential equation it is convenient to make the following change of variable $x = \alpha r,$ with the parameter $\alpha = \frac{Z^{1/3}} { \left(3\pi ^{2}\varepsilon _{0}\right)^{2/3}}\frac{2me^{2}}{\hbar ^{2}}$ whose unit is C$^{2/3}\, {\cdot}\,$s$^{2/3}\, {\cdot}\,$kg$^{-1/3}\, {\cdot}\,$m $^{-2}$. With this change of variable, the TF differential equation reduces to:

$\frac{d^{2}u(x)}{dx^{2}} = \left[\frac{u^{3}(x)}{x}\right]^{1/2}.$

(20)

The solution of this equation is obtained by numerical integration with the boundary conditions of $u(x),$ i.e., $\underset{x\rightarrow 0}{\text{lim}}% u(x) = 1$ since $\varphi (r)\rightarrow \frac{Ze}{r}$ near the nucleus, and $% \underset{x\rightarrow \infty }{\text{lim}}u(x)\rightarrow 0$. The solution $% u(x)$ is a monotonous decreasing function of $x.$ It has been shown that the curve starts from the point $u(x = 0) = 1$ with the slope $a = -1.588$ at the origin and approaches the $x$-axis asymptotically. For the numerical integration it is advantageous to introduce the change of variable $% x = w^{2}/2$. This makes the interval for each step of the numerical integration small near the origin where $u(x)$ changes appreciably, and larger where $u(x)$ has a slow variation. The detailed procedure of integration has been described by Feynman et al. ^[13]. Besides, it is convenient to use the series $u(x) = 1-1.588x+4/3x^{3/2}$ valid for small values of $x$ to initiate the numerical integration.

It is interesting to see the effect of the atomic number $Z$ on the electron density $n(r)$. To do this, consider the changes of variables $x = \alpha r$ and $u(r) = \frac{r}{Ze}\varphi (r)$. At small values of $r$, it follows that $% \varphi (r)\, {\simeq }\, \frac{Ze}{r} = \alpha \frac{Ze}{x} = \alpha \varphi (x)$ and $u(r) = \frac{x}{Ze}\varphi (x) = u(x)$, hence $u(r)$ is independent of $Z$ and $\varphi (x)$ is directly proportional to $Z$. In addition, since $% \alpha$ is proportional to $Z^{1/3},$ it results that $\varphi (r)\, \ {% \propto }\, \ Z^{4/3}$ and $n(r)\, \ {\propto }\, \ \varphi ^{3/2}(r)\, \ {% \propto }\, \ Z^{2}$.

Incidentally, note that the substitution of the electron density (Eq. (16) into the relation of electron density normalization, $% \int_{0}^{\infty }4\pi r^{2}n(r)dr = Z$, allows us to find the following normalization:

$\int_{0}^{\infty }x^{1/2}u^{3/2}(x)dx = 4\pi \varepsilon _{0},$

(21)

which becomes $\int_{0}^{\infty }x^{1/2}u^{3/2}(x)dx = 1$ in the unrationazed Gaussian unit system.

The TF method has been found to give a rough description of the electron density and the electrostatic potential in the atom, which is better when the number of electrons is large. Since the TF model does not agree well with the experiment, some authors tried to improve the model. For instance, Dirac ^[7] introduced an exchange term that modifies Eq. (20) as follows:

$\frac{d^{2}u(x)}{dx^{2}} = x\left[\epsilon +\frac{u^{1/2}(x)}{x^{1/2}}\right] ^{3},$

(22)

where $x$ has the same meaning as before and $\epsilon$ is defined as

$\epsilon = \left( \frac{3}{32\pi ^{2}}\right) ^{\!\!1/3}Z^{-2/3}.$

(23)

Schwinger ^[14] provided some justification for the addition of two corrections to the ground-state energy $E$ of a neutral atom composed of $Z$ electrons. One correction, $E_{ex}$, lies in the exchange effects, and the other, $E_{qu}$, comes from the bulk electrons (the electrons near neither the nucleus nor the surface of the atom) that experience a small varying potential. Both corrections are proportional to $Z^{5/3}$, i.e., the total correction is $E_{ex}+E_{qu} = \frac{11}{9}E_{ex} = -0.2699\, Z^{5/3}$ Hartree. In addition, it has been argued ^[15] that $E$ or its derivative is a discontinuous function of $Z$, and that $E$ has a small correction ^[16,17] varying sinusoidally with $Z$, which have something to do with the filling of atomic shells.

---

3. Summary of the Density-Functional Theory

The rigorous developments of the density-functional theory were posed by Hohenberg, Kohn and Sham ^[8,9] that legitimized the model intuitively established by Thomas, Fermi and Dirac. For practical reasons, they replace the term $U_{ei}$ with $V_{ext}$ to indicate that the electrons move in an arbitrary external potential including the Coulomb energy due to point nuclei, so that the Schrödinger equation (Eq. (2)) for $N$ electrons reads:

$\left[T+U_{ee}+V_{ext}\right] \psi _{el} = E_{el}\psi _{el}\, ,$

(24)

and the electronic energy (Eq. (11)) is a functional of the electron density profile $n(\mathbf{r})$:

$E_{el}\left[n\right] = \left\langle \psi _{el}\left\vert H\right\vert \psi _{el}\right\rangle = F\left[n\right] +\int V_{ext}(\mathbf{r})n(\mathbf{r})d% \mathbf{r}\, ,$

(25)

where

$F\left[n\right] = \left\langle \psi _{el}\left\vert T\right\vert \psi _{el}\right\rangle +\left\langle \psi _{el}\left\vert U_{ee}\right\vert \psi _{el}\right\rangle = T\left[n\right] +U_{ee}\left[n\right]$

(26)

is the so-called universal energy functional, in the sense that it does not contain the external potential $V_{ext}(\mathbf{r})$ and can be determined independently of $\int V_{ext}(\mathbf{r})n(\mathbf{r})d\mathbf{r.% }$

---

3.1. Formulation of Hohenberg and Kohn

The Hohenberg and Kohn (HK) approach may be understood in considering the ground state ^[18] of a system of electrons at $0$ K with the wave function $\psi ^{(1)}$, the Hamiltonian $H^{(1)}$ and the external potential $V_{ext}^{(1)}(\mathbf{r})$. The HK approach asserts that the external potential can be uniquely determined by an electron density, $n(% \mathbf{r})$, that minimizes the ground state energy:

$E^{(1)} = \left\langle \psi ^{(1)}\left\vert H^{(1)}\right\vert \psi ^{(1)}\right\rangle = T[n]+U_{ee}[n]+\int V_{ext}^{(1)}(\mathbf{r})n(\mathbf{r% })d\mathbf{r}\, .$

(27)

Suppose that another Hamiltonian $H^{(2)}$, only different from $H^{(1)}$ in the external potential $V_{ext}^{(2)}(\mathbf{r}){\neq }V_{ext}^{(1)}(% \mathbf{r})$, gives rise to the wave function $\psi ^{(2)}$ that yields the same density $n(\mathbf{r})$. Unless $V_{ext}^{(1)}(\mathbf{r}%)-V_{ext}^{(2)}(\mathbf{r}) =$ const, $\psi ^{(2)}$ can not be equal to $\psi ^{(1)}$ since they satisfy different Schr{ö}dinger equations. The wave function $\psi ^{(2)}$ minimizes the ground-state energy for the Hamiltonian $H^{(2)}$:

$E^{(2)} = \left\langle \psi ^{(2)}\left\vert H^{(2)}\right\vert \psi ^{(2)}\right\rangle \, .$

(28)

By force of the variational principle, the expectation energy of the Hamiltonian is the lowest with its "correct" wave function, hence

$E^{(1)} = \left\langle \psi ^{(1)}\left\vert H^{(1)}\right\vert \psi ^{(1)}\right\rangle \text{ }<\quad \left\langle \psi ^{(2)}\left\vert H^{(1)}\right\vert \psi ^{(2)}\right\rangle .$

(29)

However, the last term of this inequality reads:

$\left\langle \psi ^{(2)}\left\vert H^{(1)}\right\vert \psi ^{(2)}\right\rangle = \left\langle \psi ^{(2)}\left\vert H^{(2)}\right\vert \psi ^{(2)}\right\rangle +\left\langle \psi ^{(2)}\left\vert H^{(1)}-H^{(2)}\right\vert \psi ^{(2)}\right\rangle \, .$

(30)

The operators $T$ and $U_{ee}$ being the same in the two Hamiltonians, it follows that $\left\langle \psi ^{(2)}\left\vert H^{(1)}-H^{(2)}\right\vert \psi ^{(2)}\right\rangle = \int [V_{ext}^{(1)}(\mathbf{r})-V_{ext}^{(2)}(% \mathbf{r})]\, n(\mathbf{r})d\mathbf{r}$, hence Eq. (29) becomes:

$E^{(1)}<E^{(2)}+\int \left[V_{ext}^{(1)}(\mathbf{r})-V_{ext}^{(2)}(\mathbf{r% })\right] n(\mathbf{r})d\mathbf{r}\, .$

(31)

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

$E^{(2)}<E^{(1)}+\int \left[V_{ext}^{(2)}(\mathbf{r})-V_{ext}^{(1)}(\mathbf{r% })\right] n(\mathbf{r})d\mathbf{r}\, ,$

(32)

and the sum of the two previous relations leads to the inconsistency:

$E^{(1)}+E^{(2)}\quad \mathbf{<}\quad E^{(1)}+E^{(2)}\, ,$

(33)

indicating that the starting hypothesis is wrong. Thus the HK approach, based upon the variational principle and reductio ad absurdum, implies that it is sufficient to use the electron density (and not the full wave function) as the variational property, from which, once it is found, all ground-state properties can be deduced. The considerations above are restricted to nondegenerate ground states; however, a generalization towards degenerate states has been suggested ^[19] within a different and more systematic scheme.

While the universal energy functional $F\left[n\right]$ (composed of kinetic and potential energies of the electron gas) remains unchanged whatever the external potential $V_{ext}(\mathbf{r})$ (see Eq. (26)), two different external potentials can not yield the same ground state energy. The HK approach asserts that there is a one-to-one mapping between external potential $V_{ext}(\mathbf{r})$, electron density $% n(\mathbf{r})$ and wave function $\psi$ of the system in its ground state. It implies that, given the electron density, only a unique external potential and, consequently, its unique wave function can be determined. Conversely, given the external potential, the corresponding electron density and wave function are uniquely defined. It should also be stressed that the term $\int \!V_{ext}(\mathbf{r})\, n(\mathbf{r})d\mathbf{r}$ is relatively easy to calculate for a given external potential, whereas the expressions of the functionals $T\left[n\right]$ and $U_{ee}\left[n\right]$ are more difficult to evaluate and some approximations have to be made. In the initial density-functional theory, the main difficulty was to determine the functionals of the kinetic and potential energies for the interacting electron gas.

---

3.2. Formulation of Kohn and Sham

The Kohn and Sham (KS) approach consists in prescribing a convenient recipe to calculate the universal energy functional $F\left[n\right] = T\left[n% \right] +U_{ee}\left[n\right]$ for interacting electronic systems. The authors assumed that it was reasonable to evaluate the kinetic energy $T[n]$ of any interacting electronic system by $T_{0}[n]$, which is the kinetic energy of the corresponding non homogeneous and non interacting electron gas of the same density, plus the exchange-correlation energy $% E_{xc}[n]$, to be further specified (see Eq. (36)). The potential energy $U_{ee}[n]$ which, within the Hartree approximation, is that of electrons in the field of all other electrons (Eq. (5)), is approximated by the electron density interacting with itself, including self-interaction:

$U_{ee}\left[n\right] \sim W^{H}\left[n\right] = \frac{1}{2}\frac{e^{2}}{% 4\pi \varepsilon _{0}}\int \!\!\int \frac{n(\mathbf{r})n(\mathbf{r}^{\prime })}{\left\vert \mathbf{r}-\mathbf{r}^{\prime }\right\vert }d\mathbf{r}d% \mathbf{r}^{\prime }.$

(34)

With the KS prescription, the total energy functional (Eq. (25)) reads:

$E_{el}\left[n\right] = T_{0}\left[n\right] +E_{xc}\left[n\right] +W^{H}% \left[n\right] +\int V_{ext}(\mathbf{r})n(\mathbf{r})d\mathbf{r}\, .$

(35)

In the KS approach, the kinetic energy $T_{0}\left[n\right]$ is not an explicit functional of density, but of orbitals of which this density is constructed. The exchange-correlation energy functional $E_{xc}% \left[n\right]$ incorporates the difference between real and approximated universal energy functionals:

$E_{xc}\left[n\right] = \left( T\left[n\right] +U_{ee}\left[n\right] \right) -\left( T_{0}\left[n\right] +W^{H}\left[n\right] \right) \, .$

(36)

It should be noted that $E_{xc}\left[n\right]$ is the only unknown functional in the right-hand side of Eq. (35), since the other terms can be calculated exactly. In practical terms, $E_{xc}\left[n\right]$ is a small part of the total energy because the magnitudes of $T\left[n% \right]$ and $T_{0}\left[n\right]$ are similar, and the same is true of $% U_{ee}\left[n\right]$ and $W^{H}\left[n\right]$. Note also that, if $% E_{xc}\left[n\right]$ vanishes, the physical content of the KS procedure becomes identical to that of the Hartree approximation given by Eq. (4). The simplest and reasonable approximation for $E_{xc}\left[n% \right]$ proposed by Kohn and Sham is the so-called local density approximation (LDA):

$E_{xc}^{LDA}\left[n(\mathbf{r})\right] = \int e_{xc}\left( n(\mathbf{r}% )\right) n(\mathbf{r})d\mathbf{r}\, ,$

(37)

where $e_{xc}\left(n(\mathbf{r})\right)$ is the exchange-correlation energy, per electron, of the uniform electron gas of local electron density $% n(\mathbf{r}).$ A fairly good result for $E_{xc}^{LDA}\left[n(\mathbf{r})% \right]$ cannot be overstated when $n(\mathbf{r})$ has a very low spatial variation.

Considering the characteristics of the functionals in Eq. (35), it is not surprising that the many-electron Schrödinger equation in the KS procedure reduces to a set of one-electron equations. Indeed, the electron density would minimize the total energy, i.e., $\delta E_{el}\left[n\right] = 0$, under the normalization condition $\int n(\mathbf{r})d\mathbf{r} = N$. Therefore, the method of Lagrange with the multiplier $\mu$ may be used to obtain the equation requiring that the variation of the functional $\left\{ E_{el}\left[n\right] -\mu N\right\}$ with respect to $n(\mathbf{r})$ is equal to zero, i.e., $\delta \left\{ E_{el}\left[n\right] -\mu N\right\} = 0$. This variation allows us to define the functional derivative, denoted by $\frac{\delta }{\delta n(\mathbf{r})}\left\{ E_{el}\left[n\right] -\mu N\right\}$, by the relation^*:

$\delta \left\{ E_{el}\left[n\right] -\mu N\right\} = \int \frac{\delta }{% \delta n(\mathbf{r})}\left\{ E_{el}\left[n\right] -\mu N\right\} \delta n(% \mathbf{r})d\mathbf{r},$

(38)

^*The functional derivative is the natural extension of the standard derivative. If $F\left[f(\mathbf{r})\right]$ is the functional associated with the function $f(\mathbf{r})$, its variation $\delta F\left[f(\mathbf{r})\right]$ generated by the variation $\delta f(\mathbf{r})$ defines the functional derivative, denoted by $\frac{\delta F\left[f(\mathbf{r})\right] }{\delta f(x)}$, with the general relation:

$\delta F\left[f(\mathbf{r})\right] = \int \frac{\delta F\left[f(\mathbf{r})% \right] }{\delta f(\mathbf{r})}\delta f(\mathbf{r})d\mathbf{r}.$

where, according to Eqs. (34, 35 and 37), the functional $\left\{ E_{el}\left[n\right] -\mu N\right\}$ reads:

$\left\{ T_{0}\left[n\right] +\frac{1}{2}\frac{e^{2}}{4\pi \varepsilon _{0}}% \int \int \frac{n(\mathbf{r})n(\mathbf{r}^{\prime })}{\left\vert \mathbf{r}-% \mathbf{r}^{\prime }\right\vert }d\mathbf{r}d\mathbf{r}^{\prime }+\int e_{xc}\left( n\right) n(\mathbf{r})dr+\int V_{ext}(\mathbf{r})n(\mathbf{r})d% \mathbf{r-}\mu \int n(\mathbf{r})d\mathbf{r}\right\} .$

(39)

Now, with Eq. (38), the variational principle $\delta \left\{ E_{el}% \left[n\right] -\mu N\right\} = 0$ yields:

$\begin{eqnarray} 0 & = &\int \frac{\delta }{\delta n}\left\{ T_{0}\left[n\right] \right\} \delta nd\mathbf{r}+\frac{1}{2}\frac{e^{2}}{4\pi \varepsilon _{0}}\int \int \frac{\delta }{\delta n}\left\{ \frac{n(\mathbf{r})n(\mathbf{r}^{\prime })}{% \left\vert \mathbf{r}-\mathbf{r}^{\prime }\right\vert }\right\} \delta nd% \mathbf{r}d\mathbf{r}^{\prime }+\int \frac{\delta }{\delta n}\left\{ e_{xc}\left( n\right) n(\mathbf{r})\right\} \delta ndr \notag \\ &&+\int \frac{\delta }{\delta n}\left\{ V_{ext}(\mathbf{r})n(\mathbf{r}% )\right\} \delta nd\mathbf{r-}\mu \int \frac{\delta }{\delta n}\left\{ n(% \mathbf{r})\right\} \delta nd\mathbf{r}, \label{18} \end{eqnarray}%$

(40)

$\int \left\{ \frac{\delta }{\delta n}\left\{ T_{0}\left[n\right] \right\} +% \frac{e^{2}}{4\pi \varepsilon _{0}}\int \frac{n(\mathbf{r}^{\prime })}{% \left\vert \mathbf{r}-\mathbf{r}^{\prime }\right\vert }d\mathbf{r}^{\prime }+e_{xc}\left[n\right] +n(\mathbf{r})\frac{\delta }{\delta n}\left\{ e_{xc}% \left[n\right] \right\} +V_{ext}(\mathbf{r})-\mu \right\} \delta n(\mathbf{r% })d\mathbf{r} = 0.$

Kohn and Sham postulated that the many-electron system described by the functional derivative in brackets of the previous equation can be equivalently represented by a non interacting reference system described by the following one-electron Schrödinger equation:

$\left( -\frac{\hbar ^{2}}{2m}\nabla _{i}^{2}+v_{KS}\right) \psi _{i} = E_{i}\psi _{i},$

(41)

where each electron moves independently in the effective potential $v_{KS}(% \mathbf{r})$ consisting of the Hartree potential $w^{H}(\mathbf{r})$, the exchange-correlation potential $v_{xc}(\mathbf{r})$ and the external potential $V_{ext}(\mathbf{r}),$ namely:

$v_{KS}(\mathbf{r}) = w^{H}(\mathbf{r})+v_{xc}(\mathbf{r})+V_{ext}(\mathbf{r}),$

(42)

by setting:

$\begin{eqnarray} -\frac{\hbar ^{2}}{2m}\nabla _{i}^{2} & = &\frac{\delta }{\delta n}\left\{ T_{0}\left[n(\mathbf{r})\right] \right\}, \notag \\ w^{H}(\mathbf{r}) & = &\frac{e^{2}}{4\pi \varepsilon _{0}}\int \frac{n(\mathbf{% r}^{\prime })}{\left\vert \mathbf{r}-\mathbf{r}^{\prime }\right\vert }d% \mathbf{r}^{\prime }, \label{wH(r)} \\ v_{xc}(\mathbf{r}) & = &e_{xc}\left[n(\mathbf{r})\right] +n(\mathbf{r})\frac{% \delta }{\delta n}\left\{ e_{xc}\left[n(\mathbf{r})\right] \right\}, \notag \\ E_{i} & = &\mu, \notag \end{eqnarray}$

(43)

$\text{and}\quad \quad n(\mathbf{r}) = \sum\limits_{i}\left\vert \psi _{i}\right\vert ^{2}. \label{eldens}$

(44)

Thus, the electron density in atoms, molecules and solids can be regarded as a sum of densities due to non-interacting quasiparticles moving in the effective potential $v_{KS}(\mathbf{r}).$ In the KS equation (Eq. (41)), the potential $v_{KS}(\mathbf{r})$ is a function of the electron density $n(\mathbf{r}),$ which can be calculated from the wave functions (Eq. (44)) of independent electrons (Kohn-Sham orbitals). Provided that $e_{xc}\left[n(\mathbf{r})\right]$ is known, the KS equation may be solved in a self-consistent procedure as follows: from an initial arbitrary electron density $n_{1}(\mathbf{r}),$ the effective potential $v_{KS}(\mathbf{r})$ is determined and Eq. (41) is solved. Then, a new electron density $n_{2}(\mathbf{r})$ is deduced with Eq. (44), and the procedure is repeated until convergence of the potential $v_{KS}(\mathbf{r})$ and the electron density $n(\mathbf{r}).$

It should be noted that the KS orbitals are calculated at each iteration, in terms of a set of basis functions, and that the coefficients of linear combination of basis functions are determined as in the Hartree-Fock calculations, but the computational time is shorter. The choice of the set of basis functions is of great importance in the KS calculations.

The KS approach gives useful results for most physical and chemical applications ^[20], but a large number of subsequent approximations must be used for correcting some of the defects of the local density approximation $E_{xc}^{LDA}\left[n(\mathbf{r})\right].$

---

4. The TF Equation from the KS Approach

In this section, we rederive the Thomas-Fermi model as an example, to illustrate the use of the KS approach. All the components of Eq. (35) are written as a function of the electron density for the interacting electron gas enclosed in a box and, at the end, the integral form of the generalized Thomas-Fermi equation is recovered.

We start by determining a functional form for the kinetic energy of a non interacting electron gas, which is exact for the homogeneous electron gas $\left(n = \text{const.}\right)$. In that case, the calculation of the kinetic energy of the electron gas, $T_{0},$ may be performed in $% \mathbf{k}$-space by doing the sum of the energy states contained in the Fermi sphere of radius $k_{F}$:

$\begin{eqnarray} T_{0} & = &\sum\limits_{k<k_{F}}E(k) = \frac{2V}{\left( 2\pi \right) ^{3}}% \int_{0}^{k_{F}}\frac{\hslash ^{2}k^{2}}{2m}d\mathbf{k}\, , \notag \\ & = &\frac{2V}{\left( 2\pi \right) ^{3}}4\pi \int_{0}^{k_{F}}\frac{\hslash ^{2}k^{4}}{2m}dk = \frac{V}{\pi ^{2}}\frac{\hslash ^{2}}{10m}k_{F}^{5}\, . \label{19} \end{eqnarray}$

(45)

The total number of electrons in the Fermi sphere is obtained similarly as:

$N = \frac{2V}{\left( 2\pi \right) ^{3}}\int_{0}^{k_{F}}4\pi k^{2}dk = \frac{V}{% \left( 3\pi ^{2}\right) }k_{F}^{3}\, ,$

(46)

so that the electron density is $n = \frac{N}{V} = \frac{k_{F}^{3}}{3\pi ^{2}}$. Then, the kinetic energy density $t$ (kinetic energy per unit volume) of the homogeneous electron gas may be expressed as a function of the electron density as:

$t = \frac{T_{0}}{V} = \frac{3\hslash ^{2}}{10m}\left( 3\pi ^{2}\right) ^{2/3}n^{5/3}\, .$

(47)

Incidentally, note that the kinetic energy density may be written as a function of the Seitz radius, $r_{s} = \left(\frac{3}{4\pi n}\right) ^{1/3},$ instead of the electron density, so that the kinetic energy per electron becomes:

$\begin{eqnarray} \frac{T_{0}}{N} & = &\frac{T_{0}}{V}\frac{V}{N} = \frac{t}{n}, \notag \\ \frac{T_{0}}{N} & = &\frac{3\hslash ^{2}}{10m}\left( 3\pi ^{2}\right) ^{2/3}n^{2/3} = \frac{1.105}{r_{s}^{2}}, \end{eqnarray}$

(48)

where $r_{s}$ is expressed in atomic unit of length and $\frac{T_{0}}{N}$ in atomic unit of energy (Hartree)^*.

^* Atomic units are used throughout this article. With $m = \hslash ^{2} = e^{2} = 4\pi \varepsilon _{0} = 1$, it is found that 1 a.u. of length $% = 0.529{\times }10^{-10}$ m and $1$ Hartree $= 2$ Rydberg $= 27.21$ eV $= 4.36{% \times }10^{-18}$ J.

For a non homogeneous electron gas where the electron density $n(% \mathbf{r})$ depends on position, the same functional form is assumed for the kinetic energy density, hence the kinetic energy functional is determined by the relation $T_{0}\left[n\right] = \int t\left[n(\mathbf{r}) % \right] d\mathbf{r, }$ which becomes:

$T_{0}\left[n\right] = C_{K}\int \left\{ n(\mathbf{r})\right\}^{5/3}d\mathbf{% r, }$

(49)

where $C_{K} = \frac{3\hslash ^{2}}{10m}\left(3\pi ^{2}\right) ^{2/3}.$

For an interacting electron gas, it is necessary to evaluate the other terms of Eq. (35). In the case of the kinetic energy functional, one must add the self-interaction of electrons $W^{H}\left[n% \right]$ given by Eq. (34) in the Hartree approximation and, if the electron gas is subject to the external potential $V_{ext}(\mathbf{r}),$ the functional $\int V_{ext}(\mathbf{r})n(\mathbf{r})d\mathbf{r}$ must be also taken into account.

As far as the unknown functional $E_{xc}\left[n\right]$ is concerned, it can be approximated by the LDA (Eq. (37)). It is generally divided into two separate terms,

$E_{xc}\left[n(\mathbf{r})\right] = E_{x}\left[n(\mathbf{r})\right] +E_{c}% \left[n(\mathbf{r})\right],$

(50)

where the exchange functional $E_{x}[n(\mathbf{r})]$ follows from the Pauli exclusion principle and affects the electrons of same spin, whereas the correlation functional $E_{c}[n(\mathbf{r})]$ manifests the mutual avoiding of electrons as charged particles and, consequently, does not depend on spin. While $E_{c}[n(\mathbf{r})]$ can be obtained through computer simulations, $E_{x}[n(\mathbf{r})]$ is exactly known as:

$E_{x}\left[n\right] = C_{x}\int \left\{ n(\mathbf{r})\right\} ^{4/3}d\mathbf{% r},$

(51)

where $C_{x} = -\frac{3}{4}\frac{e^{2}}{4\pi \varepsilon _{0}}\left(\frac{3}{% \pi }\right) ^{\!1/3}$ has been specified by Dirac ^[7]. This is the LDA for $E_{x}[n]$. A dimensional argument may be used to know the power of the electron density: the dimensional homogeneity of Eq. (35) says that each term has dimension of energy $\left[\text{M}\!\cdot \!\text{L} ^{2}\!\cdot \!\text{T}^{-2}\right]$, so that the energy density has dimension $\left[\text{M}\!\cdot \!\text{L}^{-1}\!\cdot \!\text{T}^{-2} \right]$, and the electron density $n(r)$ has dimension $\left[\text{L} ^{-3}\right]$. On the other hand, the term $E_{x}\left[n\right]$ arising from the self-interaction of electrons $W^{H}\left[n\right]$ (Eq. (35)) requires that the energy density scales $\frac{dE_{x}\left[n\right] }{d\mathbf{r}}\sim \frac{e^{2}nn}{r}r^{3}$, where $e^{2}$ has dimension $% \left[\text{M}\!\cdot \!\text{L}^{3}\!\cdot \!\text{T}^{-2}\right]$. Then, if the power law assumed for the energy density is $\frac{dE_{x}\left[n% \right] }{d\mathbf{r}}\propto e^{2}n(r)^{\nu }$, the dimensional analysis of this relation yields $\left(\text{M}\!\cdot \!\text{L}^{-1}\!\cdot \!\text{T }^{-2}\right) = \left(\text{M}\!\cdot \!\text{L}^{3}\!\cdot \!\text{T} ^{-2}\right) \!\cdot \!\text{L}^{-3\nu }$, leading to the value of the exponent $\nu = 4/3$. Thus, in the absence of correlation, the total energy functional for a non homogeneous and interacting electron gas reads explicitly:

$E_{el}\left[n\right] = C_{K}\int \left\{ n(\mathbf{r})\right\} ^{5/3}d% \mathbf{r}+\int w^{H}(r)n(\mathbf{r})d\mathbf{r}+C_{x}\int \left\{ n(\mathbf{% r})\right\} ^{4/3}d\mathbf{r}+\int V_{ext}(\mathbf{r})n(\mathbf{r})d\mathbf{% r.}$

(52)

In order to recover the generalized Thomas-Fermi equation, let's seek the distribution of the electron density that would minimize the total energy, i.e., $\delta E_{el}\left[n\right] = 0,$ under the normalization condition $% \int n(\mathbf{r})d\mathbf{r} = N$. Using the method of Lagrange with the multiplier $\mu,$ write the equation requiring that $\frac{\delta }{\delta n(\mathbf{r})}\left\{ E_{el}\left[n\right] -\mu N\right\} = 0$, as in Eq. (38):

$\int \left\{ \frac{5}{3}C_{K}\left\{ n(\mathbf{r})\right\} ^{2/3}+\frac{e^{2}% }{4\pi \varepsilon _{0}}\int \frac{n(\mathbf{r}^{\prime })}{\left\vert \mathbf{r}-\mathbf{r}^{\prime }\right\vert }d\mathbf{r}^{\prime }+\frac{4}{3}% C_{x}\left\{ n(\mathbf{r})\right\} ^{1/3}+V_{ext}(\mathbf{r})-\mu \right\} \delta n(\mathbf{r})d\mathbf{r} = 0.$

(53)

In brackets is the functional derivative that must be equal to zero to fulfill the minimization condition. Note that this functional derivative corresponds to the generalized TF equation:

$\frac{5}{3}C_{K}\left\{ n(\mathbf{r})\right\} ^{2/3}+\frac{e^{2}}{4\pi \varepsilon _{0}}\int \frac{n(\mathbf{r}^{\prime })}{\left\vert \mathbf{r}-% \mathbf{r}^{\prime }\right\vert }d\mathbf{r}^{\prime }+\frac{4}{3}% C_{x}\left\{ n(\mathbf{r})\right\} ^{1/3}+V_{ext}(\mathbf{r}) = \mu .$

(54)

Comparing this equation to Eq. (15) shows that the Lagrange multiplier $\mu$ represents the classical energy of an electron on the Fermi sphere. It could be identified with the Fermi energy and it really has the meaning of the chemical potential since it is defined by $\frac{\partial E_{el}\left[n\right] }{\partial N} = \mu.$

---

5. Improved Functional Approximations in DFT

With the KS equation, the DFT is an exact theory to the extent that the electron density minimizing the total energy is the true electron density of the interacting electron gas. However, the theory is unworkable as long as the exchange-correlation energy functional $E_{xc} \left[n\right]$ is unknown.

In this section, we describe some improvements of the DFT for an electron gas with slow spatial variations of the electron density such as:

$n(\mathbf{r}) = n_{0}+n_{1}(\mathbf{r}),$

(55)

subject to the conditions $n_{1}(\mathbf{r})\ll$ $n_{0}$ and $\int n_{1}(% \mathbf{r})d\mathbf{r} = 0.$ This assumption is appropriate for the intermediate region in atoms or molecules where the electron density is almost constant, but it is no longer valid near the atomic radius where the electron density is very dense and far away where the electron density is very low.

The merit of Hohenberg and Kohn ^[8] has been to show that the universal energy functional containing the kinetic energy $T_{0}\left[n% \right],$ the Hartree potential $W^{H}\left[n\right]$ and the exchange-correlation energy $E_{xc}\left[n\right]$ in Eq. (35), i.e., $F\left[n\right] = T_{0}\left[n\right] +W^{H}\left[n\right] +E_{xc}% \left[n\right],$ has an expansion of the following form^*:

$F\left[n\right] = F\left[n_{0}\right] -\frac{1}{2}\int\!\!\int K(\left\vert \mathbf{r}-\mathbf{r}^{\prime }\right\vert )n_{1}(\mathbf{r})n_{1}(\mathbf{r}% ^{\prime })d\mathbf{r}d\mathbf{r}^{\prime }+...$

(56)

^*see Appendix Ⅱ of Ref. ^[9].

As expected, this expansion has no linear term in $n_{1}(\mathbf{r})$ because of the condition $\int n_{1}(\mathbf{r})d\mathbf{r} = 0.$

---

5.1. Without Exchange-correlation Effects

In the absence of exchange-correlation effects, HK suggest to expand the total energy (Eq. (35)), free from the exchange-correlation energy, according to the expression:

$E_{el}\left[n\right] = T_{0}\left[n_{0}\right] -\frac{1}{2}\int \!\!\int K(\left\vert \mathbf{r}-\mathbf{r}^{\prime }\right\vert )n_{1}(\mathbf{r}% )n_{1}(\mathbf{r}^{\prime })d\mathbf{r}d\mathbf{r}^{\prime }+\int \left[ w^{H}(\mathbf{r})+V_{ext}(\mathbf{r})\right] n(\mathbf{r})d\mathbf{r, }$

(57)

where the first term $T_{0}\left[n_{0}\right]$ is given by Eq. (49), i.e., $T_{0}\left[n_{0}\right] = C_{K}\int \left(n_{0}\right) ^{5/3}d% \mathbf{r}$, whereas the unknown function $K(\left\vert \mathbf{r}-\mathbf{r}% ^{\prime }\right\vert)$ remains to be determined. To define it, it is convenient to express the total energy in the $\mathbf{q}$-space owing to the Fourier transform^*, whereby the functional becomes:

$E_{el}\left[n\right] = T_{0}\left[n_{0}\right] +\frac{1}{2}\frac{1}{(2\pi )^{3}}\int K(\mathbf{q})\left\vert n(\mathbf{q})\right\vert ^{2}d\mathbf{q}+% \frac{1}{(2\pi )^{3}}\int \left[w^{H}(\mathbf{q})+V_{ext}(\mathbf{q})\right] n(\mathbf{q})d\mathbf{q}\, .$

(58)

^*The Fourier transforms of the different functions look like that of the electron density:

$\begin{eqnarray*} n(\mathbf{q}) & = &\int n_{1}(\mathbf{r})\exp \left( -i\mathbf{q\cdot r}% \right) d\mathbf{r}, \\ n_{1}(\mathbf{r}) & = &\frac{1}{(2\pi )^{3}}\int n(\mathbf{q})\exp \left( i% \mathbf{q\cdot r}\right) d\mathbf{q.} \end{eqnarray*}$

As usual, the electron density would minimize the total energy, i.e., $% \delta E_{el}\left[n\right] = 0.$ In $\mathbf{q}$-space, the minimization of $E_{el}\left[n\right]$ leads to the functional derivative:

$\frac{\delta }{\delta n(\mathbf{q})}\left\{ E_{el}\left[n\right] \right\} = % \frac{1}{(2\pi )^{3}}\int \left[K(\mathbf{q})n(\mathbf{q})+w^{H}(\mathbf{q}% )+V_{ext}(\mathbf{q})\right] d\mathbf{q} = 0,$

(59)

so that the FT of the electron density reads:

$n(\mathbf{q}) = -\frac{w^{H}(\mathbf{q})+V_{ext}(\mathbf{q})}{K(\mathbf{q})}.$

(60)

Now, to determine the function $K(\mathbf{q})$, we proceed to define the linear response function of the uniform interacting electron gas. Recall briefly the implementation of the method. According to the Poisson equation, the Hartree potential $w^{H}(\mathbf{r})$ is related to the varying part of the electron density $n_{1}(\mathbf{r})$ by the relation $\nabla ^{2}w^{H}(% \mathbf{r}) = -\frac{e^{2}}{\varepsilon _{0}}n_{1}(\mathbf{r})$ in real space and by $w^{H}(\mathbf{q}) = \frac{e^{2}}{\varepsilon _{0}q^{2}}n(\mathbf{q})$ in reciprocal space. But the Hartree potential $w^{H}(\mathbf{q})$ is also related to the total potential $\left[w^{H}(\mathbf{q})+V_{ext}(\mathbf{q})% \right]$ by the relation:

$w^{H}(\mathbf{q}) = \frac{e^{2}}{\varepsilon _{0}q^{2}}\left[w^{H}(\mathbf{q}% )+V_{ext}(\mathbf{q})\right] X_{0}(\mathbf{q})\, ,$

(61)

where $X_{0}(q)$ is the linear response function of the uniform electron gas defined by the expression ^[21]:

$X_{0}(q) = -\frac{mk_{F}}{2\pi ^{2}\hbar ^{2}}\left[1+\left( \frac{4-\eta ^{2}% }{4\eta }\right) \ln \left\vert \frac{\eta +2}{\eta-2}\right\vert \right] <0,$

(62)

with $\eta = q/k_{F}.$ Hence, the FT of the electron density $n(\mathbf{q})$ is connected to the total potential by the following relation:

$n(\mathbf{q}) = \left[w^{H}(\mathbf{q})+V_{ext}(\mathbf{q})\right] X_{0}(% \mathbf{q}).$

(63)

From the comparison of Eqs. (60) and (63), it is found that $K(q)$ identifies with $-\left[1/X_{0}(q)\right].$ As expected, one of the most significant feature of $K(q)$ is the logarithm singularity at $% q = 2k_{F}.$ There is a slight inflection in the curve in the neighborhood of $% q = 2k_{F},$ not visible on any scale. Despite this smallness, the singularity may have a large effect on properties depending upon the linear response function $X_{0}(q).$ The function $K(q)$ is expanded well in a power series of $q^{2},$ but the following approximation ^[22] may also be employed to approximate $K(q)$ in the limit of low and large regions of $q$:

$K(q) = \frac{\pi ^{2}\hbar ^{2}}{mk_{F}}\left[\gamma +3\lambda \left( \frac{q% }{2k_{F}}\right) ^{2}\right] >0\, .$

(64)

It is easy to check that this empirical relation fits the low-$q$ region (slow density variations) with the parameters $\gamma = 1$ and $\lambda = 1/9,$ and the large-$q$ region (rapid density variations) with the parameters $% \gamma = -1$ and $\lambda = 1.$

Lastly, we calculate the inverse Fourier transform of $K(q)$ in order to determine the correction $-\frac{1}{2}\int \!\!\int K(\left\vert \mathbf{r}-% \mathbf{r}^{\prime }\right\vert)n_{1}(\mathbf{r})n_{1}(\mathbf{r}^{\prime })d\mathbf{r}d\mathbf{r}^{\prime }$ to bring in the kinetic energy functional after the unknown function $K(\left\vert \mathbf{r}-\mathbf{r}% ^{\prime }\right\vert)$ is determined. In order to facilitate this, we use Eq. (64) to approximate the function $K(q)$. With the definitions:

$\delta (\mathbf{r}) = \frac{1}{(2\pi )^{3}}\int \exp (i\mathbf{q\cdot r})d% \mathbf{q}$

(65)

$\text{and}\quad \quad \nabla _{\mathbf{r}}^{2}\delta (\mathbf{r}) = \frac{1}{% (2\pi )^{3}}\int (-\mathbf{q}^{2})\exp (i\mathbf{q\cdot r})d\mathbf{q, }$

(66)

the inverse FT of $\frac{\pi ^{2}\hbar ^{2}}{mk_{F}(\mathbf{r})}\gamma$ and $\frac{\pi ^{2}\hbar ^{2}}{mk_{F}}3\lambda \left(\frac{q}{2k_{F}}\right) ^{2}$ are $\frac{\pi ^{2}\hbar ^{2}}{mk_{F}(\mathbf{r})}\gamma \delta (% \mathbf{r})$ and $-\frac{3}{4}\frac{\pi ^{2}\hbar ^{2}}{mk_{F}^{3}(\mathbf{r}%)}\lambda \nabla _{\mathbf{r}}^{2}\delta (\mathbf{r}),$ respectively, so that:

$K(\mathbf{r}) = \text{FT}^{-1}\left[K(\mathbf{q})\right] = \left[\frac{\pi ^{2}\hbar ^{2}}{mk_{F}(\mathbf{r})}\gamma-\frac{3}{4}\frac{\pi ^{2}\hbar ^{2}}{mk_{F}^{3}(\mathbf{r})}\lambda \nabla _{\mathbf{r}}^{2}\right] \delta (% \mathbf{r}).$

(67)

Taking into account that $3\pi ^{2}n(\mathbf{r}) = k_{F}^{3}(\mathbf{r}),$ the correction to the kinetic energy becomes:

$-\frac{1}{2}\int \int K(\left\vert \mathbf{r}-\mathbf{r}^{\prime }\right\vert )n_{1}(\mathbf{r})n_{1}(\mathbf{r}^{\prime })d\mathbf{r}d% \mathbf{r}^{\prime } = -\frac{\gamma }{2}\frac{\pi ^{2}\hbar ^{2}}{m}\int \frac{n_{1}^{2}(\mathbf{r})}{k_{F}(\mathbf{r})}d\mathbf{r}+\frac{\lambda }{8}% \frac{\hbar ^{2}}{m}\int \frac{\left\vert \nabla _{\mathbf{r}}n(\mathbf{r}% )\right\vert ^{2}}{n(\mathbf{r})}d\mathbf{r.}$

(68)

It should be mentioned that the gradient correction in the previous relation coincides with the von Weizsäcker ^[23] correction based on intuitive and semiclassical arguments. With a positive contribution, it improves the result of the kinetic energy, since the Thomas-Fermi contribution, $T_{0}\left[n_{0}\right] = C_{K}\int \left(n_{0}\right) ^{5/3}d\mathbf{r, }$ is always smaller than the Hartree-Fock result. Consequently, the expansion proposed by Hohenberg and Kohn for the kinetic energy is valid for both slow and rapid spatial variations of the electron density ^[24].

---

5.2. With Exchange-correlation Effects

To treat the exchange-correlation effects, Hohenberg and Kohn proposed to use the same approximation as in Eq. (56) to expand the exchange-correlation energy, namely:

$E_{xc}\left[n\right] = E_{xc}^{LDA}\left[n_{0}\right] -\frac{1}{2}\int \!\!\int K_{xc}(\left\vert \mathbf{r}-\mathbf{r}^{\prime }\right\vert )n_{1}(% \mathbf{r})n_{1}(\mathbf{r}^{\prime })d\mathbf{r}d\mathbf{r}^{\prime }\, ,$

(69)

where the leading term $E_{xc}^{LDA}\left[n_{0}\right]$ is the so-called local density approximation (LDA) founded on the uniform electron gas, whose the general form is given by Eq. (37), i.e., $E_{xc}^{LDA}\left[n(% \mathbf{r})\right] = \int e_{xc}\left(n(\mathbf{r})\right) n(\mathbf{r})d% \mathbf{r}.$ Usually, the exchange-correlation energy per electron of an interacting electron gas, $e_{xc}\left(n(\mathbf{r})\right),$ splits into two components: the exchange contribution $e_{x}\left(n(\mathbf{r})\right)$ and the correlation contribution $e_{c}\left(n(\mathbf{r})\right).$

For estimating the exchange energy $e_{x}\left(n\right),$ Slater ^[25] suggested the use of the uniform electron gas at constant electron density $n$, known exactly. This contribution may be calculated with the Hartree-Fock approximation via Eq. (9) representing the exchange energy of the electron $i$ in the field of all other electrons:

$E_{x}\left( \mathbf{k}\right) = -\sum\limits_{j\neq i}\int \frac{e^{2}}{4\pi \varepsilon _{0}}\frac{\psi _{j}^{\ast }(\mathbf{r}_{j})\psi _{i}(\mathbf{r}% _{j})\psi _{j}(\mathbf{r}_{i})}{\mathbf{r}_{ij}\psi _{i}(\mathbf{r}_{i})}d% \mathbf{r}_{j}.$

(70)

From here the subscripts $i$ and $j$ will be replaced by $\mathbf{k}$ and $% \mathbf{k}^{\prime },$ so that $E_{x}\left(\mathbf{k}\right)$ is the exchange energy for the state $\mathbf{k}.$ Applied to the uniform electron gas with the wave function $V^{-1/2}\exp (i\mathbf{k\cdot r}),$ the exchange energy for the state $\mathbf{k}$ becomes:

$\begin{eqnarray} E_{x}\left( \mathbf{k}\right) & = &-\frac{e^{2}}{4\pi \varepsilon _{0}V}\sum\limits_{% \mathbf{k}^{\prime }}\left\{ \int \frac{\exp \left[i\left( \mathbf{k}-% \mathbf{k}^{\prime }\right) \cdot \left( \mathbf{r}_{2}-\mathbf{r}% _{1}\right) \right] }{r_{12}}d\mathbf{r}_{2}\right\}, \notag \\ & = &-\frac{e^{2}}{4\pi \varepsilon _{0}V}\sum\limits_{\mathbf{k}^{\prime }}\left\{ \exp \left[-i\left( \mathbf{k}-\mathbf{k}^{\prime }\right) \cdot \mathbf{r}% _{1}\right] \int \frac{\exp \left[i\left( \mathbf{k}-\mathbf{k}^{\prime }\right) \cdot \mathbf{r}_{2}\right] }{r_{12}}d\mathbf{r}_{2}\right\} . \label{HF3} \end{eqnarray}$

(71)

To calculate the integral in brackets, the standard method consists in assuming, for the sake of calculation, that the electrostatic potential $% \varphi (\mathbf{r})$ at point $\mathbf{r}$ is produced by the formal distribution of charge $\exp \left[i\left(\mathbf{k}-\mathbf{k}^{\prime }\right) \cdot \mathbf{r}\right],$ even if it is not the real distribution. Hence, the electrostatic potential at point $\mathbf{r}_{1}$ is $\varphi (% \mathbf{r}_{1}) = \int \frac{\exp \left[i\left(\mathbf{k}-\mathbf{k}^{\prime }\right) \cdot \mathbf{r}_{2}\right] }{r_{12}}d\mathbf{r}_{2}.$ Besides, $% \varphi (\mathbf{r})$ must satisfy the Poisson equation $\nabla ^{2}\varphi (% \mathbf{r}) = -4\pi \exp \left[i\left(\mathbf{k}-\mathbf{k}^{\prime }\right) \cdot \mathbf{r}\right]$ (see footnote in the next page) whose solution is $% \varphi (\mathbf{r}) = 4\pi \frac{\exp \left[i\left(\mathbf{k}-\mathbf{k}% ^{\prime }\right) \cdot \mathbf{r}\right] }{\left\vert \mathbf{k}-\mathbf{k}% ^{\prime }\right\vert ^{2}}.$ Consequently, Eq. (71) simplifies as:

$E_{x}\left( \mathbf{k}\right) = -\frac{e^{2}}{4\pi \varepsilon _{0}V}\sum\limits_{% \mathbf{k}^{\prime }}\frac{4\pi }{\left\vert \mathbf{k}-\mathbf{k}^{\prime }\right\vert ^{2}} = -\frac{e^{2}}{\varepsilon _{0}}\frac{1}{\left( 2\pi \right) ^{3}}\int \frac{d\mathbf{k}^{\prime }}{\left\vert \mathbf{k}-\mathbf{% k}^{\prime }\right\vert ^{2}},$

(72)

where the use is made of the volume element $d\mathbf{k}^{\prime }$ containing $\frac{V}{\left(2\pi \right) ^{3}}$ $d\mathbf{k}^{\prime }$ states, which allows one to transform the sum $\sum_{\mathbf{k}^{\prime }}$ into the integral over the Fermi sphere. The previous integral on $\mathbf{k}% ^{\prime }$ may be carried out with the spherical coordinates ($k^{\prime }\in \left[0, k_{F}\right];$ $\theta \in \left[0, \pi \right];$ $\varphi \in \left[0, 2\pi \right])$, namely, ^*

$\int \frac{d\mathbf{k}^{\prime }}{\left\vert \mathbf{k}-\mathbf{k}^{\prime }\right\vert ^{2}} = 2\pi \int_{0}^{k_{F}}\frac{k^{\prime }}{k}\ln \frac{% \left\vert k+k^{\prime }\right\vert }{\left\vert k-k^{\prime }\right\vert }% dk^{\prime }\, .$

(73)

^*This integration is performed as follows:

$\begin{eqnarray*} \int \frac{d\mathbf{k}^{\prime }}{\left\vert \mathbf{k}-\mathbf{k}^{\prime }\right\vert ^{2}} & = &\int_{0}^{k_{F}}dk^{\prime }\int_{0}^{2\pi }k^{\prime }\sin \theta d\varphi \int_{0}^{\pi }\frac{k^{\prime }d\theta }{% k^{2}+k^{\prime 2}-2kk^{\prime }\cos \theta }, \\ & = &2\pi \int_{0}^{k_{F}}k^{\prime 2}dk^{\prime }\int_{-1}^{+1}\frac{dx}{% k^{2}+k^{\prime 2}-2kk^{\prime }x}, \\ & = &2\pi \int_{0}^{k_{F}}\frac{k^{\prime }}{k}\ln \left\vert \frac{% k+k^{\prime }}{k-k^{\prime }}\right\vert dk^{\prime }. \end{eqnarray*}$

Care must be taken to evaluate this integral because the integrand becomes infinite when $k^{\prime } = k$. Fortunately, the two integrals from $0$ to $k$ and from $k$ to $k_{F}$ converge, so that the integration by parts may be carried out. The final result is^*:

$\int \frac{d\mathbf{k}^{\prime }}{\left\vert \mathbf{k}-\mathbf{k}^{\prime }\right\vert ^{2}} = \pi k_{F}\left[2+\frac{k_{F}^{2}-k^{2}}{k_{F}k}\ln \left\vert \frac{k+k_{F}}{k-k_{F}}\right\vert \right] \, ,$

(74)

^*Use is made of the integrals:

$\begin{eqnarray*} \int k^{\prime }\ln \left\vert k+k^{\prime }\right\vert dk^{\prime } & = &% \frac{1}{2}\left[\left( k^{\prime 2}-k^{2}\right) \ln \left\vert k+k^{\prime }\right\vert +\left( k+k^{\prime }\right) \left( \frac{3}{2}k-% \frac{1}{2}k^{\prime }\right) \right] \, , \\ \int k^{\prime }\ln \left\vert k-k^{\prime }\right\vert dk^{\prime } & = &% \frac{1}{2}\left[\left( k^{\prime 2}-k^{2}\right) \ln \left\vert k-k^{\prime }\right\vert +\left( k-k^{\prime }\right) \left( \frac{3}{2}k+% \frac{1}{2}k^{\prime }\right) \right] \, . \end{eqnarray*}$

and the exchange energy associated with the state $\mathbf{k}$ reads:

$E_{x}\left( \mathbf{k}\right) = -\frac{e^{2}}{4\pi \varepsilon _{0}}\frac{% k_{F}}{2\pi }\left[2+\frac{k_{F}^{2}-k^{2}}{k_{F}k}\ln \left\vert \frac{% k_{F}+k}{k_{F}-k}\right\vert \right] \, .$

(75)

Now, to determine the total exchange energy of the uniform electron gas containing $N$ electrons, it is necessary to sum over the states $\mathbf{k.}$ As usual, the replacement of the sum by an integral over the Fermi sphere yields the expression:

$\sum\limits_{\mathbf{k}}E_{x}\left( \mathbf{k}\right) = -\frac{e^{2}}{4\pi \varepsilon _{0}}\frac{k_{F}}{2\pi }\frac{V}{\left( 2\pi \right) ^{3}}% \int_{0}^{k_{F}}\left[2+\frac{k_{F}^{2}-k^{2}}{k_{F}k}\ln \left\vert \frac{% k+k^{\prime }}{k-k^{\prime }}\right\vert \right] 4\pi k^{2}dk\, .$

(76)

It is found that the integral above is equal to $4\pi k_{F}^{3}.$ Hence, with $k_{F}^{3} = 3\pi ^{2}n,$ the exchange energy per electron of the uniform electron gas is:

$\begin{eqnarray} e_{x}\left( n\right) & = &\frac{1}{N}\sum\limits_{\mathbf{k}}E_{x}\left( \mathbf{k}% \right) = -\frac{e^{2}}{4\pi \varepsilon _{0}}\frac{2V}{\left( 2\pi \right) ^{3}}\frac{k_{F}^{4}}{N} \notag \\ & = &-\frac{e^{2}}{4\pi \varepsilon _{0}}\frac{1}{4\pi ^{3}}\frac{V}{N}\left( 3\pi ^{2}n\right) ^{4/3} = -\frac{e^{2}}{4\pi \varepsilon _{0}}\frac{3}{4}% \left( \frac{3}{\pi }\right) ^{1/3}n^{1/3}. \label{38} \end{eqnarray}$

(77)

Incidentally, note that the exchange energy per electron may be written as a function of the Seitz radius $r_{s} = \left(\frac{3}{4\pi n}\right) ^{1/3} = \left(\frac{9\pi }{4}\right) ^{1/3}\frac{1}{k_{F}}$, instead of the Fermi wavevector, as:

$e_{x}\left( n\right) = -\frac{e^{2}}{4\pi \varepsilon _{0}}\frac{3}{4}\left( \frac{3}{\pi }\right) ^{1/3}\left( \frac{3}{4\pi }\right) ^{1/3}\frac{1}{% r_{s}} = -\frac{0.458}{r_{s}}\, ,$

(78)

where $r_{s}$ is expressed in atomic unit of length and $e_{x}\left(n\right)$ in Hartree, with $\frac{e^{2}}{4\pi \varepsilon _{0}} = 1$. Note also that many investigations ^[7,25,26] of the exchange energy have lead to the expression obtained by multiplying the Slater exchange energy with a factor $\alpha$:

$e_{x}\left( n\right) = -\alpha \frac{e^{2}}{4\pi \varepsilon _{0}}\frac{3}{4}% \left( \frac{3}{\pi }\right) ^{1/3}n^{1/3}.$

(79)

As an example, by calculating the effective exchange potential $e_{x}\left(n\right)$ just for a state at the Fermi level, Kohn and Sham ^[9] have shown that $\alpha = 2/3$ instead of $1$ for Slater.

Finally, the dependence of the exchange energy on $\mathbf{r, }$ aiming at taking account of the inhomogeneity of the electron gas, may be introduced as in Eq. (51). As a result, the LDA exchange energy becomes:

$E_{x}^{LDA}\left[n(\mathbf{r})\right] = -\frac{e^{2}}{4\pi \varepsilon _{0}}% \frac{3}{4}\left( \frac{3}{\pi }\right) ^{1/3}\int n^{4/3}(\mathbf{r})d% \mathbf{r.}$

(80)

For the correlation energy $e_{c}\left(n\right)$, which is the difference between the exact energy of the electron gas and the kinetic and exchange energies, the situation is more complicated since $% e_{c}\left(n\right)$ of the uniform electron gas is not known exactly. Among the old estimations of $e_{c}\left(n\right),$ that of Wigner ^[27], valuable at intermediate densities ($5\geqslant r_{s}\geq 2)$ such as those in metals, is:

$e_{c}\left( n\right) = -\frac{0.44}{r_{s}+7.8}\, ,$

(81)

and that of Gell-Mann and Bruckner ^[28], for high densities ($% r_{s}\leq 2),$ is:

$e_{c}\left( n\right) = +\text{ }0.031\ln r_{s}-0.048\, .$

(82)

By using quantum Monte Carlo (QMC) calculations for the uniform electron gas, Ceperley and Alder ^[29,30] obtained precise numerical results for $e_{c}\left(n\right),$ from which many approximate expressions have been derived. For instance, Perdew and Wang ^[31] proposed the following expression, which fits well the QMC results for a large domain of density:

$e_{c}\left( n\right) = -2c_{0}\left( 1+\alpha r_{s}\right) \ln \left[1+\frac{% 1}{2c_{0}\left( \beta _{1}r_{s}^{1/2}+\beta _{2}r_{s}+\beta _{3}r_{s}^{3/2}+\beta _{4}r_{s}^{2}\right) }\right] \, ,$

(83)

with $c_{0} = 0.03109,$ $\alpha = 0.21370$, $\beta _{1} = 7.5957$, $\beta _{2} = 3.5876$, $\beta _{3} = 1.6382$, $\beta _{4} = 0.49294$.

It should be stressed that the second term of Eq. (69) resulting in a von Weizsäcker gradient correction, as in Eq. (68), almost never improves on the local density approximation. Therefore, a good level of approximation for the LDA is obtained just in doing the sum of the kinetic energy (Eq. (48)), the exchange energy (Eq. (78)) and the correlation energy (Eq. (83)), per electron, of the uniform electron gas at the electron density $n,$ namely:

$e(n) = \frac{1.105}{r_{s}^{2}}-\frac{0.458}{r_{s}}-2c_{0}\left( 1+\alpha r_{s}\right) \ln \left[1+\frac{1}{2c_{0}\left( \beta _{1}r_{s}^{1/2}+\beta _{2}r_{s}+\beta _{3}r_{s}^{3/2}+\beta _{4}r_{s}^{2}\right) }\right] \, .$

(84)

The LDA has been found to give quite good results for most applications with large negative values for $e_{x}\left(n\right)$ and small negative values for $e_{c}\left(n\right)$. The reason of this relative success is in the compensation of the errors, since the LDA underestimates the exchange energy $e_{x}\left(n\right)$ and overestimates the correlation energy $% e_{c}\left(n\right)$, and also in the fact that the exact boundary conditions concerning the exchange-correlation are satisfied in the LDA. Levy et al. ^[32,33] have addressed (on the basis of mathematical requirements) the problems of analysis that arise in discussing Coulomb systems through the density functional approach. Rigorous necessary requisites upon the exchange-correlation functionals have been found in the form of inequalities ^[34], which are satisfied within the LDA and are useful in building and improving functional approximations.

---

6. Beyond the Local Density Approximation

The LDA approach, which is based on the uniform electron gas at a given electron density $n,$ has a formal justification for systems with slow and weak variations. But the atomic and molecular systems, spatially inhomogeneous, differ from the uniform electron gas, therefore more accurate approximations have been introduced making use of the electron density and its gradients.

Since the exchange-correlation energy functional $E_{xc}\left[n\right]$ is not known exactly, the corrections to the LDA approach should ideally contain terms (ⅰ) of well-defined physical origin, (ⅱ) easily usable for practical calculations in physics and chemistry and (ⅲ) allowing transferability from one system to another. Roughly speaking, these approximations are classified into two main groups. One is based on an exact equation for $E_{xc}\left[n\right]$ derived from quantum mechanics including scaling relations, and the other of a more empirical nature is based on the gradient expansion of $E_{xc}\left[n\right]$ in power of small deviation of the density.

---

6.1. Average Exchange Hole

Before discussing these approximations, it is important to introduce the concept of average exchange hole around a given point $\mathbf{r.}$ In any electronic system, the electrons do not move independently but repel one another by the Coulomb repulsion. The coupling (correlation) among the electronic motions are usually called the Coulomb correlations. It should be stressed that Coulomb correlations are ignored in the Hartree and Hartree-Fock methods, because each electron is supposed to move in the average charge distribution of electron system. However, in the Hartree-Fock method where use is made of a determinantal wave function, there are correlations of another type among the positions of electrons with parallel spins due to the Pauli principle. This type of correlation means that (ⅰ) there is zero probability to find two such electrons at same position and (ⅱ) there is a small probability to find an electron in the immediate vicinity of each other. To take account of this aspect of Pauli's principle, it is usually assumed that each electron is surrounded by a hole in the distribution of electrons with parallel spins, called the exchange hole of limited extent.

It can be recalled that the Hartree potential energy (Eq. (5)) of an electron at position $\mathbf{r}_{1}$ is due to a charge distribution whose the density at $\mathbf{r}_{2}$ is :

$-\sum\limits_{j}\frac{e}{4\pi \varepsilon _{0}}\left\vert \psi _{j}(r_{2})\right\vert ^{2}.$

(85)

In addition, the exchange energy (Eq. (9)) of an electron at position $\mathbf{r}_{1}$ is due to an exchange charge distribution whose the density at $\mathbf{r}_{2}$ is:

$+\sum\limits_{j}\frac{e}{4\pi \varepsilon _{0}}\frac{\psi _{j}^{\ast }(\mathbf{r}% _{2})\psi _{i}(\mathbf{r}_{2})\psi _{j}(\mathbf{r}_{1})}{\psi _{i}(\mathbf{r}% _{1})}\, ,$

(86)

so that the total exchange energy reads:

$+\sum\limits_{j}\int \frac{e}{4\pi \varepsilon _{0}}\frac{\psi _{j}^{\ast }(\mathbf{% r}_{2})\psi _{i}(\mathbf{r}_{2})\psi _{j}(\mathbf{r}_{1})}{\psi _{i}(\mathbf{% r}_{1})}d\mathbf{r}_{2} = \frac{e}{4\pi \varepsilon _{0}}\, ,$

(87)

because of the orthogonality and normalization of wave functions, $\int \psi _{j}^{\ast }(\mathbf{r}_{2})\psi _{i}(\mathbf{r}_{2})d\mathbf{r}_{j} = \delta _{ij}.$ Thus, this result suggests that an electron may be viewed as surrounded by a positive cloud of total charge $e,$ such as the electron with its surrounding hole together form a neutral object.

Taking advantage of the definition of exchange charge density given by Eq. (86), we calculate the exchange charge density for the uniform electron gas by using the free electron wave function. The exchange charge density produced at position $\mathbf{r}_{2}$ by an electron in state $% \mathbf{k}$ at $\mathbf{r}_{1}$ reads:

$\sum\limits_{\mathbf{k}^{\prime }}\frac{e}{4\pi \varepsilon _{0}}\frac{\psi _{% \mathbf{k}^{\prime }}^{\ast }(\mathbf{r}_{2})\psi _{\mathbf{k}}(\mathbf{r}% _{2})\psi _{\mathbf{k}^{\prime }}(\mathbf{r}_{1})}{\psi _{\mathbf{k}}(% \mathbf{r}_{1})}\, .$

(88)

After introducing the free electron wave function and changing the sum by an integral over the Fermi sphere, the exchange charge density at position $% \mathbf{r}_{2}$ becomes:

$\frac{e}{4\pi \varepsilon _{0}}\frac{1}{V}\frac{V}{\left( 2\pi \right) ^{3}}% \!\int \!\exp \left[i\left( \mathbf{k}\!-\!\mathbf{k}^{\prime }\right) \!\cdot \!\left( \mathbf{r}_{2}\!-\!\mathbf{r}_{1}\right) \right] d\mathbf{k}% ^{\prime } = \frac{e}{4\pi \varepsilon _{0}}\frac{1}{\left( 2\pi \right) ^{3}}% \exp \left[i\mathbf{k}\!\cdot \!\left( \mathbf{r}_{2}\!-\!\mathbf{r}% _{1}\right) \right] \!\int \!\exp \left[-i\mathbf{k}^{\prime }\!\cdot \!\left( \mathbf{r}_{2}\!-\!\mathbf{r}_{1}\right) \right] d\mathbf{k}% ^{\prime }\, .$

(89)

As usual, the integral in this relation is calculated with spherical coordinates. The exchange charge density produced at position $\mathbf{r}% _{2}$ by an electron in state $\mathbf{k}$ at $\mathbf{r}_{1}$ reads^*:

$\frac{e}{4\pi \varepsilon _{0}}\frac{1}{2\pi ^{2}}\exp \left[i\mathbf{k}% \cdot \left( \mathbf{r}_{2}-\mathbf{r}_{1}\right) \right] \left( \frac{\sin \left( k_{F}r_{12}\right) -\left( k_{F}r_{12}\right) \cos \left( k_{F}r_{12}\right) }{r_{12}^{3}}\right) .$

(90)

^*If $\theta$ is the angle between the vectors $\left(\mathbf{r}_{1}-\mathbf{r}_{2}\right)$ and $\mathbf{k}^{\prime },$ the integral becomes:

$\begin{eqnarray*} \int \exp \left[-i\mathbf{k}^{\prime }\cdot \left( \mathbf{r}_{2}-\mathbf{r}% _{1}\right) \right] d\mathbf{k}^{\prime } & = &\int_{0}^{k_{F}}dk^{\prime }\int_{0}^{2\pi }k^{\prime }\sin \theta d\varphi \int_{0}^{\pi }\exp \left( ik^{\prime }r_{12}\cos \theta \right) k^{\prime }d\theta, \\ & = &2\pi \int_{0}^{k_{F}}k^{\prime 2}dk^{\prime }\int_{0}^{\pi }\exp \left( ik^{\prime }r_{12}\cos \theta \right) \sin \theta d\theta, \\ & = &4\pi \int_{0}^{k_{F}}\frac{k^{\prime }\sin \left( k^{\prime }r_{12}\right) }{r_{12}}dk^{\prime }, \\ & = &4\pi \left( \frac{\sin \left( k_{F}r_{12}\right) -\left( k_{F}r_{12}\right) \cos \left( k_{F}r_{12}\right) }{r_{12}^{3}}\right) . \end{eqnarray*}$

To determine the average exchange charge density, Eq. (90) could be used for a state at the top of Fermi distribution, with the value at $k = k_{F}.$ But it seems better to average Eq. (90) over all the occupied states $\mathbf{k, }$ containing $N/2$ electrons with same spin ^[4]. In this case, the average exchange charge density at position $\mathbf{r}_{2}$ for an electron at $\mathbf{r}_{1}\ $is:

$\frac{e}{4\pi \varepsilon _{0}}\frac{1}{2\pi ^{2}}\left( \frac{\sin \left( k_{F}r_{12}\right) -\left( k_{F}r_{12}\right) \cos \left( k_{F}r_{12}\right) }{r_{12}^{3}}\right) \frac{2}{N}\frac{V}{\left( 2\pi \right) ^{3}}\int \exp % \left[i\mathbf{k}\cdot \left( \mathbf{r}_{2}-\mathbf{r}_{1}\right) \right] d% \mathbf{k}.$

(91)

Owing to the result of the previous footnote, this expression reads:

$\frac{e}{4\pi \varepsilon _{0}}\frac{k_{F}^{6}}{2\pi ^{4}n}\left( \frac{\sin \left( k_{F}r_{12}\right) -\left( k_{F}r_{12}\right) \cos \left( k_{F}r_{12}\right) }{k_{F}^{3}r_{12}^{3}}\right) ^{\!\!2}.$

(92)

In terms of the electron density ($k_{F}^{3} = 3\pi ^{2}n)$, the average exchange charge density at position $\mathbf{r}_{2}$ for an electron at $% \mathbf{r}_{1}\ $becomes:

$\frac{e}{4\pi \varepsilon _{0}}\frac{9n}{2}\left( \frac{\sin \left( k_{F}r_{12}\right) -\left( k_{F}r_{12}\right) \cos \left( k_{F}r_{12}\right) }{k_{F}^{3}r_{12}^{3}}\right) ^{\!\!2}.$

(93)

This function decreases rapidly when the distance $r_{12}$ increases. Besides, when $r_{12}$ tends to zero its value is $\frac{e}{4\pi \varepsilon _{0}}\frac{n}{2}$. This value is opposite to the average density of electrons with same spin. It means that, if the electronic motions were uncorrelated, the average density of electrons with same spin would be ($-% \frac{e}{4\pi \varepsilon _{0}}\frac{n}{2}),$ in accordance with the fact that an electron plus its surrounding hole is a neutral object. Since the electronic motions are correlated, the average density charge distribution of electrons with same spin is given by adding Eq. (93) to ($-% \frac{e}{4\pi \varepsilon _{0}}\frac{n}{2}),$ i.e.:

$-\frac{e}{4\pi \varepsilon _{0}}\frac{n}{2}\left[1-9\left( \frac{\sin \left( k_{F}r_{12}\right)-\left( k_{F}r_{12}\right) \cos \left( k_{F}r_{12}\right) }{k_{F}^{3}r_{12}^{3}}\right) ^{2}\right] = -\frac{e}{4\pi \varepsilon _{0}}\frac{n}{2}g(r_{1}, r_{2}).$

(94)

It can be seen that the term in brackets vanishes for $r_{12} = 0$ and approaches $1$ as $r_{12}$ is increasing. It has been shown ^[35] to be the pair correlation function $g(r_{1}, r_{2})$ of the uniform electron gas with density $n.$ For uncorrelated electronic motions, $% g(r_{1}, r_{2}) = 1.$ It describes the effect of the hole dug into the electron density at $\mathbf{r}_{2}$ by the electron at $\mathbf{r}_{1},$ and it allows us to determined the probability of finding simultaneously two electrons of same spin at given points being a distance $r_{12}$ apart.

---

6.2. Exchange-correlation Energy Functional

Now we are in position to determine the exact expression for the exchange-correlation energy functional $E_{xc}\left[n(\mathbf{r})% \right]$ in interacting systems resulting from the interaction between an electron and its exchange-correlation hole $n_{xc}\left(\mathbf{r}_{1}, % \mathbf{r}_{2}\right)$. As already mentioned, the formal definition of the exchange-correlation energy $E_{xc}\left[n(\mathbf{r})\right]$ is given by Eq. (36),

$E_{xc}\left[n\right] = \left( T\left[n\right] +U_{ee}\left[n\right] \right) -T_{0}\left[n\right] -W^{H}\left[n\right],$

(95)

but, in the DFT context, it has been proved useful to consider another expression, known as the adiabatic connection formula, differing from the approximate form of Hohenberg and Kohn (Eq. (69)) by the presence of the exchange-correlation hole $n_{xc}\left(\mathbf{r}_{1}, \mathbf{r}_{2}\right).$ To achieve this, it is assumed that a system in a stationary state remains in its initial state when the parameter $\lambda$ describing the system changes very slowly. Using this parameter $\lambda,$ the exchange-correlation energy functional $E_{xc}\left[n(\mathbf{r})\right]$ containing the kinetic energy of non interacting electron gas may be written as:

$\begin{eqnarray} E_{xc}\left[n\right] & = &\left\langle \psi _{\lambda }\left\vert T\left[n% \right] +\lambda U_{ee}\left[n\right] \right\vert \psi _{\lambda }\right\rangle _{\lambda = 1}-\left\langle \psi _{\lambda }\left\vert T\left[ n\right] +\lambda U_{ee}\left[n\right] \right\vert \psi _{\lambda }\right\rangle _{\lambda = 0}-W^{H}\left[n\right], \notag \\ & = &\int_{0}^{1}\frac{d}{d\lambda }\left\langle \psi _{\lambda }\left\vert T% \left[n\right] +\lambda U_{ee}\left[n\right] \right\vert \psi _{\lambda }\right\rangle d\lambda -W^{H}\left[n\right] . \label{Exc1} \end{eqnarray}$

(96)

The adiabatic connection permits us to pass from the interacting system ($% \lambda = 1)$ to the noninteracting system ($\lambda = 0)$ by infinite number of possibilities. In the Hamiltonian of KS equation (Eq. (41)), the external potential $\lambda V_{ext}(\mathbf{r})$ and the Coulomb interaction $\frac{e^{2}}{4\pi \varepsilon _{0}}\frac{\lambda }{r_{12}}$ of partial strength $\lambda \ $are taken in order that the density $n(\mathbf{r})$ is the same at any $\lambda$ between $0$ and $1.$

Eq. (96) may be transformed with the Hellman-Feynman theorem ^[36,37,38] under the following form^*:

$E_{xc}\left[n\right] = \int_{0}^{1}\left\langle \psi _{\lambda }\left\vert U_{ee}\left[n\right] \right\vert \psi _{\lambda }\right\rangle d\lambda -W^{H}\left[n\right] .$

(97)

^*The Hellman-Feynman theorem allows us to know how the energy $E_{\lambda } = \left\langle \psi _{\lambda }\left\vert H_{\lambda }\right\vert \psi _{\lambda }\right\rangle$ varies as a function of $\lambda.$ Due to the properties of the Hamiltonian and the wavefunction, it is found that $\frac{% dE_{\lambda }}{d\lambda } = \left\langle \psi _{\lambda }^{\ast }\left\vert \frac{\partial H_{\lambda }}{\partial \lambda }\right\vert \psi _{\lambda }\right\rangle.$ The theorem has been first used to determine the forces in a molecular system.

Of course, $W^{H}\left[n\right]$ is given by Eq. (34).

$W^{H}\left[n\right] = \frac{1}{2}\frac{e^{2}}{4\pi \varepsilon _{0}}\int \int \frac{n(\mathbf{r}_{1})n(\mathbf{r}_{2})}{\left\vert \mathbf{r}_{1}-% \mathbf{r}_{2}\right\vert }d\mathbf{r}_{1}d\mathbf{r}_{2}.$

(98)

On the other hand, the first term of the right-hand side of Eq. (97) has been proved to be equal to ^[39]:

$\int_{0}^{1}\left\langle \psi _{\lambda }\left\vert U_{ee}\left[n\right] \right\vert \psi _{\lambda }\right\rangle d\lambda = \frac{1}{2}\frac{e^{2}}{% 4\pi \varepsilon _{0}}\int \int \frac{n(\mathbf{r}_{1})n_{2}\left( \mathbf{r}% _{1}, \mathbf{r}_{2}\right) }{\left\vert \mathbf{r}_{1}-\mathbf{r}% _{2}\right\vert }d\mathbf{r}_{1}d\mathbf{r}_{2},$

(99)

with:

$n_{2}\left( \mathbf{r}_{1}, \mathbf{r}_{2}\right) = n(\mathbf{r}% _{2})+n_{xc}\left( \mathbf{r}_{1}, \mathbf{r}_{2}\right),$

(100)

so that the exchange-correlation energy functional $E_{xc}\left[n(\mathbf{r})\right]$ simplifies as:

$E_{xc}\left[n(\mathbf{r})\right] = \frac{1}{2}\frac{e^{2}}{4\pi \varepsilon _{0}}\int n(\mathbf{r}_{1})d\mathbf{r}_{1}\int \frac{n_{xc}\left( \mathbf{r}% _{1}, \mathbf{r}_{2}\right) }{\left\vert \mathbf{r}_{1}-\mathbf{r}% _{2}\right\vert }d\mathbf{r}_{2}.$

(101)

As a generalization of Eq. (94), the exchange-correlation hole $% n_{xc}\left(\mathbf{r}_{1}, \mathbf{r}_{2}\right)$ is the following integral ^[40] over the coupling parameter $\lambda$:

$n_{xc}\left( \mathbf{r}_{1}, \mathbf{r}_{2}\right) = n(\mathbf{r}% _{2})\int_{0}^{1}d\lambda \left[g(r_{1}, r_{2};\lambda )-1\right] .$

(102)

Here, $g(r_{1}, r_{2};\lambda)$ is the pair correlation function of the system with the physical density $n(\mathbf{r})$ calculated from the Coulomb interaction potential $\frac{e}{4\pi \varepsilon _{0}}\frac{\lambda }{r_{12}}$. By its definition, the exchange-correlation hole $n_{xc}\left(\mathbf{r}_{1}, \mathbf{r}_{2}\right)$ must satisfy the sum rule requiring that the hole contains one electron, for all $\mathbf{r}$, i.e.:

$\int n_{xc}\left( \mathbf{r}_{1}, \mathbf{r}_{2}\right) d\mathbf{r}_{2} = -1,$

(103)

reflecting total screening of the electron at $\mathbf{r}_{1}.$ The pair correlation function $g(r_{1}, r_{2};\lambda)$ may be obtained from different theoretical approaches and computer simulations. However, for inhomogeneous systems little is known about it.

Several approximations based on Eq. (101) have been proposed to improve the LDA approach, bypassing the pair correlation function, since they only need the exchange-correlation energy per electron $e_{xc}\left(n(% \mathbf{r})\right)$ to determine $E_{xc}\left[n(\mathbf{r})\right] = \int e_{xc}\left(n(\mathbf{r})\right) n(\mathbf{r})d\mathbf{r, }$ where

$e_{xc}\left( n(\mathbf{r})\right) = \frac{1}{2}\frac{e^{2}}{4\pi \varepsilon _{0}}\int \frac{n_{xc}\left( \mathbf{r}, \mathbf{r}_{2}\right) }{\left\vert \mathbf{r}-\mathbf{r}_{2}\right\vert }d\mathbf{r}_{2}.$

(104)

Gunnarsson et al. ^[41] proposed two nonlocal approximations: (ⅰ) the averaged-density (AD) approximation, which replaces $n(\mathbf{r}% _{2})$ in Eq. (102) by the averaged density $\overline{n}(\mathbf{r}% _{1})$ obtained with a weight function to be determined, and (ⅱ) the weighted-density (WD) approximation that uses the pair correlation function of the homogeneous electron gas in Eq. (102). Langreth and Perdew ^[42,43] developed a similar theory for the exchange-correlation energy by investigating Eq. (101) in reciprocal space. The AD and WD approximations have attracted much attention as a potential tool for investigating solid state structures, metallic surfaces and quantum chemistry ^[44,45]. During the same period of time, some simple parametrized expressions have been introduced to describe the exchange-correlation energy ^[46,47,48].

Despite the use of the adiabatic connection formula making the explicit link between exchange-correlation energy $E_{xc}\left[n(\mathbf{r})\right]$ and exchange-correlation hole $n_{xc}\left(\mathbf{r}_{1}, \mathbf{r}_{2}\right)$, the results may deviate substantially from the simulation ones. Therefore, the adiabatic connection formula is much less used than the so-called gradient expansion approximations. In this class of approximations, one tries to systematically calculate the gradient corrections by a generalization of the famous lowest-order von Weizsäcker gradient correction.

Since any real system has a spatially varying density $n(\mathbf{r})$, it would clearly be useful to include information on the rate of variation in the exchange-correlation energy $E_{xc}\left[n(\mathbf{r})\right]$. A first attempt at doing this was to use Eq. (69) even if, as already mentioned, the lowest order of the expansion is often disappointing compared with the LDA results. Then, higher order corrections have been included. Using the new coordinates $\mathbf{r = r}^{\prime }-\mathbf{r}^{\prime \prime }$ and $\mathbf{R = }\left(\mathbf{r}^{\prime }+\mathbf{r}^{\prime \prime }\right) /2,$ for a slowly varying density, the electron density reads ^[49]:

$n(\mathbf{r}^{\prime }) = n(\mathbf{R})+(\frac{1}{2}\mathbf{r\cdot \nabla )}n(% \mathbf{R})+\frac{1}{2}(\frac{1}{2}\mathbf{r\cdot \nabla )}^{2}n(\mathbf{R}% )+...$

(105)

and the exchange-correlation energy functional (Eq. (69)) may be written in terms of the following gradient expansion:

$\begin{eqnarray} E_{xc}\left[n\right] & = &\int e_{xc}\left( n(\mathbf{R})\right) n(\mathbf{R}% )d\mathbf{R}-\frac{\pi }{3}\int \left\vert \mathbf{\nabla }n(\mathbf{R}% )\right\vert ^{2}d\mathbf{R}\int_{0}^{\infty }r^{4}K_{xc}(r;n(\mathbf{R}))dr \notag \\ &&-\frac{\pi }{60}\int \mathbf{\nabla }n(\mathbf{R})\mathbf{\cdot \nabla }% \left[\mathbf{\nabla }^{2}n(\mathbf{R})\right] d\mathbf{R}\int_{0}^{\infty }r^{6}K_{xc}(r;n(\mathbf{R}))dr+..., \label{Exc4} \end{eqnarray}$

(106)

where the leading term corresponds to the LDA contribution $E_{xc}^{LDA} \left[n(\mathbf{R})\right].$ The other terms are very difficult to calculate and little is known about them.

---

7. Separation of the Exchange and Correlation Functionals

The failures of the LDA for strongly inhomogeneous systems have lead to going beyond the lowest order gradient correction with different sequences of the series given in Eq.(106). Important progress has been made in deriving gradient expansion approximations for the exchange-correlation functional, under the general form:

$E_{xc}\left[n\right] = \int e_{xc}\left( n(\mathbf{r})\right) f_{xc}\left( n(% \mathbf{r}), \mathbf{\nabla }n(\mathbf{r})\right) n(\mathbf{r})d\mathbf{r.}$

(107)

This relation is known as the generalized-gradient approximation (GGA); the function $f_{xc}\left(n(\mathbf{r}), \mathbf{\nabla }n(\mathbf{r})\right)$ is the enhancement factor depending on the density and the density gradient. The construction of $f_{xc}\left(n(\mathbf{r}), \mathbf{\nabla }n(\mathbf{r})\right)$ has been made by considering dimensional analysis, sums rules as well as asymptotic behavior of density and effective potentials around atoms or molecules. It is a matter of fact that different GGAs differ greatly in the choice of the function $% f_{xc}\left(n(\mathbf{r}), \mathbf{\nabla }n(\mathbf{r})\right).$

Since the exchange energy $E_{x}$ is given by an exact expression, contrary to the correlation energy $E_{c}$ (Section 5), and since $% \left\vert E_{x}\right\vert$ is greater than $\left\vert E_{c}\right\vert,$ it seems natural to treat separately the exchange energy and the correlation energy. In this section, two approximations for pure exchange functional and correlation functional will be briefly presented, accounting for the most popular DFT investigations.

---

7.1. Exchange Energy Functionals

For the exchange energy functionals, Becke ^[11] showed that $E_{x}$ can be constructed with one parameter fitted to exchange energies of atoms. This corresponding semiempirical functional, commonly denoted $% E_{x}^{B88},$ reproduces the exact asymptotic behavior of $E_{x}$ and the effective potential. For spin polarized systems, its expression reads:

$E_{x}^{B88} = E_{x}^{LDA}-\beta \sum\limits_{\alpha }\int n_{\alpha }^{4/3}\frac{% x_{\alpha }^{2}}{\left( 1+6\beta x_{\alpha }\sinh ^{-1}x_{\alpha }\right) }d% \mathbf{r},$

(108)

where $x_{\alpha } = \frac{\left\vert \nabla n_{\alpha }\right\vert }{% n_{\alpha }^{4/3}}$ and $n_{\alpha }$ stands for spin up and spin down densities. It should be noted that this functional contains the single parameter $\beta = 0.0042$. If $\sinh ^{-1}x_{\alpha }$ is expanded into a power series, the $E_{x}^{B88}$ functional may be written under the gradient expansion:

$E_{x}^{B88} = -\frac{e^{2}}{4\pi \varepsilon _{0}}\frac{3}{4}\left( \frac{3}{% \pi }\right) ^{1/3}\int n^{4/3}(\mathbf{r})\left[1+\mu x_{\alpha }^{2}+...% \right] d\mathbf{r},$

(109)

where the leading term is $E_{x}^{LDA}\left[n(\mathbf{r})\right]$ of Eq. (80). The $E_{x}^{B88}$ functional is particularly well adapted to the calculation of the atomization energies of molecules, i.e., the required energy to fully dissociate a molecule to free atoms.

Other very convenient exchange energy functionals have been developed by the Perdew group. These functionals are based on rational principles of quantum mechanics and theoretical constraints. In particular, for spin unpolarized systems, the so-called PBE (Perdew-Burke-Ernzerhof ^[50,51]) exchange is often used under the form:

$E_{x}^{PBE} = -\frac{e^{2}}{4\pi \varepsilon _{0}}\frac{3}{4}\left( \frac{3}{% \pi }\right) ^{1/3}\int n^{4/3}(\mathbf{r})\left[1+\kappa-\frac{\kappa }{% (1+\frac{\mu s^{2}}{\kappa })}\right] d\mathbf{r},$

(110)

where $s = \frac{\left\vert \nabla n\right\vert }{2k_{F}n}$ (with $% k_{F}^{3} = 3\pi ^{2}n)$ and the two universal parameters are $\mu = 0.2195$ and $\kappa = 0.804.$ The PBE functional, bringing minor modifications to the PW (Perdew-Wang ^[31]) functional, increases the precision of the results of total energies of atoms and atomization energies of molecules. In addition, its recent version ^[52] improves equilibrium properties of packed solids and their surfaces.

---

7.2. Correlation Energy Functionals

Correlation energy functionals have been also constructed by the Perdew group. The first such expression, known as PW (Perdew-Wang ^[31]) correlation, has been refined and accurately represented by the expression:

$E_{c}^{PW} = \int \left[e_{c}(r_{s}, \zeta )+H(r_{s}, \zeta, t)\right] n(% \mathbf{r})d\mathbf{r.}$

(111)

In this relation, the variables are the Seitz radius $r_{s} = \left(\frac{3}{% 4\pi n}\right) ^{1/3}$ introduced in place of the Fermi wavevector $k_{F}$ (with $r_{s}^{3}k_{F}^{3} = 9\pi /4),$ the relative spin polarization $\zeta = % \frac{n_{1}-n_{2}}{n_{1}+n_{2}}$, which ranges from $0$ for unpolarized system to $\pm 1$ for a fully spin polarized system, and the reduced density gradient $t$ measuring the rapidity of variation of $n(r)$ on the scale of the screening length $1/k_{s}.$ The reduced density gradient $t$ is defined as:

$t = \frac{\left\vert \nabla n\right\vert }{2\phi k_{s}n},$

(112)

where $k_{s} = \left(\frac{4k_{F}}{\pi }\right) ^{1/2}$ is the Thomas-Fermi wavevector and $\phi (\zeta) = \frac{1}{2}\left[(1+\zeta)^{2/3}+(1-\zeta)^{2/3}\right]$ is the spin interpolation formula taking account of the polarization dependence.

Accurate analytic representations of $e_{c}(r_{s}, \zeta)$ are available under the general form ^[31]:

$e_{c}(r_{s}, \zeta ) = e_{c}(r_{s}, 0)+\phi (\zeta )\left[ e_{c}(r_{s}, 1)-e_{c}(r_{s}, 0)\right],$

(113)

where $e_{c}(r_{s}, 0)$ is the correlation energy of the uniform electron gas given by Eq. (83). As far as $H(r_{s}, \zeta, t)$ in Eq. (111) is concerned, it reads:

$H(r_{s}, \zeta, t) = \gamma \phi ^{3}\ln \left\{ 1+\frac{\beta }{\gamma }t^{2}% \left[\frac{1+At^{2}}{1+At^{2}+A^{2}t^{4}}\right] \right\},$

(114)

where the more recent values of the universal parameters ^[50] are $\gamma = 0.031091,$ $\beta = 0.066725$ and

$A = \frac{\beta }{\gamma }\left[\exp \left(-\frac{e_{c}(r_{s}, \zeta )}{% \gamma \phi ^{3}}\right)-1\right] ^{-1}.$

(115)

The PW correlation energy functional (Eq. (111)) has been derived from the slowly varying limit ($t\rightarrow 0)$ and the rapidly varying limit ($t\rightarrow \infty)$, as well as from the uniform scaling to the high density limit ($\gamma \rightarrow \infty),$ which arise from numerical calculations of the local density and the gradient contributions ^[50]. It is found that $H(r_{s}, \zeta, t)\sim \beta \phi ^{3}t^{2}$ when $t\rightarrow 0$ and $H(r_{s}, \zeta, t)\sim -e_{c}(r_{s}, \zeta)$ when $% t\rightarrow \infty,$ so that the correlation energy per electron vanishes in systems with rapid spatial variations of the electron density. Since the second-order gradient contribution $H(r_{s}, \zeta, t)$ is a growing function of $t$, going from $0$ to $\left\vert e_{c}(r_{s}, \zeta)\right\vert,$ $% E_{c}^{PW}$ tends to a small negative constant.

This is a much more difficult problem to test the correlation energy than the exchange energy, which is accurately defined. Therefore, a useful way to think through gradient approximations is to compare one GGA with another. There is an interesting correlation energy functional derived by Lee, Yang and Parr ^[12] from the early work of Colle and Salvetti ^[53] and known as LYP correlation energy, whose results compare favorably with those by $E_{c}^{PW}.$ For the unpolarized system, the LYP correlation energy may be written in terms of electron density as:

$E_{c}^{LYP} = -a\int \left\{ n(\mathbf{r})+bn(\mathbf{r})^{-2/3}\left[t_{HF}(% \mathbf{r})-2t_{W}(\mathbf{r})\right] \exp \left[-cn(\mathbf{r})^{-1/3}% \right] \right\} \frac{n(\mathbf{r})^{1/3}}{n(\mathbf{r})^{1/3}+f}d\mathbf{r, % }$

(116)

where the constants are $a = 0.049$, $b = 0.132$, $c = 0.2533$ and $f = 0.349$, $% t_{W}(\mathbf{r})$ is the local Weizsacker kinetic energy density, and $% t_{HF}(\mathbf{r})$ the Hartree-Fock kinetic energy to second order defined as:

$t_{W}(\mathbf{r}) = \frac{\left\vert \nabla n(\mathbf{r})\right\vert ^{2}}{% 8n(\mathbf{r})}-\frac{\nabla ^{2}n(\mathbf{r})}{8},$

(117)

$t_{HF}(\mathbf{r}) = C_{K}n(\mathbf{r})^{5/3}+\left[\frac{t_{W}(\mathbf{r})% }{9}+\frac{\nabla ^{2}n(\mathbf{r})}{18}\right],$

(118)

with $C_{K} = \frac{3}{10}\left(3\pi ^{2}\right) ^{2/3}.$

Note that the combination of $E_{x}^{B88},$ for the exchange energy, with $% E_{c}^{LYP},$ for the correlation energy, forms the GAUSSIAN program package widely used by the chemistry community. Obviously, many other combinations of the exchange energy and the correlation energy are possible. The most popular exchange-correlation functionals are composed of $% E_{x}^{B88}$ for exchange, and $E_{c}^{PW}$ or $E_{c}^{LYP}$ for correlation. These functionals provide a good compromise between computational efficiency and numerical accuracy, but are not sufficient for a correct description of some chemical and physical properties. Particularly, they do not yield significantly better results than LDA for solid state properties ^{[54,55,56,57]} and they fail for the van der Waals dispersion interactions ^[58,59,60]. However, promising exchange and correlation functionals have been developed by including additional information depending explicitly on higher order density gradients ^[61,62], i.e., the second derivative of the density.

---

8. A Classification of the Approximations and Main DFT Programs

To rank the approximated functionals, Perdew ^[63] introduced a DFT Jacob's ladder, as a guideline from the most approximate to the more sophisticated approaches for the functionals. The ladder has five rungs corresponding to different choices in which way exactly do the energy density $n(r)$ and its derivatives enter. In principle, when climbing the DFT ladder progressively, higher overall accuracy can be expected until exact results are reached. However, as with the ladder, it does not necessarily mean higher accuracy for every particular system and property but rather a smaller probability to fail.

The first rung stands for the local spin density (LSD) approximation, proposed initially by Kohn and Sham ^[9] and using only the energy density $n_{\alpha}$ ($\alpha$ being used to label the spin components).

The second rung corresponds to the GGA (generalized gradient approximation) using $\nabla n_{\alpha}$. It is an extension of the LSD that is widely used in quantum chemistry while the LSD is the most popular approximation for calculating the electronic structure in solid state. Two different philosophies may be adopted. The first one searches for exchange-correlation functionals that obey physical constraints, scaling properties and sum rules. The second one seeks for a parametric form where the parameters are fitted to a set of experiments. They fall into the category of semiempirical xc functionals.

The third rung is the meta-GGA including the higher order density gradients $% \nabla ^{2}n_{\alpha }$. It requires an explicit dependence on the kinetic energy density, and a combination of theoretical constraints satisfying the exact function and functionals that use few parameters fitting properties measured or calculated by other quantum chemistry methods. Several meta-GGA functionals have been constructed using the orbital kinetic energy densities, giving greater accuracy over a wide range of systems and properties ^[64,65,66]. Concise overview of these functionals are given by Perdew et al. ^[52].

The fourth rung is representative of the hyper-GGA making use of the exact exchange energy $E_{x}^{exact}$ and of a correlation energy functional $% E_{c}$. A notable improvement of the exchange-correlation GGAs appeared when Becke ^[67,68] proposed to mix some fraction of exchange and correlation GGA functionals with the exact exchange energy $E_{x}^{exact}$. This approach allowed to generate a new class of functionals, known as the hybrid ones, the simplest among which may read:

$E_{xc}^{hyb} = aE_{x}^{exact}+(1-a)E_{x}^{GGA}+E_{c}^{GGA},$

(119)

where $E_{x}^{exact}$ is the exact exchange energy obtained from Hartree-Fock method, $E_{x}^{GGA}$ the exchange GGA functional and $% E_{c}^{GGA}$ the correlation GGA functional. Note that $a(\sim 0.25)$ differs from $1$ because the exact exchange energy has a highly nonlocal character compared to the GGA exchange. It should be mentioned that neither the Hartree method nor the Hartree-Fock method are able to calculate the exact correlation energy. Expressions such as Eq. (119) are implemented in many DFT code packages and give almost identical results in specific applications. These explicitly nonlocal, exact exchange based functionals look very promising in use for chemical physics.

The fifth rung or fully nonlocal description refers to exact exchange combined with exact partial correlation, making use of the occupied and unoccupied Kohn-Sham orbitals. The strategy amounts to identify the part of the exchange-correlation energy for which the GGA is accurate and to treat the remainder exactly ^[69,70]. As an alternative to the standard Kohn-Sham procedure, other exact realizations of density functional theory may be carried out. Such generalized Kohn-Sham schemes minimize not only the kinetic energy for fixed density, as in the standard Kohn-Sham method, but also part of the electron-electron energy ^[71,72] with the result of gaining more flexibility.

While the performances of hybrid GGA and hybrid meta-GGA functionals are traditionally better than just plain GGA, and particularly better than LDA, results vary from functional to functional. The performances of each functional lined up on the DFT ladder depend on the type of system under consideration and on the property under evaluation. The number of density functionals currently available is very large and, among the plethora of density functionals, the difficulty lies in selecting the best one for each specific application. There is a great number of density functionals implemented in the DFT codes^*. Limited comparisons of the performance for density functionals can be found in the literature. For instance, the performances of a large number of DFT functionals have been evaluated across about fifty different density functionals ^[73] in the calculation of several properties of chemical interest against experimental data (structure, kinetics, thermochemistry, nonbounded interactions). A rigorous extension of the theory to these classes of properties is possible, but it requires the introduction of new functionals that depend not only on density but also on properties.

^* Here are a few official sites of distributed code packages that the reader may consult for calculations of quantum chemistry: ADF (http://www.scl.com/); CPMD (http://www.cpmd.org/); GAUSSIAN (http://www.gaussian.com/) and quantum physics: ABINIT (http://www.abinit.org/); VASP (http://www.vasp.at/); WIEN (http://www.wien2k.at/); Quantum ESPRESSO (Giannozzi P, et al. (2009) J Phys-Condens Mat 21: 395502.)

In this paragraph, we briefly present few ubiquitous hybrid functionals employed in the currently available DFT program packages, for the calculation of a large variety of molecular and chemical properties. The most popular hybrid functionals in computational chemistry are the so-called $E_{xc}^{B3PW91},$ $E_{xc}^{B3LYP}$ and $E_{xc}^{B3PBE}$ functionals.

(ⅰ) In the $E_{xc}^{B3PW91}$ functional, Eq. (108) is used for exchange and Eq. (111) for correlation, both mixed via the 3 parameters $a,$ $b$ and $c,$ according to the relationship:

$E_{xc}^{B3PW91} = E_{xc}^{LDA}+a(E_{x}^{exact}-E_{x}^{LDA})+b(E_{x}^{B88}-E_{x}^{LDA})+c(E_{c}^{PW91}-E_{c}^{LDA}),$

(120)

with $a = 0.20,$ $b = 0.72$ and $c = 0.81.$ The weights of the various contributions to the exchange and correlation energies are obtained by a fit either of the atomization energy of reference data ^[74] or of theoretical constraints ^[75].

(ⅱ) The $E_{xc}^{B3LYP}$ functional employs the same expression as Eq. (120) with the same parameters $a,$ $b$ and $c$, whereas the $% E_{c}^{PW91}$ correlation functional is replaced by the $E_{c}^{LYP}$ one. It should be noted that this type of functional is less successful to predict the solid state properties than the chemical properties, because of the difficulties in computing the exact (Hartree-Fock) exchange with a plane waves basis set.

(ⅲ) The $E_{xc}^{B3BPE}$ functional is the combination of the $E_{x}^{B88}$ functional for the exchange with the $E_{c}^{PBE}$ one for the correlation, with $a = 0.25,$ $b = 0.75$ and $c = 1.$ The value of 25% of the exact exchange energy, instead of 20%, is supported by a strong theoretical basis ^[76,77]. These code packages are flexible enough to generate new density functionals from any combination of exchange and correlation functionals with a fraction of the Hartree-Fock exchange energy (HF/DFT).

Standard DFT is particularly efficient to study systems where the origin of interactions is essentially electrostatic. But the interest of DFT is also to treat systems, such as biological systems, complex materials, surfaces, ..., where the dispersion interactions play a crucial role. It is a matter of fact, that the van der Waals dispersion interaction is highly problematic for conventional functionals, because no known correlation functional could generate the London dispersion interactions $(\sim -r^{-6})$ in contrast to the asymptotic interactions falling off exponentially. The treatment of the dispersion interactions within DFT is in full expansion according two principal routes. The first one, generally known as DFT-D, consists in adding a more or less empirical dispersion correction, $-C_{6}/r^{-6}$, that would account for the missing long range attraction. In this approach, the dispersion coefficient $C_{6}$ is either fitted on experimental ionization potentials and polarizabilities of isolated atoms ^[78] or obtained with a non empirical approach making $C_{6}$ environment dependent ^[79,80,81]. The second route, known as vdW-DF, describes the dispersion interaction through a DFT functional by introducing a nonlocal correlation functional ^[82,83]. It is a promising computational tool for weakly bound systems ^[84]. Some versions of the vdW-DF are implemented in widely distributed DFT software package. Compared to the first approach, the second one improves the precision of the results but it is offset by the increase of computational time by about 50%.

DFT is a powerful tool to investigate the static properties of electronic systems (geometrical structures, relative energies, ...). It is also a convenient tool to be used in conjunction with ab initio molecular dynamics (AIMD) simulations for extended systems, as described by Car and Parrinello ^[85]. Whereas classic MD describes trajectories of atoms as objects within an empirical interaction potential, the AIMD computations are simulations that involve the motions of both nuclei and electrons. For the latter, the DFT model is used, and the evolution of the kinetic energies of both electrons and nuclei is observed in an identical way, through the use of fictitious masses. Plane waves are usually employed for the description of valence orbitals, whereas pseudo-potential approximations are used for taking into account the core electrons. The combination of Car-Parrinello quantum mechanics/molecular mechanics (QM/MM) is certainly one of the most promising theoretical tools available for theoretical chemistry. Its implementation for ab initio molecular dynamics simulation of large systems has also been successfully applied to solid state and liquid state physics, as well as to materials science.

For molecular properties, the semiempirical $E_{xc}^{B3LYP}$ functional is very successful in describing a wide range of molecular properties. For periodic systems and crystalline calculations, nonempirical HF/DFT functionals are necessary but they are computationally expensive, especially for systems with metallic characteristics. The situation improves with the HSE functionals, developed by Heyd et al. ^[86,87], which use a screening technique to take advantage of the fast spatial decay of the short-range Hartree-Fock exchange. Such functionals are based on the PBE exchange-correlation functional ^[88]: the exchange energy term is split into short-range and long-range components and the HF long-range is replaced by the PBE long-range. The results obtained with HSE applied to metals exhibit significantly smaller errors than pure DFT calculations. It has been shown that the $E_{xc}^{B3LYP}$ functional yields unsatisfactory atomization energies of solids underestimated by approximately 17% and lattice parameters overestimated by about 1%, when compared with those of PBE and HSE which perform reasonably well for materials with localized electrons and free electron metals ^[89].

The theoretical description of matter as well as of many chemical, physical and biological processes requires accurate methods for investigating atomic and molecular-scale interactions, whose origins and perspectives have been widely discussed recently ^[20,90,91].

---

Acknowledgments

I am grateful to Andrei Postnikov for stimulating discussions and its very careful reading of the manuscript.

---

Conflict of Interest

There is no conflict of interest regarding the publication of this article.

---