An Overview of Self-Consistent Field Calculations Within Finite Basis Sets ^†

Abstract

A uniform derivation of the self-consistent field equations in a finite basis set is presented. Both restricted and unrestricted Hartree–Fock (HF) theory as well as various density functional approximations are considered. The unitary invariance of the HF and density functional models is discussed, paving the way for the use of localized molecular orbitals. The self-consistent field equations are derived in a non-orthogonal basis set, and their solution is discussed also in the presence of linear dependencies in the basis. It is argued why iterative diagonalization of the Kohn–Sham–Fock matrix leads to the minimization of the total energy. Alternative methods for the solution of the self-consistent field equations via direct minimization as well as stability analysis are briefly discussed. Explicit expressions are given for the contributions to the Kohn–Sham–Fock matrix up to meta-GGA functionals. Range-separated hybrids and non-local correlation functionals are summarily reviewed.

1. Introduction

Electronic structure calculations have become a cornerstone of modern-day research in chemistry and materials physics, allowing in silico modeling of chemical reactions and the first principles design of novel catalysts [1]. Electronic structure calculations on molecular systems most often employ the linear combination of atomic orbitals (LCAO) approach, where the molecular orbitals (MOs) are expanded in terms of atomic orbitals (AOs). Several possible alternatives for the form of the AOs are commonly used—Gaussian-type orbitals (GTOs), Slater-type orbitals (STOs), as well as numerical atomic orbitals (NAOs); see [2] for details. LCAO electronic structure calculations involve a variational minimization of the total energy with respect to the AO expansion coefficients of the MOs. Importantly, the formalism used in the LCAO approach is not restricted to AOs which are atom-centered basis functions; it can also be used, e.g., in combination with fully numerical basis functions as in the finite element approach, as has been recently demonstrated in [3,4]. Once the energy has been minimized and the corresponding wave function has been obtained, it is possible to compute a number of properties either directly from the electronic wave function (e.g., electron densities, orbital energies, molecular dipole moment), or from its response to external perturbations (nuclear magnetic shieldings, vibrational frequencies, etc.).

The mathematical foundations for spin-restricted Hartree–Fock (HF) theory within the LCAO approach were laid out independently by Roothaan and Hall [5,6]. In their seminal papers, Roothaan and Hall derived matrix equations that can be conveniently implemented on a computer as an iterative procedure. As will be seen later in Section 6, the Roothaan–Hall equations turn out to yield a generalized eigenvalue problem $\mathbf{F}\mathbf{C}=\mathbf{S}\mathbf{C}\mathbf{E}$ in the non-orthogonal AO basis set, which had been solved some years before by Löwdin in the context of Heitler–London theory [7].

Subsequently to the work by Roothaan and Hall, Pople and Nesbet [8] and Berthier [9] independently published the corresponding equations for an unrestricted (open-shell) HF description by an analogous scheme, without providing an explicit derivation. The Pople–Nesbet–Berthier equations assume a form similar to the Roothaan–Hall equations—constituting a coupled set of general eigenvalue equations—as will also be seen later on in the manuscript (Section 6). Restricted open-shell HF was then described by Roothaan [10]; restricted open-shell calculations will not be considered in the present work as they have been extensively reviewed by Krebs in [11] to which we refer for further details.

Density functional theory [12,13] (DFT; see also [14,15]) became popular in chemistry through the efforts of Pople and coworkers in making the method generally available to quantum chemists [16] and showing that atomization energies from DFT may agree well with experiment [17,18]. Also DFT turns out to yield self-consistent field (SCF) equations that assume the same form as in HF but with a different expression for the Fock matrix $\mathbf{F}$. Pople and coworkers reported the equations necessary for solving SCF for DFT in the LCAO context up to generalized gradient approximation (GGA) functionals in [16]; an analogous derivation was also presented by Kobayashi et al in [19]. The self-consistent implementation of meta-GGA functionals was later described by Neumann, Nobes and Handy in [20]. Density functional calculations sometimes include also non-local correlation contributions; self-consistent LCAO implementations thereof have been reported by Vydrov and coworkers [21,22,23,24].

Despite the progress in and widespread success of DFT, to our knowledge, a uniform derivation of the SCF equations for HF and DFT including all the necessary expressions for the elements of the Kohn–Sham–Fock matrix up to the level of meta-GGA functionals has, up to now, not been explicitly published in the literature. This has likely contributed to the lack of complete support for meta-GGA functionals in popular quantum chemistry programs; for instance, PSI4 [25] and PySCF [26] lack support for meta-GGAs that depend on the Laplacian of the density such as the Becke–Roussel exchange functional [27], for example. This paper, therefore, presents such a derivation, yielding expressions of the DFT contributions to the Kohn–Sham–Fock matrix up to the level of meta-GGA functionals in a consistent way, facilitating the implementation of DFT in new programs.

The present derivation also has an obvious educational value. Indeed, in what follows, HF and various flavors of DFT belonging to different rungs of Jacob’s Ladder [28]—the local spin density approximation (LDA), the GGA and meta-GGA approximations—will be explicitly described in a uniform notation, making the similarities and dissimilarities between the approaches crystal clear. Facilitated by the uniform derivation, we will discuss key issues and features in the HF and DFT methodologies that arise from the mathematical formulation.

First, the basis set expansion of the molecular orbitals and the electron density is written out in Section 2. Then, the energy expression for HF and DFT is presented in Section 3, with a brief explanation of their physical content. The HF and DFT energy is shown to be invariant to rotations of the occupied and of the virtual orbitals in Section 4, allowing the construction of localized orbitals. The possibilities and drawbacks of spin-restricted calculations are discussed in Section 5. The finite-basis SCF equations are derived as generalized eigenvalue equations in Section 6. It is shown that the general eigenvalue equations can be reduced into normal eigenvalue equations by a transformation to an orthonormal basis in Section 7, and that linear dependencies in the basis can be eliminated on the way. The reason why the solution of the SCF equations amounts to a minimization of the total energy is rationalized in Section 8. Direct minimization methods are briefly introduced and stability analysis discussed in Section 9. The SCF method and direct minimization are contrasted in Section 10. Finally, the contributions to the Kohn–Sham–Fock matrix arising from various-rung DFT functionals are listed in Section 11. The article concludes with a brief summary and discussion in Section 12. Atomic units are used throughout the text.

2. Basis Set Expansion

In the HF and DFT approaches, the electronic wave function is written as a Slater determinant, in which the electrons occupy a set of MOs $\phi \left(\overrightarrow{r}\right)$. The MOs are expanded in terms of normalized expansion functions $\chi \left(\overrightarrow{r}\right)$, which are typically AOs. Both $\phi$ and $\chi$, as well as the LCAO coefficients $C_{\mu i}$, are typically chosen to be real in the lack of magnetic interactions that would generally make the Hamiltonian operator complex. Note, however, that the use of complex coefficients has been shown to be sometimes beneficial to describe challenging systems even in the lack of magnetic fields in Hartree–Fock (as reviewed in [29]) or DFT (as shown in [30]); complex instabilities may also arise in specialized methods beyond the SCF level, see [31] for instance. The expansion functions are generally not orthonormal

$$\int \mathrm{d}\overrightarrow{r}{\chi }_{\mu }\left(\overrightarrow{r}\right){\chi }_{\nu }\left(\overrightarrow{r}\right)=S_{\mu \nu }\ne {\delta }_{\mu \nu }$$

where ${\delta }_{\mu \nu }$ is the Kronecker delta: ${\delta }_{\mu \nu }=1$ if $\mu =\nu$ and ${\delta }_{\mu \nu }=0$ otherwise. Greek letters, $\mu ,\nu ,\lambda ,\sigma ,\eta$, $\zeta$ and $\theta$ will be used to identify the expansion functions $\chi \left(\overrightarrow{r}\right)$, whereas Roman letters, i, j and k will be used to identify the MOs $\phi$. The $\alpha$ (spin-up) and $\beta$ (spin-down) MOs are expanded separately as

$$\begin{array}{c}\hfill {\phi }_i^{\alpha }\left(\overrightarrow{r}\right)=\sum _{\mu =1}^M{C}_{\mu i}^{\alpha }{\chi }_{\mu }\left(\overrightarrow{r}\right)\end{array}$$

(1)

$$\begin{array}{c}\hfill {\phi }_i^{\beta }\left(\overrightarrow{r}\right)=\sum _{\mu =1}^M{C}_{\mu i}^{\beta }{\chi }_{\mu }\left(\overrightarrow{r}\right)\end{array}$$

(2)

Both the $\alpha$ and $\beta$ MOs are orthonormal to themselves

$$\int \mathrm{d}\overrightarrow{r}{\phi }_i^{\alpha }\left(\overrightarrow{r}\right){\phi }_j^{\alpha }\left(\overrightarrow{r}\right)={\delta }_{ij}\mathrm{and}\int \mathrm{d}\overrightarrow{r}{\phi }_i^{\beta }\left(\overrightarrow{r}\right){\phi }_j^{\beta }\left(\overrightarrow{r}\right)={\delta }_{ij}$$

or equivalently within the basis set

$${\left({\mathbf{C}}^{\alpha }\right)}^{\mathrm{T}}\mathbf{S}{\mathbf{C}}^{\alpha }=\mathbf{1}\mathrm{and}{\left({\mathbf{C}}^{\beta }\right)}^{\mathrm{T}}\mathbf{S}{\mathbf{C}}^{\beta }=\mathbf{1}$$

However, the $\alpha$ orbitals are generally not orthonormal to the $\beta$ orbitals:

$$\int \mathrm{d}\overrightarrow{r}{\phi }_i^{\alpha }\left(\overrightarrow{r}\right){\phi }_j^{\beta }\left(\overrightarrow{r}\right)\ne {\delta }_{ij}$$

or

$${\left({\mathbf{C}}^{\alpha }\right)}^{\mathrm{T}}\mathbf{S}{\mathbf{C}}^{\beta }\ne \mathbf{1}$$

(3)

The electron density plays a pivotal role in quantum chemistry. In line with chemistry literature, $\rho \left(\overrightarrow{r}\right)$ will be used to denote the electron density at the point $\overrightarrow{r}$ in contrast to the physics notation $n\left(\overrightarrow{r}\right)$ which is customary in the DFT literature. The total electron density is formed from the $\alpha$ and $\beta$ densities, ${\rho }_{\alpha }$ and ${\rho }_{\beta }$, as $\rho \left(\overrightarrow{r}\right)={\rho }_{\alpha }\left(\overrightarrow{r}\right)+{\rho }_{\beta }\left(\overrightarrow{r}\right)$. The spin-$\sigma$ electron density can be evaluated as

