## The Kohn-Sham equation

The Kohn-Sham idea is the following: since we know that the energy is a functional of the density, we should isolate the parts of the potential that we do know from those that we don't know. Kohn and Sham proposed to separate the electron interaction in to two parts (the first one being like the classical Coulomb-type interaction, called the Hartree potential):

\begin{align} \hat{v}_{ee} = \int _V \frac{\rho (\vec{r})}{|\vec{r}-\vec{r}'|} \mathrm{d}\vec{r} + v_{unknown} \end{align}and then calculate the kinetic energy from single-particle SchrÃ¶dinger equations:

\begin{align} \bigg( -\frac{1}{2}\nabla ^2 + \int _V \frac{\rho (\vec{r})}{|\vec{r}-\vec{r}'|} \mathrm{d}\vec{r} + v_{unknown} + t_{unknown}\bigg) \phi _i (\vec{r}) = \epsilon _i \phi _i (\vec{r}). \end{align}The terms \(v_{unknown}\) and \(t_{unknown}\) are the unknown part of the electron-electron interaction and the unknown part of the kinetic energy, respectively. There is an unknown part about the kinetic energy since we have set up single-particle equations, and the real system has many-body wave functions instead of the single-particle ones. Now we call this unknown term \(v_{xc}(\vec{r}) = v_{unknown}+t_{unknown}\). The technical name is "exchange-correlation functional".

Now, since \(v_{xc}\) as well as the Hartree potential require knowledge of the electron density, but the electron density has to be solved from the wave functions, we seem to have something of a chicken-egg problem.

Fortunately, there is a method called "iterative solution"\ that can help us here. The idea is quite simple. First, we guess some electron density \(\rho (\vec{r})\) (a good choice would be the density of the non-interacting system). Then, we solve the Kohn-Sham equations, one for each electron pair. We then calculate the new electron density by

\begin{align} \rho (\vec{r} ) = \sum _{i=1}^N 2|\phi _i (\vec{r})|^2 \end{align}The factor 2 comes from the fact that we can occupy each orbital with two electrons. Two electrons can't occupy the same state, but, assuming there are no spin-interactions, you can put one electron spin up and one spin down in the orbital.

How does the iteration work? In some instances, it is possible to prove that an iterative method will converge to the correct solution. In practical DFT calculations, though, such a proof is generally impossible. In practice, the methods for aiding convergence are understood well enough that success is guaranteed in most systems. To give you an idea about how this works, we will solve the equation \(x^2+2x=5\) iteratively, starting from the initial guess \(x_0 = 1\). The new value is then computed from \(x_{n+1} = \sqrt{5-2x_n}\):

\begin{align} x_2 = \sqrt{5-2\cdot 1} &= 1.732 \\ x_3 = \sqrt{5-2\cdot 1.732} &= 1.239 \\ x_4 = \sqrt{5-2\cdot 1.239} &= 1.588 \\ x_5 = \sqrt{5-2\cdot 1.588} &= 1.351 \\ x_6 = \sqrt{5-2\cdot 1.351} &= 1.516 \\ x_7 = \sqrt{5-2\cdot 1.516} &= 1.403 \\ \vdots \end{align}Continuing this process, we see that the iteration converges towards the right answer, \(x \approx 1.449\). Notice that the value first overshoots, then goes too low, and then overshoots again. We can dampen this "oscillation"\ about the correct value by using the formula \(x_{n+1} = \alpha (x_n) + (1-\alpha ) \sqrt{5-2\cdot x_n }\) with \(\alpha \in [0,1]\), which takes as the new value some sort of combination of the current and previous iterations. This is a trick that is almost always used in real DFT programs.

In this book, we will again first deal with the one-dimensional version of the KS equation; we will once again return to the particle-in-a-box problem. And instead of using the Hartree potential for electron-electron interaction, we'll use

\begin{align} v_{ee} = \frac{\lambda}{4}\rho (x), \end{align}because the Hartree potential isn't really well-defined for 1D systems. This potential is supposed to also contain, in an approximate way, the exchange-correlation of the system. \(\lambda \) is a parameter we can take to be any value we like.

It doesn't really matter for the structure of our program what the potential is, in principle. When the systems get more complicated, the potential also becomes more laborious to calculate, but the basics remain the same.