Particle mesh Ewald: An N⋅log(N) method for Ewald sums in large systems

An N⋅log(N) method for evaluating electrostatic energies and forces of large periodic systems is presented. The method is based on interpolation of the reciprocal space Ewald sums and evaluation of the resulting convolutions using fast Fourier transforms. Timings and accuracies are presented for three large crystalline ionic systems.

We define the functions <I>dir(r;{3) and <I>rec(r;f3). where r is a point in U and {3 is a positive number. by and where erfc(x) is the complementary error function, V=al . a2Xa3 is the volume of the unit cell U, and nand m are given by n=nlal +n2a2+n3a3' and m=mlaT +m2a~+m3a~, for integers nk and mk (k=I,2,3). The effect of {3 on the Ewald pair potential, 'I/J(r;{3), defined by 'I/J(r) = <I>dir(r;{3) + <I>rec(r;{3), is that of an additive constant. Hence. for a neutral system. the total electrostatic energy (and its derivatives) are invariant to {3.
The infinite series defining <I>dir and <I>rec are both rapidly convergent. Their rates of convergence are controlled by adjusting the value of {3. If {3 is chosen so that only the minimum image terms in the direct space sum <I>dir(r;f3) are retained, the total electrostatic energy of a neutral unit cell U, containing N point charges Ql>Q2, ... ,qN, located at positions rl'r2 •...• rN' is given byl-3 use of Wigner potentials, 7 multiple time step ("twin range") methods,S particle-mesh techniques,9 and efficient Taylor and/or multipole expansions.lO-IS The particle-mesh Ewald (PME) method presented here involves choosing {3 sufficiently large that atom pairs for which rij exceeds a specified cutoff (e.g., 9 A..) are negligible in the direct space sum in Eq. (2) which reduces this term to order N. The reciprocal space sum in Eq. (2) is then approximated by a multidimensional piecewiseinterpolation approach inspired by the particle-mesh method of Hockney and Eastwood. 9 The approximate reciprocal energy and forces are expressed as convolutions and can thus be evaluated quickly using fast Fourier transforms (FFTs). The resulting algorithm is of order N In (N), is easily programmed, and is shown below to be efficient and accurate for macromolecular systems.

METHOD
Given positive integers K I , K 2 , andK 3 , we compute <I>rec as well as its (Cartesian) gradient on the grid of fractional coordinates (IlK!, lz/K z , 13/K3), for Ik=I, ... ,K k , k= 1,2,3, and store these quantities into arrays at the beginning of the simulation. Given a pair of atoms i and j, with fractional coordinates fik and fh' k= 1,2,3, the inwhere the functions (}p,K are obtained from the barycentric form of the weights used in (2p-l)th order Lagrangian interpolation using the grid of points (kIK) , k E z.16 More specifically, for real argument x, let [x] denote the largest integer less than or equal to x, and define the in- and finally let (}p,K(x,k) Note that (}p,K(x,k) is nonzero for 2p values of k, and that the sum of (}p,K(x,k) over all k is unity. Since (}p,K(kIK,k) = 1 and (}p,K(iIK,k) =0 for l=l=k, we see that (}p,K(x,k) are continuous in x. Note also that x is always in the central grid interval, with respect to the weights, so that we avoid the ill-conditioning associated with high order polynomial interpolation. The approximate gradient is obtained by replacing <P ree in Eq. (4) by its gradient, obtained by term-by-term differentiation, at each grid point.
Defining Q at the grid points by the expression for the approximate reciprocal sum energy Eree,pU) at atom i due to the atoms of the unit cell U is terpolation approximation (order p) to the reciprocal space pair potential, <l>ree,p, at each point is given in terms of the precomputed array values of <P ree , where <Pree*Q denotes the discrete convolution operator, and we have used the fact that <p ree (r;f3) = <pree( -r;f3).
Since (}p,K(x,k) is continuous in x, the approximate energies and forces will be continuous with respect to particle position, which is not true of the cell multipole methods.
The right hand side of Eq. (9) can be estimated by replacing the sum by an integral. We can change variables =max{[(f i -f i )2/f;jl12; l.;;;i.;;N}; where E is the total electrostatic energy and fl is the force vector at atom i, evaluated using the exact Ewald pair potential for all atom pairs. The" A" symbol indicates the approximate values evaluated with the PME method using a 9 A atom-based list for the direct space contribution. The percent overhead (% overhead) is: timePME/timeCoulombX 100%; where time PME is the time for evaluation of the total nonbond potential energy and forces usiitg the PME method for the electrostatics, and timeCoulomb is the time (2.28 s for B-DNA) required for evaluation of the total nonbond energy and forces using conventional Coulombic interactions.
A similar result can be obtained for the gradient. The error in the interpolation can be made arbitrarily small by fixing al/K I , a2/K2, and a3/K3 to be less than one, and then choosing p sufficiently large. Thus, the quantity KIK2K3 is of order the system size ala2a3 and hence of order N, for any desired tolerance, and so the proposed algorithm is of order N log (N).

RESULTS
The PME method was implemented by modifying the AMBER3.0 (Rev. A) molecular dynamics code in the following way. The direct space pair potential, cI>dir> was computed along with the van der Waals and H-bond terms from an atom-based nonbond list. The erfc(x) function and its derivative were obtained by table look up (relative error <1.0x 10-8 ) for increased efficiency. The approximate reciprocal space pair potential <i>rec was evaluated as described above in a separate routine. An additional cor- 0.5 (48,81,128) 2.0XIo-3 7.5xlO-3 0.17 39.9% 1.1 X 10-4 2.6XlO-4 6.1 X 10-3 45.2% 1.1 X 10-5 1.9X 10-5 4.0XIO-4 58.3% 1.7 X 10-6 5.7XIO-6 1.5 X 10-4 83.8% rection was needed to account for exclusion of nearest images of masked atom pairs (bonded pairs, etc.). All tests were run as single processor jobs on the Cray YMP at the Frederick Cancer Research and Development Center.
The accuracy and computational efficiency of the PME method were tested on several macromolecular crystals. Unit cells were constructed from the space groups using the crystallographic coordinates for the heavy atom positions, as in previous work. 18 ,19 Point charges were assigned using the AMBER force field. 2o Electrostatic energies and forces for each atom are compared to the corresponding values calculated using the exact Ewald pair potential (relative accuracy < 1 X 10 -7 ) . Table I shows the results for an ionic B-DNA crystal. 21 Results are given for grid densities (a l /K 1 xa2/K2 Xa3/K3)1/3 of approximately 1.0, 0.75, and 0.5 A, and interpolation order p= 1,2,,3,4. The computational cost of the PME method is compared to that of normal evaluation of nonbond interactions using a 9 A atom-based Verlet list on a single processor Cray YMP. The overhead above normal nonbond interactions ranges from approximately 16% to 84% corresponding to root mean square (rms) relative force errors of 2.9 X 10-2 and 5.7 X 10-6 , respectively. We conclude that reasonable relative accuracy (about 2 X 10-4 rms force error) can be obtained with approximately a 40% overhead, by using triquintic interpolation (p = 3) on a 0.75 A grid. Similar accuracies and computational costs are obtained for large protein crystals 22 ,23 using comparable grid densities and interpolation orders (Table II). Note the PME method is completely general to nonorthogonal unit cells. Traditional particle-mesh techniques have been criticized as not being able to attain high accuracy efficiently, especially for nonuniform particle distributions. 11 Our results indicate that high precision (rms relative force error < 1.0 X 10-5 ) is easily obtained for macromolecular systems by using higher order interpolation with only a modest increase in computational time. Although the PME method has increased memory requirements over conventional nonbond list-based methods, the cost of memory appears to be rapidly decreasing, and hence may not be an issue in the future.
The PME method offers several advantages as a method for the treatment of long-range forces in macromolecular systems. These include (1) High accuracy: High accuracy (z 5 X 10-6 relative force) can be obtained with relatively little increase in computational effort.
(2) East o/implementation: The PME method can be efficiently implemented into conventional MD algorithms such as AMBER which use a Verlet list.
(3) Continuity: The PME pair potential and its derivatives are continuous functions of position, regardless of the accuracy required, and thus avoid problems involved with integration of discontinuous functions.
( 4) Efficiency: The PME method is fast. For large macromolecular systems, the PME method requires only about a 40% overhead over conventional truncated listbased methods to obtain relative force accuracies of z2 X 10-4 • The FORTRAN subroutines for performing the PME approximate energies and forces are available upon request from the authors.