Electrostatic Energy

The electrostatic energy of a system of nuclear charges $Z_{I}$ at positions ${\bf R}_{I}$ and an electronic charge distribution $n({\bf r})$ consists of three parts: the Hartree energy of the electrons, the interaction energy of the electrons with the nuclei and the internuclear interactions

$$E_{\rm ES} = {1 \over 2} \int \! \! \int \! d{\bf r} d{\bf r'} { n({\bf r}) n({\bf r'}) \over \vert{\bf r} - {\bf r'}\vert}$$

$$+ \sum_{I} \int \! d{\bf r} V^{I}_{\rm core}({\bf r}) n({\bf r... ...m_{I \neq J} {Z_{I} Z_{J} \over \vert{\bf R}_{I} - {\bf R}_{J}\vert} \enspace .$$ (82)

The Ewald method (see e.g. Ref. [12]) can be used to avoid singularities in the individual terms when the system size is infinite. In order to achieve this a Gaussian core charge distribution associated with each nuclei is defined

$$n^{I}_{\rm c}({\bf r}) = - {Z_{I} \over \left( {\rm R}_{I}^{\rm c... ... {\bf r} - {\bf R}_{I} \over {\rm R}_{I}^{\rm c}} \right)^2 \right] \enspace .$$ (83)

It is convenient at this point to use a special definition for the core potential and define it to be the potential of the Gaussian charge distribution of Eq. (83)

$$V^{I}_{\rm core}({\bf r}) = \int \! d{\bf r'} {n^{I}_{\rm c}({\bf... ...vert{\bf r} - {\bf r'}\vert} = - {Z_{I} \over \vert{\bf r} - {\bf R}_{I}\vert} \left[ { \vert{\bf r} - {\bf R}_{I}\vert \over {\rm R}_{I}^{\rm c}} \right] \enspace ,$$ (84)

where erf is the error function. This potential has the correct long range behavior but we will have to add a correction potential for the short range part. The interaction energy of this Gaussian charge distributions is now added and subtracted from the total electrostatic energy

$$E_{\rm ES} = {1 \over 2} \int \! \! \int \! d{\bf r} d{\bf r'} { n({\bf ... ... { n_{\rm c}({\bf r}) n_{\rm c}({\bf r'}) \over \vert{\bf r} - {\bf r'}\vert}$$

$$+ \int \! \! \int \! d{\bf r} d{\bf r'} { n_{\rm c}({\bf r}... ... c}({\bf r}) n_{\rm c}({\bf r'}) \over \vert{\bf r} - {\bf r'}\vert} \enspace ,$$ (85)

where $n_{\rm c}({\bf r}) = \sum_{I} n^{I}_{\rm c}({\bf r})$ . The first three terms can be combined to the electrostatic energy of a total charge distribution $n_{\rm tot}({\bf r}) = n({\bf r}) + n_{\rm c}({\bf r})$ . The remaining terms are rewritten as a double sum over nuclei and a sum over self-energy terms of the Gaussian charge distributions

$$E_{\rm ES} = {1 \over 2} \int \! \! \int \! d{\bf r} d{\bf r'} { n_{\rm tot}({\bf r}) n_{\rm tot}({\bf r'}) \over \vert{\bf r} - {\bf r'}\vert}$$

$$+ {1 \over 2} \sum_{I \neq J} {Z_{I} Z_{J} \over \vert{\bf R}_{I} - {\bf R}_{J}\vert} \left[{ \vert{\bf R}_{I} - {\bf R}_{J}\vert \over \sqrt{{{\rm R}_... ...ght] - \sum_{I} {1 \over \sqrt{2 \pi}} {Z_{I}^2 \over R_{I}^{\rm c}} \enspace ,$$ (86)

where erfc denotes the complementary error function.

For a periodically repeated system the total energy per unit cell is derived from the above expression by using the solution to Poisson's equation in Fourier space for the first term and make use of the rapid convergence of the second term in real space. The total charge is expanded in plane waves with expansion coefficients

$$n_{\rm tot}({\bf G}) = n({\bf G}) + \sum_{I} n_{\rm c}^{I}({\bf G}) S_{I}({\bf G})$$ (87)

$$= n({\bf G}) - {1 \over \Omega} \sum_{I} {Z_{I} \over \sqrt{4 \pi}}... ...\! \left[ - {1 \over 2} G^2 {R_{I}^{\rm c}}^2 \right] S_{I}({\bf G}) \enspace .$$ (88)

The structure factor of an atom is defined by

$$S_I({\bf G}) = \exp [-i {\bf G}\cdot {\bf R}_I] \enspace .$$ (89)

This leads to the electrostatic energy for a periodic system

$$E_{\rm ES} = 2 \pi \Omega \sum_{{\bf G} \neq 0} { \vert n_{\rm tot}({\bf G})\vert^2 \over G^2} + E_{\rm ovrl} - E_{\rm self} \enspace ,$$ (90)

where

$$E_{\rm ovrl} = \mathop{{\sum}'}_{I,J} \sum_{\bf L} {Z_{I} Z_{J} \over \vert{\bf R}_{I} - {\bf R}_{J} - {\bf L}\vert} \left[{ \vert{\bf R}_{I} - {\bf R}_{J} - {\bf L}\vert \over \sqrt{{{\rm R}_{I}^{\rm c}}^2+{{\rm R}_{J}^{\rm c}}^2} } \right]$$ (91)

and

$$E_{\rm self} = \sum_{I} {1 \over \sqrt{2 \pi}} {Z_{I}^2 \over R_{I}^{\rm c}} \enspace .$$ (92)

Here, the sums expand over all atoms in the simulation cell, all direct lattice vectors L, and the prime in the first sum indicates that ${I} < {J}$ is imposed for ${\bf L} = {\bf0}$ .