3.1 Density-functional theory

In this section we will describe the remarkable theorems of density-functional theory (DFT) which allow us to find ground-state properties of a system without dealing directly with the many-electron state $\vert \Psi \rangle$. We deal with a system of $N$ electrons moving in a static potential, and adopt a conventional normalisation in which $\langle \Psi \vert \Psi \rangle = N$.

3.1.1 The Hohenberg-Kohn theorems

As a result of the Born-Oppenheimer approximation, the Coulomb potential arising from the nuclei is treated as a static external potential $V_{\mathrm{ext}}({\bf r})$:

$$V_{\mathrm{ext}}({\bf r}) = - \sum_{\alpha} \frac{Z_{\alpha}}{\left\vert {\bf r} - {\bf r}_{\alpha} \right\vert} .$$ (3.1)

We define the remainder of the electronic Hamiltonian given in (2.19) as ${\hat F}$:

$${\hat F} = -\frac{1}{2} \sum_i \nabla_i^2 + \frac{1}{2} \su... ...not= i} \frac{1}{\left\vert {\bf r}_i - {\bf r}_j \right\vert}$$ (3.2)

such that ${\hat H} = {\hat F} + {\hat V}_{\mathrm{ext}}$ where

$${\hat V}_{\mathrm{ext}} = \sum_i V_{\mathrm{ext}}({\bf r}_i) .$$ (3.3)

${\hat F}$ is the same for all $N$-electron systems, so that the Hamiltonian, and hence the ground-state $\vert \Psi_0 \rangle$, are completely determined by $N$ and $V_{\mathrm{ext}}({\bf r})$. The ground-state $\vert \Psi_0 \rangle$ for this Hamiltonian gives rise to a ground-state electronic density $n_0({\bf r})$

$$n_0({\bf r}) = \langle \Psi_0 \vert {\hat n} \vert \Psi_0 \r... ...{\bf r}_2 , {\bf r}_3 \ldots {\bf r}_N \right) \right\vert^2 .$$ (3.4)

Thus the ground-state $\vert \Psi_0 \rangle$ and density $n_0({\bf r})$ are both functionals of the number of electrons $N$ and the external potential $V_{\mathrm{ext}}({\bf r})$. Density-functional theory, introduced in 1964 by Hohenberg and Kohn [16], makes two remarkable statements.

