Kohn–Sham method

Kohn–Sham equation

Early attempts to employ density functional theory for electronic structure investigations of atoms led to relatively poor results due to inaccurate approximations of the kinetic energy term in the functional \(F\). In 1965, Kohn and Sham [KS65] proposed a solution to this problem by introducing the concept of Kohn–Sham (KS) orbitals and evaluating the kinetic energy according to

\[\begin{equation*} T_s[n(\mathbf{r})] = - \frac{\hbar^2}{2 m_\mathrm{e}} \sum_i \int \psi_i^*(\mathbf{r}) \nabla^2 \psi_i(\mathbf{r}) \, \mathrm{d}^3\mathbf{r} \end{equation*}\]

where the electron density is defined as

\[\begin{equation*} n(\mathbf{r}) = \sum_{i} | \psi_i(\mathbf{r})|^2 \end{equation*}\]

Underlying this idea, there is a representation of the electronic ground state in terms of a Slater determinant known as the Kohn–Sham reference state

\[\begin{equation*} | \Psi_s \rangle = \frac{1}{\sqrt{N!}} \det \big[\psi_1 \psi_2\cdots\psi_N\big] \end{equation*}\]

which ensures that the density is \(N\)-representable.

To within a small residual correction, the functional \(T_s\) amounts to the true kinetic energy functional \(T\). Correspondingly, the energy functional can be written

\[\begin{equation*} E[n(\mathbf{r})] = T_s[n(\mathbf{r})] + J[n(\mathbf{r})] + E_\mathrm{xc}[n(\mathbf{r})] + \int v(\mathbf{r}) n(\mathbf{r}) \, \mathrm{d}^3\mathbf{r} \end{equation*}\]

where the electron–electron interaction contribution is partitioned into the classical Coulomb interaction energy

