Contents   Index

Ewald Sum

The Ewald sum [12] is the best technique for calculating electrostatic interactions in a periodic (or pseudo-periodic) system.

The basic model for a neutral periodic system is a system of charged point ions mutually interacting via the Coulomb potential. The Ewald method makes two amendments to this simple model. Firstly each ion is effectively neutralised (at long range) by the superposition of a spherical gaussian cloud of opposite charge centred on the ion. The combined assembly of point ions and gaussian charges becomes the Real Space part of the Ewald sum, which is now short ranged and treatable by the methods described above (section 2.1). ^2.3 The second modification is to superimpose a second set of gaussian charges, this time with the same charges as the original point ions and again centred on the point ions (so nullifying the effect of the first set of gaussians). The potential due to these gaussians is obtained from Poisson's equation and is solved as a Fourier series in Reciprocal Space. The complete Ewald sum requires an additional correction, known as the self energy correction, which arises from a gaussian acting on its own site, and is constant. Ewald's method therefore replaces a potentially infinite sum in real space by two finite sums: one in real space and one in reciprocal space; and the self energy correction.

For molecular systems, as opposed to systems comprised simply of point ions, additional modifications are necessary to correct for the excluded (intra-molecular) Coulombic interactions. In the real space sum these are simply omitted. In reciprocal space however, the effects of individual gaussian charges cannot easily be extracted, and the correction is made in real space. It amounts to removing terms corresponding to the potential energy of an ion $\ell$ due to the gaussian charge on a neighbouring charge $m$ (or vice versa). This correction appears as the final term in the full Ewald formula below. The distinction between the error function $erf$ and the more usual complementary error function $erfc$ found in the real space sum, should be noted.

The total electrostatic energy is given by the following formula.

$$U_{c} = \frac{1}{2V_{o}\epsilon_{0}} \sum_{\mbox{$\underline{k}$}\neq\mbo... ...\pi\epsilon_{0}}\sum_{n<j}^{N^{*}}\frac{q_{j}q_{n}}{r_{nj}} erfc(\alpha r_{nj})$$

$$\phantom{=xxxxx}-\frac{1}{4\pi\epsilon_{0}} \sum_{molecules}\sum_... ...rd \pi}+\frac{erf(\alpha r_{\ell m})}{r_{\ell m}^{1-\delta_{\ell m}}}\right \},$$ (2.131)

where $N$ is the number of ions in the system and $N^{*}$ the same number discounting any excluded (intramolecular) interactions. $M^{*}$ represents the number of excluded atoms in a given molecule and includes the atomic self correction. $V_{o}$ is the simulation cell volume and $\mbox{$\underline{k}$}$ is a reciprocal lattice vector defined by

$$\mbox{$\underline{k}$}=\ell \mbox{$\underline{u}$} + m \mbox{$\underline{v}$} + n \mbox{$\underline{w}$}$$ (2.132)

where $\ell,m,n$ are integers and $\mbox{$\underline{u}$},\mbox{$\underline{v}$},\mbox{$\underline{w}$}$ are the reciprocal space basis vectors. Both $V_{o}$ and $\mbox{$\underline{u}$},\mbox{$\underline{v}$},\mbox{$\underline{w}$}$ are derived from the vectors ( $\mbox{$\underline{a}$},\mbox{$\underline{b}$},\mbox{$\underline{c}$}$) defining the simulation cell. Thus

$$V_{o}=\vert\mbox{$\underline{a}$}\cdot\mbox{$\underline{b}$}\times\mbox{$\underline{c}$}\vert$$ (2.133)

and

$$\mbox{$\underline{u}$} = 2\pi\frac{\mbox{$\underline{b}$}\times \mbox{$\underline{c}$}}{\mbox{$\underline{a}$}\cdot\mbox{$\underline{b}$}\times\mbox{$\underline{c}$}}$$

$$\mbox{$\underline{v}$} = 2\pi\frac{\mbox{$\underline{c}$}\times \mbox{$\underline{a}$}}{\mbox{$\underline{a}$}\cdot\mbox{$\underline{b}$}\times\mbox{$\underline{c}$}}$$ (2.134)

