# Ewald summation

This notebook provides a detailed look at classical electrostatics, demonstrating problems related to convergence and how to overcome those issues.  The techniques discussed are ubiquitous in periodic simulations, both classical and quantum.


```python
%matplotlib notebook
from mpl_toolkits.mplot3d import Axes3D
import matplotlib
import numpy as np
import functools
import scipy.special
import itertools
import matplotlib.pyplot as plt
import matplotlib.animation as animation
from ipywidgets import interact, interactive, fixed, FloatSlider, widgets
```

## Classical electrostatics

Molecular simulations of periodic systems and crystals are an extremely important area of study.  To demonstrate the difficulties, and some popular methods to overcome them, we'll investigate a simple system of unit charges.  Central to our discussion is the potential $\phi$ at a point $\mathbf{R}$ generated by a Gaussian charge density centered at $\mathbf{R}_0$ with exponent $\alpha$

\begin{equation}
    q G(\alpha, \mathbf{R}_0, \mathbf{R}) \equiv q \left(\frac{\alpha}{\sqrt{\pi}}\right)^3
        e^{-\alpha^2(\mathbf{R}-\mathbf{R}_0)^2},
\end{equation}

which can be obtained as the solution of Poisson's equation:

\begin{equation}
    \nabla^2 \phi(\alpha, \mathbf{R}) = -\frac{q G(\alpha, \mathbf{0}, \mathbf{R})}{\epsilon_0}.
\end{equation}

The above involves the permittivity of free space $\epsilon_0$ and we have assumed the Gaussian is placed at the origin.  The Laplacian in spherical coordinates is

\begin{align}
    \nabla^2 & = \frac{1}{r^2}\frac{\partial}{\partial r}\left(r^2\frac{\partial}{\partial r}\right)
    +\frac{1}{r^2 \sin^2\phi} \frac{\partial^2}{\partial \theta^2}
    +\frac{1}{r^2 \sin \theta}\frac{\partial}{\partial \phi}\left(\sin \phi \frac{\partial}{\partial \phi}\right),
\end{align}

but only the leading term will survive when applied to a spherically symmetric density.  Recognizing that $\frac{1}{r^2}\frac{\partial}{\partial r}\left(r^2\frac{\partial}{\partial r}\right)f = \frac{1}{r}\frac{\partial^2}{\partial r^2}(rf)$ Poisson's equation for a Gaussian potential is

\begin{align}
    \frac{1}{r}\frac{\partial^2}{\partial r^2}\left[r \phi(\alpha, \mathbf{R})\right] 
        &= -\frac{q G(\alpha, \mathbf{0}, \mathbf{R})}{\epsilon_0} \\
    \frac{\partial^2}{\partial r^2}\left[r \phi(\alpha, \mathbf{R})\right] 
        &= -\frac{r q G(\alpha, \mathbf{0}, \mathbf{R})}{\epsilon_0} \\
    \frac{\partial}{\partial r}\left[r \phi(\alpha, \mathbf{R})\right] 
        &= -\int_r^\infty \frac{r q G(\alpha, \mathbf{0}, \mathbf{R})}{\epsilon_0}\mathrm{d}r 
            = \frac{q G(\alpha, \mathbf{0}, \mathbf{R})}{2\alpha^2\epsilon_0} \\
    r \phi(\alpha, \mathbf{R}) 
        &= \int_0^r \frac{q G(\alpha, \mathbf{0}, \mathbf{R})}{2\alpha^2\epsilon_0}\mathrm{d}r \\
    \phi(\alpha, \mathbf{R}) 
        &= \frac{q \mathrm{Erf}(\alpha r)}{4 \pi \epsilon_0 r},\\
\end{align}

which introduces the definition of the error function:
\begin{equation}
    \mathrm{Erf}(z) \equiv \frac{2}{\sqrt{\pi}}\int_0^z{-t^2}\mathrm{d}t
\end{equation}

and the scalar distance $r=|\mathbf{R}-\mathbf{R}_0|$.
Recognizing that the infinite exponent limit of a Gaussian is simply a delta function, we get a very familiar result

\begin{equation}
    \lim_{\alpha\rightarrow\infty}\frac{q \mathrm{Erf}(\alpha r)}{4 \pi \epsilon_0r} = \frac{q}{4 \pi \epsilon_0 r}
\end{equation}

*i.e.* the potential at a distance $r$ from a point charge $q$ is $\phi(r) = \frac{q}{4 \pi \epsilon_0 r}$.  Placing a test charge $Q$ that same distance away then gets us to Coulomb's law:

\begin{equation}
    U = Q\phi(r) = \frac{Qq}{4 \pi \epsilon_0 r}.
\end{equation}

If we want to embed a unit cell in an infinite number of replica ("image") unit cells in all directions, we just need to sum over all interactions within the primary unit cell, and over all primary-image cell interactions.  The unit cell is described by a set of vectors, which for the orthorhombic case with lengths $a,b,c$ are simply 

$$
    \mathbf{a} = \left(\begin{eqnarray} &a&0&0\\ &0&b&0\\ &0&0&c\\ \end{eqnarray}\right).
$$

With this notation, and labeling atoms $\{i,j\}$ possessing charges $\{q_i,q_j\}$, image cells can be compactly denoted by a vector $\mathbf{n} = \{n_1 \mathbf{a}_1, n_2 \mathbf{a}_2, n_3 \mathbf{a}_3\}$, where $n_1, n_2$ and $n_3$ are integers:

