5.2 Kohn-Sham Density Functional Theory

The density functional theory by Hohenberg, Kohn, and Sham[Hohenberg and Kohn(1964), Kohn and Sham(1965)] stems from earlier work by Dirac,[Dirac(1930)] who showed that the exchange energy of a uniform electron gas can be computed exactly from the charge density along. However, while this traditional density functional approach, nowadays called “orbital-free” DFT, makes a direct connection to the density alone, in practice it is constitutes a direct approach where the necessary equations contain only the electron density, difficult to obtain decent approximations for the kinetic energy functional. Kohn and Sham sidestepped this difficulty via an indirect approach in which the kinetic energy is computed exactly for a noninteracting reference system, namely, the Kohn-Sham determinant.[Kohn and Sham(1965)] It is the Kohn-Sham approach that first made DFT into a practical tool for calculations.

Within the Kohn-Sham formalism,[Kohn and Sham(1965)] the ground state electronic energy, $E$ , can be written as

$\begin{equation} \label{eq:DFT_ energy_ components} E = E_{\text {T}}^{} + E_{\text {V}}^{} + E_{\text {J}}^{} + E_{\text {XC}}^{} \end{equation}$

(5.1)

where $E_{\text {T}}^{}$ is the kinetic energy, $E_{\mathrm{V}}^{}$ is the electron–nuclear interaction energy, $E_{\text {J}}^{}$ is the Coulomb self-interaction of the electron density, $\rho (\ensuremath{\mathbf{r}})$ and $E_{\text {XC}}^{}$ is the exchange-correlation energy. Adopting an unrestricted format, the $\alpha$ and $\beta$ total electron densities can be written as

$\displaystyle \label{eq429} \begin{aligned} \rho _{\alpha }(\ensuremath{\mathbf{r}}) & = \sum _{i=1}^{n_\alpha } |\psi _ i^{\alpha }|^2 \\ \rho _{\beta }(\ensuremath{\mathbf{r}}) & = \sum _{i=1}^{n_\beta } |\psi _ i^{\beta } |^2 \end{aligned}$

(5.2)

where $n_{\alpha }$ and $n_{\beta }$ are the number of alpha and beta electron respectively, and $\psi _ i$ are the Kohn-Sham orbitals. Thus, the total electron density is

$\begin{equation} \label{eq430} \rho (\ensuremath{\mathbf{r}}) =\rho _\alpha (\ensuremath{\mathbf{r}}) + \rho _\beta (\ensuremath{\mathbf{r}}) \end{equation}$

(5.5)

Within a finite basis set, the density is represented by[Pople et al.(1992)Pople, Gill, and Johnson]

$\begin{equation} \label{eq431} \rho (\ensuremath{\mathbf{r}}) = \sum _{\mu \nu } P_{\mu \nu } \phi _\mu (\ensuremath{\mathbf{r}}) \phi _\nu (\ensuremath{\mathbf{r}}) \; , \end{equation}$

(5.6)

where the $P_{\mu \nu }$ are the elements of the one-electron density matrix; see Eq. eq:P(mu,nu) in the discussion of Hartree-Fock theory. The various energy components in Eq. eq:DFT_energy_components can now be written

$\displaystyle \label{eq432} E_{\mathrm{T}}^{}$	$\displaystyle =$	$\displaystyle \sum _{i=1}^{n_\alpha } \left\langle \psi _ i^\alpha \left\|-\frac{1}{2}\hat{\nabla }^2\right\|\psi _ i^{\alpha } \right\rangle + \sum _{i=1}^{n_\beta } \left\langle \psi _ i^\beta \left\|-\frac{1}{2}\hat{\nabla }^2\right\|\psi _ i^{\beta } \right\rangle$	(5.7)
$\displaystyle$	$\displaystyle =$	$\displaystyle \sum _{\mu \nu } P_{\mu \nu } \left\langle \phi _\mu (\ensuremath{\mathbf{r}}) \left\|-\frac{1}{2}\hat{\nabla }^2\right\| \phi _\nu (\ensuremath{\mathbf{r}}) \right\rangle$	(5.8)
$\displaystyle E_{\mathrm{V}}^{}$	$\displaystyle =$	$\displaystyle -\sum _{A=1}^ M Z_ A\int \frac{\rho (\ensuremath{\mathbf{r}})}{\|\ensuremath{\mathbf{r}}-\ensuremath{\mathbf{R}}_ A\|} d\ensuremath{\mathbf{r}}$	(5.9)
$\displaystyle$	$\displaystyle =$	$\displaystyle -\sum _{\mu \nu } P_{\mu \nu } \sum _ A \left\langle \phi _\mu (\ensuremath{\mathbf{r}}) \left\| \frac{Z_ A}{\|\ensuremath{\mathbf{r}}-\ensuremath{\mathbf{R}}_ A\|} \right\| \phi _\nu (\ensuremath{\mathbf{r}}) \right\rangle$	(5.10)
$\displaystyle E_{\mathrm{J}}^{}$	$\displaystyle =$	$\displaystyle \frac{1}{2}\left\langle \rho (\ensuremath{\mathbf{r}}_1) \left\| \frac{1}{\|\ensuremath{\mathbf{r}}_1 - \ensuremath{\mathbf{r}}_2\|} \right\| \rho (\ensuremath{\mathbf{r}}_2) \right\rangle$	(5.11)
$\displaystyle$	$\displaystyle =$	$\displaystyle \frac{1}{2}\sum _{\mu \nu } \sum _{\lambda \sigma } P_{\mu \nu } P_{\lambda \sigma } \left(\mu \nu \vert \lambda \sigma \right)$	(5.12)
$\displaystyle \label{E_ XC} E_{\mathrm{XC}}^{}$	$\displaystyle =$	$\displaystyle \int f\bigl [\rho (\ensuremath{\mathbf{r}}),\hat{\bm@general \boldmath \m@ne \mv@bold \bm@command {\nabla }}\rho (\ensuremath{\mathbf{r}}),\ldots \bigr ] \, \rho (\mathbf{r}) \; d\ensuremath{\mathbf{r}} \; .$	(5.13)

Minimizing $E$ with respect to the unknown Kohn-Sham orbital coefficients yields a set of matrix equations exactly analogous to Pople-Nesbet equations of the UHF case, Eq. eq:Pople-Nesbet, but with modified Fock matrix elements [cf. Eq. eq:F(mu,nu)]

$\displaystyle \label{eq437} \begin{aligned} F_{\mu \nu }^\alpha & = H_{\mu \nu }^{\mathrm{core}} + J_{\mu \nu } - F_{\mu \nu }^{\mathrm{XC}\alpha } \\ F_{\mu \nu }^\beta & = H_{\mu \nu }^{\mathrm{core}} + J_{\mu \nu } - F_{\mu \nu }^{\mathrm{XC}\beta } \; . \end{aligned}$

(5.14)

Here, $\ensuremath{\mathbf{F}}^{\mathrm{XC}\alpha }$ and $\ensuremath{\mathbf{F}}^{\mathrm{XC}\beta }$ are the exchange-correlation parts of the Fock matrices and depend on the exchange-correlation functional used. UHF theory is recovered as a special case simply by taking $F_{\mu \nu }^{\mathrm{XC}\alpha } = K_{\mu \nu }^\alpha$ , and similarly for $\beta$ . Thus, the density and energy are obtained in a manner analogous to that for the HF method. Initial guesses are made for the MO coefficients and an iterative process is applied until self-consistency is achieved.