• The external potential $V_{\mathrm{ext}}({\bf r})$ is uniquely determined by the corresponding ground-state electronic density, to within an additive constant.

  Proof by reductio ad absurdum: assume that a second different external potential $V_{\mathrm{ext}}'({\bf r})$ with ground-state $\vert \Psi_0' \rangle$ gives rise to the same density $n_0({\bf r})$. The ground-state energies are $E_0 = \langle \Psi_0 \vert {\hat H} \vert \Psi_0 \rangle$ and $E_0' = \langle \Psi'_0 \vert {\hat H}' \vert \Psi'_0 \rangle$ where ${\hat H} = {\hat F} + {\hat V}_{\mathrm{ext}}$ and ${\hat H}' = {\hat F} + {\hat V}_{\mathrm{ext}}'$. Taking $\vert \Psi_0' \rangle$ as a trial wave-function for the Hamiltonian ${\hat H}$, we obtain the strict inequality
  

  $$E_0 < \langle \Psi_0' \vert {\hat H} \vert \Psi_0' \rangle = \langle \Psi_0' \vert {\hat H}' \vert \Psi_0' \rangle + \langle \Psi_0' \vert \left( {\hat H} - {\hat H}' \right) \vert \Psi_0' \rangle$$

  $$= E_0' + \int {\mathrm d}{\bf r} ~n_0({\bf r}) \left[ V_{\mathrm{ext}}({\bf r}) - V_{\mathrm{ext}}'({\bf r}) \right] ,$$ (3.5)

  
  
  whereas taking $\vert \Psi_0 \rangle$ as a trial wave-function for ${\hat H}'$ gives
  

  $$E_0' < \langle \Psi_0 \vert {\hat H}' \vert \Psi_0 \rangle = \langle \Psi_0 \vert {\hat H} \vert \Psi_0 \rangle + \langle \Psi_0 \vert \left( {\hat H}' - {\hat H} \right) \vert \Psi_0 \rangle$$

  $$= E_0 - \int {\mathrm d}{\bf r} ~n_0({\bf r}) \left[ V_{\mathrm{ext}}({\bf r}) - V_{\mathrm{ext}}'({\bf r}) \right]$$ (3.6)

  
  
  and adding these two equations together results in the contradiction
  
  $$E_0 + E_0' < E_0 + E_0' .$$
  
  
  Thus, at least in principle, the ground-state density determines (to within a constant) the external potential of the Schrödinger equation of which it is a solution. The external potential and number of electrons $${N = \int {\mathrm d}{\bf r} ~ n_0({\bf r})}$$ determine all the ground-state properties of the system since the Hamiltonian and ground-state wave-function are determined by them.

  So for all densities $n({\bf r})$ which are ground-state densities for some external potential ($v$-representable) the functional $F[n] = \langle \Psi \vert {\hat F} \vert \Psi \rangle$ is unique and well-defined, since $n({\bf r})$ determines the external potential and $N$ (and therefore ${\hat F}$) and thence $\vert \Psi \rangle$. Now a functional for an arbitrary external potential $V({\bf r})$ unrelated to the $V_{\mathrm{ext}}({\bf r})$ determined by $n({\bf r})$ can be defined:
  

  $$E_V[n] = F[n] + \int {\mathrm d}{\bf r} ~V({\bf r}) n({\bf r}) .$$ (3.7)

  
  

• For all $v$-representable densities $n({\bf r})$, $E_V[n] \geq E_0$ where $E_0$ is now the ground-state energy for $N$ electrons in the external potential $V({\bf r})$.

  Proof of this energy variational principle: by the first theorem, a given $n({\bf r})$ determines its own external potential $V_{\mathrm{ext}}({\bf r})$ and ground-state $\vert \Psi \rangle$. If this state is used as a trial state for the Hamiltonian with external potential $V({\bf r})$, we have
  

  $$\langle \Psi \vert {\hat H} \vert \Psi \rangle = \langle \Ps... ...int {\mathrm d}{\bf r}~V({\bf r}) n({\bf r}) = E_V[n] \geq E_0$$ (3.8)

  
  
  by the variational principle. For non-degenerate ground-states, equality only holds if $\vert \Psi \rangle$ is the ground-state for potential $V({\bf r})$.

Thus the problem of solving the Schrödinger equation for non-degenerate ground-states can be recast into a variational problem of minimising the functional $E_V[n]$ with respect to $v$-representable densities. It should be noted that simple counter-examples of $v$-representable densities have been found [17,18,19], but this restriction and the non-degeneracy requirement are overcome by the constrained search formulation.

3.1.2 The constrained search formulation

Following Levy [20,21] we define a functional of the density $n({\bf r})$ for the operator ${\hat F}$ (defined above) as:

$$F[n] = \mathop{\rm min}\limits _{\vert \Psi \rangle \rightarrow n} \langle \Psi \vert {\hat F} \vert \Psi \rangle$$ (3.9)

i.e. the functional takes the minimum value of the expectation value with respect to all states $\vert \Psi \rangle$ which give the density $n({\bf r})$. For a system with external potential $V({\bf r})$ and ground-state $\vert \Psi_0 \rangle$ with energy $E_0$, consider a state $\vert \Psi_{[n]} \rangle$, an $N$-electron state which yields density $n({\bf r})$ and minimises $F[n]$. Define $E_V[n]$ as:

$$E_V[n] = F[n] + \int {\mathrm d}{\bf r}~n({\bf r}) V({\bf r}... ...Psi_{[n]} \vert ({\hat F} + {\hat V}) \vert \Psi_{[n]} \rangle$$ (3.10)

but since ${\hat H} = {\hat F} + {\hat V}$, by the variational principle we obtain

$$E_V[n] \geq E_0$$ (3.11)

with equality only if $\vert \Psi_{[n]} \rangle = \vert \Psi_0 \rangle$. This holds for all densities which can be obtained from an $N$-electron wave-function ($N$-representable). But from the definition of $F[n]$ (3.9) we must also have

$$F[n_0] \leq \langle \Psi_0 \vert {\hat F} \vert \Psi_0 \rangle$$ (3.12)

since $\vert \Psi_0 \rangle$ must be one of states which yields $n_0({\bf r})$. Adding $${\int {\mathrm d}{\bf r}~n_0({\bf r}) V({\bf r})}$$ gives

$$E_V[n_0] \leq E_0$$ (3.13)

which when combined with (3.11) gives the desired result that $E_V[n] \geq E_V[n_0] = E_0$.

Thus the ground-state density $n_0({\bf r})$ minimises the functional $E_V[n]$ and the minimum value is the ground-state electronic energy. Note that the requirement for non-degeneracy of the ground-state has disappeared, and further that instead of considering only $v$-representable densities, we can now consider $N$-representable densities. The requirements of $N$-representability are much weaker and satisfied by any well-behaved density, indeed the only condition [22] is proper differentiability i.e. that the quantity

$$\int {\mathrm d}{\bf r} \left\vert \nabla n^{1 \over 2}({\bf r}) \right\vert^2$$

is real and finite.

3.1.3 Exchange and correlation

The remarkable results of density-functional theory are the existence of the universal functional $F[n]$, which is independent of the external potential, and that instead of dealing with a function of $3 N$ variables (the many-electron wave-function) we can instead deal with a function of only three variables (the density). The complexity of the problem has thus been much reduced, and we note here that this complexity now scales linearly with system-size $N$, so that quantum-mechanical calculations based on density-functional theory can in principle be performed with an effort which scales linearly with system-size.

The exact form of the universal functional $F[n]$ is unknown. The Thomas-Fermi functional [23,24,25]

$$F_{\mathrm{TF}}[n] = \frac{3}{10} \left( 3 \pi^2 \right)^{ 2... ...bf r}) n({\bf r'})}{\left\vert {\bf r} - {\bf r'} \right\vert}$$ (3.14)

can, with hindsight, be viewed as a tentative approximation to this universal functional, but fails to provide even qualitatively correct predictions for systems other than isolated atoms [26,27] although recent, more accurate developments [28,29,30,31,32] have led to the implementation of linear-scaling orbital-free methods for nearly-free electron metals.

The failure to find accurate expressions for the density-functional is a result of the complexity of the many-body problem which is at the heart of the definition of the universal functional. For the electron gas, a system of many interacting particles, the effects of exchange and correlation are crucial to an accurate description of its behaviour. In a non-interacting system, the antisymmetry of the wave-function requires that particles with the same spin occupy distinct orthogonal orbitals, and this results in the particles becoming spatially separated. In an interacting system such as the electron gas in which all the particles repel each other, exchange will thus lead to a lowering of the energy. Moreover, the interactions cause the motion of the particles to become correlated to further reduce the energy of interaction. Thus it is impossible to treat the electrons as independent particles. These effects are completely neglected by the Thomas-Fermi model, and must in part account for its failure, the other source of error being the local approximation for the kinetic energy.

3.1.4 The Kohn-Sham equations

In order to take advantage of the power of DFT without sacrificing accuracy (i.e. including exchange and correlation effects) we follow the method of Kohn and Sham [33] to map the problem of the system of interacting electrons onto a fictitious system of non-interacting ``electrons''. We write the variational problem for the Hohenberg-Kohn density-functional, introducing a Lagrange multiplier $\mu$ to constrain the number of electrons to be $N$:

$$\delta \left[ F[n] + \int {\mathrm d}{\bf r}~V_{\mathrm{ext}... ...( \int {\mathrm d}{\bf r}~n({\bf r}) - N \right) \right] = 0 .$$ (3.15)

Kohn and Sham separated $F[n]$ into three parts

$$F[n] = T_{\mathrm s}[n] + \frac{1}{2} \int {\mathrm d}{\bf r... ...left\vert {\bf r} - {\bf r'} \right\vert} + E_{\mathrm{xc}}[n]$$ (3.16)

in which $T_{\mathrm s}[n]$ is defined as the kinetic energy of a non-interacting gas with density $n({\bf r})$ (not the same as that of the interacting system, although we might hope that the two quantities were of the same order of magnitude), the second term is the classical electrostatic (Hartree) energy and the final term is an implicit definition of the exchange-correlation energy which contains the non-classical electrostatic interaction energy and the difference between the kinetic energies of the interacting and non-interacting systems. The aim of this separation is that the first two terms can be dealt with simply, and the last term, which contains the effects of the complex behaviour, is a small fraction of the total energy and can be approximated surprisingly well.

Using this separation, equation 3.15 can be rewritten:

$$\frac{\delta T_{\mathrm s}[n]}{\delta n({\bf r})} + V_{\mathrm{KS}}({\bf r}) = \mu$$ (3.17)

in which the Kohn-Sham potential $V_{\mathrm{KS}}({\bf r})$ is given by

$$V_{\mathrm{KS}}({\bf r}) = \int {\mathrm d}{\bf r'} \frac{n(... ...\vert } + V_{\mathrm{xc}}({\bf r}) + V_{\mathrm{ext}}({\bf r})$$ (3.18)

and the exchange-correlation potential $V_{\mathrm{xc}}({\bf r})$ is

$$V_{\mathrm{xc}}({\bf r}) = \frac{ \delta E_{\mathrm{xc}}[n] } { \delta n({\bf r}) } .$$ (3.19)

The crucial point to note here is that equation 3.17 is precisely the same equation which would be obtained for a non-interacting system of particles moving in an external potential $V_{\mathrm{KS}}({\bf r})$. To find the ground-state density $n_0({\bf r})$ for this non-interacting system we simply solve the one-electron Schrödinger equations;

$$\left[ -\textstyle{1 \over 2} \nabla^2 + V_{\mathrm{KS}}({\bf r}) \right] \psi_i({\bf r}) = \varepsilon_i \psi_i({\bf r})$$ (3.20)

for ${\textstyle{1 \over 2}N}$ single-particle states^3.1 $\vert \psi_i \rangle$ with energies $\varepsilon_i$, constructing the density from

$$n({\bf r}) = 2 \sum_{i=1}^{N/2} \left\vert \psi_i ({\bf r}) \right\vert^2$$ (3.21)

(the factor 2 is for spin degeneracy - we assume the orbitals are singly-occupied) and the non-interacting kinetic energy $T_{\mathrm s}[n]$ from

$$T_{\mathrm s}[n] = - \sum_{i=1}^{N/2} \int {\mathrm d}{\bf r}~\psi_i^{\ast}({\bf r}) \nabla^2 \psi_i({\bf r}) .$$ (3.22)

Since the Kohn-Sham potential $V_{\mathrm{KS}}({\bf r})$ depends upon the density $n({\bf r})$ it is necessary to solve these equations self-consistently i.e. having made a guess for the form of the density, the Schrödinger equation is solved to obtain a set of orbitals $\{ \psi_i ({\bf r}) \}$ from which a new density is constructed, and the process repeated until the input and output densities are the same. In practice there is no problem converging to the ground-state minimum because of the convex nature of the density-functional [34].

The energy of the non-interacting system, the sum of one-electron eigenvalues, is

$$2 \sum_{i=1}^{N/2} \varepsilon_i = T_{\mathrm s}[n] + \int {\mathrm d}{\bf r}~ n({\bf r}) V_{\mathrm{KS}}({\bf r})$$

$$= T_{\mathrm s}[n] + \int {\mathrm d}{\bf r}~{\mathrm d}{\bf r'}~ \... ...rm{xc}}({\bf r}) + \int {\mathrm d}{\bf r}~n({\bf r}) V_{\mathrm{ext}}({\bf r})$$ (3.23)

which, compared to the interacting system, double-counts the Hartree energy and over-counts the exchange-correlation energy so that the interacting energy is

$$E = 2 \sum_{i=1}^{N/2} \varepsilon_i - \frac{1}{2} \int {\ma... ... r}~n({\bf r}) V_{\mathrm{xc}}({\bf r}) + E_{\mathrm{xc}}[n] .$$ (3.24)

Direct solution of the Schrödinger equation for the extended non-interacting orbitals $\{ \psi_i ({\bf r}) \}$ requires a computational effort which scales as the cube of the system-size $N$, due to the cost of diagonalising the Hamiltonian or orthogonalising the orbitals, whereas the original complexity of finding a minimum of the Hohenberg-Kohn functional only required an effort which scaled linearly with $N$. Thus a linear-scaling method must modify this Kohn-Sham scheme.

3.1.5 The local density approximation

The results so far are exact, provided that the functional form of $E_{\mathrm{xc}}[n]$ is known. The problem of determining the functional form of the universal Hohenberg-Kohn density functional has now been transferred to this one term, and therefore this term is not known exactly. Remarkably, it is possible to make simple approximations for the exchange-correlation energy which work extremely well, and the simplest of these, which is the approximation adopted in this work, is the local density approximation (LDA).

In the LDA, the contribution to the exchange-correlation energy from each infinitesimal volume in space, ${\mathrm d}{\bf r}$, is taken to be the value it would have if the whole of space were filled with a homogeneous electron gas with the same density as is found in ${\mathrm d}{\bf r}$ i.e.

$$E_{\mathrm{xc}}[n] = \int {\mathrm d}{\bf r}~\epsilon_{\mathrm{xc}}\left( n({\bf r})\right) n({\bf r})$$ (3.25)

where $\epsilon_{\mathrm{xc}}\left(n({\bf r})\right)$ is the exchange-correlation energy per electron in a homogeneous electron gas of density $n({\bf r})$. The exchange-correlation potential $V_{\mathrm{xc}}({\bf r})$ then takes the form

$$V_{\mathrm{xc}}({\bf r}) = \frac{\delta E_{\mathrm{xc}}[n]}{... ...c}}\left(n\right)}{{\mathrm d}n} \right\vert _{n=n({\bf r})} .$$ (3.26)

The exchange-correlation energy for the homogeneous electron gas has been calculated by Ceperley and Alder [35] using Monte Carlo methods and in this work we use a parameterisation by Perdew and Zunger [36]. The LDA is exact in the limit of slowly-varying densities, however, the density in systems of interest is generally rapidly varying, and the LDA would appear to be a crude approximation in these cases. Its use is justified a posteriori by its surprising success at predicting physical properties in real systems. This success may be due in part to the fact that the sum rule for the exchange-correlation hole, which must be obeyed by the real functional, is reproduced by the LDA [37]. We can connect the interacting and non-interacting systems using a variable coupling constant $\lambda$ which varies between 0 and 1. We replace the Coulomb interaction by

$$\frac{\lambda}{\left\vert {\bf r} - {\bf r'} \right\vert}$$

and vary $\lambda$ in the presence of an external potential $V_{\lambda}({\bf r})$ so that the ground-state density for all values of $\lambda$ is the same [38]. The Hamiltonian is therefore

$${\hat H}_{\lambda} = -\frac{1}{2} \sum_i \nabla_i^2 + \frac{... ... \right\vert} + {\hat V}_{\mathrm{ext}} + {\hat V}_{\lambda} .$$ (3.27)

The exchange-correlation hole $n_{\mathrm{xc}}({\bf r},{\bf r'})$ is then defined in terms of a coupling-constant integration of the pair correlation function $g({\bf r},{\bf r'};\lambda)$ of the system with density $n({\bf r})$ and scaled Coulomb interaction [39,40];

$$n_{\mathrm{xc}}({\bf r},{\bf r'}) = n({\bf r'}) \int_0^1 {\mathrm d}\lambda \left[ g({\bf r},{\bf r'};\lambda) - 1 \right] .$$ (3.28)

The exchange-correlation energy can then be expressed in the form of a classical electrostatic interaction between the density $n({\bf r})$ and the hole density $n_{\mathrm{xc}}({\bf r},{\bf r'})$;

$$E_{\mathrm{xc}}[n] = \frac{1}{2} \int {\mathrm d}{\bf r}~{\m... ...bf r},{\bf r'})} {\left\vert {\bf r} - {\bf r'} \right\vert} .$$ (3.29)

The sum rule follows from the definition of the pair correlation function [41]

$$\int {\mathrm d}{\bf r'}~n_{\mathrm{xc}}({\bf r},{\bf r'}) = -1 ,$$ (3.30)

which is interpreted by saying that the exchange-correlation hole excludes one electron as expected. It can also be shown that the exchange-correlation energy depends only weakly on the detailed shape of the exchange-correlation hole [42], and these two facts account, at least in part, for the success of the LDA. This view is supported by the fact that improvements to the LDA involving gradient expansions show no consistent improvement unless they enforce the sum rule obeyed by the LDA [43,44].