# Ion-ion interaction potential in Kohn-Sham DFT

The Kohn-Sham equation as described in "Density Functional Theory: A Practical Introduction" by Dr David Sholl is:

$$\tag{1}\left[-\frac{\hbar^2}{2m}\nabla^2+V({\bf r})+V_H({\bf r})+V_{XC}({\bf r})\right]\psi_i({\bf r})=\varepsilon_i\psi_i({\bf r}).$$ The first term on the left-hand side represents the kinetic energy of the electrons. The second term represents electron-ion interaction and the third represents electron-electron interaction. The last term is the exchange-correlation potential.

With the Born-Oppenheimer approximation we neglect the kinetic energy of ions. What about the ion-ion interaction potential?

In lattice dynamics studies, we account specifically for the ion-ion interaction but with minute oscillations of the ions. I believe DFT would still account for the interaction among stationary ions in the static lattice viewpoint.

• Could'nt the ion-ion interaction being a constant, be absorbed into the $\epsilon_{i}$ term. Dec 29 '20 at 21:31
• Hi @AnoopANair, that is the first reason I got from a friend. More in the lines of it being constant, hence neglected. I don't see why it should be neglected because it's constant. What exactly do you mean when you say, "absorbed in the $\epsilon_i$ term"? Dec 29 '20 at 23:09
• The eigenvalue of a constant is just itself, so a constant term $C$ in a Hamiltonian, will introduce a term $C\psi$ on the left, and an identical term $C\psi$ on the right. The right side then becomes $(C + \epsilon = \tilde\epsilon)\psi$ where $\tilde\epsilon$ is just $\epsilon$ but with $C$ absorbed into it (additively). Dec 29 '20 at 23:47
• Hi I hope you've got your answer :). You can check this : [mattermodeling.stackexchange.com/questions/683/… ], its like a rephrased form of your question. Dec 30 '20 at 7:35
• @AnoopANair, understood, thank you. I checked the VASP output and they write: "FREE ENERGIE OF THE ION-ELECTRON SYSTEM (eV)" before they report the total energy. Does this mean they don't add the constant ion-ion interaction at the end? Dec 30 '20 at 10:19

If the ion-ion interaction contributes a constant term to the Hamiltonian $$H$$, then our new Hamiltonian is $$H+C$$. The eigenvalue of a constant is just itself, so we have:

$$\tag{1} (H + C )\psi = (\epsilon + C)\psi$$

So if your DFT code calculates only $$\epsilon$$ (the energy if you neglect the ion-ion interaction), it is easy to obtain the energy with the ion-ion interaction by simply adding the constant $$C$$, which is something that does not need a complicated DFT code. The DFT code can easily add the energy coming from the ion-ion interaction at the end of the calculation in the same way that things like the nuclear-nuclear repulsion energy might be added in a molecular quantum chemistry software.

The matter can be viewed as a set of ions and electrons. The Kohn-Sham equation listed in your post aims to solve the electronic part. As for the ionic part, which is usually treated classically in the framework of Newton's mechanics. The ion-ion potential or force can be calculated with the empirical method (classical molecular dynamics) or the first-principles method (ab-initio molecular dynamics).

Within the first-principles method, the total energy of the system is calculated with density functional theory then the force is calculated by energy derivative.

I would like to emphasise a few aspects that seem to be a little bit between the lines in the other answers.

Density functional theory is based on the fact that observables of an interacting-electron system can in principle be obtained from its ground-state electron density. The Kohn-Sham system is a means of obtaining this density (and a few other objects that make certain calculations more reasonable). Obviously the interaction between the nuclei does not directly affect the ground-state electron density and therefore it is not required to include this interaction directly in the Kohn-Sham system$$^1$$.

Nevertheless this interaction is very important when calculating the total energy of a system. For a system with a unit cell $$\Omega$$ containing atoms with core charges $$Z_\alpha$$ at $$\mathbf{\tau}_\alpha$$ and featuring a spin-dependent ground-state electron density $$\rho^\sigma$$ and Kohn-Sham eigenvalues $$E_{\nu,\sigma}$$ the total energy functional is