$${\rho }_{\sigma }\left(\overrightarrow{r}\right)=\sum _{i=1}^{N_{\sigma }}{\left|{\phi }_i^{\sigma }\left(\overrightarrow{r}\right)\right|}^2=\sum _{i=1}^{N_{\sigma }}\sum _{\mu \nu }{C}_{\mu i}^{\sigma }{C}_{\nu i}^{\sigma }{\chi }_{\mu }\left(\overrightarrow{r}\right){\chi }_{\nu }\left(\overrightarrow{r}\right)=\sum _{\mu \nu }{P}_{\mu \nu }^{\sigma }{\chi }_{\mu }\left(\overrightarrow{r}\right){\chi }_{\nu }\left(\overrightarrow{r}\right)$$

(4)

in which $N_{\sigma }$ is the number spin-$\sigma$ electrons in the system, and the sums over basis functions ${\sum }_{\mu }$ and ${\sum }_{\nu }$ indicate ${\sum }_{\mu =1}^M$ and ${\sum }_{\nu =1}^M$, respectively; this convention for easier readability of sums over basis functions is used throughout the rest of this work.

The density matrix ${\mathbf{P}}^{\sigma }$ has been defined in Equation (4) as

$$P_{\mu \nu }^{\sigma }=\sum _{i=1}^{N_{\sigma }}{C}_{\mu i}^{\sigma }{C}_{\nu i}^{\sigma }$$

(5)

As is evident from the form of Equation (5), the density matrices are symmetric, $P_{\mu \nu }^{\sigma }=P_{\nu \mu }^{\sigma }$. As was already mentioned above, the total electron density is obtained from the sum of the $\alpha$ and $\beta$ densities. Correspondingly, a total density matrix is given by

$$\mathbf{P}={\mathbf{P}}^{\alpha }+{\mathbf{P}}^{\beta }$$

(6)

from which the total density can be evaluated using a relation analogous to Equation (4).

3. Energy Expression

The starting point for the derivation is the non-relativistic energy expression [5,8,13,16],

$$E=\sum _{\mu \nu }{P}_{\mu \nu }{H}_{\mu \nu }+\frac{1}{2}\sum _{\mu \nu \lambda \sigma }{P}_{\mu \nu }{P}_{\lambda \sigma }\left(\mu \nu \right|\lambda \sigma )-\frac{a}{2}\sum _{\mu \nu \lambda \sigma }(P_{\mu \lambda }^{\alpha }{P}_{\nu \sigma }^{\alpha }+P_{\mu \lambda }^{\beta }{P}_{\nu \sigma }^{\beta })\left(\mu \nu \right|\lambda \sigma )+b\int f\left(\overrightarrow{r}\right)\mathrm{d}\overrightarrow{r}$$

(7)

where the electron repulsion integral $\left(\mu \nu \right|\lambda \sigma )$ is defined as

$$\left(\mu \nu \right|\lambda \sigma )=\int \mathrm{d}{\overrightarrow{r}}_1\int \mathrm{d}{\overrightarrow{r}}_2{\chi }_{\mu }\left(\overrightarrow{r_1}\right){\chi }_{\nu }\left(\overrightarrow{r_1}\right)\frac{1}{r_{12}}{\chi }_{\lambda }\left(\overrightarrow{r_2}\right){\chi }_{\sigma }\left(\overrightarrow{r_2}\right)$$

(8)

and a and b are constants that define the fraction of HF exchange and the weight of the density functional approximation, respectively. The choice $a=1$ and $b=0$ corresponds to HF, whereas $a=0$ and $b=1$ yields a “pure” density functional without exact exchange such as the Perdew–Burke–Ernzerhof functional [32]. The choice $a\ne 0$ and $b\ne 0$ is the most general one, which corresponds to a hybrid functional [33] that are popular in quantum chemistry; perhaps the most famous example being the historical B3LYP functional [34].

The first term in Equation (7), which will be referred to as $E^H$, describes the kinetic energy of the electrons and the Coulombic attraction of the N nuclei in the system, with the matrix elements

$$H_{\mu \nu }=\int dr{\chi }_{\mu }\left(\overrightarrow{r}\right)\left(-\frac{1}{2}{\nabla }^2+\sum _N\frac{Z_N}{r_N}\right){\chi }_{\nu }\left(\overrightarrow{r}\right)$$

(9)

The one-electron operator in Equation (9) is commonly known as the core Hamiltonian, and the resulting $E^H$ is the dominating contribution to the total energy.

However, the core Hamiltonian lacks electronic interactions. These are described by the second and third terms in Equation (7), which describe the classical Coulomb and the quantum mechanical “exchange” energy, and are referred to as $E^J$ and $E^K$, respectively. The $E^J$ contribution to the total energy can be straightforwardly derived from the expression for the Coulomb repulsion between the electrons described by the electron density $\rho \left(\overrightarrow{r}\right)$

$$E^J=\frac{1}{2}\int d\overrightarrow{r_1}\int d\overrightarrow{r_2}\rho \left(\overrightarrow{r_1}\right)\frac{1}{r_{12}}\rho \left(\overrightarrow{r_2}\right)$$

whereas the expression for the exchange energy contribution $E^K$ can be obtained, for instance, using Slater’s rules for a HF wave function ($a=1$).

The final term in Equation (7), referred to as $E^{\mathrm{XC}}$, describes the DFT exchange-correlation contributions which, alike $E^J$ and $E^K$, arise from electronic interactions. The exchange-correlation term is commonly written as

$$E^{\mathrm{XC}}=\int \mathrm{d}\overrightarrow{r}f\left(\overrightarrow{r}\right)=\int \mathrm{d}\overrightarrow{r}\rho \left(\overrightarrow{r}\right){ϵ}^{\mathrm{XC}}\left(\overrightarrow{r}\right)$$

where ${ϵ}^{\mathrm{XC}}$ is the exchange-correlation energy density per electron. Usually $f\left(\overrightarrow{r}\right)$ is a function of the electron density $\rho \left(\overrightarrow{r}\right)$; it may also depend on the derivatives of $\rho \left(\overrightarrow{r}\right)$ and the kinetic energy density $\tau \left(\overrightarrow{r}\right)$, depending on which rung of Jacob’s Ladder [28] is used to the describe the exchange-correlation effects. The various rungs are discussed in Section 11.

4. Unitary Invariance

The ${\mathbf{P}}^{\alpha }$ and ${\mathbf{P}}^{\beta }$ matrices turn out to be invariant to rotations of the occupied orbitals among themselves. Rotating the occupied subset of the molecular orbitals $\phi$ by a orthogonal matrix $\mathbf{U}$ defines a new set of occupied orbitals

$${\phi }_i^{{\alpha }^{\prime }}=\sum _{k=1}^{N_{\alpha }}{\phi }_k^{\alpha }{U}_{ki}$$

the MO coefficients of which can be obtained as

$$C_{\mu i}^{{\alpha }^{\prime }}=\sum _{k=1}^{N_{\alpha }}{C}_{\mu k}^{\alpha }{U}_{ki}$$

This can also be written in matrix notation as

$${\mathbf{C}}^{{\alpha }^{\prime }}={\mathbf{C}}^{\alpha }\mathbf{U}\mathrm{or}{\mathbf{C}}^{{\alpha }^{\prime }}{\mathbf{U}}^{\mathrm{T}}={\mathbf{C}}^{\alpha }$$

(10)

The invariance to rotations in the occupied-occupied block is easy to prove, as

$$P_{\mu \nu }^{{\alpha }^{\prime }}=\sum _{i=1}^{N_{\alpha }}{C}_{\mu i}^{{\alpha }^{\prime }}{C}_{\nu i}^{{\alpha }^{\prime }}=\sum _{ikl=1}^{N_{\alpha }}{U}_{ik}{C}_{\mu k}^{\alpha }{U}_{il}{C}_{\nu l}^{\alpha }=\sum _{kl=1}^{N_{\alpha }}{\delta }_{kl}{C}_{\mu k}^{\alpha }{C}_{\nu l}^{\alpha }=\sum _{k=1}^{N_{\alpha }}{C}_{\mu k}^{\alpha }{C}_{\nu k}^{\alpha }=P_{\mu \nu }^{\alpha }$$

(11)

where we have used the orthogonality of $\mathbf{U}$, ${\mathbf{U}}^{\mathrm{T}}\mathbf{U}=\mathbf{1}=\mathbf{U}{\mathbf{U}}^{\mathrm{T}}$.

The invariance to rotations in the occupied-occupied block can be used to fashion localized orbitals, for instance using an unitary optimization procedure [35]. Although localized orbitals are not strictly speaking observables—due to which several localization criteria have been suggested in the literature [36,37,38,39,40]—they have been shown to offer an effective way to study chemical reactions with ab initio calculations [41,42,43].

In addition to the occupied orbitals, in general there are also a number of unoccupied orbitals, which are commonly known as virtual orbitals. The number of virtual orbitals in any given calculation depends on the size of the basis set: the bigger the basis is, the more virtual orbitals there are. Because the virtual orbitals do not enter into the density matrix, the HF and DFT energy expression, Equation (7), is also invariant to rotations in the virtual-virtual block, similarly allowing their localization. However, as will be seen below in Section 9, the energy can be changed by mixing virtual orbitals into the occupied orbitals [44]. This approach provides another way to optimize the orbitals directly with, e.g., a gradient descent method, such as the geometric direct minimization method described in [45].

5. Spin-Restriction vs. Unrestriction

The molecular orbitals are obtained from the requirement that they minimize the total energy according to Equation (7). However, one must first choose the used spin formalism. The general choice is to use different spatial orbitals for the $\alpha$ and $\beta$ electrons, in which case a spin-unrestricted approach is obtained. The unrestricted approach is often used even in systems in which there are an equal number of alpha and beta electrons, $N_{\alpha }=N_{\beta }$: although the spin-restricted and unrestricted descriptions often reproduce matching results for such systems near the equilibrium, only the unrestricted formalism is able to break bonds in general. The reason for this is that when molecules are stretched past the Coulson–Fischer point [46], the optimal orbitals spontaneously break spin symmetry, which can only be described in the unrestricted formalism. At variance, in the spin-restricted case the electrons occupy a common set of $N=N_{\alpha }=N_{\beta }=(N_{\alpha }+N_{\beta })/2$ spatial orbitals. The limitation of the spatial orbitals to be the same for both spins, $C_{\mu i}=C_{\mu i}^{\alpha }=C_{\mu i}^{\beta }$, yields less variational freedom, and prevents the correct dissociation of e.g., the H${}_2$ molecule. As a flip side, the spin-restricted formalism affords computational savings over the unrestricted approach. The spin-restricted density matrices, Equation (5), reduce to

