Useful links:

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 Short-Ranged Potential
Mapping Particles to and from the Mesh
Reciprocal Space Convolution of the Long-Ranged Potential

Introduction: What PME is, and Four Myths that It is Not

The (Smooth) Particle-Mesh Ewald (PME) method found in Amber is a special case of the Particle-Particle / Particle-Mesh (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 long-ranged, slowly-decaying potentials in periodic systems. Moreover, the P3M methods make it cheaper to solve O(N²) problems.

The central concept is to split the potential function into short- and long-ranged components. The short-ranged, rapidly decaying part of the potential is readily solved as a cutoff-based Particle-Particle problem while the long-ranged part of the potential can be solved by representing each particle, which can still have any real-valued position in three-dimensional space, as a weighted stencil on a regular mesh grid. The potential function for particle-particle 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 system--default 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ų of space. But, wait--with ten times the density, doesn't that make an O(N²) 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 element-wise at each mesh point, and then back-transforming the product. It is then, finally, the Fast Fourier Transform that turns the O(N²) problem into one which is merely O(N log N). It was therefore the Cooley-Tukey FFT algorithm, a collaboration between computer scientists at IBM and Princeton that re-discovered mathematics originally credited to Carl Friedrich Gauss, which propelled the development of particle-mesh methods in the following years.

There are some myths about Particle-Mesh Ewald electrostatics that should be busted. First is the notion that the interactions of particles are not pairwise-decomposable 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 three-dimensional all-to-all problem to a series of one-dimensional all-to-all 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 all-to-all 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 Particle-Mesh 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 this--the short-ranged 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 short-ranged potential is handled by the complement, 1 - erf(αr), which for the Coulomb potential translates to kq_iq_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 short-ranged 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.0e-5 (the interaction of two point charges i and j according to kq_iq_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_iq_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 short-ranged 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 128-bit 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 true--the same quantity will be calculated every time. A longer cutoff with the same direct sum tolerance will still neglect about the same amount of short-range 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 long-ranged potential can be calculated. Similarly, a denser grid does not change the overall result‐again, it increases the precision of the long-ranged part of the potential. To increase the precision of the short-ranged 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 long-ranged 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 short-ranged part of the potential is the most important. Some approaches neglect the long-ranged part entirely. Another code has even referred to the long-ranged 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 short-ranged 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 short-ranged cutoff, will be equal to the direct sum tolerance.) With a denser mesh, the short-ranged 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 short-ranged sum. Similarly, if the charges in the simulation are spherical Gaussians to begin with, the short-ranged part of the potential can be attenuated or removed entirely.

Hopefully this has removed some misconceptions about Ewald sums in long-ranged 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