\begin{align} E_\text{total}[\rho^\uparrow,\rho^\downarrow] &= \underbrace{\left[\sum\limits_\sigma \left(\sum\limits_{\nu=1}^{N_\text{occ}^\sigma} E_{\nu,\sigma}\right) - \int\limits_{\Omega} \rho^\sigma(\mathbf{r}) V_{\text{eff},\sigma}(\mathbf{r}) d^3 r \right]}_{E_\text{kin}}\nonumber \\ &\phantom{=} + \underbrace{\frac{1}{2}\int\limits_{\Omega}\int\limits_{\Omega}\frac{\rho(\mathbf{r})\rho(\mathbf{r}')}{\vert\mathbf{r}-\mathbf{r}'\vert} d^3rd^3r' + \int\limits_{\mathbb{R}^3\backslash \Omega}\int\limits_{\Omega}\frac{\rho(\mathbf{r})\rho(\mathbf{r}')}{\vert\mathbf{r}-\mathbf{r}'\vert} d^3rd^3r'}_{E_\text{H}} \\ &\phantom{=} + \underbrace{\int\limits_{\Omega} V_\text{ext}(\mathbf{r}) \rho(\mathbf{r})d^3r \nonumber}_{E_\text{ext}} + E_\text{xc}[\rho^\uparrow,\rho^\downarrow] \\ &\phantom{=} + \underbrace{\frac{1}{2}\sum\limits_{\alpha \in \Omega}^{N_\text{atom}} \sum\limits_{\substack{\beta \in \Omega \\ \alpha\neq \beta}}^{N_\text{atom}} \frac{Z_\alpha Z_\beta}{\vert\mathbf{\tau}_\alpha - \mathbf{\tau}_\beta\vert} + \sum\limits_{\alpha \not\in \Omega} \sum\limits_{\beta \in \Omega}^{N_\text{atom}} \frac{Z_\alpha Z_\beta}{\vert\mathbf{\tau}_\alpha - \mathbf{\tau}_\beta\vert}}_{E_\text{II}}. \end{align}

In this expression $$E_\text{kin}$$ denotes the kinetic energy of the occupied Kohn-Sham orbitals, $$E_\text{H}$$ the Hartree energy, $$E_\text{ext}$$ the energy due to the interaction between the electrons and the external potential, $$E_\text{XC}$$ the exchange-correlation energy, and $$E_\text{II}$$ the energy due to the Coulomb interaction between the ionized atomic nuclei.

By having a look at this expression two properties directly become obvious:

1. $$E_\text{II}$$ gives an energy contribution that depends on the coordinates of the atomic nuclei relative to each other. This term therefore is important when calculating forces $$\mathbf{F}_\alpha = -\frac{\delta E_\text{total}}{\delta \mathbf{\tau}_\alpha}$$ and also when only relating different structures to each other that have slightly different atom distances, e.g., when calculating a lattice constant.
2. For periodic systems like crystals $$E_\text{H}$$, $$E_\text{ext}$$, and $$E_\text{II}$$ each are divergent. This is because of the long range of the Coulomb interaction together with the inclusion of contributions from the whole space outside the unit cell. These energy contributions only become finite when combined. For such systems neglecting $$E_\text{II}$$ therefore would result in a divergent total energy for the unit cell. Also care has to be taken to evaluate these contributions such that intermediate results don't diverge. A similar divergence arises if the periodically repeated unit cell is not charge neutral. Such a situation would lead to an infinite charge in the whole crystal implying an infinite electrostatic energy.

Taking into account the ion-ion interaction within a DFT procedure therefore is essential, not optional. But you will not see it explicitly in the Kohn-Sham equations.

[1] Of course, the issue of divergent contributions for infinite setups also has to be taken care of in the Kohn-Sham system.

• +1. I'm very glad you wrote this answer, and you did a phenomenal job with the equation. There is indeed a "nuclear DFT" too, but it's not so mainstream and it's probably not needed for this user's application. Dec 30 '20 at 18:36
• Thank you for the detailed answer. I couldn't understand the 2nd point about each of those energies being divergent. Can you explain that in a bit more detail? Dec 30 '20 at 19:02
• @Hitanshu Sachania: Actually I had to correct the equation slightly, but this change also makes the problem more obvious. I also added a few explaining sentences. In crystals there are also interactions between the charges in the unit cell with charges outside the unit cell. Since the Coulomb interaction has a very long range, this can lead to such a divergent behavior. One has to make sure that in summing everything up the contributions from outside the unit cell are charge neutral, otherwise the result will not be finite. If charge neutrality is present the divergent parts cancel each other. Dec 30 '20 at 19:46
• Thank you, @GregorMichalicek. Periodicity makes crystals quite interesting but also complicated, haha. Jan 2 '21 at 16:27