7 Long-Range Interactions: The Ewald Summation

So far, we have considered interparticle interactions that are short-ranged by construction. Because the Lennard-Jones potential decays so strongly with distance (as \(r^{-6}\)), it is acceptable to cut off this interaction at moderate distances and, if desired, add a correction factor which is the result of integrating the potential over a uniform particle density out to \(r = \infty \). However, Coulomb interactions, common in molecular simulation, decay relatively much more slowly (as \(r^{-1}\)) and as a consequence, we cannot compute a correction factor; the integral diverges. There are several ways to handle long-ranged interactions, but the most popular is the Ewald summation [14], which we discuss here. This discussion is drawn primarily from F&S chapter 12 [1], and the excellent paper by Markus Deserno and Christian Holm [15, 16].

7.1 The Ewald Coulombic energy

First, assume we have a collection of charged particles in a cubic box with side length \(L\), with periodic boundary conditions. The collection is assumed neutral; there is an equal number of positive and negative charges. The total Coulombic energy in this system is given by:

\(\phi (r_i)\) is the electrostatic potential at position \(r_i\): \begin{equation} \phi \left (r_i\right ) = {\sum _{j=1}^{N}}^\prime {\sum _{{\bf n}\in \mathbb {Z}^3}} \frac {q_j}{\left |{\bf r}_{ij} + {\bf n}L\right |} \end{equation} \(\bf n\) is a three dimensional integer vector. The prime on the first summation indicates that we do not admit the term for which \(j=i\) if \({\bf n} = (0,0,0)\). That is, we allow each particle to interact with its periodic images, but not with itself.

To evaluate \(\mathscr {U}\) efficiently, we break it into two parts:

• A short-ranged potential treated with a simple cutoff;
• A long-ranged potential which is periodic and slowly varying, which can therefore be represented to an acceptable level of accuracy by a finite Fourier series.

How can we do this? The idea of Ewald is to do two things: first, screen each point charge using a diffuse cloud of opposite charge around each point charge, and then compensate for these screening charges using a smoothly varying, periodic charge density. The screening charge is constructed to make the electrostatic potential due to a charge at position \(r_j\) decay rapidly to near zero at a prescribed distance. These interactions are treated in real space. The compensating charge density, which is the sum of all screening densities except with opposite charges, is treated using a Fourier series.

The standard choice for a screening potential is Gaussian: \begin{equation} \rho _{\rm s}(r) = -q_i(\alpha /\pi )^{3/2} e^{-\alpha r^2} \end{equation}

So for each charge, we add such a screening potential. Now, to evaluate \(\mathscr {U}_{\rm Coul}\), we have to evaluate the potential of a charge density that compensates for the screening charge densities at each particle. This is done in Fourier space.

The potential of a given charge distribution is given by Poisson’s equation: \begin{equation} \label {eq:ft_pois} -\nabla ^2\phi (r) = 4\pi \rho (r) \end{equation}

Now, the compensating charge distribution, denoted \(\rho _1\), can be written: \begin{equation} \rho _1(r) = \sum _{j=1}^{N}\sum _{{\bf n}\in \mathbb {Z}^3} q_j (\alpha /\pi )^{3/2} \exp \left [-\alpha \left |{\bf r} - \left ({\bf r}_j + {\bf n}L\right )\right |\right ] \end{equation} Notice that the sum over \(j\) includes the self-interaction when we include the potential due to this charge density in the calcuation of the total Coulombic energy (i.e., we have omitted the prime on the outer summation over particle indices).

Now, consider the Fourier transform of Poisson’s equation: \begin{equation} k^2 \tilde {\phi }(k) = 4\pi \tilde {\rho }(k) \end{equation} The Fourier transform of \(\rho _1(r)\) is given by

(The math required to show this involves noting that the the Fourier transform of a Gaussian is another Gaussian, and that the integral over all space of a normalized Gaussian is unity.) The \(\bf k\)-vectors are given by \begin{equation} {\bf k} = \frac {2\pi }{L} {\bf l} \ \ \ \ {\bf l} \in \mathbb {Z}^3 \end{equation}