\begin{equation}
    U = \frac{1}{8 \pi \epsilon_0}\sum_{\mathbf{n}'}\sum_{i,j}\frac{q_i q_j}
            {\left|\mathbf{r}_i - \mathbf{r}_j + \mathbf{n}\right|} 
      = \frac{1}{8 \pi \epsilon_0}\sum_{\mathbf{n}'}\sum_{i,j}\frac{q_i q_j}{r_{ij}}
\end{equation}

The prime on the summation of over image cells is a reminder that $\mathbf{n}=\mathbf{0}$ terms where $i=j$ should be omited and we have introduced the shorthand notation $r_{ij}=\left|\mathbf{r}_i - \mathbf{r}_j + \mathbf{n}\right|$.  Implementing electrostatics by brute force summation over many lattice vectors is not a good idea; to see why we'll take a brief diversion.

## Conditional convergence
To get a feel for conditional convergence consider the classic example of the alternating harmonic series:

\begin{equation}
    \sum_{n=0}^\infty \frac{(-1)^{n+1}}{n} = \left(1 - \frac{1}{2} + \frac{1}{3}
        - \frac{1}{4} + \frac{1}{5} - \frac{1}{6} + \frac{1}{7} - \frac{1}{8} + \ldots\right) = \ln(2).
\end{equation}

Summing the terms in ascending order yields the result show.  However, we can also reorder the terms in parentheses before summing

\begin{equation}
    \underbrace{1 - \frac{1}{2}}_{\frac{1}{2}} -
    \underbrace{\frac{1}{4}}_{\frac{1}{4} } +
    \underbrace{\frac{1}{3} - \frac{1}{6} }_{\frac{1}{6}} -
    \underbrace{\frac{1}{8}}_{\frac{1}{8} } +
    \underbrace{\frac{1}{5} - \frac{1}{10} }_{\frac{1}{10}} -
    \underbrace{\frac{1}{12}}_{\frac{1}{12} } +
    \underbrace{\frac{1}{7} - \frac{1}{14} }_{\frac{1}{14}} -
    \underbrace{\frac{1}{16}}_{\frac{1}{16} } +
    \ldots,
\end{equation}

which yields exactly half of what we got above:

\begin{equation}
    \frac{1}{2}\left(1 - \frac{1}{2} + \frac{1}{3} - \frac{1}{4} + \frac{1}{5} - \frac{1}{6}
        + \frac{1}{7} - \frac{1}{8} + \ldots\right) = \frac{1}{2} \ln(2).
\end{equation}

There exists some specific permutation of the summands that will generate any given result when the sum is conditionally convergent.  A sum over some series $a_n$ is conditionally convergent if $\sum_{n=1}^\infty |a_n| = \infty$.  Returning to electrostatics, we introduce the integral analog of the potential expression above, avoiding the singularity by neglecting the infinitesimal point $\mathbf{r}=\mathbf{0}$:
$$\int_{r>0} \frac{1}{|\mathbf{r}|} \mathrm{d}^3\mathbf{r} = 4\pi \int_{\delta}^\infty \frac{r^2}{r} \mathrm{d}r = \infty$$
The factor of $r^2$ comes from the conversion to polar coordinates.  We find that

\begin{equation}
    \int_{\delta}^\infty \frac{1}{r^{n-2}} \mathrm{d}r = 
    \begin{cases}
    \frac{r_c^{3-n}}{n-3}, & \text{if}\ n>3 \\
    \infty, & \text{otherwise}
    \end{cases}
\end{equation}

so, by analogy to the integrals, the Coulomb summation is conditionally convergent in three dimensional space.  We don't arrive at an *absolutely* convergent series until we deal with an operator that decays as $\frac{1}{r^4}$ in three dimensions.  Not only is the Coulombic summation conditionally convergent, its convergence is very slow w.r.t. increasing cutoffs due to the slow decay of the operator.

The purpose of this notebook is to demonstrate a clever trick to decompose the slow, conditionally convergent summation into two rapid, absolutely convergent terms.  We use kJ/mol for energies, nm for length and atomic units for charge throughout.


```python
# Setup up a cubic simulation box
atoms_per_dim = 5
natoms = atoms_per_dim**3
epsilon0 = 5.72765E-4 # (e^2 mol) / (kJ nm)
coulomb = 1/(4*np.pi*epsilon0)
boxlength = 0.8 #nm
# Initial coordinates; regular lattice
spacing = np.linspace(-boxlength/2.0, boxlength/2.0, atoms_per_dim, endpoint=False)
spacing += boxlength/(2.0*atoms_per_dim)
# Center the array
coords = np.array(list(itertools.product(spacing, spacing, spacing)))
# Move the atoms very slightly away from their regular grid positions
np.random.seed(0)
coords += 0.2*(np.random.rand(natoms,3)-0.5)
coords[coords>boxlength/2.0] = boxlength/2.0
coords[coords<-boxlength/2.0] = -boxlength/2.0
charges = np.zeros(len(coords))
if 1:
    # Generate a large dipole
    charges[:natoms//2] = 1.0
    charges[(natoms+1)//2:] = -1.0
else:
    # Generate a more homogeneous distribution
    charges[1:natoms:2] = -1.0
    charges[0:natoms-1:2] = 1.0
poscrd = coords[charges > 0]
negcrd = coords[charges < 0]
neucrd = coords[charges == 0]
print("Sum of charges", np.sum(charges))
# Visualize the coordinates and charges in a rotatable 3D plot.  Atoms are colored according to their charge.
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.scatter(poscrd[:,0], poscrd[:,1], poscrd[:,2], 'bo')
ax.scatter(negcrd[:,0], negcrd[:,1], negcrd[:,2], 'ro')
ax.scatter(neucrd[:,0], neucrd[:,1], neucrd[:,2], 'ko')
ax.set_xlim(-boxlength/2.0, boxlength/2.0)
ax.set_ylim(-boxlength/2.0, boxlength/2.0)
ax.set_zlim(-boxlength/2.0, boxlength/2.0)
ax.set_xticks([])
ax.set_yticks([])
ax.set_zticks([])
plt.show()
```


```python
# Define the lattice energy using brute force summation, and evaluate over a range of lattice cutoffs.

def coulomb_energy(qij, R):
    """ Evaluate the conventional Coulomb energy, given products of charges and
        distances between the centers bearing those charges. """
    return np.sum(qij/R)


def lattice_sum(charges, coords, nboxes, efun):
    """ Compute the Coulomb energy of a system of charges, by considering interactions
        between the primary box and all image boxes that are within nboxes of the
        primary cell. The energy is evaluated by the provided function, efun. """
    natoms = coords.shape[0]
    boxlist = range(-nboxes,nboxes+1)

    # Build some lists to describe interactions 
    # products of charges between two boxes
    chargeij = np.prod(list(itertools.product(charges, charges)), axis=1)
    dRij = np.repeat(coords, natoms, axis=0) - np.tile(coords, (natoms,1))
    energy = 0.0
    truncated_energy = 0.0
    for nvecs in itertools.product(boxlist, boxlist, boxlist):
        nvecs = np.array(nvecs)
        shiftvec = nvecs*boxlength
        absvecs = np.abs(nvecs)
        if nvecs.any():
            # This is an image box
            shiftedRij = np.linalg.norm(dRij+shiftvec, axis=1)
            eterm = 0.5*efun(qij=chargeij, R=shiftedRij)
        else:
            # Handle the 'home' box
            eterm = 0.0
            for i in range(natoms):
                for j in range(i):
                    dR = np.linalg.norm(coords[i]-coords[j])
                    eterm += efun(qij=charges[i]*charges[j], R=dR)
        eterm *= coulomb
        if np.linalg.norm(absvecs) <= nboxes:
            truncated_energy += eterm
        energy += eterm        
    return energy,truncated_energy
    
vals = []
for n in range(7):
    vals.append(lattice_sum(charges, coords, n, efun=coulomb_energy))
vals = np.array(vals)

# Plot the convergence of the energy w.r.t. cutoff using the two truncation schemes
fig = plt.figure()
ax = fig.add_subplot(111)
xvals = range(1,len(vals))
ax.plot(xvals, vals[1:,0], 'bo', label ='Cubic')
ax.plot(xvals, vals[1:,1], 'ro', label = 'Spherical')
ax.set_xticks(xvals)
ax.set_xlabel('Number of image box layers')
ax.set_ylabel('Energy (kJ/mol)')
ax.legend(loc='lower right')
plt.show()

# Check the reference values
assert np.allclose(vals,
[[ 361515.2359571,   361515.2359571 ],
 [ 528282.46725449,  557057.25818972],
 [ 534335.79047581,  536496.90616012],
 [ 535962.70789396,  536005.76078745],
 [ 536633.65606731,  537475.71261986],
 [ 536973.60561882,  537518.58787992],
 [ 537169.21808044,  537527.68088663]])
```

The above code sets up a cell and performs the summation in two different ways 1) adding a complete layer of image boxes, which corresponds to a cubic macroscopic crystal and 2) adding image boxes in such a way that the macroscopic system is roughly spherical.  The convergence is both slow and conditional; for this contrived example of a large net dipole the energies for the two summation schemes converge to slightly different results.  When the second layer of vectors is added, we are introducing interaction up to three times the box length and the results are still not fully converged!