$$\mbox{$\underline{w}$} = 2\pi\frac{\mbox{$\underline{a}$}\times \mbox{$\underline{b}$}}{\mbox{$\underline{a}$}\cdot\mbox{$\underline{b}$}\times\mbox{$\underline{c}$}}.$$

With these definitions, the Ewald formula above is applicable to general periodic systems. (A small additional modification is necessary for rhombic dodecahedral and truncated octahedral simulation cells [30].)

In practice the convergence of the Ewald sum is controlled by three variables: the real space cutoff $r_{cut}$; the convergence parameter $\alpha$ and the largest reciprocal space vector $\mbox{$\underline{k}$}_{max}$ used in the reciprocal space sum. These are discussed more fully in section 3.3.5. DL_POLY_2 can provide estimates if requested (see CONTROL file description 4.1.1.

The force on an atom $j$ is obtained by differentiation and is

$$\mbox{$\underline{f}$}_{j} = -\frac{q_{j}}{V_{o}\epsilon_{0}} \sum_{\mbox{$\underline{k}$}\neq... ... \sum_{n}^{N}q_{n}\exp(-i\mbox{$\underline{k}$}\cdot\mbox{$\underline{r}$}_{n})$$

$$\phantom{xxxxx} +\frac{q_{j}}{4\pi\epsilon_{0}}\sum_{n}^{N^{*}}\f... ...{nj}}{\surd\pi}\exp(-\alpha^{2}r_{nj}^{2})\right \} \mbox{$\underline{r}$}_{nj}$$ (2.135)

$$\phantom{xxxxxxxxxx} -\frac{q_{j}}{4\pi\epsilon_{0}}\sum_{\ell}^{... ...surd\pi}\exp(-\alpha^{2}r_{\ell j}^{2})\right \}\mbox{$\underline{r}$}_{\ell j}$$

The electrostatic contribution to the system virial can be obtained as the negative of the Coulombic energy. However in DL_POLY_2 this formal equality can be used as a check on the convergence of the Ewald sum. The actual electrostatic virial is obtained during the calculation of the diagonal of the stress tensor.

The electrostatic contribution to the stress tensor is given by

$$\mbox{$\underline{\underline{\bf\sigma}}$} = \frac{1}{2V_{o}\epsilon_{0}} \sum_{\mbox{$\underline{k}$}\neq\mbo... ...^{N}q_{j}\exp(-i\mbox{$\underline{k}$}\cdot\mbox{$\underline{r}$}_{j})\vert^{2}$$

$$\phantom{xxxxx}+\frac{1}{4\pi\epsilon_{0}}\sum_{j<n}^{N^{*}}\frac... ...exp(-\alpha^{2}r_{nj}^{2})\right \} \mbox{$\underline{\underline{\bf R_{nj}}}$}$$ (2.136)

$$\phantom{xxxxxxxxxx} -\frac{1}{4\pi\epsilon_{0}}\sum_{j<\ell}^{M^... ...pha^{2}r_{\ell j}^{2})\right \}\mbox{$\underline{\underline{\bf R_{\ell j}}}$},$$

where matrices $\mbox{$\underline{\underline{\bf K}}$}$ and $\mbox{$\underline{\underline{\bf R_{\ell j}}}$}$ are defined as follows.

$$K^{\alpha\beta} = k^{\alpha}k^{\beta}$$ (2.137)

$$R_{\ell j}^{\alpha\beta} = r_{\ell j}^{\alpha}r_{\ell j}^{\beta}$$ (2.138)

In DL_POLY_2 the full Ewald sum is handled by several routines: EWALD1 and EWALD1A handle the reciprocal space terms; EWALD2, EWALD2_3PT, EWALD2_RSQ and EWALD4, EWALD4_3PT handle the real space terms (with the same Verlet neighbour list routines that are used to calculate the short range forces); and EWALD3 calculates the self interaction corrections. It should be noted that the Ewald potential and force interpolation arrays in DL_POLY_2 are erc and fer respectively.

  Contents   Index