$$P_{\mu \nu }^{\alpha }=P_{\mu \nu }^{\beta }=\frac{1}{2}{P}_{\mu \nu }=\sum _i^N{C}_{\mu i}{C}_{\nu i}$$

(12)

meaning, e.g., that the $\alpha$ and $\beta$ exchange terms in Equation (7) coincide and can be simplified.

Spin-restriction is also possible in the case in which $N_{\alpha }\ne N_{\beta }$. In this case, a restricted open-shell method is obtained. Restricted open-shell methods are more involved than the spin-restricted and spin-unrestricted methods discussed in the present work. Restricted open-shell methods have been extensively discussed in [11] to which we refer for further discussion.

6. Self-Consistent Field Equations

Having chosen to use either spin-restricted or spin-unrestricted orbitals, one can proceed to the minimization of the energy expression in Equation (7). The energy expression depends only on the $\alpha$ and $\beta$ density matrices ${\mathbf{P}}^{\alpha }$ and ${\mathbf{P}}^{\beta }$ and their sum $\mathbf{P}$. The density matrices, in turn, are determined by the lowest $N_{\alpha }$ and $N_{\beta }$ molecular orbitals according to Equation (5). Because the energy expression in Equation (7) thus only depends on the density matrices ${\mathbf{P}}^{\alpha }$ and ${\mathbf{P}}^{\beta }$, it is expedient to use the chain rule to write, e.g.,

$$\frac{\partial }{\partial C_{\theta k}^{\alpha }}=\sum _{\eta \zeta }\frac{\partial P_{\eta \zeta }^{\alpha }}{\partial C_{\theta k}^{\alpha }}\frac{\partial }{\partial P_{\eta \zeta }^{\alpha }}$$

(13)

where the partial derivative of the density matrix element $P_{\eta \zeta }^{\alpha }$ is

$$\begin{array}{lll}\frac{\partial P_{\eta \zeta }^{\alpha }}{\partial C_{\theta k}^{\alpha }}& =& \frac{\partial }{\partial C_{\theta k}^{\alpha }}{\displaystyle \sum _{i=1}^{N_{\alpha }}}{C}_{\eta i}^{\alpha }{C}_{\zeta i}^{\alpha }\\ & =& {\displaystyle \sum _{i=1}^{N_{\alpha }}}{\delta }_{\theta \eta }{\delta }_{ki}{C}_{\zeta i}^{\alpha }+{\displaystyle \sum _{i=1}^{N_{\alpha }}}{C}_{\eta i}^{\alpha }{\delta }_{\theta \zeta }{\delta }_{ki}\hfill \\ & =& {\delta }_{\theta \eta }{C}_{\zeta k}^{\alpha }+C_{\eta k}^{\alpha }{\delta }_{\theta \zeta }\end{array}$$

(14)

The $\beta$ orbital derivative of the total density matrix has the same form as Equation (14), where all $\alpha$ are replaced with $\beta$. Note that these findings hold even when the same orbitals are used for both $\alpha$ and $\beta$ in a spin-restricted formalism, since the $\alpha$ and $\beta$ orbitals are formally independent.

Due to the chain rule, Equation (13), all we need are the density matrix derivatives of the energy expression. We only have to calculate the derivatives of the energy expression of Equation (7) for one spin, as the energy expression is symmetric with respect to the $\alpha$ and $\beta$ densities. It does not matter which spin we choose to be “up”; the expressions for the other spin will follow by symmetry by interchanging $\alpha \leftrightarrow \beta$. The first term of Equation (7) yields simply

$$\frac{\partial E^H}{\partial P_{\eta \zeta }^{\alpha }}=\frac{\partial }{\partial P_{\eta \zeta }^{\alpha }}\sum _{\mu \nu }{P}_{\mu \nu }{H}_{\mu \nu }=H_{\eta \zeta }$$

(15)

Next, taking the partial derivative with respect to $P_{\eta \zeta }^{\alpha }$ of the Coulomb and exchange terms in Equation (7) results in

$$\begin{array}{ccc}\hfill \frac{\partial E^J}{\partial P_{\eta \zeta }^{\alpha }}& =& \sum _{\mu \nu }(P_{\mu \nu }^{\alpha }+P_{\mu \nu }^{\beta })\left(\mu \nu \right|\eta \zeta )=\sum _{\mu \nu }{P}_{\mu \nu }\left(\mu \nu \right|\eta \zeta )=J_{\eta \zeta }\hfill \end{array}$$

(16)

where $\mathbf{J}$ is known as the Coulomb matrix, and

$$\begin{array}{ccc}\hfill \frac{\partial E^K}{\partial P_{\eta \zeta }^{\alpha }}& =& -a\sum _{\mu \nu }{P}_{\mu \nu }^{\alpha }\left(\mu \eta \right|\nu \zeta )=-a{K}_{\eta \zeta }^{\alpha }\hfill \end{array}$$

(17)

where ${\mathbf{K}}^{\alpha }$ is the spin-$\alpha$ exchange matrix, respectively. The Coulomb and exchange matrices can be used to rewrite the energy expression in Equation (7) as

$$E=\sum _{\mu \nu }{P}_{\mu \nu }{H}_{\mu \nu }+\frac{1}{2}\sum _{\mu \nu }{P}_{\mu \nu }{J}_{\mu \nu }-\frac{a}{2}\sum _{\mu \nu }(P_{\mu \nu }^{\alpha }{K}_{\mu \nu }^{\alpha }+P_{\mu \nu }^{\beta }{K}_{\mu \nu }^{\beta })+b\int f\left(\overrightarrow{r}\right)\mathrm{d}\overrightarrow{r}$$

(18)

Note that in contrast to the Coulomb and exact exchange terms, the exchange-correlation term does not undergo simplifications, because the exchange-correlation term is not quadratic in the density matrix, as will be seen later in Section 11. For the time being, we will denote the partial derivative of $E^{\mathrm{XC}}$ with respect to $P_{\eta \zeta }^{\alpha }$ as

$$\frac{\partial E^{\mathrm{XC}}}{\partial P_{\eta \zeta }^{\alpha }}=b{K}_{\eta \zeta }^{\mathrm{XC};\alpha }$$

(19)

as the full expressions for ${\mathbf{K}}^{\mathrm{XC};\alpha }$ will be presented in Section 11. Now, collecting the partial derivatives in Equations (15)–(19) gives us the density matrix derivatives of the energy expression as

$$\frac{\partial E}{\partial P_{\eta \zeta }^{\sigma }}=H_{\eta \zeta }+J_{\eta \zeta }-a{K}_{\eta \zeta }^{\sigma }+b{K}_{\eta \zeta }^{\mathrm{XC};\sigma }=F_{\eta \zeta }^{\sigma }$$

(20)

where we have identified the Kohn–Sham–Fock matrices ${\mathbf{F}}^{\sigma }$, where $\sigma$ denotes $\alpha$ or $\beta$. Because the density matrices defined by Equation (5) are symmetric, also the Fock matrices are symmetric, $F_{\eta \zeta }^{\sigma }=F_{\zeta \eta }^{\sigma }$. Note that since the Fock matrices only depend on the density matrices, also they are invariant to occupied-occupied and virtual-virtual rotations, ${\mathbf{F}}^{{\sigma }^{\prime }}={\mathbf{F}}^{\sigma }$.

Naïvely, one would obtain the orbital derivative of the full energy expression in Equation (7) with Equations (13), (14) and (20) and then set it to zero to yield an equation for the unknown expansion coefficients ${\mathbf{C}}^{\alpha }$. However, the molecular orbitals cannot be varied freely: one must make sure that the orbitals stay orthonormal during the variation, as otherwise the Pauli exclusion principle would be violated. For instance, the orthonormality condition for the $\alpha$ electrons is

$$\int {\phi }_i^{\alpha }\left(\overrightarrow{r}\right){\phi }_j^{\alpha }\left(\overrightarrow{r}\right)dr={\delta }_{ij}$$

(21)

The way to enforce these conditions is to use Lagrangian multipliers ${\epsilon }_{ij}$. That is, instead of the bare energy expression E, we will optimize the Lagrangian

$$L=E-\sum _{ij}{\epsilon }_{ij}^{\alpha }\left[\int \mathrm{d}\overrightarrow{r}{\phi }_i^{\alpha }\left(\overrightarrow{r}\right){\phi }_j^{\alpha }\left(\overrightarrow{r}\right)-{\delta }_{ij}\right]-\sum _{ij}{\epsilon }_{ij}^{\beta }\left[\int \mathrm{d}\overrightarrow{r}{\phi }_i^{\beta }\left(\overrightarrow{r}\right){\phi }_j^{\beta }\left(\overrightarrow{r}\right)-{\delta }_{ij}\right]$$

(22)

where the sums over i and j run over all orbitals; that is, both the occupied and the virtual ones. We can see from Equation (22) that the matrices of Lagrangian multipliers ${\epsilon }^{\alpha }$ and ${\epsilon }^{\beta }$ can be chosen to be symmetric. For instance, if ${\epsilon }^{\alpha }$ contained a symmetric part ${\epsilon }_s^{\alpha }$ and an antisymmetric part ${\epsilon }_a^{\alpha }$, ${\epsilon }^{\alpha }={\epsilon }_s^{\alpha }+{\epsilon }_a^{\alpha }$, the contribution from the antisymmetric part would vanish because it is multiplied with the orbital overlap that is symmetric.

Next, we can calculate $\partial L/\partial C_{\theta k}^{\alpha }$, where $\partial E/\partial C_{\theta k}^{\alpha }$ is given by Equations (13), (14), (20) and the derivative of the constraint term is given by

$$\begin{array}{ccc}\frac{\partial }{\partial C_{\theta k}^{\alpha }}\sum _{ij}{\epsilon }_{ij}^{\alpha }\left[\int \mathrm{d}\overrightarrow{r}{\phi }_i^{\alpha }\left(\overrightarrow{r}\right){\phi }_j^{\alpha }\left(\overrightarrow{r}\right)-{\delta }_{ij}\right]\hfill & =& \sum _{ij}{\epsilon }_{ij}^{\alpha }\frac{\partial }{\partial C_{\theta k}^{\alpha }}\sum _{\eta \zeta }\left[C_{\eta i}^{\alpha }{S}_{\eta \zeta }{C}_{\zeta j}^{\alpha }-{\delta }_{ij}\right]\hfill \\ & =& \sum _{ij}{\epsilon }_{ij}^{\alpha }\sum _{\eta \zeta }\left[{\delta }_{\eta \theta }{\delta }_{ki}{C}_{\zeta j}^{\alpha }+C_{\eta i}^{\alpha }{\delta }_{\zeta \theta }{\delta }_{kj}\right]S_{\eta \zeta }\hfill \\ & =& \sum _{\eta }\sum _i{\epsilon }_{ki}^{\alpha }{C}_{\eta i}^{\alpha }{S}_{\theta \eta }+\sum _{\eta }\sum _i{\epsilon }_{ik}^{\alpha }{C}_{\eta i}^{\alpha }{S}_{\eta \theta }\hfill \end{array}$$