To speed up the convergence, Ewald summation introduces Gaussian functions at the atomic positions, with the opposite sign to a given atom's charge, to screen the point charges, making them short range.  To compensate, a Gaussian with the same sign as the charge is added in, but is handled separately.  From the above arguments, this affects the potential by splitting the Coulomb operator into short- and long-range components:

\begin{align}
    \frac{1}{r} &= \frac{1-\mathrm{Erf}(\alpha r)}{r} + \frac{\mathrm{Erf}(\alpha r)}{r}\\
                &= \frac{\mathrm{Erfc}(\alpha r)}{r} + \frac{\mathrm{Erf}(\alpha r)}{r}
\end{align}

In the plot below we show the density (left) where point charges and their screening Gaussians are shown in blue; the resulting $\mathrm{Erfc}$ potential is shown in the right panel.  Increasing the attenuation parameter $\alpha$, which is inversely proportional to the width of the Gaussians, screens this short range potential more effectively, but also increases the overall contribution of the long-range potential (red) that results from the compensating Gaussians.  The rapid decay of the short-range term means we can evaluate it in the standard real-space treatment using pairwise loops with short cutoffs. The red part is still long-range so will still decay slowly; however it is now free from the singularity that plagues the full (black) potential, and we can develop a treatment for it.


```python
# Define some point charges to demonstrate the density and potential
positions = [1, 3, 5]
qvals  = [1, -1, 1]

def plot_gaussians(gausplot, erfplot, alpha=3):
    """ Generate a plot that shows point charges in a 1D coordinates system, whose
        density is represented by blue delta functions, with blue screening Gaussians
        introduced to attenuate these point charges.  The Gaussian is also subtracted
        out, which is shown in red, to make the method exact. The potential corresponding
        to the regular Coulomb operator (black) is also shown, decomposed into its short
        range (blue) component and its long range (red) term. """
    prefac = alpha/np.sqrt(np.pi) # 1D Gaussian normalization
    def plot_gaussian(axis, xvals, origin, exponent, prefac, col):
        axis.plot(xvals, prefac*np.exp(-np.square(exponent*(xvals-origin))), color=col)
    npts = 150
    xvals = np.linspace(0.001, 6, npts)
    gausplot.cla()
    erfplot.cla()
    for p,c in zip(positions, qvals):
        plot_gaussian(gausplot, xvals, p, alpha,  c*prefac, 'r')
        plot_gaussian(gausplot, xvals, p, alpha, -c*prefac, 'b')
        gausplot.plot([p, p], [0, c], color='b')
    gausplot.plot([0, 40], [0, 0], color='k', linestyle='-', linewidth=2)
    erfvals = scipy.special.erf(alpha*xvals)
    inv = 1/xvals
    erfplot.set_title('Potential')
    erfplot.plot(xvals, erfvals*inv, color='r', label='Erf(alpha R) / R')
    erfplot.plot(xvals, (1-erfvals)*inv, color='b', label='Erfc(alpha R) / R')
    erfplot.plot(xvals, inv, '--', color='k', label='1 / R')
    erfplot.set_xlabel('r (nm)')
    erfplot.set_ylabel('$\\phi(r)$')
    erfplot.legend()
    gausplot.set_title('Density')
    gausplot.set_xlim(0, 6)
    gausplot.set_ylim(-5,5)
    gausplot.set_ylabel('$\\rho(r)$')
    gausplot.set_xlabel('Position')
    erfplot.set_xlim(0, 1.5)
    erfplot.set_ylim(0, 10)
    fig.canvas.draw()

# Draw the plots with a default alpha value.  Run the box below
# to get a slider to vary alpha and update the plots in real time.
fig, (gausplot, erfplot) = plt.subplots(1,2,figsize=(8,4))
plot_gaussians(gausplot, erfplot)
```