\[ J[n(\mathbf{r})] = \frac{1}{2} \frac{e^2}{4\pi\varepsilon_0} \iint \frac{n(\mathbf{r}) n(\mathbf{r}')}{|\mathbf{r} - \mathbf{r}'|} \, \mathrm{d}^3\mathbf{r} \, \mathrm{d}^3\mathbf{r}' \]

and the exchange–correlation (xc) functional, \(E_\mathrm{xc}[n(\mathbf{r})]\), accounting for the non-classical electron–electron interaction and the residual kinetic energy correction, \((T[n(\mathbf{r})] - T_s[n(\mathbf{r})])\).

Within the Kohn–Sham formulation of DFT, the act of minimizing the energy under the constraint of particle conservation closely follows the derivation of the Hartree–Fock equation and the resulting equation from which the best orbitals are found is known as the Kohn–Sham equation. In its canonical form, it reads

\[\begin{equation*} \label{eq:ks} \hat f \psi_i(\mathbf{r}) = \left[ -\frac{\hbar^2}{2 m_\mathrm{e}} \nabla^2 + v_\mathrm{eff}(\mathbf{r}) \right] \psi_i(\mathbf{r}) = \varepsilon_i \psi_i(\mathbf{r}) \end{equation*}\]

where \(\varepsilon_i\) is the KS orbital energy and \(v_\mathrm{eff}(\mathbf{r})\) is the effective potential equal to

\[ v_\mathrm{eff}(\mathbf{r}) = \frac{e^2}{4\pi\varepsilon_0} \int \frac{n(\mathbf{r'})}{|\mathbf{r} - \mathbf{r}'|} \, \mathrm{d}^3\mathbf{r}' + v_\mathrm{xc}(\mathbf{r}) + v(\mathbf{r}) \]

with

\[ v_\mathrm{xc}(\mathbf{r}) = \frac{\delta E_{\text{xc}}}{\delta n(\mathbf{r})} \]

Exchange–correlation functionals

The accuracy of the Kohn-Sham method is determined by the capability of the xc-functional to capture the non-classical electron–electron interactions in the molecular system. As the exact form of xc-functional is unknown, one introduces physically justified approximations. The development of new xc-functionals is an on-going process with new functionals being developed to address the inaccuracies of previous generations.

The xc-functionals can be classified by various criteria such as functional parameterization, arguments, etc. Here, we will describe four categories of functionals sorted by argument types and highlight the most prominent representatives of each group. We will assume the xc-functional to be written in terms on an energy density \(\varepsilon_\mathrm{xc}\) on the form

\[ E_\mathrm{xc}[n(\mathbf{r})] = \int n(\mathbf{r}) \, \varepsilon_\mathrm{xc}(n, \nabla n, \nabla^2 n) \, d^3\mathbf{r} \]

where it is implicitly understood that the electron density can be decomposed into \(\alpha\)- and \(\beta\)-spin densities that in general may not be equal.

Fig. 1 Jacob’s ladder of density functionals.

Local density approximation (LDA)

Functionals of this group depend solely on \(n(\mathbf{r})\). In general, they provide a quite poor description of electron-electron interactions in regions of rapidly varying density, like the bonding regions in molecules. Nevertheless, until the late 1980s, almost all DFT calculations of the electronic structure of molecules were carried out using such functionals due to a lack of alternatives. Within LDA, the combination of the Dirac exchange [Dir30] and the Vosko-Wilk-Nussair correlation [VWN80] functionals provides the most accurate results in a wide range of applications and it has become recognized as the de facto standard LDA functional.

In the LDA, the exchange–correlation potential takes the form

\[ v_\mathrm{xc}(\mathbf{r}) = \varepsilon_\mathrm{xc}(n(\mathbf{r})) + n(\mathbf{r}) \frac{\partial \varepsilon_\mathrm{xc}}{\partial n} \]

Generalized gradient approximation (GGA)

Functionals of this group depend on \(n(\mathbf{r})\) and \(\nabla n(\mathbf{r})\). The inclusion of the electron density gradient enables a more accurate description of electron-electron interactions in regions of rapidly varying density. The first attempts to develop xc.functionals of this type were undertaken by Sham [Sha71] and Herman and co-workers [HVDO69] in the 1970s but from a performance point of view they were unsuccessful. Only in 1988, the first reliable exchange functional (B88) of this type was proposed by Becke [Bec88]. Adopting the Becke formulation, several other xc-functionals have since been developed. Among these, most notable are the Lee-Yang-Parr (LYP) correlation functional [LYP88] and the Perdew-Burke-Ernzerhof (PBE) xc-functional [ES99]. The former LYP functional is an integral component of the most widely used hybrid functional (B3LYP) [Bec93, SDCF94], and the latter PBE functional is the most popular GGA functional in solid-state calculations.

Meta generalized gradient approximation (mGGA)

Functionals of this group depend on \(n(\mathbf{r})\), \(\nabla n(\mathbf{r})\), and higher-order derivatives of the electron or/and kinetic energy densities. The development of xc-functionals of this type has been actively pursued only during the last two decades and after Tao and co-workers introduced the first reliable mGGA functional in 2003 namely Tao-Perdew-Staroverov-Scuseria (TPSS) [TPSS03]. Among the mGGA xc-functionals, the functionals developed by the Truhlar group are the most widely used [ZT11], and the revised M06-L functional [WJY+17] offers accuracy comparable to hybrid xc-functional in many situations.

Hybrid exchange-correlation functionals

Functionals of this group explicitly incorporate a fraction (\(0 < c_{\text{x}} < 1\)) of the exact (“Hartree–Fock”) exchange in the energy expression. Since the introduction by Becke [Bec93], these functionals dominate in calculations of molecular electronic structure. In particular, the B3LYP functional, which combines Dirac, Becke exchange functionals with VWN, LYP correlation functionals via the three-parameters scheme [Bec93], has become the gold standard for DFT calculations of electronic structure and properties of molecules in their electronic ground states. Besides the B3LYP functional, the various hybrid functionals from the Minnesota functional family [ZT11] have reached considerable popularity.

Range separated exchange-correlation functionals

Functionals of this group explicitly partition the xc-contribution to the energy functional into short- and long-range parts. Typically, the same partitioning is also applied to the classical Coulomb term to ensure correct asymptotic behavior of \((J[n(\mathbf{r})] + E_{\text{xc}}[n(\mathbf{r})])\). The functionals of this type are frequently employed in studies of excited states and optical properties of molecules in connection with time-dependent DFT methods. Among these functionals, the CAM-B3LYP xc-functional [YTH04] is the most popular one, and frequently provides a valued alternative to ab initio methods in studies of molecular systems with charge-transfer character.

Known deficiencies

After introducing the most popular classes of xc-functionals, it is important to mention some of their known deficiencies.

First, we recognize the inherent inability of standard xc-functionals to describe van der Waals interactions. This issue has been addressed in recent years by adding empirical dispersion corrections or by incorporating dispersion corrections via a kinetic-energy density dependent contribution. The former strategy has been adopted by Grimme and co-workers [CEH+19, GAEK10], who introduced empirical D3 and D4 models of dispersion interactions and successfully combined these models with several xc-functionals. As a result, the B3LYP+D3 xc-functional became a new standard for structure optimizations of both weakly and strongly bonded molecular systems. Truhlar and co-workers have adopted the latter strategy to develop the Minnesota functional family [WJY+17, ZT11], and these functionals provide a viable alternative to the D3- or D4-corrected standard functionals.

Second, we draw attention to the fact that there is an amount of self interaction built into the treatment of the electrons as due to the non-exact form of the xc-functional. In Hartree–Fock theory, there is an exact cancellation of the self-interaction terms in between the Coulomb and exchange contributions to the energy, but in DFT this symmetry is broken.