(23)

where on the third line dummy summation indices have been renamed from j to i and $\zeta$ to $\eta$. The derivative can be evaluated as

$$\begin{array}{ccc}\frac{\partial L}{\partial C_{\theta k}^{\alpha }}\hfill & =& \sum _{\eta \zeta }\frac{\partial E}{\partial P_{\eta \zeta }^{\alpha }}\frac{\partial P_{\eta \zeta }^{\alpha }}{\partial C_{\theta k}^{\alpha }}-\sum _{ij}{\epsilon }_{ij}^{\alpha }\frac{\partial }{\partial C_{\theta k}^{\alpha }}\int \mathrm{d}\overrightarrow{r}{\phi }_i^{\alpha }\left(\overrightarrow{r}\right){\phi }_j^{\alpha }\left(\overrightarrow{r}\right)\hfill \\ & =& \sum _{\eta \zeta }{F}_{\eta \zeta }^{\alpha }\left[{\delta }_{\theta \eta }{C}_{\zeta k}^{\alpha }+C_{\eta k}^{\alpha }{\delta }_{\theta \zeta }\right]-\sum _{\zeta }\sum _i{\epsilon }_{ki}^{\alpha }{C}_{\zeta i}^{\alpha }{S}_{\theta \zeta }-\sum _{\eta }\sum _i{\epsilon }_{ik}^{\alpha }{C}_{\eta i}^{\alpha }{S}_{\eta \theta }\hfill \\ & =& 2\sum _{\eta }{F}_{\theta \eta }^{\alpha }{C}_{\eta k}^{\alpha }-2\sum _{\eta }\sum _i{\epsilon }_{ki}^{\alpha }{S}_{\theta \eta }{C}_{\eta i}^{\alpha }\hfill \end{array}$$

(24)

because $\mathbf{F}$ and $\epsilon$ are symmetric, and dummy summation indices can be renamed.

The optimal orbitals satisfy the stationary condition $\partial L/\partial C_{\theta k}^{\alpha }=0$ from which

$$\begin{array}{c}\hfill \sum _{\eta }{F}_{\theta \eta }^{\alpha }{C}_{\eta k}^{\alpha }=\sum _{\eta }\sum _i{S}_{\theta \eta }{C}_{\eta i}^{\alpha }{\epsilon }_{ik}^{\alpha }\end{array}$$

(25)

Equation (25) can thus be written in matrix form as

$${\mathbf{F}}^{\alpha }{\mathbf{C}}^{\alpha }=\mathbf{S}{\mathbf{C}}^{\alpha }{\mathbf{E}}^{\alpha }$$

(26)

where ${\mathbf{E}}^{\alpha }$ is the (symmetric) matrix of Lagrangian multipliers.

Because ${\mathbf{E}}^{\alpha }$ is symmetric, it can be diagonalized and it has real eigenvalues. Let us now assume that ${\mathbf{U}}^{\alpha }$ is an orthogonal matrix that diagonalizes ${\mathbf{E}}^{\alpha }$

$$E_{ij}^{\alpha }\to E_{ij}^{{\alpha }^{\prime }}=\sum _{kl}{U}_{ki}^{\alpha }{E}_{kl}^{\alpha }{U}_{lj}^{\alpha }={\delta }_{ij}{\epsilon }_i^{\alpha }$$

(27)

where ${\epsilon }^{\alpha }$ are the eigenvalues. Re-expressing the orbital coefficients ${\mathbf{C}}^{\alpha }$ in terms of a new set of orbitals rotated by ${\mathbf{U}}^{\alpha }$ with Equation (10), ${\mathbf{C}}^{\alpha }={\mathbf{C}}^{{\alpha }^{\prime }}{\left({\mathbf{U}}^{\alpha }\right)}^{\mathrm{T}}$, rewrites Equation (26) in the form

$${\mathbf{F}}^{\alpha }{\mathbf{C}}^{{\alpha }^{\prime }}{\left({\mathbf{U}}^{\alpha }\right)}^{\mathrm{T}}=\mathbf{S}{\mathbf{C}}^{{\alpha }^{\prime }}{\left({\mathbf{U}}^{\alpha }\right)}^{\mathrm{T}}{\mathbf{E}}^{\alpha }$$

that can be multiplied from the right by ${\mathbf{U}}^{\alpha }$ producing

$${\mathbf{F}}^{\alpha }{\mathbf{C}}^{{\alpha }^{\prime }}=\mathbf{S}{\mathbf{C}}^{{\alpha }^{\prime }}{\mathbf{E}}^{{\alpha }^{\prime }}$$

(28)

where, according to Equation (27), ${\mathbf{E}}^{{\alpha }^{\prime }}={\left({\mathbf{U}}^{\alpha }\right)}^{\mathrm{T}}{\mathbf{E}}^{\alpha }{\mathbf{U}}^{\alpha }$ is a diagonal matrix with elements ${\epsilon }_i^{\alpha }$.

Equation (28) is almost what we want—an equation in the rotated basis that looks like Equation (26) with a diagonal ${\mathbf{E}}^{{\alpha }^{\prime }}$—but one problem remains: the Kohn–Sham–Fock matrix is still the one corresponding to the original orbitals ${\mathbf{C}}^{\alpha }$ instead of the transformed orbitals ${\mathbf{C}}^{{\alpha }^{\prime }}$, while the orbital rotation by $\mathbf{U}$ that takes us from ${\mathbf{C}}^{\alpha }$ to ${\mathbf{C}}^{{\alpha }^{\prime }}$ may lead to a different Kohn–Sham–Fock matrix ${\mathbf{F}}^{{\alpha }^{\prime }}\ne {\mathbf{F}}^{\alpha }$. However, if we choose the form of $\mathbf{U}$ such that the occupied-virtual (ov) and virtual-occupied (vo) blocks vanish

$$\mathbf{U}=\left(\begin{array}{cc}{\mathbf{U}}_{\mathrm{oo}}& {\mathbf{U}}_{\mathrm{ov}}\\ {\mathbf{U}}_{\mathrm{vo}}& {\mathbf{U}}_{\mathrm{vv}}\end{array}\right)=\left(\begin{array}{cc}{\mathbf{U}}_{\mathrm{oo}}& \mathbf{0}\\ \mathbf{0}& {\mathbf{U}}_{\mathrm{vv}}\end{array}\right)$$

(29)

then $\mathbf{U}$ only rotates occupied orbitals with occupied orbitals and virtual orbitals with virtual orbitals, meaning that the orbital rotation does not change the density matrix given in Equation (11). Then, the Fock matrix corresponding to the rotated orbitals coincides with that of the original orbitals, ${\mathbf{F}}^{\prime }=\mathbf{F}$, completing the proof that ${\mathbf{E}}^{\alpha }$ can be chosen to be diagonal. (Occupied-virtual rotations ${\mathbf{U}}_{\mathrm{ov}}$ or ${\mathbf{U}}_{\mathrm{ov}}$, discussed in more detail in Section 9, are in fact here forbidden: the SCF equations were derived with the assumption that the energy is stationary, but this condition would instantly be violated by such rotations.)

We have thus obtained the Berthier–Pople–Nesbet [8,9,13,16] equations for the orbital coefficients