```python
# Evaluate this cell, then play with the sliders to see the effect of
# varying alpha on the long and short range densities and potentials.
interact(plot_gaussians, gausplot=fixed(gausplot), erfplot=fixed(erfplot),
         alpha=widgets.FloatSlider(min=0,max=9,step=.2,value=3),);
```

Looking at the form of the long-range density in the left panel of the above plot, it is a periodic, smoothly oscillating function so any solution is likely to involve trigonometric functions somehow.  Returning to the discussion at the top of the notebook, we remember that the potential comes from solving the Poisson equation:

\begin{equation}
    \nabla^2 \phi(\alpha, \mathbf{R}) = -\frac{q G(\alpha, \mathbf{0}, \mathbf{R})}{\epsilon_0}.
\end{equation}

Combining these considerations, it's easy to imagine that the trigonometric functions will be introduced via complex exponentials; because $\frac{\partial}{\partial x}e^{-imx} = -im e^{-imx}$, introducing an exponential will effectively turn the derivative operator in the Poisson equation into a constant, making it easy to handle.  As a reminder, the connection between complex exponentials and trigonometric functions is provided by Euler's formula: $e^{ix} = \mathrm{cos}(x)+i\mathrm{sin}(x)$.  The transformation that will introduce the complex exponential is the Fourier transform.

First, we normalize the coordinates into the range $[0,1]$, appropriate for integration, by introducing the reciprocal lattice vectors $\mathbf{a}^*$, which are defined in terms of the unit cell vectors by the orthonormality relation $\mathbf{a}_i^*\cdot\mathbf{a}_j=\delta_{ij}$ and for the orthorhombic unit cell are

$$ \mathbf{a}^* = \left(\begin{eqnarray} &\frac1a&0&0\\ &0&\frac1b&0\\ &0&0&\frac1c\\ \end{eqnarray}\right). $$

The resulting reciprocal space coordinate system has orthogonal axes, even for non-orthorhombic crystals.  Analogous to $\mathbf{n}$ indexing summation over image unit cells, we use $\mathbf{m}$ to index summation over *reciprocal* space lattice vectors, $\mathbf{m} = \{m_1 \mathbf{a}_1^*, m_2 \mathbf{a}_2^*, m_3 \mathbf{a}_3^*\}$; because we're evaluating the energy using these vectors, we call the long range energy the *reciprocal* energy, $U_\mathbf{rec}$.

The long-range density $\rho^L(\mathbf{r})$ and potential $\phi^L(\mathbf{r})$ have Fourier transforms defined by

\begin{align}
    \tilde\rho^L(\mathbf{m}) &= \int_V \rho^L(\mathbf{r}) e^{-2\pi i\mathbf{m}\cdot\mathbf{r}} \mathrm{d}^3\mathbf{r}\\
    \tilde\phi^L(\mathbf{m}) &= \int_V \phi^L(\mathbf{r}) e^{-2\pi i\mathbf{m}\cdot\mathbf{r}} \mathrm{d}^3\mathbf{r}
\end{align}

and the inverse transforms are defined by summation over reciprocal lattice vectors $\mathbf{m}$

\begin{align}
    \rho^L(\mathbf{r}) &= \frac1V \sum_{\mathbf{m}}\tilde\rho^L(\mathbf{m}) e^{2\pi i\mathbf{r}\cdot\mathbf{m}}\\
    \phi^L(\mathbf{r}) &= \frac1V \sum_{\mathbf{m}}\tilde\phi^L(\mathbf{m}) e^{2\pi i\mathbf{r}\cdot\mathbf{m}}
\end{align}

Inserting $\tilde\rho^L(\mathbf{m})$ into Poisson's equation yields the Fourier space expression

