The Kohn-Sham method

Let us introduce a reference system of non-interacting electrons ^8.5 with Kohn-Sham (KS) hamiltonian:

$$\hat{H}_{KS} ({\bf r}_1,\dots, {\bf r}_M)= \sum_{n=1}^M \hat{h}_{KS}[\rho_0] ({\bf r}_n)$$

$$\hat{h}_{KS}[\rho_0] ({\bf r}_n) = -\frac{1}{2} \nabla_n^2 + v_{KS}[\rho_0] ({\bf r}_n),$$ (8.21)

where the KS potential $v_{KS}[\rho_0] ({\bf r})$ will be obtained below, using the condition that the ground state density of the non-interacting system must be the same of that of the fully interacting system, $\rho_0$. The ground state wavefunction of this system is:

$$\Psi_{KS} = \frac{1}{\sqrt{M!}} det[\psi_1, \dots, \psi_M],$$ (8.22)

where the single particle wavefunctions $\psi_n$ are the lowest $M$ eigenstates of the single particle KS hamiltonian $\hat{h}_{KS}[\rho_0]$. If there are $N$ electrons in the system then $M = N/2$, because we can accommodate two electrons in each orbital ^8.6. Since the orbitals ${\psi_n}$ are ortho-normal, the electron density of this system is given by:

$$\rho_0({\bf r}) = 2 \sum_{n=1}^{N/2} \vert\psi_n({\bf r})\vert^2,$$ (8.23)

and the kinetic energy:

$$T_{ne}[\rho_0] = -\sum_{n=1}^{N/2} \int_V \psi_n^*({\bf r}) \nabla^2 \psi_n({\bf r}) d^3{\bf r}.$$ (8.24)

Consider now the electrostatic (Coulomb) repulsive energy $E_H[\rho]$ due to a generic charge density $\rho$:

$$E_H[\rho] = \frac{1}{2} \int_V \frac{\rho({\bf r})\rho({\bf r'})}{\vert{\bf r} - {\bf r'}\vert} d^3{\bf r}d^3{\bf r'}.$$ (8.25)

We can write the energy of the fully interacting system in terms of $T_{ne}[\rho]$ and $E_H[\rho]$:

$$E[\rho] = T_{ne}[\rho] + V_{ext}[\rho] + E_H[\rho] + E_{xc}[\rho],$$ (8.26)

which defines the exchange-correlation functional

$$E_{xc}[\rho] = T[\rho] - T_{ne}[\rho] + V_{ee}[\rho] - E_H[\rho] = F[\rho] - T_{ne}[\rho] - E_H[\rho].$$ (8.27)

Since we do not know $F[\rho]$, we do not know $E_{ex}[\rho]$ either, but the hope is that a large part of the kinetic energy of the fully interacting system is captured by $T_{ne}[\rho]$ and, similarly, a large parte of the electron-electron interaction energies are captured by $E_H[\rho]$. If this is the case, then we are left with the need to approximate a smaller quantity. We shall come back later to an approximation of $E_{ex}[\rho]$, but now we want to see how we can exploit the introduction of the non-interacting system to find the GS energy of the system. Coming back to the constrained minimisation introduced in , it is useful to recast the variation of the density in terms of variations of the single particle wavefunctions $\psi_n$. Such a variation has to be performed while keeping the wavefunctions ortho-normal, which gives a set of $M^2$ constraints:

$$\int_V d^3 {\bf r} \psi_i^*({\bf r}) \psi_j({\bf r}) = \delta_{ij}.$$ (8.28)

To do that, we define the functional:

$$\Omega[\{\psi_n\},\{\epsilon_{ij}\}] = E[\rho] - \sum_{i,j=1}^{N/... ... ( \int_V d^3 {\bf r} \psi_i^*({\bf r}) \psi_j({\bf r}) - \delta_{ij} \right ),$$ (8.29)

and impose the condition:

$$\delta \Omega[\{\psi_n\},\{\epsilon_{ij}\}] = 0.$$ (8.30)

Since the single particle wavefunctions are complex we need to consider variations of $\Omega$ w.r.t. the real and the imaginary part of $\psi_n$, and so the minimum condition requires $\left ( \frac{\delta \Omega}{\delta \psi_n^r}\right )_{\psi_n^i}=0$ and $\left ( \frac{\delta \Omega}{\delta \psi_n^i}\right )_{\psi_n^r}= 0$, with each partial derivatives taken while holding all other wavefunctions constant. However, the functional $T_{ne}[\rho]$ is written in terms of $\psi_n$ and $\psi^*_n$, and so it would be more useful to be able to consider variations w.r.t to $\psi_n^*$ and $\psi_n$ rather than $\psi_n^r$ and $\psi_n^i$. Using the usual chain rule for the derivative, we can write:

$$\left(\frac{\delta \Omega}{\delta \psi_n^*}\right )_{\psi_n} = \l... ...t )_{\psi_n^r} \left(\frac{\delta \psi_n^i}{\delta \psi_n^*} \right )_{\psi_n},$$ (8.31)

and since

$$\psi_n^r = \frac{1}{2} (\psi_n + \psi_n^*); \quad \quad \quad \psi_n^i = \frac{1}{2 i} (\psi_n - \psi_n^*),$$ (8.32)

we have

$$\left(\frac{\delta \Omega}{\delta \psi_n^*}\right )_{\psi_n} = \f... ...i} + i \left(\frac{\delta \Omega}{\delta \psi_n^i}\right )_{\psi_n^r} \right ].$$ (8.33)

