by Dave Cerutti
Download the library of
Matlab or
GNU Octave functions for
this demonstration. Aside from the Introduction, this demonstration is
divided into multiple sections:
Unwrapping the library, and Coordinate
Transformations
Computing the ShortRanged Potential
Mapping Particles to and from the Mesh
Reciprocal Space Convolution of the LongRanged
Potential
Introduction: What PME is, and Four Myths that It is Not
The (Smooth) ParticleMesh Ewald (PME) method found in Amber is a special
case of the ParticleParticle / ParticleMesh (P3M) methods that have been used
in physics simulations since the mid 1970s. A good reference on the general
principles can be found in Computer Simulation Using Particles by R.W. Hockney
and J.W. Eastwood. The Amber implementation is based on
Ulrich
Essmann's work. The problem that these techniques solve is the
indeterminate nature of longranged, slowlydecaying potentials in periodic
systems. Moreover, the P3M methods make it cheaper to solve O(N^{2})
problems.
The central concept is to split the potential function into short
and longranged components. The shortranged, rapidly decaying part of the
potential is readily solved as a cutoffbased ParticleParticle problem while
the longranged part of the potential can be solved by representing each
particle, which can still have any realvalued position in threedimensional
space, as a weighted stencil on a regular mesh grid. The potential function
for particleparticle interactions is then expressed as an interaction between
mesh points: each mesh point becomes its own particle, in a sense,
participating in a convolution with all other mesh point "particles." In Amber,
there are generally about ten times as many mesh points as atoms in the
systemdefault settings set a mesh spacing of as close to 1Å as
possible without going over, and the typical density of water places three
atoms in about 29Å^{3} of space. But, waitwith ten times the
density, doesn't that make an O(N^{2}) problem a hundred times
worse? Yes, but the regularity of the grid enables us to make use of the
convolution theorem, which states that the convolution of every mesh point with
every other can be performed by taking the Fourier transform of the original
density mesh, plus the Fourier transform of the potential function, multiplying
the two results elementwise at each mesh point, and then backtransforming the
product. It is then, finally, the Fast Fourier Transform that turns the
O(N^{2}) problem into one which is merely O(N log N). It was therefore
the CooleyTukey
FFT algorithm, a collaboration between computer scientists at IBM and
Princeton that rediscovered mathematics originally credited to Carl Friedrich
Gauss, which propelled the development of particlemesh methods in the
following years.
There are some myths about ParticleMesh Ewald electrostatics that should be
busted. First is the notion that the interactions of particles are not
pairwisedecomposable in such a scheme. They are, it's just rather expensive
to do so! Using the FFT to simultaneously convolve each mesh point with every
other‐more or less reducing the threedimensional alltoall problem to a
series of onedimensional alltoall problems which, in turn, reduce to
(N log N) thanks to Cooley and Tukey‐may lead one to believe that the
way to extract the influence of individual particles on the rest of the system
is to delete them, repeat the entire calculation, and take the difference. But
a smarter way to go about it is to isolate individual particles, separately
compute the potential grids projected by each of them throughout the simulation
box, and then take the (charge) density mapping of any other particle in the
system in the context of one of those fields to get a potential energy of any
desired pair interaction. Even more efficient would be to compute the
influence (electrostatic potential) due to a unit charge density at exactly one
mesh point (not a particle sitting directly on a mesh point and projecting
density over a stencil of 64 or more points, just a unit of charge at only one
mesh point) over the rest of the box. That influence map constitutes a lookup
table by which the density maps of any two system particles interact‐if
each particle maps to a stencil of 64 points, as occurs in the default Amber
settings, then the alltoall interaction of these mapped points, the sum of
4096 table lookups from the influence map of the unit charge density at a
single mesh point, is the interaction of the particles. Voila, pairwise
decomposable interactions in ParticleMesh Ewald electrostatics (really painful
to calculate, but still better than repeating the entire FFT if all that's
needed is a few of the pairs). All of this is feasible with the
PME library and some custom scripting.
A second myth about PME is that the primary image is treated specially.
Technically, there is some truth to thisthe shortranged part of the
potential may yet extend indefinitely even though it is negligible beyond the
cutoff. This is the case in PME, where the splitting function is the error
function erf(αr) (&alpha, the Ewald coefficient, controls the strength of
the splitting as will be seen later in this lecture). The shortranged
potential is handled by the complement, 1  erf(αr), which for the
Coulomb potential translates to
kq_{i}q_{j}(1 
erf(αr))/r_{ij} for particles
i and j (k is Coulomb's constant). This quantity is
computed up to a cutoff, at which the another setting known as the "direct sum
tolerance," our term for the point at which neglecting the shortranged part of
the potential can be no worse than a given amount, actually determines the
Ewald coefficient α. For the default settings, a cutoff of 8Å and
a direct sum tolerance of 1.0e5 (the interaction of two point charges
i and j according to
kq_{i}q_{j}(1 
erf(αr_{ij}))/r_{ij}
beyond 8Å must be not more than
0.001&0025; of their interaction under the original potential
kq_{i}q_{j}/r_{ij}. This
determines an α of about 0.35, which can also be thought of as a Gaussian
charge spreading of about 1.4Å. In this sense, the primary image is
treated differently, but even if the shortranged sum were carried over to
other images, the contributions from the nearest images in most simulations
would be beneath machine precision, even if computed in an insanely long
128bit floating point format.
A third myth about Ewald electrostatics, perhaps more precisely a rookie
mistake, is that using longer cutoffs on the direct space part will change the
results. Because of the splitting and the direct sum tolerance, this is not
truethe same quantity will be calculated every time. A longer cutoff with
the same direct sum tolerance will still neglect about the same amount of
shortrange potential contributions (although the sources of those neglected
contributions may shift around). What longer cutoffs will change is the
Ewald coefficient‐which is more aptly thought of as the Gaussian charge
spreading factor‐increasing the spread more or less proportionally to the
cutoff for a given direct sum tolerance and increasing the precision with which
the longranged potential can be calculated. Similarly, a denser grid does not
change the overall result‐again, it increases the precision of the
longranged part of the potential. To increase the precision of the
shortranged part of the potential, the direct sum tolerance must be tightened,
but this necessarily implies a smaller Gaussian charge spreading factor, which
reduces precision in the longranged potential sum. In summary, the cutoff,
direct sum tolerance, and charge spreading factor are all linked, whereas the
mesh grid density can be adjusted to compensate. The cutoff and mesh grid
density are the only two parameters that affect the cost of the calculation.
Most commonly used simulation settings are chosen to balance the precision in
both parts while minimizing the overall effort.
The final myth to dispel about Ewald electrostatics is that the shortranged
part of the potential is the most important. Some approaches
neglect the longranged
part entirely. Another code has even referred to the longranged potential
gradients as the "Dark Side of the Force." However, Coulomb's law is really
implemented on the mesh. The only caveat in Ewald electrostatics is that the
mesh facilitates computation of the electrostatic potential due to a system of
spherical Gaussian charges in exactly the same positions as the point
charges that the original system contains. The shortranged part of the
potential, which is often substantially more expensive to compute, is a
correction for the difference in interactions between point and spherical
Gaussian charges. (The relative difference between the interaction of two
spherical Gaussian charges and two point charges, either pair separated by the
shortranged cutoff, will be equal to the direct sum tolerance.) With a denser
mesh, the shortranged potential cutoff can be decreased without sacrificing
accuracy, and with a very dense mesh the cutoff can be driven so low that no
two atoms in the simulation will ever come close enough to contribute to the
shortranged sum. Similarly, if the charges in the simulation are spherical
Gaussians to begin with, the shortranged part of the potential can be
attenuated or removed entirely.
Hopefully this has removed some misconceptions about Ewald sums in
longranged electrostatics potentials. They are a special case of a much
broader class of methods that make use of the Fast Fourier Transform to compute
infinite potentials using mesh grids. In the following pages we will explore
the basic operations of PME using the Matlab library. Readers that do not have
Matlab can download GNU Octave for free
here.
Go to the next section
