Ion-ion interaction potential in Kohn-Sham DFT

Question

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.

Answer 1

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.

Answer 2

Add more information to @Nike Dattani's answer:

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.

Answer 3

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.