# Long Range Electrostatics¶

## Ewald summation¶

The total electrostatic energy of $$N$$ particles and their periodic images is given by

(300)$V=\frac{f}{2}\sum_{n_x}\sum_{n_y} \sum_{n_{z}*} \sum_{i}^{N} \sum_{j}^{N} \frac{q_i q_j}{{\bf r}_{ij,{\bf n}}}.$

$$(n_x,n_y,n_z)={\bf n}$$ is the box index vector, and the star indicates that terms with $$i=j$$ should be omitted when $$(n_x,n_y,n_z)=(0,0,0)$$. The distance $${\bf r}_{ij,{\bf n}}$$ is the real distance between the charges and not the minimum-image. This sum is conditionally convergent, but very slow.

Ewald summation was first introduced as a method to calculate long-range interactions of the periodic images in crystals 105. The idea is to convert the single slowly-converging sum (300) into two quickly-converging terms and a constant term:

(301)\begin{split}\begin{aligned} V &=& V_{\mathrm{dir}} + V_{\mathrm{rec}} + V_{0} \\[0.5ex] V_{\mathrm{dir}} &=& \frac{f}{2} \sum_{i,j}^{N} \sum_{n_x}\sum_{n_y} \sum_{n_{z}*} q_i q_j \frac{\mbox{erfc}(\beta {r}_{ij,{\bf n}} )}{{r}_{ij,{\bf n}}} \\[0.5ex] V_{\mathrm{rec}} &=& \frac{f}{2 \pi V} \sum_{i,j}^{N} q_i q_j \sum_{m_x}\sum_{m_y} \sum_{m_{z}*} \frac{\exp{\left( -(\pi {\bf m}/\beta)^2 + 2 \pi i {\bf m} \cdot ({\bf r}_i - {\bf r}_j)\right)}}{{\bf m}^2} \\[0.5ex] V_{0} &=& -\frac{f \beta}{\sqrt{\pi}}\sum_{i}^{N} q_i^2,\end{aligned}\end{split}

where $$\beta$$ is a parameter that determines the relative weight of the direct and reciprocal sums and $${\bf m}=(m_x,m_y,m_z)$$. In this way we can use a short cut-off (of the order of $$1$$  nm) in the direct space sum and a short cut-off in the reciprocal space sum (e.g. 10 wave vectors in each direction). Unfortunately, the computational cost of the reciprocal part of the sum increases as $$N^2$$ (or $$N^{3/2}$$ with a slightly better algorithm) and it is therefore not realistic for use in large systems.

### Using Ewald¶

Don’t use Ewald unless you are absolutely sure this is what you want - for almost all cases the PME method below will perform much better. If you still want to employ classical Ewald summation enter this in your mdp file, if the side of your box is about $$3$$  nm:

coulombtype     = Ewald
rvdw            = 0.9
rlist           = 0.9
rcoulomb        = 0.9
fourierspacing  = 0.6
ewald-rtol      = 1e-5


The ratio of the box dimensions and the fourierspacing parameter determines the highest magnitude of wave vectors $$m_x,m_y,m_z$$ to use in each direction. With a 3-nm cubic box this example would use $$11$$ wave vectors (from $$-5$$ to $$5$$) in each direction. The ewald-rtol parameter is the relative strength of the electrostatic interaction at the cut-off. Decreasing this gives you a more accurate direct sum, but a less accurate reciprocal sum.

## PME¶

Particle-mesh Ewald is a method proposed by Tom Darden 14 to improve the performance of the reciprocal sum. Instead of directly summing wave vectors, the charges are assigned to a grid using interpolation. The implementation in GROMACS uses cardinal B-spline interpolation 15, which is referred to as smooth PME (SPME). The grid is then Fourier transformed with a 3D FFT algorithm and the reciprocal energy term obtained by a single sum over the grid in k-space.

The potential at the grid points is calculated by inverse transformation, and by using the interpolation factors we get the forces on each atom.

The PME algorithm scales as $$N \log(N)$$, and is substantially faster than ordinary Ewald summation on medium to large systems. On very small systems it might still be better to use Ewald to avoid the overhead in setting up grids and transforms. For the parallelization of PME see the section on MPMD PME (Multiple-Program, Multiple-Data PME parallelization).

With the Verlet cut-off scheme, the PME direct space potential is shifted by a constant such that the potential is zero at the cut-off. This shift is small and since the net system charge is close to zero, the total shift is very small, unlike in the case of the Lennard-Jones potential where all shifts add up. We apply the shift anyhow, such that the potential is the exact integral of the force.

### Using PME¶

As an example for using Particle-mesh Ewald summation in GROMACS, specify the following lines in your mdp file:

coulombtype     = PME
rvdw            = 0.9
rlist           = 0.9
rcoulomb        = 0.9
fourierspacing  = 0.12
pme-order       = 4
ewald-rtol      = 1e-5


In this case the fourierspacing parameter determines the maximum spacing for the FFT grid (i.e. minimum number of grid points), and pme-order controls the interpolation order. Using fourth-order (cubic) interpolation and this spacing should give electrostatic energies accurate to about $$5\cdot10^{-3}$$. Since the Lennard-Jones energies are not this accurate it might even be possible to increase this spacing slightly.

Pressure scaling works with PME, but be aware of the fact that anisotropic scaling can introduce artificial ordering in some systems.

coulombtype     = P3M-AD