We can use Eq. 165 to solve for \(\tilde {\phi }(k)\): \begin{equation} \label {eq:ft_phi1} \tilde {\phi }(k) = \frac {4\pi }{k^2} \sum _{j=1}^{N}q_j e^{-i{\bf k}\cdot {\bf r}_j}e^{-k^2/4\alpha } \end{equation} Note that this solution is not defined for \(k = 0\). In fact, we have to assume that \(\tilde \phi (0) = 0\), which is consistent with the notion that our system and all its periodic images is embedded in a medium of infinite dielectric constant (a perfect conductor; the “tinfoil” boundary condition).

Fourier inverting \(\tilde \phi (k)\) gives \begin{equation} \phi _i(r) = \frac {1}{V}\sum _{{\bf k}\ne 0}\tilde \phi (k)e^{i{\bf k}\cdot {\bf r}} \end{equation} which, when we substitute for \(\tilde {\phi }(k)\) from Eq. 169 yields \begin{equation} \phi _i(r) = \sum _{{\bf k}\ne 0}\sum _{j=1}^{N} \frac {4\pi q_j}{Vk^2} e^{i{\bf k}\cdot \left ({\bf r}-{\bf r}_j\right )}e^{-k^2/4\alpha } \end{equation}

So, the total Coulombic energy due to the compensating charge distribution is

where \begin{equation} \rho (k) = \sum _{i=1}^{N}q_ie^{i{\bf k}\cdot {\bf r}} \end{equation}

Notice that this does indeed include a spurious self-self interaction, because the point charge at \({\bf r}_i\) interacts with the compensating charge cloud also at \({\bf r}_i\). This self-interaction is the potential at the center of a Gaussian charge distribution. First, we solve Poisson’s equation for the potential due to a Gaussian charge distribution (details in F&S): \begin{equation} -\frac {1}{r}\frac {\partial ^2r\phi _{\rm Gauss}}{\partial r^2} = 4\pi \rho _{\rm Gauss}(r) \end{equation} yielding \begin{equation} \phi _{\rm Gauss}(r) = \frac {q_i}{r}{\rm erf}(\sqrt {a}r) \end{equation} where \(\rm erf\) is the error function: \begin{equation} {\rm erf}(x) = \frac {2}{\sqrt {\pi }}\int _0^x e^{-r^2}dr \end{equation} At \(r=0\), we have \begin{equation} \phi _{\rm self} = \phi _{\rm Gauss}(0) = 2\left (\frac {\alpha }{\pi }\right )^{1/2}q_i \end{equation}

So the total self-interaction energy becomes \begin{equation} \mathscr {U}_{\rm self} = \left (\frac {\alpha }{\pi }\right )^{1/2}\sum _{i=1}^{N}q_i^2 \end{equation} which must be subtracted from the total Coulombic energy.

Finally, the real-space contribution of the point charge at \({\bf r}_i\) is the screened potential: \begin{equation} \phi _{\rm short}(r) = \frac {q_i}{r} - \frac {q_i}{r}{\rm erf}\left (\sqrt {a}r\right ) \equiv \frac {q_i}{r}{\rm erfc}\left (\sqrt {\alpha }r\right ) \end{equation} where \(\rm erfc\) is the complementary error function. The total real-space Coulombic potential energy is therefore \begin{equation} \mathscr {U}_{\rm short} = \frac {1}{2}\sum _{i\ne j}^{N} \frac {q_i q_j}{r_{ij}} {\rm erfc}\left (\sqrt {\alpha }r_{ij}\right ) \end{equation}

Putting it all together:

where \begin{equation} \rho (k) = \sum _{i=1}^{N}q_ie^{i{\bf k}\cdot {\bf r}}. \end{equation}

Now, the arbitrariness left to us at this point is in a choice for the parameter \(\alpha \). Clearly, very small alphas make the Gaussians tighter and therefore the compensating charge distribution less smoothly varying. This means a Fourier series representation of \(\mathscr {U}_1\) with a given number of terms is more accurate for larger \(\alpha \). We’ll evaluate choice of \(\alpha \) in Sec. 7.3.

7.2 Ewald Forces

Now, we can use Eq. ?? in a Monte Carlo simulation of a system of charges, provided that periodic boundary conditions are used and the domain is cubic. (Extensions to non-cubic boxes and slab geometries are discussed to a limited extent in F&S.) We can also use the Ewald technique to calculate forces for use in molecular dynamics simulations.

The force on particle \(i\) due to the charges in the system is given by \begin{equation} {\bf F}_i = -\frac {\partial }{\partial {\bf r}_i}\mathscr {U}_{\rm Coul} \end{equation}

