Amber masthead
Filler image AmberTools23 Amber22 Manuals Tutorials Force Fields Contacts History
Filler image

Useful links:

Amber Home
Download Amber
Amber Citations
GPU Support
Mailing Lists
For Educators
File Formats
Materials for Educators → Particle Mesh Ewald Electrostatics

Computing the Short-Ranged Potential

The short-ranged potential, also known as the "direct space sum," is conceptually the simplest part of PME. Simply loop over all particles i and j, compute the minimum image distance between them, and if the distance rij is within the cutoff, apply kqiqj(1 - erf(αrij))/rij to get the energy contribution (differentiate to get the force--this part of the potential is often splined for efficiency).

This is an ideal time to look at the Ewald coefficient, which to reiterate is the inverse of the Gaussian charge spreading factor. The function getalpha takes two arguments, the short-ranged potential cutoff and the direct sum tolerance:

getalpha(8.0, 1.0e-5)

ans = 0.3486

Running that in a loop to generate multiple Gaussian charge spreading factors shows a nearly linear relationship with the cutoff:

i = 0;
gspread = zeros(1,61);
for Lcut = 6.0:0.1:12.0, i = i+1; gspread(i) = 0.5/getalpha(Lcut, 1.0e-5); end
plot(6.0:0.1:12.0, gspread, 'r.', 'markersize', 24);
hold on;
i = 0;
for Lcut = 6.0:0.1:12.0, i = i+1; gspread(i) = 0.5/getalpha(Lcut, 5.0e-6); end
plot(6.0:0.1:12.0, gspread, 'k.', 'markersize', 24);
i = 0;
for Lcut = 6.0:0.1:12.0, i = i+1; gspread(i) = 0.5/getalpha(Lcut, 1.0e-6); end
plot(6.0:0.1:12.0, gspread, '.', 'color', [1.0 0.4 0.0], 'markersize', 24);
axis([6 12 1.0 2.4]);
daspect([4 1 1]);
Gaussian Sigma Vs. Lcut plot

This is a useful result--for all practical values of the short-ranged cutoff and direct sum tolerance, the Gaussian charge spreading factor is nearly proportional to the cutoff. The ratio of the charge spreading factor to mesh grid spacing is essential to precision in the reciprocal space part of the calculation. For the default Amber interpolation method (4th order, this will be explained in detail later), the charge spreading factor needs to be about 1.5 times the mesh spacing to meet the short-ranged potential precision imparted by a direct sum tolerance of 1.0e-5 (note how this matches a spreading factor of 1.434Å to mesh spacings of close to but not exceeding 1.0Å). For a more expensive interpolation scheme (6th order), the spreading factor can be as small as the 1.0 times the mesh spacing and still provide a precise long-ranged potential calculation. That relationship breaks down a bit for really large mesh spacings (⋝ 1.8Å), but it is seldom worthwhile to extend the cutoff that far just to get a coarser mesh.

While the short-ranged, "direct space" potential is conceptually simple, the O(N2) problem arises again, and the code can be much more efficient by testing only interactions that are likely to yield success on the cutoff criterion. This is the origin of neighbor lists, hash tables, neutral territory methods, and multiple time steps to avoid recalculating the tail of the electrostatic short-ranged potential, which is small and changes slowly. The complexities introduced by all of these factors create the preponderance of code devoted to direct space electrostatics in most modern MD packages. The PME library contains pmeDirectBasicLoop.m for the basic all-to-all particle-particle computation, pmeDirectMidpointOrthog.m for the midpoint method, and pmeDirectTowerPlate.m for the tower-plate method. While these codes illustrate various distinct methods of performing the direct space decomposition, none of them is highly optimized.

To run the direct space calculation as a standalone operation, we first need a system: a four-column matrix containing coordinates and charges, [ X Y Z Q ], plus a box size. Take the cubeions system in the Tests/ folder. The dimensions array has only three elements, but the box angles are then assumed to all be 90°. Any of the three functions can then be used to calculate the forces on the particles:

load Tests/cubeions.mat
Lcut = 8.0;
dtol = 1.0e-5;
ewcoeff = getalpha(Lcut, dtol);
fbasic = pmeDirectBasicLoop(crdq, 0, gdim, ewcoeff, Lcut);
ftower = pmeDirectTowerPlate(crdq, 0, gdim, ewcoeff, Lcut);
fmidpt = pmeDirectMidpointOrthog(crdq, 0, gdim, ewcoeff, Lcut);

