How is the tight binding model derived from the Kohn-Sham DFT energy?

Question

I am following along a publication titled "Density-functional tight-binding for beginners" (P. Koskinen, V. Mäkinen, Computational Materials Science 2009, 47, 237–253) and I get stomped by a step in their derivation.

The authors start with the following expression for the Kohn-Sham DFT energy: $$ E[n]=\sum_{a}f_{a}\langle\psi_{a}|\left(-\frac{1}{2}\,\nabla^{2}+\int V_{\mathrm{ext}}(r)n(r)dr\right)|\psi_{a}\rangle + \frac{1}{2}\int\int\frac{n n^{\prime}}{\left|r-r^{\prime}\right|}dr^{\prime}dr+E_{x c}[n]+E_{I I} $$ where $f_a \in [0, 2]$ is the occupation of a single-particle state $\psi_a$, $r$ is a position vector in 3D, $n(r)$ is the function of electron density and $V_{ext}(r)$ is an external potential. $E_{xc}[n]$ is the exchange-correlation functional which has an unknown form at this point and $E_{II}$ is the internuclear repulsion which can be taken as constant.

I believe that the second term in the first parenthesis should instead be just $V_{ext}(r)$, but I could be wrong.

The authors then write the following:

Consider a system with density $n_0(r)$ that is composed of atomic densities, as if atoms in the system were free and neutral. Hence $n_0(r)$ contains (artificially) no charge transfer. The density $n_0(r)$ does not minimize the functional $E[n(r)]$, but neighbors the true minimizing density $n_{min}(r) = n_0(r) + \delta n_0(r)$, where $\delta n_0(r)$ is supposed to be small. Expanding $E[n]$ at $n_0(r)$ to second order in fluctuation $\delta n(r)$ the energy reads:

$$ E[\delta n]\approx\sum_{a}f_{a}\langle\psi_{a}|-\frac{1}{2}\nabla^{2}+V_{\mathrm{ext}}+V_{H}[n_{0}]+V_{x c}[n_{0}]|\psi_{a}\rangle \\ +\frac{1}{2}\int\int\biggl(\frac{\delta^{2}E_{x c}[n_{0}]}{\delta n\delta n^{\prime}}+\frac{1}{|r-r^{\prime}|}\biggr)\delta n\delta n^{\prime}dr^{\prime}dr-\frac{1}{2}\int V_{H}[n_{0}](r)n_{0}(r)dr \\ +\,E_{x c}[n_{0}]+E_{I I}-\int V_{x c}[n_{0}]n_{0}(r)dr $$

They note that the linear terms in $\delta n$ vanish, which I suppose is because we are considering the ground-state density, and thus the derivative of the energy with regards to the density is 0 since we are at a minimum.

The first term and $E_{II}$ form the zeroth order term, simple enough. I also managed to derive the second order term (the double integral, the second term). However, I cannot figure out where the three other terms come from and what they mean. Isn't $E_{xc}[n_0] = \int V_{xc}[n_0]n_0(r)dr$? And why is there another term for $E_{H}$ which is equivalent (?) to the one in the zeroth order term, after some manipulations?

Answer 1

The additional terms come from the Kohn-Sham equations. Everything is detailed quite nicely in W. M. C. Foulkes, R. Haydock, Phys. Rev. B 1989, 39, 12520–12536. Here is a summary of mine:

We start by considering a fictitious system of non-interacting electrons with the same groundstate density $n(r)$ as the real system. From the Schrödinger equation, we know the energy of each electron in the fictitious system: $$ \left( \frac{1}{2}\nabla^2 + V(r) \right)\psi_i(r) = \epsilon_i \psi_i(r) $$ where $V(r)$ is the potential which yields the groundstate density $n(r)$ in this non-interacting system. The total energy is then simply the sum of one-electron energies $\epsilon_i$. We now consider this case from a DFT perspective: $$ E[n] = T_s[n] + F[n] $$ where $T_s[n]$ is the kinetic energy functional for non-interacting electrons and $F[n]$ is the energy associated with the potential $V(r)$. In the first equation, we can left-multiply both sides by $\psi_i^*(r)$ and derive the following relationship:

$$ \int \psi_i^*(r)\frac{1}{2}\nabla^2\psi_id^3r = T_s[n] = \sum_{i=1}^N\epsilon_i - \int V(r)n(r) d^3r $$

Then, we consider the condition that we want a groundstate density, meaning that $\delta E[n]/\delta n(r) = 0$. Starting from the second equation:

$$ \begin{align} \int \frac{\delta E[n]}{\delta n(r)}\delta n(r) d^3r &= \int\frac{\delta T_s[n]}{\delta n(r)}\delta n(r) d^3r + \int\frac{\delta F[n]}{\delta n(r)}\delta n(r) d^3r \approx 0 \\ &= \int \left(-V(r) + \frac{\delta F[n]}{\delta n(r)}\right)\delta n(r) d^3r \approx 0 \\ \end{align} $$

We conclude that $V(r) = \delta F[n]/\delta n(r) = V_H[n](r) + V_{ext}(r) + \mu_{xc}[n](r)$ for this to be true. Putting it all together:

$$ \begin{align} E_0[n_0] &= \sum_{i=1}^N\epsilon_i - \int V(r)n_0(r)d^3r + F[n_0(r)] \\ &= \sum_{i=1}^N\epsilon_i - \int \left.\frac{\delta F}{\delta n}\right|_{n_0(r)} n_0(r)d^3r + F[n_0(r)] \\ &= \sum_{i=1}^N\epsilon_i - \int \left(V_{ext}(r) + V_H[n_0(r)] + \mu_{xc}[n_0(r)]\right)n_0(r)d^3r + E_{ext}[n_0(r)] + E_H[n_0(r)] + E_{xc}[n_0(r)] \\ &= \sum_{i=1}^N\epsilon_i - \int V_H[n_0(r)]n_0(r)d^3r - \int\mu_{xc}[n_0(r)]n_0(r)d^3r + E_H[n_0(r)] + E_{xc}[n_0(r)] \\ &= \sum_{i=1}^N\epsilon_i - 2E_H[n_0(r)] - \int\mu_{xc}[n_0(r)]n_0(r)d^3r + E_H[n_0(r)] + E_{xc}[n_0(r)] \\ &= \sum_{i=1}^N\epsilon_i - E_H[n_0(r)] - \int\mu_{xc}[n_0(r)]n_0(r)d^3r + E_{xc}[n_0(r)] \\ \end{align} $$

And this is the true zeroth order term in the functional expansion. The second order terms are added to this to yield the final expression.