$$\left\{\begin{array}{c}{\mathbf{F}}^{\alpha }{\mathbf{C}}^{\alpha }=\mathbf{S}{\mathbf{C}}^{\alpha }{\mathbf{E}}^{\alpha }\hfill \\ {\mathbf{F}}^{\beta }{\mathbf{C}}^{\beta }=\mathbf{S}{\mathbf{C}}^{\beta }{\mathbf{E}}^{\beta }\hfill \end{array}\right.$$

(30)

where the primes have become unnecessary and have been omitted for simplicity. The elements of the Kohn–Sham–Fock matrices ${\mathbf{F}}^{\alpha }$ and ${\mathbf{F}}^{\beta }$ are given by Equation (20), and the orbital energy matrices ${\mathbf{E}}^{\alpha }$ and ${\mathbf{E}}^{\beta }$ are diagonal. In the spin-restricted case [5,16] the $\alpha$ and $\beta$ molecular orbitals coincide, leading to identical density matrices ${\mathbf{P}}^{\alpha }$ and ${\mathbf{P}}^{\beta }$, and identical Fock matrices ${\mathbf{F}}^{\alpha }$ and ${\mathbf{F}}^{\beta }$. In this case, the SCF equations simplify to the Roothaan–Hall form

$$\mathbf{F}\mathbf{C}=\mathbf{S}\mathbf{C}\mathbf{E}$$

(31)

which was already mentioned in the Introduction.

7. Solution of Self-Consistent Field Equations

The Roothaan–Hall and Berthier–Pople–Nesbet expressions take the form of a generalized eigenvalue equation. The conventional way to solve these equations is to re-express the (unknown) orbital coefficients in terms of a matrix $\mathbf{X}$ as

$$\mathbf{C}=\mathbf{X}\tilde{\mathbf{C}}$$

(32)

Inserting Equation (32) into the Roothaan–Hall equation, Equation (31), yields

$$\mathbf{F}\mathbf{X}\tilde{\mathbf{C}}=\mathbf{S}\mathbf{X}\tilde{\mathbf{C}}\mathbf{E}$$

which can be multiplied from the left with ${\mathbf{X}}^{\mathrm{T}}$ to yield

$${\mathbf{X}}^{\mathrm{T}}\mathbf{F}\mathbf{X}\tilde{\mathbf{C}}={\mathbf{X}}^{\mathrm{T}}\mathbf{S}\mathbf{X}\tilde{\mathbf{C}}\mathbf{E}$$

This means that the orbital transform of Equation (32) yields a new generalized eigenvalue equation

$$\tilde{\mathbf{F}}\tilde{\mathbf{C}}=\tilde{\mathbf{S}}\tilde{\mathbf{C}}\mathbf{E}$$

(33)

where $\tilde{\mathbf{F}}={\mathbf{X}}^{\mathrm{T}}\mathbf{F}\mathbf{X}$ and $\tilde{\mathbf{S}}={\mathbf{X}}^{\mathrm{T}}\mathbf{S}\mathbf{X}$. Now, if we choose $\mathbf{X}$ in such a way that $\tilde{\mathbf{S}}=\mathbf{1}$, Equation (33) reduces to a normal eigenvalue equation

$$\tilde{\mathbf{F}}\tilde{\mathbf{C}}=\tilde{\mathbf{C}}\mathbf{E}$$

(34)

which can be solved with standard techniques. Then, the wanted orbital coefficients $\mathbf{C}$ can be calculated from $\tilde{\mathbf{C}}$ using Equation (32).

If the basis set is well-conditioned, the matrix $\mathbf{X}$ can be chosen as

$$\mathbf{X}=\mathbf{V}{\mathsf{\Lambda }}^{-1/2}{\mathbf{V}}^{\mathrm{T}}$$

(35)

where $\mathbf{V}$ and $\mathsf{\Lambda }$ are the eigenvectors and eigenvalues of $\mathbf{S}$

$$\mathbf{S}=\mathbf{V}\mathsf{\Lambda }{\mathbf{V}}^{\mathrm{T}}$$

(36)

This procedure is known as symmetric orthogonalization [7].

However, if a large LCAO basis is used, the atomic orbital basis functions centered on different atoms may generate linear dependencies in the basis, making the basis set expansion ambiguous. These linear dependencies can be removed with the “canonical” orthonormalization procedure [47], in which

$$\mathbf{X}={\mathbf{V}}^{\prime }{{\mathsf{\Lambda }}^{\prime }}^{-1/2}$$

(37)

where only those eigenvectors ${\mathbf{V}}_i$ with large enough eigenvalues ${\lambda }_i\ge \tau$ are included. The threshold $\tau$ is typically of the order of $\tau ={10}^{-7}\dots {10}^{-5}$, and its value may have a noticeable effect on, e.g., the absolute energies that result from a SCF calculation; relative energies, however, should be less sensitive to $\tau$. If no eigenvalues fall under the threshold $\tau$, the symmetric and canonical orthogonalization approaches become equivalent for the purposes of SCF calculations in the case of a well-conditioned basis set: both yield an orthonormal basis of the same size, which will yield the same variational ground state energy.

Unnormalized basis sets can also be handled easily by the orthogonalization procedure. Although in principle it is not necessary to normalize the individual basis functions before obtaining an orthonormal basis by Equations (35) and (37), computer linear algebra packages may fail to find the eigenvalues and eigenvectors in a reliable fashion if the basis functions have pronouncedly different norms. Moreover, missing normalization of the basis set affects the eigenvalues, which has repercussions for canonical orthogonalization. These issues can be circumvented by normalizing the overlap matrix $\mathbf{S}\to {\mathbf{S}}^{\prime }=\mathbf{n}\mathbf{S}\mathbf{n}$ where $n_{ij}=S_{ii}^{-1/2}{\delta }_{ij}$ before using Equations (35) and (37) [3,4]. The orthogonalizer for the unnormalized basis set is obtained as $\mathbf{X}\to \mathbf{n}\mathbf{X}$; it is easy to see that this satisfies the necessary condition ${\mathbf{X}}^{\mathrm{T}}\mathbf{S}\mathbf{X}=\mathbf{1}$ even though the symmetricity of $\mathbf{X}$ for the case of Equation (35) will be lost.

Even if $\mathbf{S}$ has been properly normalized, the use of the symmetric or canonical orthogonalization procedures still requires that the diagonalization of $\mathbf{S}$ is numerically stable. However, whenever a large number of linear dependencies exists in the basis set (e.g., a large number of diffuse functions are used or two nuclei are close together), $\mathbf{S}$ may become so ill-conditioned it cannot be accurately diagonalized. In such cases it is possible to reduce the size of the basis set without losing a significant amount of accuracy by an automatic procedure, see [48,49] for details.

8. Why Does the Self-Consistent Field Method Minimize the Energy?

The SCF equations, Equation (30) or Equation (31), offer a way to solve for the molecular orbitals described by $\mathbf{C}$ from a Kohn–Sham–Fock matrix $\mathbf{F}$ by finding its eigenvectors from Equation (34). However, the Kohn–Sham–Fock matrix depends on the density matrices, which are built from the molecular orbitals according to the Aufbau principle. In the SCF procedure, one tries to find a self-consistent solution: $\mathbf{C}$ yields $\mathbf{F}$, whose eigenvectors are $\mathbf{C}$. The procedure starts from an initial guess for the orbitals ${\mathbf{C}}^{\sigma }$ or the density matrices ${\mathbf{P}}^{\sigma }$, which have been recently reviewed and benchmarked in [50] to which we refer for further details.

Why does the self-consistent field procedure—diagonalizing ${\mathbf{F}}^{\sigma }$ to update the orbital coefficients ${\mathbf{C}}^{\sigma }$—correspond to minimization of the Hartree–Fock/Kohn–Sham energy? For simplicity, let us examine the case of HF theory. The energy expression, Equation (18), can be written in this case ($a=1$, $b=0$) as

$$E=\mathrm{Tr}\mathbf{P}\mathbf{H}+\frac{1}{2}\mathrm{Tr}\mathbf{P}\mathbf{J}-\frac{1}{2}\mathrm{Tr}{\mathbf{P}}^{\alpha }{\mathbf{K}}^{\alpha }-\frac{1}{2}\mathrm{Tr}{\mathbf{P}}^{\beta }{\mathbf{K}}^{\beta }$$

(38)

The Fock matrix elements, Equation (20), are given by

$$\begin{array}{ccc}\hfill {\mathbf{F}}^{\alpha }& =& \mathbf{H}+\mathbf{J}-{\mathbf{K}}^{\alpha }\hfill \end{array}$$

(39)

$$\begin{array}{ccc}\hfill {\mathbf{F}}^{\beta }& =& \mathbf{H}+\mathbf{J}-{\mathbf{K}}^{\beta }\hfill \end{array}$$

(40)

Equation (38) can be rewritten with Equations (39) and (40) as

$$E=\frac{1}{2}\mathrm{Tr}{\mathbf{P}}^{\alpha }(\mathbf{H}+{\mathbf{F}}^{\alpha })+\frac{1}{2}\mathrm{Tr}{\mathbf{P}}^{\beta }(\mathbf{H}+{\mathbf{F}}^{\beta })$$

(41)

Expanding the density matrices using Equation (5) we see that Equation (41) can be written as

$$E=\frac{1}{2}\sum _{i=1}^{N_{\alpha }}(h_{ii}^{\alpha }+f_{ii}^{\alpha })+\frac{1}{2}\sum _{i=1}^{N_{\beta }}(h_{ii}^{\beta }+f_{ii}^{\beta })$$

(42)

where the core Hamiltonian and Fock matrices have been written in the molecular orbital basis, ${\mathbf{h}}^{\sigma }={\left({\mathbf{C}}^{\sigma }\right)}^{\mathrm{T}}\mathbf{H}{\mathbf{C}}^{\sigma }$ and ${\mathbf{f}}^{\sigma }={\left({\mathbf{C}}^{\sigma }\right)}^{\mathrm{T}}{\mathbf{F}}^{\sigma }{\mathbf{C}}^{\sigma }$.

If one were to start the calculation from the core guess, then ${\sum }_i{h}_{ii}^{\alpha }$ and ${\sum }_i{h}_{ii}^{\beta }$ would be minimized. However, as discussed in [50], this is a horrible choice as it completely disregards electronic repulsion effects, meaning that the ${\sum }_i{f}_{ii}^{\alpha }$ and ${\sum }_i{f}_{ii}^{\beta }$ terms are far from optimal. The Roothaan step—obtaining new molecular orbitals by diagonalization of the Fock matrix—results in a minimization of the ${\sum }_i{f}_{ii}^{\alpha }$ and ${\sum }_i{f}_{ii}^{\beta }$ terms, as after diagonalization only the lowest orbitals become populated and the sum thus runs only over the lowest eigenvalues $f_{ii}^{\sigma }$. After the update, ${\sum }_i{h}_{ii}^{\alpha }$ and ${\sum }_i{h}_{ii}^{\beta }$ no longer yield their lowest possible values. However, the increase in the value of ${\sum }_i{h}_{ii}^{\alpha }+{\sum }_i{h}_{ii}^{\beta }$ should be much smaller than the decrease in the value of ${\sum }_i{f}_{ii}^{\alpha }+{\sum }_i{f}_{ii}^{\beta }$, as the Fock matrices ${\mathbf{f}}^{\alpha }$ and ${\mathbf{f}}^{\beta }$ also contain the core Hamiltonian. It is thus seen that Roothaan’s self-consistent field method, that is, the iterative diagonalization of the Fock matrix minimizes the energy.

However, the minimization is only valid for a fixed potential ${\mathbf{f}}^{\sigma }$ in which the electrons are moving. When the orbitals are changed—as happens when $\mathbf{f}$ is made diagonal and its lowest eigenvectors occupied—a new Fock matrix $\mathbf{F}$ must be built and a new $\mathbf{f}$ constructed: the potential also changes with the electron density. If the orbitals were far from their optimal values, $\mathbf{P}$ and therefore $\mathbf{F}$ may change quite radically by the orbital update. This means that even though $\mathbf{f}$ was made diagonal in the previous iteration, it is no longer diagonally dominant after it has been updated. Indeed, the straightforward iterative diagonalization procedure often fails to converge for all but the simplest systems, because the density tends to undergo large oscillations in the naïve self-consistency cycle. To make the method usable, the convergence of the fixed-point problem of finding a $\mathbf{C}$ that generates $\mathbf{F}$ that generates $\mathbf{C}$ must be stabilized or accelerated in some way. This can be achieved, e.g., by damping [51,52], level shifts [53,54,55], or extrapolation [56,57,58,59]. Fractional occupations can also be used in the initial iterations to aid convergence [60].

The argument for why density functional calculations converge similarly to HF with the iterative Roothaan procedure is somewhat less obvious, because unlike HF the exchange-correlation functional is not generally quadratic in the density. However, the total energy expression is approximately quadratic also in DFT when one is sufficiently close to an extremal point, as is easily seen by a Taylor expansion of Equation (7). In practice the iterative procedure works well also for DFT, whose contributions to the Kohn–Sham–Fock matrix we will discuss in Section 11.

9. Direct Minimization of the Energy

Instead of solving the orbitals from the SCF equations, which were obtained in Section 6 from the stationary condition for the energy under the constraint of orthonormal orbitals, the orbitals can also be optimized by a direct minimization of the energy. As was discussed in Section 4, the energy expression of Equation (7) is invariant to occupied-occupied and virtual-virtual rotations. This means that if we have $o^{\sigma }$ occupied orbitals and $v^{\sigma }$ virtual orbitals from some initial guess (see possible choices in [50]) for spin $\sigma$, we can consider the energy as a function of a set of $o^{\sigma }{v}^{\sigma }$ rotation angles [44] by examining a rotation of the orbitals via Equation (10) by an orthogonal matrix

$${\mathbf{U}}^{\sigma }\left({\mathbf{\theta }}^{\sigma }\right)=exp\left(\begin{array}{cc}\mathbf{0}& {\mathbf{\theta }}^{\sigma }\\ -{\left({\mathbf{\theta }}^{\sigma }\right)}^{\mathrm{T}}& \mathbf{0}\end{array}\right)$$

(43)

where ${\mathbf{\theta }}^{\sigma }$ is an $o^{\sigma }\times v^{\sigma }$ matrix containing the rotation angles. The rotation matrix determined by Equation (43) reduces to an identity matrix for vanishing rotation parameters, ${\mathbf{\theta }}^{\sigma }=\mathbf{0}$. Because the rotation matrix of Equation (43) is orthogonal, it automatically preserves the orthonormality of the orbitals, and special tricks i.e. Lagrangian multipliers are not needed to enforce this behavior.

The change in the density matrix is given by

$$\frac{\partial P_{\eta \zeta }^{\sigma }}{\partial {\theta }_{ia}^{\sigma }}=\sum _{k=1}^{N_{\sigma }}\left[\frac{\partial C_{\eta k}^{\sigma }}{\partial {\theta }_{ia}^{\sigma }}{C}_{\zeta k}^{\sigma }+C_{\eta k}^{\sigma }\frac{\partial C_{\zeta k}^{\sigma }}{\partial {\theta }_{ia}^{\sigma }}\right]$$

(44)

How do the orbital coefficients change? Remembering that the first $N_{\sigma }$ orbitals are occupied, and the rest are virtual, we can write ${\mathbf{C}}^{\sigma }=(\begin{array}{cc}{\mathbf{C}}_o^{\sigma }& {\mathbf{C}}_v^{\sigma })\end{array}$. After an infinitesimal rotation $\mathbf{\theta }\approx \mathbf{0}$, the occupied orbitals change into ${\mathbf{C}}_o^{{\sigma }^{\prime }}\approx {\mathbf{C}}_o^{\sigma }-{\mathbf{C}}_v^{\sigma }{\left({\mathbf{\theta }}^{\sigma }\right)}^{\mathrm{T}}$, that is

$$C_{\eta k}^{{\sigma }^{\prime }}=C_{\eta k}^{\sigma }-\sum _{b\mathrm{virtual}}{C}_{\eta b}^{\sigma }{\theta }_{kb}^{\sigma }$$

(45)

from which $\partial C_{\eta k}^{\sigma }/\partial {\theta }_{ia}^{\sigma }=-C_{\eta a}^{\sigma }{\delta }_{ik}$. Now the gradient of the energy with respect to rotation of the current set of orbitals can be obtained as

$$\frac{\partial E}{\partial {\theta }_{ia}^{\sigma }}=\sum _{\eta \zeta }\frac{\partial E}{\partial P_{\eta \zeta }^{\sigma }}\frac{\partial P_{\eta \zeta }^{\sigma }}{\partial {\theta }_{ia}^{\sigma }}=\sum _{\eta \zeta }{F}_{\eta \zeta }^{\sigma }\sum _{k=1}^{N_{\sigma }}\left[\frac{\partial C_{\eta k}^{\sigma }}{\partial {\theta }_{ia}^{\sigma }}{C}_{\zeta k}^{\sigma }+C_{\eta k}^{\sigma }\frac{\partial C_{\zeta k}^{\sigma }}{\partial {\theta }_{ia}^{\sigma }}\right]=-\sum _{\eta \zeta }{F}_{\eta \zeta }^{\sigma }\left[C_{\eta a}^{\sigma }{C}_{\zeta i}^{\sigma }+C_{\eta i}^{\sigma }{C}_{\eta a}^{\sigma }\right]=-2f_{ia}^{\sigma }$$

(46)

where ${\mathbf{f}}^{\sigma }={\left({\mathbf{C}}^{\sigma }\right)}^{\mathrm{T}}{\mathbf{F}}^{\sigma }{\mathbf{C}}^{\sigma }$ is the Fock matrix in the MO basis. Direct minimization of Equation (7) can then be pursued using Equation (46) with, e.g., gradient descent methods. However, a proper preconditioning of the search direction is essential in order for the algorithm to be usable; see, e.g., the geometric direct minimization method described in [45]. Many other direct minimization methods for the HF or DFT energy have also been proposed, and we refer the interested reader to the vast existing literature that cannot be comprehensively cited here.

10. SCF vs. Direct Minimization

Having described two alternative ways for solving the orbitals, we can discuss their advantages and disadvantages. The self-consistent field method is hard to beat for systems where convergence is straightforward: a suitably stabilized and accelerated SCF procedure often converges within 10 to 20 iterations when a suitable initial guess (see [50]) has been provided. However, when the gap between the highest occupied and lowest unoccupied orbital is small, which commonly occurs in, e.g., first-row transition metal complexes, the SCF procedure may become extremely slow to converge, oscillate between two or more solutions, converge to a higher-lying solution, or even to a saddle-point solution. Namely, it is critically important to realize that even if the orbital gradient vanishes, or equivalently, that the SCF equations are fulfilled, this does not mean that the energy expression Equation (7) truly has been minimized. Because there are typically several occupied as well as virtual orbitals, the minimization problem involves a large number of degrees of freedom. In multivariate calculus, a vanishing gradient only means that the orbitals correspond to some kind of extremum of the energy: a local minimum, a saddle-point solution, or even a local maximum, although the lattermost is highly improbable in SCF calculations.

In contrast to the sometimes erratic behavior of the SCF method, direct minimization based on orbital rotations is guaranteed to converge onto an extremal point $f_{ia}=0$ per the theory of numerical analysis; this is of great worth when studying systems with complicated electronic structures for which conventional SCF algorithms fail. However, more predictable convergence does not come for free: the downside of direct minimization methods is that they carry a higher computational cost due to, e.g., the use of line searches in the orbital optimization. Direct minimization methods can also be formulated at the second order, yielding more robust convergence to a local minimum at the cost of more computational resources per iteration, see, e.g., [61,62,63]. Because direct minimization methods are based on an explicit rotation of the orbitals, they are able to always follow the same solution at variance to SCF methods where the orbital occupations are typically reset at every iteration according to the Aufbau principle. Because of this, direct minimization can lead to a solution where the Aufbau rule is violated, that is, the highest occupied orbital lies higher in energy than the lowest unoccupied orbital. Direct minimization methods can also be straightforwardly applied in more complicated electronic structure theories than self-consistent field theory. Such methods may especially include explicit dependence on the molecular orbitals, as discussed by one of the present authors in [31,64,65] for the Perdew–Zunger self-interaction correction [66] which depends explicitly on the $N_{\sigma }$ occupied orbitals, and [67] for the perfect quadruples [68] and perfect hextuples [69] models that also depend on the $N_{\sigma }$ corresponding virtual orbitals.

In order to check the character of the extremum found by the SCF procedure or a direct minimization method, it is necessary to continue the analysis to second-order changes in the energy with respect to the orbital rotations by finding the lowest eigenvalue of the Hessian matrix: if it is negative, rotating the orbitals in the direction of the corresponding eigenvector will result in a further decrease of the energy. Whenever post-HF calculations are performed, or benchmark-quality values are sought at the SCF level, stability analysis [70,71] should be used to guarantee that the wave function indeed corresponds to a local minimum. Alternatively, trust-region methods [72,73,74] can be employed to ensure that the orbitals converge onto a true local minimum.

As always in the minimization of multivariate functions, locating the global minimum is difficult, and typically the best one can hope for is to find a local minimum. Some systems permit several local electronic minima: for instance, charge transfer complexes may allow both a neutral (X ⋯ Y) as well as an ionic (X${}^{+}\cdots$ Y${}^{-}$) solution. Finding such physically motivated solutions is often straightforward by suitable manipulations of the initial guess, for instance, by constructing guesses via the superposition of atomic potentials [75,76] with the correct atomic charges. Sometimes it may also be interesting to locate saddle-point solutions, which have physical interpretations as excited states. Specific excited states can be explored within the SCF approach by replacing Aufbau population of the orbitals with overlap criteria [77,78] or with direct minimization by replacing the energy with the square of the gradient [79]; for instance, such an approach has been recently shown to predict highly accurate core spectra [80]. The full space of SCF solutions can be explored via, e.g., meta-dynamics [81].

11. Density Functional Contributions to Kohn–Sham–Fock Matrix

In Section 6 we derived expressions for the Kohn–Sham–Fock matrix elements for all but the density functional contribution

$$E^{\mathrm{XC}}=b\int f\left(\overrightarrow{r}\right)\mathrm{d}\overrightarrow{r}$$

(47)

which we will consider next. Hundreds of density functionals $f\left(\overrightarrow{r}\right)$ of various forms have been published in the literature in the recent decades [82], and offering a comprehensive selection thereof poses a considerable challenge to quantum chemistry software developers. This problem is further exacerbated by the need to keep track with the several new functionals still being published every year. Moreover, the density functionals $f\left(\overrightarrow{r}\right)$ typically carry extremely complicated functional forms, making their correct implementation painstaking work. The implementation is made even more difficult by the need to compute the first derivatives of $f\left(\overrightarrow{r}\right)$ for the SCF procedure, as well as several higher-order ones for, e.g., the calculation of various properties.

Fortunately, these challenges have been obviated by freely available, portable standard implementations such as LibXC [83] and XCFun [84]. The LibXC software package strives to implement all DFT functionals published in the literature, and provides a uniform interface to ∼500 functionals of various forms. At present, LibXC is used by ∼30 electronic structure programs based on various numerical representations that range from basis set approaches (Gaussian-type orbitals, Slater-type orbitals, numerical atomic orbitals, finite elements, plane waves) to finite difference procedures. New functionals only have to be added once to LibXC, meaning the library is easily kept up to date, after which they become available to all programs that support the corresponding rung of functionals on Jacob’s ladder [28]. Next, we will derive the equations necessary to implement the various rungs’ functionals in the variational basis set approach.

11.1. LDA Functionals

The simplest density functional approximations (DFAs), belonging to the first rung of Jacob’s Ladder [28], are generally referred to as local (spin) density approximations (LDAs). These are functions of only the electron density [13]

$$f\left(\overrightarrow{r}\right)=f({\rho }_{\alpha }\left(\overrightarrow{r}\right),{\rho }_{\beta }\left(\overrightarrow{r}\right))$$

(48)

such as the LDA exchange functional [85,86]

$$f\left(\overrightarrow{r}\right)=-\frac{3}{2}{\left(\frac{3}{4\pi }\right)}^{\frac{1}{3}}\left({\rho }_{\alpha }^{4/3}\left(\overrightarrow{r}\right)+{\rho }_{\beta }^{4/3}\left(\overrightarrow{r}\right)\right)$$

Assuming f has the form of Equation (48), the resulting contribution to the Kohn–Sham–Fock matrix $F_{\mu \nu }^{\alpha ;\mathrm{XC}}=\partial E^{\mathrm{XC}}/\partial P_{\mu \nu }^{\alpha }$ can be evaluated using Equation (5) for the densities at point $\overrightarrow{r}$

$${\rho }_{\sigma }\left(\overrightarrow{r}\right)=\sum _{\mu \nu }{P}_{\mu \nu }^{\sigma }{\chi }_{\mu }\left(\overrightarrow{r}\right){\chi }_{\nu }\left(\overrightarrow{r}\right)$$

(49)

as [16]

$$F_{\mu \nu }^{\alpha ;\mathrm{LDA}}=\frac{\partial E^{\mathrm{XC}}}{\partial P_{\mu \nu }^{\alpha }}=b\int \mathrm{d}\overrightarrow{r}\frac{\partial f\left(\overrightarrow{r}\right)}{\partial {\rho }_{\alpha }}\frac{\partial {\rho }_{\alpha }}{\partial P_{\mu \nu }^{\alpha }}=b\int \mathrm{d}\overrightarrow{r}\frac{\partial f\left(\overrightarrow{r}\right)}{\partial {\rho }_{\alpha }}{\chi }_{\mu }\left(\overrightarrow{r}\right){\chi }_{\nu }\left(\overrightarrow{r}\right)$$

(50)

with ${\mathbf{F}}^{\beta ;\mathrm{LDA}}$ having an analogous expression. Note that if the integral is evaluated using numerical quadrature,

$$F_{\mu \nu }^{\alpha ;\mathrm{LDA}}\approx b\sum _i{w}_i\frac{\partial f\left({\overrightarrow{r}}_i\right)}{\partial {\rho }_{\alpha }}{\chi }_{\mu }\left(\overrightarrow{r_i}\right){\chi }_{\nu }\left(\overrightarrow{r_i}\right)$$

(51)

Becke’s multigrid approach [87] and further developments thereof being the standard approach in LCAO programs, the expression of Equation (51) can be most efficiently formulated with matrix products. Storing the values of the basis functions at the quadrature points as a matrix

$$X_{\mu i}={\chi }_{\mu }\left({\overrightarrow{r}}_i\right)$$

and defining a scaled version thereof as

$$Y_{\mu i}^{\alpha }=w_i\frac{\partial f\left({\overrightarrow{r}}_i\right)}{\partial {\rho }_{\alpha }}{\chi }_{\mu }\left({\overrightarrow{r}}_i\right)$$

the Fock matrix contribution can be evaluated as simply as

$${\mathbf{F}}^{\alpha ;\mathrm{LDA}}=b\mathbf{X}{\left({\mathbf{Y}}^{\alpha }\right)}^{\mathrm{T}}$$

which is orders of magnitude faster than a simple for loop based algorithm.

11.2. GGA Functionals

The second rung of Jacob’s Ladder [28] is referred to as the Generalised Gradient Approximation [88] (GGA). Density functional approximations on this rung also depend on the derivatives of the density

$$f\left(\overrightarrow{r}\right)=f({\rho }_{\alpha },{\rho }_{\beta },{\gamma }_{\alpha \alpha },{\gamma }_{\alpha \beta },{\gamma }_{\beta \beta })$$

via the reduced gradients ${\gamma }_{\alpha \alpha }=\nabla {\rho }_{\alpha }·\nabla {\rho }_{\alpha }$, ${\gamma }_{\alpha \beta }=\nabla {\rho }_{\alpha }·\nabla {\rho }_{\beta }$, and ${\gamma }_{\beta \beta }=\nabla {\rho }_{\beta }·\nabla {\rho }_{\beta }$, with the gradient of the density $\nabla {\rho }_{\sigma }$ being determined by

$$\nabla {\rho }_{\sigma }=\sum _{\mu \nu }{P}_{\mu \nu }^{\sigma }\nabla \left[{\chi }_{\mu }\left(\overrightarrow{r}\right){\chi }_{\nu }\left(\overrightarrow{r}\right)\right]=\sum _{\mu \nu }{P}_{\mu \nu }^{\sigma }[{\chi }_{\nu }\left(\overrightarrow{r}\right)\nabla {\chi }_{\mu }\left(\overrightarrow{r}\right)+{\chi }_{\mu }\left(\overrightarrow{r}\right)\nabla {\chi }_{\nu }\left(\overrightarrow{r}\right)]$$

The GGA contribution to the Fock matrix is given by [16]

$$\begin{array}{ccc}{F}_{\mu \nu }^{\alpha ;\mathrm{GGA}}\hfill & =& \frac{\partial E^{\mathrm{XC}}}{\partial P_{\mu \nu }^{\alpha }}\hfill \\ & =& \int \mathrm{d}\overrightarrow{r}\left[\frac{\partial f\left(\overrightarrow{r}\right)}{\partial {\rho }_{\alpha }}\frac{\partial {\rho }_{\alpha }}{\partial P_{\mu \nu }^{\alpha }}+\frac{\partial f\left(\overrightarrow{r}\right)}{\partial {\gamma }_{\alpha \alpha }}\frac{\partial {\gamma }_{\alpha \alpha }}{\partial \nabla {\rho }_{\alpha }}·\frac{\partial \nabla {\rho }_{\alpha }}{\partial P_{\mu \nu }^{\alpha }}+\frac{\partial f\left(\overrightarrow{r}\right)}{\partial {\gamma }_{\alpha \beta }}\frac{\partial {\gamma }_{\alpha \beta }}{\partial \nabla {\rho }_{\alpha }}·\frac{\partial \nabla {\rho }_{\alpha }}{\partial P_{\mu \nu }^{\alpha }}\right]=F_{\mu \nu }^{\alpha ;\mathrm{LDA}}\hfill \\ & +& \int \mathrm{d}\overrightarrow{r}\left(2\frac{\partial f\left(\overrightarrow{r}\right)}{\partial {\gamma }_{\alpha \alpha }}\nabla {\rho }_{\alpha }\left(\overrightarrow{r}\right)+\frac{\partial f\left(\overrightarrow{r}\right)}{\partial {\gamma }_{\alpha \beta }}\nabla {\rho }_{\beta }\left(\overrightarrow{r}\right)\right)·[{\chi }_{\nu }\left(\overrightarrow{r}\right)\nabla {\chi }_{\mu }\left(\overrightarrow{r}\right)+{\chi }_{\mu }\left(\overrightarrow{r}\right)\nabla {\chi }_{\nu }\left(\overrightarrow{r}\right)]\hfill \end{array}$$

(52)

The $\beta$ expression can be obtained by switching $\alpha$ and $\beta$ in Equation (52). In the restricted case, $E^{\mathrm{XC}}=\int \mathrm{d}\overrightarrow{r}f(\rho ,\gamma )$ with $\gamma =\nabla \rho ·\nabla \rho$, which leads to the DFT contributions to ${\mathbf{F}}^{\mathrm{GGA}}$ given by the simpler expression

$$F_{\mu \nu }^{\mathrm{GGA}}=\int \mathrm{d}\overrightarrow{r}\left[\frac{\partial f\left(\overrightarrow{r}\right)}{\partial \rho }{\chi }_{\mu }\left(\overrightarrow{r}\right){\chi }_{\nu }\left(\overrightarrow{r}\right)+2\frac{\partial f\left(\overrightarrow{r}\right)}{\partial \gamma }\nabla \rho \left(\overrightarrow{r}\right)·\left[{\chi }_{\nu }\left(\overrightarrow{r}\right)\nabla {\chi }_{\mu }\left(\overrightarrow{r}\right)+{\chi }_{\mu }\left(\overrightarrow{r}\right)\nabla {\chi }_{\nu }\left(\overrightarrow{r}\right)\right]\right]$$

(53)

Practical implementations of Equations (52) and (53) can again be formulated using matrix products.

11.3. Meta-GGA Functionals

On the third rung on Jacob’s Ladder [28] are the meta-GGA (mGGA) approximations

$$f\left(\overrightarrow{r}\right)=f({\rho }_{\alpha },{\rho }_{\beta },{\gamma }_{\alpha \alpha },{\gamma }_{\alpha \beta },{\gamma }_{\beta \beta },{\tau }_{\alpha },{\tau }_{\beta },{\nabla }^2{\rho }_{\alpha },{\nabla }^2{\rho }_{\beta })$$

in which ${\tau }_{\sigma }$ and ${\nabla }^2{\rho }_{\sigma }$ are obtained as

$$\begin{array}{ccc}{\tau }_{\sigma }\hfill & =& \frac{1}{2}\sum _{i=1}^{N_{\sigma }}{|\nabla {\phi }_i\left(\overrightarrow{r}\right)|}^2=\frac{1}{2}\sum _{i=1}^{N_{\sigma }}\nabla {\phi }_i\left(\overrightarrow{r}\right)·\nabla {\phi }_i\left(\overrightarrow{r}\right)=\frac{1}{2}\sum _{i=1}^{N_{\sigma }}\sum _{\mu \nu }{C}_{\mu i}^{\sigma }{C}_{\nu i}^{\sigma }\nabla {\chi }_{\mu }\left(\overrightarrow{r}\right)·\nabla {\chi }_{\nu }\left(\overrightarrow{r}\right)\hfill \\ & =& \frac{1}{2}\sum _{\mu \nu }{P}_{\mu \nu }^{\sigma }\nabla {\chi }_{\mu }\left(\overrightarrow{r}\right)·\nabla {\chi }_{\nu }\left(\overrightarrow{r}\right)\hfill \end{array}$$

(54)

$$\begin{array}{ccc}{\nabla }^2{\rho }_{\sigma }\hfill & =& \sum _{\mu \nu }{P}_{\mu \nu }^{\sigma }{\nabla }^2\left[{\chi }_{\mu }\left(\overrightarrow{r}\right){\chi }_{\nu }\left(\overrightarrow{r}\right)\right]\hfill \\ & =& \sum _{\mu \nu }{P}_{\mu \nu }^{\sigma }[{\chi }_{\mu }\left(\overrightarrow{r}\right){\nabla }^2{\chi }_{\nu }\left(\overrightarrow{r}\right)+2\nabla {\chi }_{\mu }\left(\overrightarrow{r}\right)·\nabla {\chi }_{\nu }\left(\overrightarrow{r}\right)+{\chi }_{\nu }\left(\overrightarrow{r}\right){\nabla }^2{\chi }_{\mu }\left(\overrightarrow{r}\right)]\hfill \\ & =& \sum _{\mu \nu }{P}_{\mu \nu }^{\sigma }[{\chi }_{\mu }\left(\overrightarrow{r}\right){\nabla }^2{\chi }_{\nu }\left(\overrightarrow{r}\right)+{\chi }_{\nu }\left(\overrightarrow{r}\right){\nabla }^2{\chi }_{\mu }\left(\overrightarrow{r}\right)]+4{\tau }_{\sigma }\hfill \end{array}$$

(55)

The meta-GGA contributions to the Kohn–Sham–Fock matrix are straightforwardly obtained as [20]

$$\begin{array}{ccc}{F}_{\mu \nu }^{\alpha ;\mathrm{mGGA}}\hfill & =& \frac{\partial E^{\mathrm{XC}}}{\partial P_{\mu \nu }^{\alpha }}\hfill \\ & =& F_{\mu \nu }^{\alpha ;\mathrm{GGA}}+\int \mathrm{d}\overrightarrow{r}\left[\frac{\partial f\left(\overrightarrow{r}\right)}{\partial {\tau }_{\alpha }}\frac{\partial {\tau }_{\alpha }}{\partial P_{\mu \nu }^{\alpha }}+\frac{\partial f\left(\overrightarrow{r}\right)}{\partial {\nabla }^2{\rho }_{\alpha }\left(\overrightarrow{r}\right)}\frac{\partial {\nabla }^2{\rho }_{\alpha }\left(\overrightarrow{r}\right)}{\partial P_{\mu \nu }^{\alpha }}\right]\hfill \\ & =& F_{\mu \nu }^{\alpha ;\mathrm{GGA}}+\int \mathrm{d}\overrightarrow{r}\left[\left(\frac{1}{2}\frac{\partial f\left(\overrightarrow{r}\right)}{\partial {\tau }_{\alpha }}+2\frac{\partial f\left(\overrightarrow{r}\right)}{\partial {\nabla }^2{\rho }_{\alpha }\left(\overrightarrow{r}\right)}\right)\nabla {\chi }_{\mu }\left(\overrightarrow{r}\right)·\nabla {\chi }_{\nu }\left(\overrightarrow{r}\right)\right.\hfill \\ & & +\left.\frac{\partial f\left(\overrightarrow{r}\right)}{\partial {\nabla }^2{\rho }_{\alpha }\left(\overrightarrow{r}\right)}[{\chi }_{\mu }\left(\overrightarrow{r}\right){\nabla }^2{\chi }_{\nu }\left(\overrightarrow{r}\right)+{\chi }_{\nu }\left(\overrightarrow{r}\right){\nabla }^2{\chi }_{\mu }\left(\overrightarrow{r}\right)]\right]\hfill \end{array}$$

(56)

which can again be expressed in terms of matrix products to achieve faster quadrature. The expressions remain formally the same in the restricted case, but the quantities correspond to the total electron density.

11.4. Range-Separated Hybrid Functionals

As was mentioned before, the use of non-zero values for the constants a and b in Equation (7) allows the inclusion of exact exchange effects in a DFT calculation. Such functionals represent the fourth rung of Jacob’s Ladder [28], and are generally referred to as hybrids. A further development on hybrid functionals are range-separated hybrids [89,90], in which the interelectronic interaction is divided into a short-range (sr) and a long-range (lr) part with a resolution of the identity

$$\frac{1}{r_{12}}=\frac{{\varphi }^{\mathrm{sr}}\left(r_{12}\right)}{r_{12}}+\frac{{\varphi }^{\mathrm{lr}}\left(r_{12}\right)}{r_{12}}$$

(57)

where ${\varphi }^{\mathrm{sr}}\left(r_{12}\right)+{\varphi }^{\mathrm{lr}}\left(r_{12}\right)=1$. The rationale for range separation is that since density functional approximations for the exchange are based only on local information about the density, they fail to reproduce accurate estimates for, e.g., charge transfer processes. Separating the interaction by range per Equation (57) leads to a hybrid exchange functional that has four contributions

$$E^{\mathrm{X}}=a^{\mathrm{sr}}{E}^{\text{sr-HF}}+a^{\mathrm{lr}}{E}^{\text{lr-HF}}+b^{\mathrm{sr}}{E}^{\text{sr-DFT}}+b^{\mathrm{lr}}{E}^{\text{lr-DFT}}=a^{\mathrm{sr}}{E}^{\text{sr-HF}}+a^{\mathrm{lr}}{E}^{\text{lr-HF}}+E^{\mathrm{DFT}}$$

(58)

where we have stressed that since the DFT contributions are evaluated based only on the density (and possibly its derivatives), $b^{\mathrm{sr}}{E}^{\text{sr-DFT}}+b^{\mathrm{lr}}{E}^{\text{lr-DFT}}$ is nothing but a definition of a new density functional. In contrast, the HF contributions to the energy and the Kohn–Sham–Fock matrix have to be evaluated separately with range-separated ERIs

$${\left(\mu \nu \right|\lambda \sigma )}^{\mathrm{sr}/\mathrm{lr}}=\int \mathrm{d}{\overrightarrow{r}}_1\int \mathrm{d}{\overrightarrow{r}}_2{\chi }_{\mu }\left(\overrightarrow{r_1}\right){\chi }_{\nu }\left(\overrightarrow{r_1}\right)\frac{{\varphi }^{\mathrm{sr}/\mathrm{lr}}\left(r_{12}\right)}{r_{12}}{\chi }_{\lambda }\left(\overrightarrow{r_2}\right){\chi }_{\sigma }\left(\overrightarrow{r_2}\right)$$

(59)

Several kinds of range-separation kernels ${\varphi }^{\mathrm{sr}/\mathrm{lr}}\left(r\right)$ have been proposed; however, the error function based kernel ${\varphi }^{\mathrm{sr}}(r;\omega )=\mathrm{erfc}\left(\omega r\right)$, ${\varphi }^{\mathrm{lr}}(r;\omega )=\mathrm{erf}\left(\omega r\right)$, where $\omega$ is the range-separation parameter, is by far the most commonly used one because it is exceedingly simple to implement in codes employing a plane-wave or Gaussian basis set [91,92]. The error function kernel is used, for instance, in the Heyd–Scuseria–Ernzerhof (HSE) functionals for solid-state calculations [91,93], as well as in the $\omega$B97M-V [94] functional that is discussed below in Section 11.5. Some functionals based on Yukawa kernels, ${\varphi }^{\mathrm{sr}}(r;\omega )=e^{-\lambda r}$, ${\varphi }^{\mathrm{lr}}(r;\omega )=1-e^{-\lambda r}$, have also been published and are available in LibXC, for instance. It is important to check that the range-separation kernel used in the density functional implementation matches the one used in the computation of the range-separated ERIs in Equation (59), as the results will be incorrect otherwise.

11.5. Non-Local Correlation

Dispersion effects, i.e., van der Waals interactions, can be modeled in an ab initio DFT setting with non-local correlation functionals [95]

$$E^{\mathrm{nlc}}=\int \mathrm{d}{\overrightarrow{r}}_1\mathrm{d}{\overrightarrow{r}}_2\rho \left({\overrightarrow{r}}_1\right){\mathsf{\Phi }}_0({\overrightarrow{r}}_1,{\overrightarrow{r}}_2)\rho \left({\overrightarrow{r}}_2\right)$$

(60)

Because the non-local correlation energy term depends explicitly on the electron density, it also needs to be included in the SCF procedure, in principle. In contrast, empirical dispersion corrections such as Grimme’s various DFT-D approaches [96,97,98] do not depend on the electron density, and are added only as an ad hoc correction onto the electronic energy.

Perhaps the most accurate rung-3 and rung-4 functionals currently available [99,100,101], the pure B97M-V [102] mGGA as well as the range-separated $\omega$B97M-V [94] hybrid mGGA, respectively, are built on top of [103] the VV10 non-local correlation functional [23] which is controlled by two adjustable parameters, b and C, which have been trained alongside the density functional in B97M-V and $\omega$B97M-V. The results of a recent benchmark study suggest that the VV10 contributions on densities and orbital energies are negligible, and that sufficiently accurate energetics may be obtained by a one-shot evaluation of $E^{\mathrm{nlc}}$ in a post-SCF fashion [99]. Still, a rigorous minimization of the energy requires considering the effects of the non-local correlation on the wave function. Although Equation (60) does not appear to fit on the rungs of functionals discussed above, the VV10 kernel turns out to yield a GGA-type contribution to the Kohn–Sham–Fock matrix as discussed in [23], to which we refer for further details.

12. Summary and Discussion

We have presented an overview of self-consistent field calculations within a variational basis set formalism, and discussed the solution of the self-consistent field equations arising from Hartree–Fock as well as various levels of density functional approximations using either the traditional fixed-point equations or direct minimization, as well as various conceptual and numerical issues arising in their implementation. No assumptions have been made on the underlying basis set in the present work: the self-consistent field formalism is the same regardless of the form of the basis functions, which can be chosen to be, e.g., Gaussian-type orbitals (GTOs), Slater-type orbitals (STOs), numerical atomic orbitals (NAOs), or finite element shape functions. The basis set is only reflected in the molecular integrals, that is, the matrix elements of the core Hamiltonian and the two-electron integrals. The various choices for the basis set have different advantages and disadvantages, including the evaluation of the molecular integrals; see [2] for discussion. Instead, the main restriction in the present work is the implicit assumption that the basis set is compact enough so that the $N\times N$ density and Fock matrices are small enough to allow $O\left(N^2\right)$ dense matrix storage and $O\left(N^3\right)$ diagonalization. Although the present overview has focused on non-relativistic calculations on molecules, the discussion for relativistic calculations as well as crystalline systems is analogous to a large degree. We hope that the present, consistent and thorough derivation and analysis will be useful for reference as well as teaching purposes, and that the results presented herein will lead to a wider availability of density functionals in electronic structure programs.