Even in Matlab, with all its clunky loops, it is evident that the methods based on hash cells are faster than the all-to-all implementation, which is so slow for systems of any significant size that it periodically prints progress to the terminal to let users know that Matlab didn't hang. The differences in forces will be below 1.0e-4 kcal/mol-Å, but not trivial due to splined versus analytic computation of the short-ranged potential and also non-associativity in floating-point addition.

A different exercise can illustrate how the Ewald coefficient / Gaussian spreading factor controls the way in which Coulomb's law is split into short- and long-ranged parts. Try re-calculating the forces using a different cutoff, but the same Ewald coefficient.

Lcut = 8.0;
ewcoeff  = getalpha(Lcut, 1.0e-5);
fmidpt8  = pmeDirectMidpointOrthog(crdq, 0, gdim, ewcoeff, Lcut);
fmidpt9  = pmeDirectMidpointOrthog(crdq, 0, gdim, ewcoeff, Lcut + 1.0);
fmidpt10 = pmeDirectMidpointOrthog(crdq, 0, gdim, ewcoeff, Lcut + 2.0);
fmidpt12 = pmeDirectMidpointOrthog(crdq, 0, gdim, ewcoeff, Lcut + 4.0);

It is easy to check that forces in fmidpt8 and fmidpt9 may differ by up to a hundredth of a kcal/mol-Å. Not a big deal, given the magnitudes of the actual forces. As would be expected, fmidpt9 and fmidpt10 differ by much less‐about an order of magnitude‐and fmidpt10 and fmidpt12 hardly differ at all. The short-ranged potential and its gradients converge. But, this is not an exercise that can be done with the sander or pmemd programs--they calculate the Ewald coefficient given a short-ranged cutoff and a direct sum tolerance, whereas here we have kept the Ewald coefficient fixed and in effect driven down the direct sum tolerance with increasingly long cutoffs, computing the sum beyond the original limit of 8Å. In the mdgx program, a niche application, it is possible to set the Gaussian charge spreading factor explicitly, in which case the direct sum tolerance will be determined based on the short-ranged cutoff. We can emulate the behavior of the workhorse Amber MD engines as follows:

Lcut = 8.0;
ewcoeff  = getalpha(Lcut, 1.0e-5);
fmidpt8  = pmeDirectMidpointOrthog(crdq, 0, gdim, ewcoeff, Lcut);
Lcut = 9.0;
ewcoeff  = getalpha(Lcut, 1.0e-5);
fmidpt9  = pmeDirectMidpointOrthog(crdq, 0, gdim, ewcoeff, Lcut);

A simple comparison of the forces reveals that the two are markedly different‐this is not because either is wrong, but rather the splitting has changed, with the Gaussian charge spreading factor being 1.434Å in the first case and 1.625Å in the second. We can use getalpha with guess-and-check to determine what the direct sum tolerance would need to be in order to match the original splitting:

0.5/getalpha(9.0, 5.0e-6) = 1.5597
0.5/getalpha(9.0, 2.0e-6) = 1.4840
0.5/getalpha(9.0, 1.0e-6) = 1.4334
0.5/getalpha(9.0, 1.1000e-6) - 0.5/getalpha(8.0, 1.0e-5) = 0.0059
0.5/getalpha(9.0, 1.0500e-6) - 0.5/getalpha(8.0, 1.0e-5) = 0.0026
0.5/getalpha(9.0, 1.0250e-6) - 0.5/getalpha(8.0, 1.0e-5) = 0.0009
0.5/getalpha(9.0, 1.0250e-6) - 0.5/getalpha(8.0, 1.0e-5) = 0.0009
0.5/getalpha(9.0, 1.0115e-6) - 0.5/getalpha(8.0, 1.0e-5) = 3.1969e-6

Extending the cutoff by 1.0Å allows us to converge the direct space sum almost one order of magnitude more tightly while conserving the precision of the reciprocal space sum.

Go to the next section

Go to the previous section

"How's that for maxed out?"

Last modified: Aug 17, 2018