For our purposes, the two contributions to \({\bf F}_i\) are due to the \(k\)-space energy and the short-ranged, real-space energy: \begin{equation} {\bf F}_i = {\bf F}_i^{(k)} + {\bf F}_i^{(r)} \end{equation} Notice that there is no change in \(\mathscr {U}_{\rm self}\) when \({\bf r}_i\) changes, so no forces arise from \(\mathscr {U}_{\rm self}\).

The \(k\)-space contribution is given by \begin{equation} {\bf F}_i^{(k)} = q_i\sum _j q_j \frac {1}{V} \sum _{{\bf k}\ne {\bf 0}} \frac {4\pi {\bf k}}{k^2}e^{-k^2/4\alpha }\sin \left ({\bf k}\cdot {\bf r}_{ij}\right ) \end{equation}

The real-space contribution is given by \begin{equation} {\bf F}_i^{(r)} = q_i\sum _j q_j \left [2\sqrt {\frac {\alpha }{\pi }}e^{-\alpha r_{ij}^2}+\frac {1}{r_{ij}}{\rm erfc}(\sqrt {\alpha } r_{ij})\right ]\frac {{\bf r}_{ij}}{r_{ij}^2} \end{equation}

7.3 Implementation and Evaluation

We will consider an Ewald implementation which is a modified version of the ewald code written for Berend Smit’s Molecular Simulation course. (All of Prof. Smit’s codes are available in the FrenkelSmitCodes directory of the instructional-codes respository.) This code simply computes the Ewald energy for a cubic lattice, given an appropriate number of particles, and a value for \(\sqrt {\alpha }\) (which is called \(\alpha \) in the code), and a value for \(k_{\rm max}\), the maximum integer index for enumerating \(k\)-vectors. ¹

The units used in a system with electrostatics differ depending on community. So far, we have assumed that the units of electrostatic potential are charge \(\mathscr {C}\), divided by length \(\mathscr {L}\), because we write potential as \(\phi = q/|{\bf r}|\), where \(q\) is measured in units of \(\mathscr {C}\) and distance in units of \(\mathscr {L}\). Energy is therefore written in units of \(\mathscr {C}^2\) over \(\mathscr {L}\), and force in units of \(\mathscr {C}^2\) over \(\mathscr {L}^2\). If we want the final energy in more familiar units, we can choose \(\mathscr {C}\) and \(\mathscr {L}\), and use the standard prefactor \(1/4\pi \epsilon _0\) to convert from “charge squared per length” to “energy”. For example, in SI units, \(\epsilon _0 = 8.85419\times 10^{-12}\) (C\(^2\)/m)/J. In this implementation, we use a length of \(\mathscr {L} = 1\) and \(\mathscr {C} = 1\) and measure energy such that \(1/4\pi \epsilon _0 = 1\).

We will examine two configurations, both with \(N\) = 8\(^3\) = 512 particles, with alternating + and - charges. One configuration has the particle on a cubic lattice with lattice spacing \(r_0 = 1\), which is the standard NaCl crystal structure. We will call this the “crystal” configuration. The other is like the crystal, only each particle is displaced by a random amount from its Self Part lattice position with a maximum displacement of 0.3. We will call this the “liquid” configuration. We compute the total electrostatic energy via the Ewald sum technique for various values of \(1/\sqrt {\alpha }\) and maximum \(k\)-vector index. As we increase the number of \(k\)-vectors taken in the sum, we would like to show that the total energy converges to a certain value. We will measure this in terms of the Madelung constant, \(M\): \begin{equation} \mathscr {U}_{\rm Coul} = -\frac {N q^2 M}{4\pi \epsilon _0r_0} \end{equation}

Table 2 shows results of Ewald summation for the perfect lattice, and Table 3 shows resuts for the “liquid”. We see several interesting things from these example calculations:

1.

The self-energy is the dominant contributor from the long-range compensating charge smears.

2.

The Fourier-space contribution is relatively small, and is much smaller for the perfect lattice than for the liquid. This means the precision of the Fourier transform of the charge distribution (which is larger for larger \(k_{\rm max}\) is relatively unimportant to the calculation.

3.

The overall Coulomb energy is insensitive to the real-space cutoff for the perfect lattice, while it is weakly decreasing in magnitude (it is negative overall) for increasing real-space cutoff in the liquid.