Charged systems with density functional theory

On this page, we describe technical issues with computing the energies of charged systems with periodic density functional theory (DFT) calculations. We then discuss why the energies of charged systems diverge for systems with lower dimensionality, such as with surfaces (2D), nanowires (1D) and molecules (0D) while potentially providing useful information for bulk (3D) systems. Finally, we present methods implemented in VASP which allow for calculations of charged 3D, 2D and 0D systems.

Treating divergence in charge neutral calculations

VASP makes use of efficient fast Fourier transforms (FFT) to compute the electrostatic potential from the charge density using the Poisson equation,

$$V(\mathbf{r}) = 4\pi \int \frac{\rho(\mathbf{r}^\prime)}{\left | \mathbf{r} - \mathbf{r}^\prime \right|} d\mathbf{r}^\prime,$$

where ${\textstyle \mathbf{r}}$ and ${\textstyle \mathbf{r}^\prime}$ are all points in real space. Fourier transforming the Poisson equation to reciprocal space,

$$V(\mathbf{G}) = \frac{4\pi}{\mathrm{G}^2} \rho(\mathbf{G}),$$

where ${\textstyle \mathbf{G}}$ is the reciprocal lattice vector and ${\textstyle \mathrm{G}}$ is its norm.

An obvious issue with computing ${\textstyle V(\mathbf{G})}$ is that it diverges for ${\textstyle \mathrm{G}\to 0}$. This divergence is handled in charge neutral density DFT calculations by cancelling out individual divergences for the electron-electron, ion-electron and ion-ion energies as can be seen by explicitly writing out their functional forms,

$$E_{\mathrm{electron-electron}} = \frac{1}{2} \sum_{\mathrm{G}} \rho_{\mathrm{elec}}(\mathbf{G}) V_{\mathrm{elec}}(\mathbf{G}),$$

where ${\textstyle \rho_{\mathrm{elec}}}$ and ${\textstyle V_{\mathrm{elec}}}$ are the electronic charge density and potential respectively.

$$E_{\mathrm{ion-electron}} = -\sum_{\mathrm{G}} \rho_{\mathrm{ion}}(\mathbf{G}) V_{\mathrm{elec}}(\mathbf{G}),$$

where ${\textstyle \rho_{\mathrm{ion}}}$ is the ion charge density. VASP does not explicitly compute ${\textstyle E_{\mathrm{ion-electron}}}$, but the potential of the ion is reflected through the eigenvalues. The ion-ion interactions are treated by Ewald summation

$$E_{\mathrm{ion-ion}} = E_{\mathrm{long-range}} + E_{\mathrm{short-range}} + E_{\mathrm{self}} + E_{\mathrm{homogeneous}},$$

where ${\textstyle E_{\mathrm{long-range}}}$ is the only component which sums over the ${\textstyle \mathrm{G}}$ vectors and has the form,

$$E_{\mathrm{ion-ion}} = \frac{1}{2} \sum_{\mathrm{G}} \rho_{\mathrm{ion}}(\mathbf{G}) V_{\mathrm{ion}}(\mathbf{G}).$$

The ${\textstyle \mathrm{G}=0}$ terms of ${\textstyle E_{\mathrm{electron-electron}}}$, ${\textstyle E_{\mathrm{ion-electron}}}$ and ${\textstyle E_{\mathrm{ion-ion}}}$ are given by,

$$\frac{1}{2}\rho_{\mathrm{elec}}(\mathbf{G}=0) V_{\mathrm{elec}}(\mathbf{G}=0) - \rho_{\mathrm{ion}}(\mathbf{G}=0) V_{\mathrm{elec}}(\mathbf{G}=0) + \frac{1}{2}\rho_{\mathrm{ion}}(\mathbf{G}=0) V_{\mathrm{ion}}(\mathbf{G}=0),$$

which, can be expressed through the reciprocal space Poisson equation as,

$$\frac{1}{\mathrm{G}^2} \left [\frac{1}{2}\rho_{\mathrm{elec}}(\mathbf{G}=0)^2 - \rho_{\mathrm{ion}}(\mathbf{G}=0)^2 + \frac{1}{2}\rho_{\mathrm{ion}}(\mathbf{G}=0)^2 \right].$$

Under the condition that ${\textstyle \rho_{\mathrm{elec}}(\mathbf{G}=0)=\rho_{\mathrm{ion}}(\mathbf{G}=0)}$, the individual divergences would cancel out and the total energy is convergent. Recall that the <math display="inline">\mathrm{G}=0<\math> term of a Fourier transformed quantity is nothing by its average which implies that the above condition is satisfied when the average electron and ion charge density are the same, i.e. the system is charge neutral.