Since $\Omega$ is real, the condition $\left(\frac{\delta \Omega}{\delta \psi_n^*}\right )_{\psi_n} = 0$ satisfies at once both conditions $\left ( \frac{\delta \Omega}{\delta \psi_n^r}\right )_{\psi_n^i}=0$ and $\left ( \frac{\delta \Omega}{\delta \psi_n^i}\right )_{\psi_n^r}= 0$. Of course, we also need to impose the conditions $\left ( \frac{\delta \Omega}{\delta \epsilon_{ij}}\right ) = 0$, which ensures that the constraints are satisfied. To compute these variations, consider the functional $\Omega$ explicitly:

$$\Omega [\{\psi_n\},\{\epsilon_{ij}\}] = - \sum_{n=1}^{N/2} \int_V... ...psi_n({\bf r'}) d^3{\bf r'} + \int_V v({\bf r'}) \rho ({\bf r'}) d^3 {\bf r'} +$$

$$\frac{1}{2} \int_V \frac{\rho({\bf r'})\rho({\bf r''})}{\vert{\bf... ...\int_V \psi_i^*({\bf r'}) \psi_j({\bf r'}) d^3 {\bf r'} - \delta_{ij} \right ).$$ (8.34)

Considering point variations of the type $\left(\frac{\delta G}{\delta \psi_n^*({\bf r})}\right )_{\psi_n}$, which select the integrands evaluated at point ${\bf r}$ ^8.7, we obtain (see Appendix ):

$$\left [- \nabla^2 + 2v({\bf r}) + 2\int_V \frac{\rho({\bf r'})}{\... ...r})} \right ] \psi_n({\bf r}) = 2\sum_{j=1}^{N/2}\epsilon_{nj} \psi_n({\bf r}).$$ (8.35)

Dividing everything by 2 and setting:

$$v_{KS}[\rho] ({\bf r}) = v({\bf r}) + \int_V \frac{\rho({\bf r'})... ... r} - {\bf r'}\vert}d^3 {\bf r'} + \frac{\delta E_{xc}}{\delta \rho ({\bf r})},$$ (8.36)

we can rewrite Eq.  as:

$$\hat{h}_{KS}[\rho]({\bf r}) \psi_n({\bf r}) = \left [ - \frac{1}{... ...f r}) \right ] \psi_n({\bf r}) = \sum_{j=1}^{N/2}\epsilon_{nj} \psi_j({\bf r}),$$ (8.37)

which is a linear system of $N/2$ equations. Since the constraint requires $\langle \psi_i \vert \psi_j \rangle = \delta_{ij}$, the Lagrange multipliers $\epsilon_{ij}$ are obtained from

$$\epsilon_{ij} = \langle \psi_i\vert \hat{h}_{KS}[\rho] \vert \psi_j \rangle,$$ (8.38)

which is a hermitian matrix. It always possible, therefore, to find a unitary transformation of the ${\psi_n}$ such that the linear system takes the diagonal form:

$$\left [ - \frac{1}{2}\nabla^2 + v_{KS}[\rho] ({\bf r}) \right ] \psi_n({\bf r}) = \epsilon_n \psi_n({\bf r}),$$ (8.39)

which is a system of single particle Schrödinger like equations and are known as the Kohn-Sham equations. They are coupled via the effective potential $v_{KS}[\rho]$, as $\rho$ depends on all the ${\psi_n}$. The unitary transformation that takes us from  to  is the one that diagonalises the matrix $\epsilon_{ij}$. The KS equations show that the extremes of the functional $\Omega$ are obtained for any ensemble of $N/2$ eigenstates of $\hat{h}_{KS}[\rho]$, and the ground state is obtained by finding its minimum of the total energy w.r.t. any choice of $N/2$ state. In practice one always takes the lowest $N/2$ states, although this may not necessarily be the correct choice (see below). One way to do this in practice is to introduce a set of $L$ basis functions $\{\phi_m\}$, construct the matrix $\epsilon_{ij} = \langle \phi_i \vert \hat{h}_{KS}[\rho] \vert \phi_j \rangle$ and diagonalise it. Note the subtle difference with the discussion above that brought us from  to . In that case we used the diagonalisation procedure to rewrite the linear system as an ensemble of independent (though coupled via $\rho$) equations. The diagonalisation was performed in the subspace spanned by the $N/2$ eigenvectors of the hamiltonian to write the KS equations in their canonical form. Here we are using diagonalisation in the whole space spanned by the $L$ basis functions to find the lowest $N/2$ eigenstates of the hamiltonian, and with those build the ground state of the full system. The eigenvectors are of the type:

$$\psi_n = \sum_{i=1}^L c_i^n \phi_i.$$ (8.40)

If the basis set $\{\phi_l\}$ is complete the solution is exact. However, usually this means including an infinite number of elements, which cannot be done in practice, therefore the solution to the KS equations is only approximate. By increasing the number of basis functions one can drive the calculations to convergence, where the energy and other properties are obtained within some predefined threshold. Note that the usual variational principle applies, and so by including more and more elements in the basis set the ground state energy decreases monotonically. The rate of decrease of the energy can be used to judge the level of convergence.

Since the KS potential depends on $\rho$, Eqs.  have to be solved self-consistently. This is typically done by iteration, in which one starts with some initial guess for the electron density $\rho_1$ (a superposition of atomic densities, for example), constructs $v_{KS}[\rho_1]$ and solves . With those solutions construct $\rho_2$ using  ^8.8, and solve again , this time using the newly constructed $v_{KS}[\rho_2]$. Keep repeating until the difference between $\rho_{j+1}$ and $\rho_j$ is below some acceptable threshold.

If we multiply  by $\psi_n^*({\bf r})$, integrate over ${\bf r}$, sum over $n$ and compare with Eq.  and , we see that the total energy can be obtained from:

$$E = 2\sum_{n=1}^{N/2} \epsilon_n - \frac{1}{2} \int_V \frac{\rho(... ... \frac{\delta E_{xc}}{\delta \rho ({\bf r})}\rho ({\bf r})d^3 {\bf r} + E_{xc}.$$ (8.41)

The condition that guarantees that the first $N/2$ states give the minimum energy in  is known as $v-representability$, that is the ground state density of the interacting system is the same as the ground state density of $some$ non-interactive system. This also implies that the total energy  would be exact if one had the exact form of the exchange-correlation functional $E_{xc}$. It is possible, however, that the ground state density of a system is not $v-representable$, i.e. there is no non-interacting system whose ground state density is the same as that of the interacting system. In this case we need to go back to the variational principle, and the $N/2$ states in the non-interacting system that minimise the energy may not be the first $N/2$. ^8.9