\begin{equation}
    4\pi^2m^2 \tilde\phi^L(\mathbf{m}) = \frac{\tilde\rho^L(\mathbf{m})}{\epsilon_0}
\end{equation}

whose solution is trivial if we neglect the $\mathbf{m}=\mathbf{0}$ term, which is zero for neutral unit cells with $\sum_i q_i = 0$. We define the scalar $m=|\mathbf{m}|$.  The reciprocal space long-range density is obtained from the Fourier transform:

\begin{align}
    \tilde\rho(\mathbf{m}) &= \int_V \sum_\mathbf{n}\sum_j q_j 
        G(\alpha, \mathbf{r}_j, \mathbf{r+\mathbf{n}}) 
        e^{-2\pi i \mathbf{r}\cdot\mathbf{m}}\mathrm{d}\mathbf{r}^3 \\
     &=\sum_\mathbf{n}\sum_j q_j  \int_\mathrm{all\ space} 
         G(\alpha, \mathbf{r}_j, \mathbf{\mathbf{w}}) 
         e^{-2\pi i \mathbf{(w-n)}\cdot\mathbf{m}}\mathrm{d}\mathbf{w}^3 \\
     &=\sum_\mathbf{n}\sum_j q_j  \int_\mathrm{all\ space}
         G(\alpha, \mathbf{r}_j, \mathbf{\mathbf{w}})
         e^{-2\pi i \mathbf{w}\cdot\mathbf{m}}\mathrm{d}\mathbf{w}^3 \\
     &= \sum_j q_j e^{-\frac{m^2\pi^2}{\alpha^2}}
         e^{-2\pi i \mathbf{r}_j\cdot\mathbf{m}}.
\end{align}

In developing this expression we first change variables $\mathbf{w}=\mathbf{r+n}$ before taking advantage of the fact that $e^{2\pi i \mathbf{n}\cdot\mathbf{m}} = 1$. The integral on the third line is the definition of the Fourier transform of a 3D Gaussian, which is evaluated for each dimension using the standard result $\int e^{a^2x^2}e^{-2\pi i m x}\mathrm{d}x = \frac{\sqrt{\pi}}{a}e^{-\frac{\pi^2m^2}{a^2}}$.  With the Fourier space expression for the density in hand, we can now get the Fourier space potential from the simplified Poisson equation:

\begin{equation}
    \tilde\phi(\mathbf{m}) = \frac{1}{4\pi^2\epsilon_0}
    \frac{e^{-\frac{m^2\pi^2}{\alpha^2}}}{m^2} \sum_j q_j  e^{-2\pi i \mathbf{r}_j\cdot\mathbf{m}}
\end{equation}

This is inverse Fourier transformed to obtain the real space potential at some point $\mathbf{r}$:

\begin{equation}
    \phi(\mathbf{r}) = \frac{1}{4\pi^2\epsilon_0 V}
        \sum_{\mathbf{m}\ne\mathbf{0}}\frac{e^{-\frac{m^2\pi^2}{\alpha^2}}}{m^2}
        \sum_j q_j  e^{-2\pi i \mathbf{r}_j\cdot\mathbf{m}}
        e^{2\pi i \mathbf{r}\cdot\mathbf{m}}
\end{equation}

And finally, as above, we evaluate the potential at each atom, multiply by the charge of that atom, and divide by 2 to avoid double-counting:

\begin{align}
  U_\mathrm{rec} &= \frac{1}{8\pi^2\epsilon_0 V}
      \sum_{\mathbf{m}\ne\mathbf{0}}\frac{e^{-\frac{m^2\pi^2}{\alpha^2}}}{m^2}
      \sum_{i} q_i e^{2\pi i \mathbf{r}_i\cdot\mathbf{m}}
      \sum_{j} q_j  e^{-2\pi i \mathbf{r}_j\cdot\mathbf{m}} \\
   &= \frac{1}{8\pi^2\epsilon_0 V}\sum_{\mathbf{m}\ne\mathbf{0}}
       \frac{e^{-\frac{m^2\pi^2}{\alpha^2}}}{m^2}
       S(\mathbf{m})
       S(-\mathbf{m})
\end{align}

On the surface, this doesn't seem like a simplification.  However, the presence of the factor $\frac{e^{-\frac{m^2\pi^2}{\alpha^2}}}{m^2}$ leads to a rapidly convergent summation involving a limited number of reciprocal lattice vectors $\mathbf{m}$.  The quantity $S(\mathbf{m})$ is known as a structure factor.  The only remaining issue to deal with is the fact that the real space term neglects the $i=j$ term in the primary unit cell summation, but the reciprocal space energy contains this term as a result of the reciprocal space potential containing contributions from all atoms.  The spurious term that we need to remove is the interaction of the probe charge with its own contribution to the reciprocal space potential, which comes from its screening Gaussian:

\begin{split}
    U_\mathrm{self} &= -\sum_i \lim_{r\rightarrow0} \frac{q_i \mathrm{Erf}(\alpha r)q_i}{4 \pi \epsilon_0 r} \\
   &= -\frac{1}{4 \pi \epsilon_0}\frac{\alpha}{\sqrt\pi}\sum_i q_i^2
\end{split}

In the code below, we investigate the effect of introducing the screening to the short-range *direct space* (as opposed to reciprocal space, used to evaluate the long-range term) energy:

\begin{equation}
    U = \frac{1}{2}\sum_{\mathbf{n}'}\sum_{i,j}
        \frac{q_i q_j\mathrm{Erfc}(\alpha\left|\mathbf{r}_i - \mathbf{r}_j + \mathbf{n}\right|)}
                                {\left|\mathbf{r}_i - \mathbf{r}_j + \mathbf{n}\right|}.
\end{equation}

After choosing an attenuation parameter $\alpha=4 \mathrm{nm}^{-1}$ we repeat the lattice summation shown above:


```python
def ewald_dir_energy(alpha, maxn):
    """ An energy function, to be passed to the lattice sum code, that
        evaluated the direct contribution to the Ewald energy. """
    # Define the direct space energy as the attenuated direct Ewald energy.
    def ewald_dir_interaction(qij, R, alpha):
        """ Evaluate the energy of pairs of charges, separated by R,
            using the attenuated Coulomb operator. """
        return np.sum(qij*scipy.special.erfc(alpha*R)/R)
    enerfunc = functools.partial(ewald_dir_interaction, alpha=alpha)
    return lattice_sum(charges, coords, n, efun=enerfunc)
    
vals = np.array([ewald_dir_energy(alpha=4, maxn=n) for n in range(4)])

# Plot the convergence of the lattice sum, with respect to number of
# image cells included, using cubic and spherical layering schemes.
fig = plt.figure()
ax = fig.add_subplot(111)
xvals = range(1,len(vals))
ax.plot(xvals, vals[1:,0], 'bo', label ='Cubic')
ax.plot(xvals, vals[1:,1], 'ro', label = 'Spherical')
ax.set_xticks(xvals)
ax.set_xlabel('Number of image box layers')
ax.set_ylabel('Energy (kJ/mol)')
ax.legend()

plt.show()
```

Now both summation schemes converge rapidly to the same results, showing that the direct space part of the lattice sum is well-behaved.  Implementing the reciprocal and self terms, we find that the overall Ewald energy is also well-behaved:


```python
def ewald_rec_energy(alpha, maxm):
    """ Evaluate the long range (reciprocal space) Ewald energy using an
        attenuation parameter of alpha, and using up to maxm wavevectors. """
    mdim = 2*maxm+1
    twopii = np.pi*2j
    pi2_a2 = np.pi**2/alpha**2 if alpha != 0.0 else 1e50
    miter = [ 0 ]
    for m in range(1,maxm+1):
        miter.extend([m,-m])
    miter = np.array(miter)
    mvals = miter / boxlength
    xvecs = np.exp(twopii*np.outer(coords[:,0], mvals))
    yvecs = np.exp(twopii*np.outer(coords[:,1], mvals))
    zvecs = np.exp(twopii*np.outer(coords[:,2], mvals))
    Sm = np.einsum('n,nx,ny,nz->xyz', charges, xvecs, yvecs, zvecs)
    m = np.array(list(itertools.product(mvals, mvals, mvals)))
    m2 = np.einsum('nx,nx->n', m, m).reshape(mdim,mdim,mdim)
    prefac = np.exp(-pi2_a2*m2)
    m2[0,0,0] = 1e50 # avoid division by 0
    prefac /= m2
    V = boxlength**3
    return coulomb*np.einsum('xyz,xyz', prefac, np.square(Sm.real) + np.square(Sm.imag)) / (2*np.pi*V)

def ewald_self_energy(alpha):
    """ Evaluate the Ewald self interaction energy. """
    return -coulomb*(alpha/np.sqrt(np.pi))*np.sum(np.square(charges))
    
def ewald_energy(alpha, maxn, maxm):
    """ Evaluate all components of the Ewald sum, to get the total energy. """
    Udir = ewald_dir_energy(alpha, maxn)[0]
    Urec = ewald_rec_energy(alpha, maxm)
    Uslf = ewald_self_energy(alpha)
    U = Udir+Urec+Uslf
    return Udir,Urec,Uslf,U

# A list of alpha, maxn, maxm to loop over
params = np.array([ 
    [0, 4, 2],
    [1, 3, 2],
    [2, 2, 2],
    [3, 2, 2],
    [4, 1, 3],
    [5, 1, 3],
    [6, 1, 4],
])
params = np.array([ 
    [0, 4, 2],
    [1, 3, 3],
    [2, 2, 4],
    [3, 2, 4],
    [4, 1, 5],
    [5, 1, 5],
    [6, 1, 6],
])

# Generate Ewald energies, varying the attentuation paramter alpha,
# as well as the real and recpiprocal space cutoffs needed to reach
# convergence for each alpha value.
Uvals = np.array([ewald_energy(*p) for p in params])

# Plot the total energy, and its components, as a function of the
# attenuation parameter alpha.
fig = plt.figure()
ax = fig.add_subplot(111)
alphavals = params[:,0]
ax.plot(alphavals, Uvals[:,0], '--r', label = 'Udir')
ax.plot(alphavals, Uvals[:,1], '--b', label = 'Urec')
ax.plot(alphavals, Uvals[:,2], '--g', label = 'Uslf')
ax.plot(alphavals, Uvals[:,3], 'k', label = 'U')
ax.set_xticks(alphavals)
ax.set_xlabel('Alpha (1/nm)')
ax.set_ylabel('Energy (kJ/mol)')
ax.legend()

plt.show()

# Confirm that all high alpha values have the same total energy
assert np.allclose(Uvals[2,3], Uvals[3:,3])
```

For almost any nonzero screening parameter, the total energy is the same, although its composition in terms of reciprocal, direct and self terms changes; this useful fact is very handy when testing the correctness of an Ewald implementation.  For small $\alpha$ more $\mathbf{n}$ vectors are needed, as the short-range operator is not well attenuated, while a large $\alpha$ demands more $\mathbf{m}$ vectors to evaluate the reciprocal space term; this is evident from the values chosen in the `params` array above.  Because cutoffs can be used for the direct space term, it scales as $\mathcal{O}(N)$.  Being attenuated, the direct space term can be evaluated within the primary unit cell, using the minimum image convention.

Constructing the stucture factor $S(\mathbf{m})$ involves a summation over $N$ atoms as well as $\mathbf{m}$; for a given $\alpha$ value the latter scales with the volume of the box and is therefore proportional to $N$.  Overall, the Ewald reciprocal energy therefore scales as $\mathcal{O}(N^2)$, although it can be shown that optimally tuning $\alpha$ to minimize computational cost yields a $\mathcal{O}(N^{{3/2}})$ algorithm.  A clever trick introduced by Darden, York and Pedersen in 1993 reduces this cost further, as we will now demonstrate.

The key insight is that a discrete Fourier transform can be performed using a fast recursion algorithm if the atoms are arranged on a uniform mesh, hence the name particle mesh Ewald (PME).  A method that only works for selected configurations would be useless, so PME uses spline interpolation to effectively coerce the atoms to predefined, regularly spaced locations.  Consider some function $f$ we've tabulated at every integer.  If we want either $f(12)$ or $f(13)$, we're in luck and can just look up the tabulated values.  However, $f(12.5)$ would come from $\frac{f(12)+f(13)}{2}$, while $f(12.1)$ would be best approximated with two points by $0.9f(12) + 0.1f(13)$. This simple order 2 example can be generalized to order $n$, with the coefficients of each pre-tabulated data point provided by an order $n$ *spline* function.  The cardinal basis spline (B-spline) used in PME is defined for interpolation using the nearerst $n$ points by the recursion

\begin{equation}
    M_n(u) = \frac{u}{n-1}M_{n-1}(u) + \frac{n-u}{n-1} M_{n-1}(u-1)
\end{equation}

where $u$ is the fractional distance from the nearest grid point to the left of the point of interest.  Using this method to discretize the density to a regular lattice the fast Fourier transform (FFT) can then be used to perform the discrete Fourier transform step, and scales as $\mathcal{O}(N log(N))$ instead of $\mathcal{O}(N^2)$; the same is true of the inverse tranformation step.  The PME algorithm is otherwise similar to conventional Ewald method, as demonstrated below:


```python
def PME_rec_energy(alpha, maxm, splineorder):
    """ Evaluate the particle mesh Ewald long range (reciprocal space)
        energy via discretization onto a regular cubic grid. """
    
    def build_spline(order, u):
        """ Construct a B spline of given order, from a seed value u """
        def one_pass_bspline(array, val, n):
            div = 1.0 / float(n - 1)
            array[n-1] = div * val * array[n-2]
            for j in range(1, n-1): 
                array[n-j-1] = div * ((val+j) * array[n-j-2] + (n-j-val) * array[n-j-1])
            array[0] *= div * (1 - val)
            return array
        
        array = np.zeros(order)
        array[0] = 1.0 - u
        array[1] = u
        for m_n in range(1, order-1):
            array = one_pass_bspline(array, u, m_n+2)
        return array
    
    def dftmod(splineorder, nfft):
        """ Compute the norm of a B spline in Fourier space. """
        ftarr = np.zeros(nfft)
        # modulus is invariant to position, so just make a spline with 0
        ftarr[:splineorder] = build_spline(splineorder, 0.0)
        ftarr = np.fft.fft(ftarr)
        return np.square(ftarr.real) + np.square(ftarr.imag)

    mdim = 2*maxm+1
    # A convenient array to allow iteration over the grid from any given
    # starting point,  correctly wrapping around the boundaries.
    griditer = np.array([ [ ((p+n)%mdim) for p in range(splineorder) ] for n in range(mdim) ])
    # Define a simple iterator array to bring m into the range -0.5 < m <= 0.5
    shiftm = np.array([ k - mdim if k >= (mdim+1)/2 else k for k in range(mdim) ])

    splines = np.zeros((natoms, 3, splineorder))
    gridstart = np.zeros((natoms, 3), dtype=int)
    # Define coordinates in fractional coordinates
    reduced_coords = coords/boxlength + 0.5
    for atom,coord in enumerate(reduced_coords):
        # Define gridpoint to the left in each dimension, small shift
        # to handle atoms located on the rightmost grid point
        gridstart[atom,:] = list(map(int, mdim*coord-1e-8))
        wvals = mdim*coord - gridstart[atom]
        for xyz in range(3):
            splines[atom,xyz,:] = build_spline(splineorder, wvals[xyz])
            
    # Form the influence function, including spline normalization terms.  Note
    # that this is independent of particle positions and can therefore be formed
    # at the start of the simulation and cached for reuse in each time step.
    pi2_a2 = np.pi**2/alpha**2 if alpha != 0.0 else 1e50
    mvals = shiftm / float(boxlength)
    splinemod = 1/dftmod(splineorder, mdim)
    splinemods = np.einsum('x,y,z->xyz', splinemod, splinemod, splinemod)
    m = np.array(list(itertools.product(mvals, mvals, mvals)))
    m2 = np.einsum('nx,nx->n', m, m).reshape(mdim,mdim,mdim)
    gF = np.multiply(splinemods, np.exp(-pi2_a2*m2))
    m2[0,0,0] = 1e50 # avoid division by 0
    gF /= m2
        
    # Spread the charges onto the charge grid, using spline coefficients at O(N) cost
    Qgrid = np.zeros((mdim, mdim, mdim))
    for atom in range(natoms):
        q = charges[atom]
        for nx,xgrid in enumerate(griditer[gridstart[atom,0]]):
            xspline = splines[atom,0,nx]
            for ny,ygrid in enumerate(griditer[gridstart[atom,1]]):
                yspline = splines[atom,1,ny]
                for nz,zgrid in enumerate(griditer[gridstart[atom,2]]):
                    zspline = splines[atom,2,nz]
                    Qgrid[xgrid,ygrid,zgrid] += q*xspline*yspline*zspline
                    
    # Forward FFT Q to form S(m), at O(Nlog(N)) cost
    Sm = np.fft.fftn(Qgrid)
    
    # Convolution of gF with S(m), at O(N) cost
    conv = np.multiply(gF, Sm)
    
    # Inverse FFT (S(m)*gF) real potential grid, at O(Nlog(N)) cost
    phigrid = mdim**3*np.fft.ifftn(conv).real

    # Probe the potential grid using the splines, to get the energy at O(N) cost
    e = 0.0
    for atom in range(natoms):
        q = charges[atom]
        for nx,xgrid in enumerate(griditer[gridstart[atom,0]]):
            xspline = splines[atom,0,nx]
            for ny,ygrid in enumerate(griditer[gridstart[atom,1]]):
                yspline = splines[atom,1,ny]
                for nz,zgrid in enumerate(griditer[gridstart[atom,2]]):
                    zspline = splines[atom,2,nz]
                    e += q*xspline*yspline*zspline*phigrid[xgrid,ygrid,zgrid]
                
    return e*coulomb / (2*np.pi*boxlength**3)

    
def PME_energy(alpha, maxn, maxm, splineorder):
    """ Compute the total PME energy by evaluating all components. """
    Udir = ewald_dir_energy(alpha, maxn)[0]
    Urec = PME_rec_energy(alpha, maxm, splineorder)
    Uslf = ewald_self_energy(alpha)
    U = Udir+Urec+Uslf
    return Udir,Urec,Uslf,U

# Repeat the same scan performed above for Ewald, using PME with 4th order spline.
efunc = functools.partial(PME_energy, splineorder=5)
Uvals = np.array([efunc(*p) for p in params])

# Plot the PME energy and its components as a function of attenuation parameter.
fig = plt.figure()
ax = fig.add_subplot(111)
alphavals = params[:,0]
ax.plot(alphavals, Uvals[:,0], '--r', label = 'Udir')
ax.plot(alphavals, Uvals[:,1], '--b', label = 'Urec')
ax.plot(alphavals, Uvals[:,2], '--g', label = 'Uslf')
ax.plot(alphavals, Uvals[:,3], 'k', label = 'U')
ax.set_xticks(alphavals)
ax.set_xlabel('Alpha (1/nm)')
ax.set_ylabel('Energy (kJ/mol)')
ax.legend()

plt.show()

# Confirm that all high alpha values have the same total energy
assert np.allclose(Uvals[2,3], Uvals[3:,3], atol=5)
```

The PME energy performs similarly to regular Ewald summation, with significantly better asymptotic scaling.  Forces may be evaluated easily by replacing the splines used to probe the potential grid with their derivatives, which are defined by a very simple recursion relation, as outlined by Darden and co-workers.

In developing the absolutely convergent sums above, we seem to have magically removed the conditional convergence.  The reason for this is somewhat complicated and beyond the scope of this discussion.  The conditional convergence results from the mutual interaction of the image cells' dipole moments.  By thinking of the periodic system not as an infinite system, but a large macroscopic crystal with a given shape (we considered cubes and spheres in the lattice sum code above) is is possible to derive a "correction" to the Ewald energy that describes the conditionally (and therefore crystal shape) dependent part of the lattice sum.  Such corrections have been derived for various crystal shapes, and they all depend on the cell dipole.  In dynamic simulations under periodic boundary condisions, a charge particle exiting the top of a unit cell will reappear at the bottom, possibly flipping the dipole of the cell.  If the extra shape-dependent term were considered in this circumstance, the potential energy would be discontinuous, so these corrections are usually neglected.



References
======

1. M. P. Allen and D. J. Tildesley, Computer Simulation of Liquids, Oxford University Press, New York 1987.
2. D. Frenkel and B. Smit, Understanding Molecular Simulation: From Algorithms to Applications, Academic Press, San Diego 1996.
3. U. Essmann, L. Perera, M. L. Berkowitz, T. A. Darden, H. Lee, and L. G. Pedersen, J. Chem. Phys. 103, 8577 (1995).
4. C. Sagui, L. G. Pedersen, and T. A. Darden, J. Chem. Phys. 120, 73 (2004).

5dab578f1b22455e2eb7dac2a068d2bb072f36ff
ipynb
Jupyter Notebook
Example/Psi4Numpy/12-MD/12b_ewald.ipynb
yychuang/109-2-compchem-lite
cbf17e542f9447e89fb48de1b28759419ffff956
BSD-3-Clause
2017-03-01T08:04:48.000Z
2022-03-23T08:52:04.000Z
Example/Psi4Numpy/12-MD/12b_ewald.ipynb
yychuang/109-2-compchem-lite
cbf17e542f9447e89fb48de1b28759419ffff956
BSD-3-Clause
2017-03-03T13:20:20.000Z
2022-03-05T18:20:27.000Z
Example/Psi4Numpy/12-MD/12b_ewald.ipynb
yychuang/109-2-compchem-lite
cbf17e542f9447e89fb48de1b28759419ffff956
BSD-3-Clause
2017-02-17T19:44:47.000Z
2022-03-22T05:52:43.000Z
Qwen/Qwen-72B
1. YES
2. YES