Long Range Electrostatics#
Ewald summation#
The total electrostatic energy of
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 (301) into two quickly-converging terms and a constant term:
where
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
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 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
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
Pressure scaling works with PME, but be aware of the fact that anisotropic scaling can introduce artificial ordering in some systems.
P3M-AD#
The Particle-Particle Particle-Mesh methods of Hockney & Eastwood can also be applied in GROMACS for the treatment of long range electrostatic interactions 106. Although the P3M method was the first efficient long-range electrostatics method for molecular simulation, the smooth PME (SPME) method has largely replaced P3M as the method of choice in atomistic simulations. One performance disadvantage of the original P3M method was that it required 3 3D-FFT back transforms to obtain the forces on the particles. But this is not required for P3M and the forces can be derived through analytical differentiation of the potential, as done in PME. The resulting method is termed P3M-AD. The only remaining difference between P3M-AD and PME is the optimization of the lattice Green influence function for error minimization that P3M uses. However, in 2012 it has been shown that the SPME influence function can be modified to obtain P3M 107. This means that the advantage of error minimization in P3M-AD can be used at the same computational cost and with the same code as PME, just by adding a few lines to modify the influence function. However, at optimal parameter setting the effect of error minimization in P3M-AD is less than 10%. P3M-AD does show large accuracy gains with interlaced (also known as staggered) grids, but that is not supported in GROMACS (yet).
P3M is used in GROMACS with exactly the same options as used with PME by selecting the electrostatics type:
coulombtype = P3M-AD
Optimizing Fourier transforms and PME calculations#
It is recommended to optimize the parameters for calculation of electrostatic interaction such as PME grid dimensions and cut-off radii. This is particularly relevant to do before launching long production runs.
gmx mdrun will automatically do a lot of PME optimization, and GROMACS also includes a special tool, gmx tune_pme, which automates the process of selecting the optimal number of PME-only ranks.