Multiple Staggered Mesh Ewald: Boosting the Accuracy of the Smooth

Subscriber access provided by CORNELL UNIVERSITY LIBRARY

Article

Multiple Staggered Mesh Ewald: Boosting the Accuracy of the Smooth Particle Mesh Ewald Method Han Wang, Xingyu Gao, and Jun Fang J. Chem. Theory Comput., Just Accepted Manuscript • DOI: 10.1021/acs.jctc.6b00701 • Publication Date (Web): 19 Oct 2016 Downloaded from http://pubs.acs.org on October 21, 2016

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

Journal of Chemical Theory and Computation is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 40

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

Multiple Staggered Mesh Ewald: Boosting the Accuracy of the Smooth Particle Mesh Ewald Method Han Wang,∗,†,‡ Xingyu Gao,¶,†,‡ and Jun Fang†,‡ †Institute of Applied Physics and Computational Mathematics, Fenghao East Road 2, Beijing 100094, P.R. China ‡CAEP Software Center for High Performance Numerical Simulation, Huayuan Road 6, Beijing 100088, P.R. China ¶Laboratory of Computational Physics, Huayuan Road 6, Beijing 100088, P.R. China E-mail: [email protected]

Abstract The smooth particle mesh Ewald (SPME) method is the standard method for computing the electrostatic interactions in the molecular simulations. In this work, we develop the multiple staggered mesh Ewald (MSME) method, which averages the SPME forces computed on, e.g. M , staggered meshes. We prove, from a theoretical perspective, that the MSME is as accurate as the SPME, but uses M 2 times less mesh points in a certain parameter range. In the complementary parameter range, the MSME is as accurate as the SPME with twice of the interpolation order. The theoretical conclusions are numerically validated both by a uniform and uncorrelated charge system, and by a three-point-charge water system that is widely used as solvent for the biomacromolecules.

1

ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

1

Page 2 of 40

Introduction

The electrostatic interaction is one of the most important interactions, and, perhaps, the most computationally demanding part in molecular dynamics (MD) simulations. One popular way of computing the electrostatic interaction is the Ewald summation. 1 It was shown that the optimal computational complexity of this method is O(N 3/2 ) 2 (N being the number of point charges in the system), therefore, as the number of charges increases, e.g. to several hundreds, 3 the Ewald summation becomes relatively expensive. This stimulates the development of the Ewald-based fast algorithms like the particle mesh Ewald (PME) method, 4 the smooth particle mesh Ewald (SPME) method, 5 the particle-particle particlemesh (PPPM) method 6,7 and the methods based on the non-equispaced fast Fourier transforms. 8 All these fast methods reduce the computational complexity to O(N log N ) by interpolating the charges on a mesh and solving the Poisson’s equation with the fast Fourier transform (FFT). The Ewald-based fast methods were also proposed to compute the dispersion interactions in the systems that present interfaces. 9–11 They were shown to be more accurate and even faster than the traditional treatment of the dispersion interactions in these systems, viz. using a large cut-off radius. 12,13 The standard way to improve the accuracy of the Ewald-based fast methods (in the reciprocal space) is to refine the FFT mesh or to use interpolation bases of higher orders. An alternative way is to average the interactions computed on multiple copies of meshes, so that the aliasing errors in different phases could be canceled. This idea was discussed in the one dimensional case in the textbook of the PPPM by Hockney and Eastwood. 6 They also investigated the three dimensional particle simulation using two interlaced meshes in detail. Recently, the idea of two interlaced meshes was revived in the SPME method 14 (called the staggered mesh Ewald method) and the PPPM method. 15 The theoretical error estimates for the methods 15,16 showed that the staggered meshes substantially improve the accuracy compared with the non-staggered mesh counterparts.

By far, to the best of

our knowledge, only the two-mesh approaches have been investigated and proved to be 2


Page 3 of 40

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


effective. The questions on how to construct the multiple (more than two) meshes in the three dimensional particle simulations, how much the accuracy can be improved, and how many meshes should be used in practice have not been answered. In this work, we develop the multiple staggered mesh Ewald (MSME) method. It takes the average of the reciprocal forces computed by the SPME method (with a modified influence function) on M meshes, which are shifted to the M equally partitioning points of the mesh subcell diagonal. The MSME method uses M times more FFT mesh points, and achieves, in a certain range of the splitting parameter, the same accuracy that would need M 3 times more FFT mesh points in the SPME. In the complementary range of the splitting parameter, the MSME achieves the same accuracy as the SPME that uses twice of the interpolation order. These properties of the MSME are understood by an a priori theoretical error estimate. We prove that the first order part of the error is identical to that of the SPME with a mesh M times finer in each direction, while the second order part is roughly the same as that of the SPME using twice of the interpolation order. The first order error dominates the reciprocal error in the case of M = 1, so the M = 2 MSME is always more accurate than M = 1. When M > 1, either the first or the second order error may dominate the reciprocal error, depending on the splitting parameter. As the M increases, the first order error decreases while the second order error does not, thus the reciprocal error reduces but will eventually reach the lower bound that is defined by the second order error. The proposed error estimate is numerically validated both by a uniform and uncorrelated point charge system, and by a rigid three-point-charge water model that is widely used as solvent for biological macromolecules. The main difficulty of applying the Ewald-based fast methods in practical simulations is how to determine the working parameters. It is well known that the arbitrarily chosen parameters may lead to a substantial slowdown of the computation, or results that are several orders of magnitude less accurate. This problem may be solved, in a posteriori manner, by comparing the computed forces while scanning the parameter space for a representative

3



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

snapshot of the system. 17 An alternative way is by a parameter tuning algorithm that automatically determines the most efficient combination of parameters under the restraint of a prerequisite accuracy. 18 The success of this algorithm relies on the quality of the a priori error estimate that quantitatively describes the accuracy of the Ewald-based fast methods as a function of working parameters, and a large amount of work have been dedicated to this direction. 15,16,18–24 Using the error estimate developed by this work, we propose a parameter tuning algorithm that provides a practical guide on how to select the working parameters (including the number of meshes) for the MSME method. This paper is organized as follows. In Section 2, the Ewald summation and the particle mesh Ewald method are briefly introduced. The notations are setup, and the relation between our implementation of the SPME, the original SPME and the PPPM methods are discussed. In Section 3, the MSME method is introduced, and its accuracy is numerically discussed and compared to the SPME method. In Section 4, the numerical phenomena of MSME is analyzed by the error estimate. Section 5 validates the quality of the error estimate by comparing it to the numerically computed error. Section 6 proposes a parameter tuning algorithm that automatically searches for the most time-saving parameters under the constraint of a prerequisite accuracy. In Section 7, the work is concluded and the performance issue is discussed in more detail.

2

The Ewald summation and the particle mesh Ewald method

We denote the N charged particles in the system by {q1 , · · · , qN }, and their positions by {r 1 , · · · , r N }. If the system is subject to the periodic boundary condition, then the electro-

4


Page 4 of 40

Page 5 of 40

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


static interaction is given by: ∗

N

1 X X qi qj E= , 2 n i,j=1 |r ij + n|

(1)

where r ij = r i − r j , n = n1 a 1 + n2 a 2 + n3 a 3 is the lattice with (n1 , n2 , n3 ) ∈ Z3 and (a 1 , a 2 , a 3 ) are the unit cell vectors. When n = 0, the inner summation runs over all Coulomb interactions between the particle pairs in the unit cell, otherwise, it is counting the interaction between the unit cell and its periodic images. The “∗” over the outer summation means that when n = 0, i 6= j. The Ewald summation decomposes the electrostatic interaction (1) into three parts, the direct part, the reciprocal part and the correction part

E = Edir + Erec + Ecorrection ,

(2)

with the definitions ∗

N

Edir

1 X X qi qj erfc(β|r ij + n|) , = 2 n i,j=1 |r ij + n|

(3)

Erec

1 X exp(−π 2 m 2 /β 2 ) S(m)S(−m), = 2πV m 6=0 m2

(4)

N β X 2 = −√ q , π i=1 i

(5)

Ecorrection

where m = m1 a ∗1 + m2 a ∗2 + m3 a ∗3 is the reciprocal space lattice with (m1 , m2 , m3 ) ∈ Z3 and a ∗α , α = 1, 2, 3 are the conjugate reciprocal vectors that are defined by relations a α ·a ∗β = δαβ , α, β = 1, 2, 3. The direction index “β” here should be clearly distinguished with the splitting parameter β in Eqs. (3)–(5) by the context. V is the volume of the unit cell calculated by V = a 1 · (a 2 × a 3 ). The erfc(x) in Eq. (3) is the complementary error function, and the

5



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 6 of 40

S(m) in Eq. (4) is the structure factor defined by

S(m) =

N X

qj e2πim ·r j .

(6)

j=1

The “i” at the exponent is the imaginary unit that should not be confused with the particle index i. The complementary error function in the direct energy (3) converges exponentially fast w.r.t. increasing distance |r ij + n|, thus it can be cut off at a certain radius rc . As a consequence, one needs to sum a finite number of terms in (3), and the cost of this computation is O(N ) by using the standard neighbor list algorithm. 25 The summation in the reciprocal energy (4) also converges exponentially fast, so it can be truncated at the reciprocal space cut-off: −Kα /2 ≤ mα < Kα /2. To achieve a prerequisite accuracy, the number of Fourier modes in summation (4) should be taken the same order as the number of charged particles in the system, i.e. K1 K2 K3 ∼ N . Therefore, the computational complexity of the reciprocal interaction is O(N 2 ). We temporally assume that N = K1 K2 K3 , and that the charged particles locate on a K1 × K2 × K3 mesh, then the structure factor (6) is nothing but a discrete Fourier transform and can be computed at the cost of O(N log N ) by using the FFT. In MD simulations, the particles do not necessarily locate on a mesh, so the fast methods interpolate the particle charges on the uniform mesh points, and then use FFT to accelerate the computation. The computational complexity of the interpolation and the FFT are O(N ) and O(N log N ), respectively. The total computational expense of the fast methods is O(N log N ). Before developing the MSME method, we briefly introduce the SPME method and refer the readers to Ref. 5 for more details. We let the numbers of mesh points in three directions be Kα , α = 1, 2, 3. For the coordinate of a particle r , we introduce the notation uα = Kα a ∗α · r , and it is clear that 0 ≤ uα < Kα . We approximate the complex exponential e2πimα uα /Kα measured on particle position uα by a linear combination of the complex exponentials measured

6


Page 7 of 40

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


on the mesh (see Eq. (S15) in the Supporting Information): e2πimα uα /Kα ≈

1 X 1 ϕn (uα − lα ) e2πimα lα /Kα , Kα l ∈I ϕˆn (mα ) α

(7)

K

where IK = { l | l ∈ Z, −K/2 ≤ l < K/2}. ϕn denotes the n-th order cardinal B-spline that is defined by the recursive formula:

ϕ1 (x) = χ[− 1 , 1 ] (x), 2 2

ϕn (x) = Kα ϕn−1 ∗ ϕ1 (x),

(8)

where χ[− 1 , 1 ] denotes the characteristic function on interval [− 12 , 12 ]. ϕˆn (mα ) is the Fourier 2 2

transform of ϕn , given by πmα in 1 h sinc( ) . ϕˆn (mα ) = Kα Kα

(9)

By using the approximation (7) in three dimensional space, we reach

e

2πim ·r

m1 l 1 m2 l2 m3 l 3 B(m) X Pr (l1 , l2 , l3 ) exp[2πi( + + )] , ≈ K1 K2 K3 l ,l ,l K1 K2 K3

(10)

1 2 3

where 1 , ϕˆn (mα ) α Y Pr (l1 , l2 , l3 ) = ϕn (uα − lα ). B(m) =

Y

(11) (12)

α

Inserting the interpolation (10) into Eq. (4), we have the fast method of computing the reciprocal energy,

Erec ≈

X

l1 ,l2 ,l3

Q(l1 , l2 , l3 )[ Q ∗ (F B 2 )∨ ](l1 , l2 , l3 ),

7


(13)


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 8 of 40

where  exp(−π 2 m 2 /β 2 )   1 m2 F (m) = ×  2πV 0

|m| = 6 0,

(14)

|m| = 0,

and Q is the charge distribution defined on the mesh

Q(l1 , l2 , l3 ) =

X

qj Pr j (l1 , l2 , l3 ).

(15)

j

The symbol “∨” denotes the backward Fourier transform. The operator “∗” denotes the ˆ × (F B 2 ) ]∨ with “∧” denoting convolution, which can be computed by Q ∗ (F B 2 )∨ = [ Q the forward Fourier transform. The r.h.s. of the equation is roughly explained as firstly ˆ solving the Poisson’s equation converting the charge distribution to the reciprocal space (Q), ˆ × F ), and then transforming the result back ([ Q ˆ × (F B 2 ) ]∨ ). The with smeared source (Q function B appears due to the interpolation of the point charge distribution. By using the FFT, the computational expense of the convolution is of order O(N log N ). The force of a particle can be computed in two ways. The first is known as the ikdifferentiation, which takes the negative gradient w.r.t. particle coordinate on the reciprocal energy (4), and then interpolates the complex exponential e2πim ·r . This yields F ik rec,i ≈ qi

X

l1 ,l2 ,l3

Pr i (l1 , l2 , l3 )[ Q ∗ (GB 2 )∨ ](l1 , l2 , l3 ),

(16)

where

G(m) = −4πimF (m).

(17)

The alternative way is called the analytical differentiation, which takes the negative gradient

8


Page 9 of 40

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


on the approximated reciprocal energy (13), and leads to F ad rec,i ≈ qi

X

l1 ,l2 ,l3

−2∇r i Pr i (l1 , l2 , l3 )[ Q ∗ (F B 2 )∨ ](l1 , l2 , l3 ).

(18)

ˆ and three backwards Remark 1: The ik-differentiation needs four FFTs: one forward for Q, 2 ˆ for the three components of Q(GB ). The three backward FFTs are independent and of

the same size, so they can be executed as a multiple variable FFT, which usually saves time comparing with executing three single variable FFTs consecutively. Remark 2: The function B defined by Eq. (11) should be distinguished with that defined in the original SPME paper, 5 saying

BO-SPME (m) =

Y α

e2πi(n−1)mα /Kα , Pn−2 2πimα k/Kα ϕ (k + 1)e n k=0

(19)

which can be equivalently written as 26

BO-SPME (m) =

Y α

P

1 . ˆn (mα + lα Kα ) lα ϕ

(20)

In this work, we still call the reciprocal computation on one mesh by the SPME method, although the function B is different from the original one. Other options on the function B are available, for example, the PPPM method, which minimizes the error of the reciprocal force computation in homogeneous and uncorrelated charge systems, derives a function B of

BPPPM (m) =

Y α

P

ϕˆn (mα ) ˆ2n (mα + lα Kα ) lα ϕ

(21)

for the ik-differentiation. 7,26 Our implementation of the function B (Eq. (11)) is much closer to the PPPM method than to the original SPME method. Therefore, the accuracy of the reciprocal computation on one mesh in our work is higher than the original SPME, and is very close to the PPPM method. A detailed comparison is provided in the Appendix A.

9



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

3

Multiple staggered mesh Ewald method

2/3 1/3

Figure 1: The schematic plot of a triple staggered mesh Ewald. The F i with different color on the r.h.s. of the equation denotes the reciprocal SPME force computed on the mesh with the same color. The MSME method averages the reciprocal interactions computed by SPME on M identical meshes, which locate on the M equally partitioning points of the mesh subcell diagonal, as illustrated by Fig. 1. The reciprocal force can be computed by either ik- or analytical differentiation, but all meshes in MSME should use the identical force scheme. These two branches of MSME are named by IK-MSME and AD-MSME, respectively. The MSME method can be easily implemented by using the existing SPME codes (with function B slightly modified), because the computation on each mesh is the SPME and the different meshes are independent. We claim that the reciprocal accuracy of MSME that uses M > 1 meshes is either of the following two cases: I. In a certain range of the splitting parameter β, the accuracy is the same as the SPME that refines the mesh in each direction by M times, i.e. using an M K1 × M K2 × M K3 mesh. II. In the complementary range of β, the accuracy is almost the same as the SPME that uses twice order of the B-spline interpolation (2n in our notation). These claims will be carefully checked both by numerical examples and by theoretical accuracy analyses. 10


Page 10 of 40

Page 11 of 40

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


The numerical investigation of the MSME is carried out in two testing systems: • System 1: The simulation region is a cubic cell of size 3.724nm × 3.724nm × 3.724nm. In total, 5184 charged particles are uniformly and independently distributed in the region. Each of the 1728 particles carries a negative partial charge of −0.834 e, while each of the rest (3456) particles carries a positive partial charge of 0.417 e. The total charge of the system is neutral. • System 2: The simulation region is a cubic cell of size 3.724nm × 3.724nm × 3.724nm, and contains 1728 TIP3P 27 water molecules. The partial charges on the hydrogen and oxygen atoms are 0.417 e and −0.834 e, respectively. The configuration was taken from an equilibrium NPT simulation 28 at 300 K and 1 Bar. It should be noted that the number of particles and the amount of the partial charge of each particle are exactly the same as the System 1. The difference is that the charges in System 1 are randomly and independently distributed, while the charges in System 2 are bonded and correlated according to the equilibrium water system configuration. All the numerical studies in this work are carried out by using our in-house MD package MOASP. The accuracy in computing the electrostatic interactions is measured by the root mean square (RMS) force error that is defined by

E=

p h|∆F |2 i,

∆F = F − F ∗ ,

(22)

where F and F ∗ are the computed and accurate electrostatic forces, respectively. The RMS force error will be called the “error” in short. The direct and reciprocal errors are defined similarly by replacing the electrostatic force in Eq. (22) with the direct and reciprocal forces, respectively. The direct force of MSME is computed in the same way as the Ewald summation, the accuracy of which has already been studied. 19,24 Therefore, if not stated otherwise, the error estimate is developed only for the reciprocal error. 11



10

1

RMS Force Error [ kJ/mol/nm ]

100

System 1

n=4,Kα=32

-1

n=4,Kα=64

10

10-2

n=4,Kα=128

-3

10

n=4,Kα=256

-4

10

-5

10

n = 4, Kα = 32, M = 1 n = 4, Kα = 32, M = 2 n = 4, Kα = 32, M = 4 n = 4, Kα = 32, M = 8

-6

10

n=8,Kα=32

-7

10

1

2

3

4

5

-1

β [nm ]

Figure 2: The reciprocal error of the IK-MSME as a function of the splitting parameter β investigated in System 1. The chromatic solid lines are errors computed from the IK-SPME method (parameters presented along with the lines), while the red “+”, the green “×”, the blue “⊡” and the pink “⊙” denote the reciprocal errors computed from the IK-MSME with M = 1, 2, 4 and 8, respectively. The M = 1 IK-MSME is identical to the IK-SPME method.

10

1

100 RMS Force Error [ kJ/mol/nm ]

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

System 2

n=4,Kα=32

-1

10

n=4,Kα=64

10-2

n=4,Kα=128

-3

10

n=4,Kα=256

-4

10

10-5 n = 4, Kα = 32, M = 1 n = 4, Kα = 32, M = 2 n = 4, Kα = 32, M = 4 n = 4, Kα = 32, M = 8

-6

10

10-7

n=8,Kα=32 1

2

3

4

5

-1

β [nm ]

Figure 3: The reciprocal error of the IK-MSME as a function of the splitting parameter β investigated in System 2. The meaning of the symbols is the same as Fig. 2.

12


Page 12 of 40

102 10 RMS Force Error [ kJ/mol/nm ]

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


1

System 1

n=4,Kα=32

n=4,Kα=64

100 -1

10

n=4,Kα=128

-2

10

n=4,Kα=256

-3

10

10-4 n = 4, Kα = 32, M = 1 n = 4, Kα = 32, M = 2 n = 4, Kα = 32, M = 4 n = 4, Kα = 32, M = 8

10-5 10-6

n=8,Kα=32 1

2

3

4

5

β [nm-1]

Figure 4: The reciprocal error of the AD-MSME as a function of the splitting parameter β investigated in System 1. The correction of the “self-interaction” has been applied (see Sec. 4.2 for details). The chromatic solid lines are errors computed from the AD-SPME method (parameters presented along with the lines), while the red “+”, the green “×”, the blue “⊡” and the pink “⊙” denote the reciprocal errors computed from the AD-MSME with M = 1, 2, 4 and 8, respectively. The M = 1 AD-MSME is identical to the AD-SPME method.

10

1


Page 13 of 40

System 2

n=4,Kα=32 n=4,Kα=64

-1

10

n=4,Kα=128

10-2

n=4,Kα=256

10-3 -4

10

-5

10

n = 4, Kα = 32, M = 1 n = 4, Kα = 32, M = 2 n = 4, Kα = 32, M = 4 n = 4, Kα = 32, M = 8

-6

10

n=8,Kα=32

-7

10

1

2

3

4

5

-1

β [nm ]

Figure 5: The reciprocal error of the AD-MSME as a function of the splitting parameter β investigated in System 2. The correction of the “self-interaction” has been applied (see Sec. 4.2 for details). The meaning of the symbols is the same as Fig. 4.

13



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The reciprocal errors of the IK-MSME in Systems 1 and 2 are presented in Figs. 2 and 3, respectively. The reciprocal errors of AD-MSME in Systems 1 and 2 are presented in Figs. 4 and 5, respectively. For both force schemes, we investigate the number of meshes M = 1, 2, 4, 8 with the interpolation order n = 4 and the number of mesh points K1 = K2 = K3 = 32. As a comparison, we present the reciprocal error of the SPME that uses n = 8, Kα = 32 and n = 4, Kα = 32, 64, 128, 256. The M = 1 MSME is identical to the SPME that uses the same mesh (Kα = 32) and interpolation order (n = 4). The accurate forces are computed by well-converged Ewald summations. In the Figures, the points present the reciprocal error of the MSME method, and the solid lines present the reciprocal error of the SPME method. The numerical phenomena of both force schemes in both systems are similar: At relatively small β, the accuracy of the MSME that uses M > 1 meshes matches the SPME that uses a mesh refined M times in three directions, while at relatively large β, the MSME error of M > 1 roughly follows the error of SPME that uses twice of the B-spline interpolation order. These numerical results are clear evidences that support our claims I and II. It is noticed that there is a clear transition of β, at which the accuracy behavior of M > 1 MSME changes from Claim I to Claim II. The transition happens roughly at the crossover between the SPME that uses an n-th order interpolation and an M K1 × M K2 × M K3 mesh, and the SPME that uses a 2n-th order interpolation and a K1 × K2 × K3 mesh. As the number of the staggered meshes M increases, the transition moves toward the smaller β side. Therefore, the accuracy of the MSME cannot be infinitely improved by using increasingly many meshes, because any β > 0 will eventually falls in the region of Claim II, which actually defines the upper bound for the accuracy of the MSME.

14


Page 14 of 40

Page 15 of 40

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


4

Discussion of the numerical phenomena and the error estimate

The numerical phenomena of the MSME observed in the testing systems are understood by the error estimate, which is derived under the framework proposed in Ref. 16 It has been shown that the error of a pairwise interaction (electrostatic interaction is pairwise) is composed of three additive parts, which are the homogeneity, inhomogeneity and correlation errors 2 2 E 2 = Ehomo + Einhomo + Ecorr .

(23)

The homogeneity error stems from the fluctuation of the error force ∆F (see Eq. (22)), while the inhomogeneity error is due to the inhomogeneous charge distribution and is essentially the bias of ∆F . If the positions of the charges in the system are correlated (due to e.g. covalent bonds, hydrogen bonds, van der Waals interaction and so on), then the correlation error arises. The homogeneity and inhomogeneity error contributions are positive definite, while the correlation error contribution may be negative, which means that the charge correlation reduces the error in force computation. 16 In this work, for simplicity, we assume that the charges are uniformly distributed, so the inhomogeneity error vanishes. The error estimate that only includes the homogeneity error is precise for the systems, in which the charges are uniformly and independently distributed (e.g. System 1). Besides the homogeneity error, we estimate the correlation error by using the nearest neighbor approximation technique 16 for System 2 (TIP3P water). This correlation error estimate can be directly extended to other rigid water models containing three point charges, like the SPC, 29 SPCE 30 and TIP4P 27 water models. The error estimate is provided without proof, and the readers are referred to the Supporting Information for more details on the derivation of the error estimate. The Supporting Information also provides the error kernels for the IK- and AD-MSME, from which the error estimates can be easily extended

15



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 16 of 40

to systems that have non-uniform charge distributions.

4.1

ik-differentiation

Firstly we introduce the short-hand notation

Zα,l (m) =

ϕˆn (mα + lKα ) . ϕˆn (mα )

(24)

By using the definition of ϕˆn , we have the estimate |Zα,l (m)| ≤ |mα /(mα + lKα )|n . Noticing that −Kα /2 ≤ mα < Kα /2 and l 6= 0, |Zα,l (m)| is fast decaying w.r.t. growing |l|. We further denote

(25)

G α,l (m) = G(m)Zα,l (m),

(26)

G α,l1 ;β,l2 (m) = G(m)Zα,l1 (m)Zβ,l2 (m).

G α,l (m) has the first power of Zα,l (m), and G α,l1 ;β,l2 (m) has the second power of Zα,l (m). Therefore, in general, we have |G α,l (m)| ≫ |G α,l1 ;β,l2 (m)| for l 6= 0. The homogeneity error is estimated, for IK-MSME, by (see Eqs. (S52)–(S54)) ik,(1)

ik,(2)

ik |Ehomo |2 ≈ |Ehomo |2 + |Ehomo |2 .

(27)

ik,(2)

ik,(1)

On the r.h.s., Ehomo and Ehomo denote the first and second order homogeneity errors, respectively, and are defined by ik,(1)

|Ehomo |2 = 2q 2 Q2 ik,(2)

|Ehomo |2 = q 2 Q2 +

XXX α

l6=0

θM (l)|G α,l (m)|2 ,

(28)

m

X X X

θM (l1 + l2 )|G α,l1 ;β,l2 (m)|2

α6=β l1 ,l2 6=0 m

X X X α,β l1 ,l2 6=0 m

θM (l1 − l2 )|G α,l1 ;β,l2 (m)| ,

16


2

(29)

Page 17 of 40

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


where Q2 = by

P

2 i qi ,

q 2 = Q2 /N . In the estimates (28) and (29), the function θM is defined

θM (l) =

   0   1

(l mod M ) 6= 0,

(30)

(l mod M ) = 0,

and is introduced by using the multiple staggered meshes. For the IK-MSME with M = 1 ik,(1)

ik,(2)

and θM (l) ≡ 1, noticing |G α,l (m)| ≫ |G α,l1 ;β,l2 (m)|, we have |Ehomo | ≫ |Ehomo |, thus the first order error dominates the reciprocal error. This is why the error estimate of SPME (M = 1 MSME) only considered the first order contribution. 18 In the cases of non-trivial IK-MSME where M > 1, we compute the first order error (28) by ik,(1)

|Ehomo (β, n, K , M )|2 XXX = 2q 2 Q2 θM (l)|G α,l (m)|2 α

= 2q 2 Q2

l6=0

m

XXX α

l6=0

m

|G α,M l (m)|2

X X X ϕˆn (mα + M lKα ) 2 = 2q 2 Q2 G(m) ϕ ˆ (m ) n α α l6=0 m ik,(1)

= |Ehomo (β, n, M K , 1)|2 .

(31)

Here we explicitly write the dependency of the error estimate on the working parameters: β (splitting parameter), n (B-spline order), K = (K1 , K2 , K3 ) (number of mesh points) and M (number of meshes). The last equation holds because of the identity M lKα = l(M Kα ). Thus the leading order error of MSME using M identical K1 × K2 × K3 meshes is the same as the IK-SPME using an M K1 × M K2 × M K3 mesh. In the estimate of the second order error (29), only the terms with α = β contribute significantly (will be shown by numerical examples later), and in this case l1 6= l2 terms are 17



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 18 of 40

much smaller than the l1 = l2 terms due to the fast decaying of Zα,l (m) w.r.t. |l|. Therefore, we have ik,(2)

|Ehomo (β, n, K , M )|2 XXX ≈ |G α,l;α,l (m)|2 α

l6=0

m

h ϕˆ (m + lK ) i2 2 X X X n α α 2 2 =q Q G(m) ϕˆn (mα ) α l6=0 m X X X ϕˆ2n (mα + lKα ) 2 = q 2 Q2 G(m) ϕˆ2n (mα ) α l6=0 m =

1 ik,(1) (β, 2n, K , 1)|2 . |E 2 homo

(32)

The second equation holds because of the definition of ϕˆn (m), i.e. Eq. (9). Eq. (32) means √ that the second order error of the IK-MSME is approximately 1/ 2 times of the first order error of the IK-SPME method using twice order of the B-spline interpolation. As the number of staggered meshes M increases, the first order error decreases, while the second order error does not, therefore, the second order error defines the lower bound of the total reciprocal error of the MSME. The Claim II (see Sec. 3) holds in the β range that the first order error decreases to smaller than the second order error. The Claim I holds in the complementary range, in which the first order error still dominates. The boundary of the ik,(1)

ik,(2)

ranges can be determined by solving the equation Ehomo (β, n, K , M ) = Ehomo (β, n, K , M ) w.r.t. variable β. This crossover moves towards the smaller β side as M increases, because ik,(2)

the second order error increases monotonically w.r.t. β (the derivative of Ehomo w.r.t. β is positive definite). Since the M = 1 MSME is dominated by the first order error for all β, the M = 2 MSME is always more accurate than M = 1. Increasing the M larger than 2 may not be able to further improve the accuracy, depending on which part of the error dominates for a certain β. The correlation error in the water system contains bonded and non-bonded contributions. The former is easily estimated by the nearest neighbor approximation with the knowledge 18


Page 19 of 40

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


of the O-H bond length and H-O-H angle, while the latter is partially estimated by the same technique with the knowledge of the radial distribution functions (RDFs). It should be noted that the estimate for the bonded correlation error is a priori, because the values of O-H bond length and H-O-H angle are set up by the water model, and are constrained over the simulations (for the rigid water models). By contrast, the estimate for the non-bonded correlation error is a posteriori, because the RDFs are, in general, computed out of the molecular configurations sampled by the simulation. Ref. 16 showed that, for a water system, the bonded contribution accounts for the substantial part of the correlation error, and is good enough for the applications like the parameter tuning. 16 Therefore, we only consider the nearest neighbor approximation for the bonded correlation error. For the MSME using the ik-differentiation, the nearest neighbor approximation to the correlation error in a three-point-charge water system is (see Eq. (S65)) ik Ecorr = q 2 Q2

XXX α

+ q 2 Q2

l6=0

θM (l)T w (m)|G α,l (m)|2

m

XXX α

l6=0

θM (l)T w (m + lKα a ∗α )|G α,l (m)|2 ,

(33)

m

where the function T w provides the information of the bonded charge correlation in the water molecule, and is defined by T w (m) =

4qH qO 2qH2 T (m) + T (m). s O 2 2 sH 2qH2 + qO 2qH2 + qO

(34)

s O is the vector connecting the oxygen and the hydrogen atoms, and s H is the vector connecting two hydrogen atoms. The function Tb (m) is the structure factor of bond b averaged over all possible directions, therefore, it only depends on the size of the Fourier mode m = |m|, and the bond length b = |b|, and writes Tb (m) = he2πim ·b idirections = 19

sin(2πmb) . 2πmb


(35)


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 20 of 40

In the ensemble average in Eq. (35), we assumed that all directions are equally possible. For TIP3P water model, |s O | = 0.09572 nm, and |s H | = 0.15139 nm.

4.2

Analytical differentiation

An extra error due to the “self-interaction” presents in the AD-SPME force computation, and can be removed by subtracting the analytic formula of the self-interaction from the force. 16,31 For the AD-MSME, the self-interaction is given by (see Eq. (S79)) F self = qi2 i

XXX α

l6=0

θM (l)F α, (m) e2πiluα

m

+ qi2

X X X

θM (l1 + l2 ) F α,l1 ;β,l2 (m) e2πi(l1 uα +l2 uβ )

+ qi2

X X X

θM (l1 − l2 ) F α,l1 ;β,l2 (m) e2πi(l1 uα −l2 uβ ) ,

α6=β l1 ,l2 6=0 m

α,β l1 ,l2 6=0 m

(36)

where we introduce the notations

F α,l (m) = −4πilKα a ∗α F (m)Zα,l (m), F α,l1 ;β,l2 (m) = −4πil1 Kα a ∗α F (m)Zα,l1 (m)Zβ,l2 (m).

(37) (38)

We also have |F α,l (m)| ≫ |F α,l1 ;β,l2 (m)| for l 6= 0. When the self-interaction is removed, the homogeneity error of the AD-MSME is estimated as (see Eqs. (S93)–(S95)) ad,(1)

ad,(2)

ad |Ehomo |2 ≈ |Ehomo |2 + |Ehomo |2 ,

20


(39)

Page 21 of 40

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


where, similar to the IK-MSME, the homogeneity error of the AD-MSME is also composed of the first and second order contributions, which are defined by ad,(1)

|Ehomo |2 = q 2 Q2 ad,(2)

|Ehomo |2 = q 2 Q2

XXX α

l6=0

m

i h θM (l) |G α,l (m)|2 + |G α,l (m) + F α,l (m)|2 .

X X X

θM (l1 + l2 )

α6=β l1 ,l2 6=0 m

h1 2

(40)

|G α,l1 ;β,l2 (m)|2

1 + | G α,l1 ;β,l2 (m) + F α,l1 ;β,l2 (m)|2 2 1 1 i + G α,l1 ;β,l2 (m) + F α,l1 ;β,l2 (m) · G α,l1 ;β,l2 (m) + F β,l2 ;α,l1 (m) 2 2 X X X 2 + (41) θM (l1 − l2 )|G α,l1 ;β,l2 (m) + F α,l1 ;β,l2 (m)| . α,β l1 ,l2 6=0 m

The first order error in the estimate is ad,(1)

|Ehomo (β, n, K , M )|2 i h XXX = q 2 Q2 θM (l) |G α,l (m)|2 + |G α,l (m) + F α,l (m)|2 α

l6=0

m

= q 2 Q2

XXXh

2

XXXh

α

=q Q

2

α

l6=0

l6=0

m

m

|G α,M l (m)|2 + |G α,M l (m) + F α,M l (m)|2 2

|G(m)| + |G(m) −

i

4πiM lKα a ∗α F (m)|2

ad,(1)

= |Ehomo (β, n, M K , 1)|2 .

i

ϕˆ (m + M lK ) 2 n α α × ϕˆn (mα )

(42)

Thus the first order error of AD-MSME using M identical K1 × K2 × K3 meshes is the same as the AD-SPME using a mesh M times finer in each direction, i.e. an M K1 × M K2 × M K3 mesh. In the second order error, only the terms with α = β, l1 = l2 contribute significantly, so

21



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 22 of 40

we have the approximation ad,(2)

|Ehomo (β, n, K , M )|2 XXX ≈ |G α,l;α,l (m) + F α,l;α,l (m)|2 α

2

=q Q

m

l6=0

2

XXX

|G(m) −

XXX

|G(m) − 4πilKα a ∗α F (m)|2

α

= q 2 Q2

α

l6=0

l6=0

m

m

ad,(1)

= |Ehomo (β, 2n, K , 1)|2 −

4πilKα a ∗α F (m)|2

h ϕˆ (m + lK ) i4 n α α ϕˆn (mα ) h ϕˆ (m + lK ) i2 2n

α

α

ϕˆ2n (mα )

1 ik,(1) |Ehomo (β, 2n, K , 1)|2 2

ad,(1)

≈ |Ehomo (β, 2n, K , 1)|2 .

(43)

The last approximation holds because the IK-SPME is usually much more accurate than the AD-SPME when the parameters are the same (see, e.g. Figs. 2 and 4). Eq. (43) means that the second order error of AD-MSME is roughly the same as the error of the AD-SPME with twice order of the B-spline interpolation. Finally, the nearest neighbor approximation to the correlation error of AD-MSME is given by (see Eq. (S101)) ad Ecorr = q 2 Q2

XXX α

+ q 2 Q2

l6=0

m

XXX α

5

θM (l) T w (m)|G α,l (m) + F α,l (m)|2

l6=0

θM (l) T w (m + lKα a ∗α )|G α,l (m)|2 .

(44)

m

Numerical validation of the error estimate

In this section we numerically check the quality of the error estimate. The reciprocal error of the MSME and the corresponding error estimate are shown in Figs. 6–9. The reciprocal error is obtained by comparing the MSME reciprocal force with a well converged Ewald reciprocal force using the same splitting parameter β. For System 1 (see Figs. 6 and 8), the

22


101 100 RMS Force Error [ kJ/mol/nm ]

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


System 1

10-1 10-2 -3

10

M=1 M=2 M=4 M=8 M=1, estimated M=2, estimated M=4, estimated M=8, estimated

-4

10

10-5 -6

10

10-7

1

2

3

4

5

β [nm-1]

Figure 6: The reciprocal error (points) of the IK-MSME and the error estimate (solid lines) as a function of the splitting parameter β investigated in System 1. The number of mesh points is Kα = 32 and the order of B-spline interpolation is n = 4. The IK-MSME with M = 1, 2, 4 and 8 are plotted with color red, green, blue and pink, respectively.

10

1


Page 23 of 40

System 2

10-1 10-2 -3

10

M=1 M=2 M=4 M=8 M=1, estimated M=2, estimated M=4, estimated M=8, estimated M=1, esti no water corr M=8, esti no water corr

-4

10

10-5 -6

10

10-7

1

2

3

4

5

β [nm-1]

Figure 7: The reciprocal error (points) of the IK-MSME and its error estimate (solid lines) as a function of the splitting parameter β investigated System 2. The meaning of the symbols is the same as Fig. 6. The error estimate without taking into account the correlation error is shown by the dashed lines for M = 1 and 8.

23



102


101

System 1

100 10-1 -2

10

M=1 M=2 M=4 M=8 M=1, estimated M=2, estimated M=4, estimated M=8, estimated

-3

10

10-4 -5

10

10-6

1

2

3

4

5

β [nm-1]

Figure 8: The reciprocal error (points) of the AD-MSME and the error estimate (lines) as a function of the splitting parameter β investigated in System 1. The number of mesh points is Kα = 32 and the order of B-spline interpolation is n = 4. The AD-MSME with M = 1, 2, 4 and 8 are plotted with color red, green, blue and pink, respectively.

10

1


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

System 2

10-1 10-2 -3

10

M=1 M=2 M=4 M=8 M=1, estimated M=2, estimated M=4, estimated M=8, estimated M=1, esti no water corr M=8, esti no water corr

-4

10

10-5 -6

10

10-7

1

2

3

4

5

β [nm-1]

Figure 9: The reciprocal error (points) of the AD-MSME and its error estimate (lines) as a function of the splitting parameter β investigated in System 2. The meaning of the symbols is the same as Fig. 8. The error estimate without taking into account the correlation error is shown by the dashed lines for M = 1 and 8.

24


Page 24 of 40

Page 25 of 40

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


error estimate is sharp for both ik and analytical differentiations in most of the β range. The agreement between the actual and estimated errors is expected, because the estimate catches all error contributions in a system that has a uniform and uncorrelated charge distribution. Deviation of the error estimate is observed when β > 5.0 nm−1 , and the error in this range is larger than 1 kJ/mol/nm. The deviation may stem from the truncated terms in the Ewald summation, which become increasingly significant with larger parameter β. 18 In this work, we do not consider this contribution, because the error estimate would become too complicated. The error estimates of the IK- and AD-MSME in System 2 are presented and compared with the actual error in Figs. 7 and 9, respectively. Both the error estimates of IK- and ADMSME are sharp, but they are less precise when β is smaller than 1.6 nm−1 . The maximum deviation at β = 0.6 nm−1 is 87%, which is acceptable for the applications like the parameter tuning. 16,18 The error estimates without counting the correlation error are also presented in the Figures by the dashed lines. For clarity, we only plot the M = 1 and 8 cases. It is noticed that the error of IK-MSME is not sensitive to the charge correlation in the system, while the error of the AD-MSME is very sensitive to the charge correlation: At relatively small β, the charge correlation reduces the error by more than one order of magnitude. This reminds us that, if the correlation error is important but difficult to estimate, simply using the error estimates without considering the charge correlation may be more reliable for the IK-MSME than for the AD-MSME. In Sec. 4, the Claims I and II were explained by the estimates of the first and second order homogeneity errors. We plot those of the IK- and AD-MSME for System 1 in Figs. 10 and 11, respectively. When M = 1 (MSME reduces to the SPME), the first order error dominates. This explains why only the first order error was estimated for the SPME. For the non-trivial MSME with M > 1, first order error dominates at the relatively small β, while the second order error dominates at the relatively large β. Therefore, we observe that the actual error firstly matches the first order error estimate at relatively small β, and then

25



10

1


100

System 1

-1

10

10-2 10-3 -4

10

M=1, 1st order err M=2, 1st order err M=4, 1st order err M=8, 1st order err 2nd order err 2nd order err, α ≠ β terms 2nd order err, α = β terms

-5

10

-6

10

-7

10

1

2

3

4

5

-1

β [nm ]

ik,(2)

ik,(1)

Figure 10: The estimates of the first (Ehomo ) and second order (Ehomo ) homogeneity errors of the IK-MSME in System 1. The points plot the reciprocal errors of IK-MSME with M = 1 (red “+”), 2 (green “×”), 4 (blue “⊡”) and 8 (pink “⊙”). The first order errors of M = 1, 2, 4 and 8 are presented by red, green, blue and pink lines, respectively. The second order error does not depends on M , therefore, only one line is plotted. The dashed cyan line plots the α 6= β contribution in the second order error. The dotted cyan line (overlapping with the solid line) plots the α = β contribution in the second order error.


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

10

2

10

1

10

0

System 1

-1

10

10-2 10-3

M=1, 1st order err M=2, 1st order err M=4, 1st order err M=8, 1st order err 2nd order err 2nd order err, α ≠ β terms 2nd order err, α = β terms

-4

10

10-5 -6

10

1

2

3

4

5

-1

β [nm ]

ad,(2)

ad,(1)

Figure 11: The estimates of the first (Ehomo ) and second order (Ehomo ) homogeneity errors of AD-MSME in System 1. The meaning of the symbols is the same as Fig. 10.

26


Page 26 of 40

Page 27 of 40

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


follows the second order error estimate at relatively large β. When the M increases, the first order error decreases, but the second order error does not, therefore, the range that the first order error dominates shrinks. The dashed and dotted cyan lines in Figs. 10 and 11 denote the α 6= β and α = β contributions (see Eqs. (29) and (41)) in the second order error, respectively. The α = β contribution overlaps with the second order error, while the α 6= β contribution is at least one order of magnitude smaller, therefore, it is demonstrated that the α = β contribution dominates the second order error.

6

The optimal working parameters for MSME in terms of efficiency

In the applications of the MSME, one needs to know how to tune the working parameters so that the time-to-solution is minimized while a desired accuracy is achieved. In this work, we improve the parameter tuning algorithm proposed by Ref., 18 which suggested to answer the question by solving the following optimization problem:

min T (β, rc , n, K , M ),

(45)

E(β, rc , n, K , M ) ≤ E ∗ ,

(46)

subject to the constraint of

where E ∗ is the target accuracy, E is the error estimate, and T is the performance model that estimates the time-to-solution as a function of the working parameters. The performance model T depends on the studied system, the hardware architecture and the software implementation, and can be obtained by fitting the model parameters to systematical timing tests using varying working parameters. It should be noticed that the variables β and rc are continuous, while n, K and M are discretized. Moreover, the components of K , i.e. Kα , 27



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 28 of 40

should be able to decomposed into small prime numbers to achieve the best performance in the FFTs. In this work, we require that Kα is even, and only prime factors 2, 3, 5 and 7 are allowed. We utilize the independency of the direct and reciprocal computations, and the monotony of the function T and E to simplify the problem (45)–(46). Firstly, for each fixed β0 , the direct and reciprocal computations are decoupled, so the original problem (45)–(46) can be approximately solved by two sub-problems min T dir (β0 , rc ), 1 s.t. E dir (β0 , rc ) ≤ √ E ∗ 2

(47) (48)

and min T rec (β0 , n, K , M ), 1 s.t. E rec (β0 , n, K , M ) ≤ √ E ∗ . 2

(49) (50)

Where T dir and T rec are the time-to-solution of the direct and reciprocal computations, respectively. E dir is the error estimate for the direct part, 19 while E rec is the reciprocal error estimate that is one of the main results of the current work. Taking the IK-MSME for ik ik ik ik example, it is given by E rec = (|Ehomo |2 + Ecorrection )1/2 , where Ehomo and Ecorrection are de-

fined by Eq. (27) and Eq. (33), respectively. Remembering that the splitting parameter is fixed to β0 , the direct time-to-solution T dir (β0 , rc ) increases monotonically w.r.t. rc , while the direct error E dir (β0 , rc ) decreases monotonically w.r.t. rc . The reciprocal time-to-solution T rec (β0 , n, K , M ) increases monotonically w.r.t. n, Kα and M , while the reciprocal error E rec (β0 , n, K , M ) decrease monotonically w.r.t. n, Kα and M . Therefore, the solutions to (47)–(48) and (49)–(50) must satisfy the equations in the constraints (48) and (50), respec√ tively. The direct space constraint E dir (β0 , rc ) = E ∗ / 2 can be efficiently solved by the bisection algorithm. In the reciprocal space, when n and M are fixed, the mesh size is also 28


Page 29 of 40

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


efficiently computed by the bisection algorithm. For all n and M pairs in a finite set, the K is solved, then the solution to (49)–(50) is the most time saving combination of {n, K , M }. The order of interpolation n was suggested to take less than 10. 18 The choice of M is upper bounded, because for an increasingly large M the reciprocal error is lower bounded. When the sub-problems (47)–(48) and (49)–(50) are solved, we set up a mapping from β to the total computational time T = T (β, rc (β), n(β), K (β), M (β)), which can be minimized by the 0.618 searching method. The resulting set of parameters is the solution to the original problem (45)–(46). We take the System 1 for instance and use the aforementioned process to tune the working parameters for the IK-MSME. The target accuracy varies from E ∗ = 1 kJ/(mol nm) to 10−6 kJ/(mol nm). For each target accuracy, the optimization problem (45)–(46) is solved with M fixed to {1, 2, 3, 4}, respectively. Here we do not take M as a variable because we intend to check the effect of using different numbers of meshes on the efficiency of the computation. The M value that achieves the shortest time-to-solution is the solution to the optimization problem (45)–(46), and should be used for simulating the System 1. The time model T = T dir + T rec is built for MOASP running on one core of an Intel Xeon E5-2692 v2 CPU. We take T dir (β, rc ) = C dir rc3 and T rec (β, n, K , M ) = C3intpl n3 + C1intpl n + C0intpl + C fft K1 K2 K3 log(K1 K2 K3 ). Model parameters C dir , C0intpl , C1intpl , C3intpl and C fft are fit to short profiling simulations taking different values of rc , n, K for each M . It should be noted that C0intpl , C1intpl and C3intpl do not increase linearly w.r.t. M , because the single instruction multiple data (SIMD) architecture of the CPU is used to compute the particlemesh interpolations on different meshes simultaneously. We call the third-party package P3DFFT 32 for FFTs, and find that the C fft almost increases linearly w.r.t. M in the single core situation, although the multivariable FFT interface is used. One of the difficulties in building a performance model for the parallel MSME is the parallel three dimensional FFT. It needs all-to-all data exchanges among the processors, therefore, the communication cost may contribute significantly to the time-to-solution be-

29



sides the floating point operations that scale as K1 K2 K3 log(K1 K2 K3 ). The modeling of the communication cost is not a trivial task, and depends on the computer architecture and the data partitioning. For example, Richards et. al. worked out the communication cost of the FFT for computers with multidimensional torus networks, such as the BlueGene/L computer. 33 Jung et. al. studied the communication overhead for different volumetric partitions, and determined the model parameters for the K computer. 34 The authors noticed that the communication cost is sensitive to the network latency, and can be effectively reduced by improving the data partitioning on a domestic computer. 35 Due to this complexity, the performance model for the parallel MSME could be an independent topic that goes beyond the scope of the current work. 0


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

10

System 1

-2

10-4 β = 4.39 nm 10

-6

10

-8

10-10

-1

n = 6, Kα = 40 M=1 M=2 M=3 M=4 1

2

3

4

5

β [nm-1]

Figure 12: The reciprocal error of IK-MSME as a function of the splitting parameter in the System 1. The order of interpolation is n = 6, and the number of mesh points is Kα = 40. M = 1, 2, 3 and 4 are presented in the Figure. The tuned working parameters for each M and the actual error and time-to-solution are listed in Table 1. Several observations are made from the Table. Firstly, the optimized parameters achieve the target accuracy. In several cases, the actual error is slightly (10% at most) above the required accuracy, because the error is slightly underestimated by the error estimate in these cases. Secondly, for each M , in general, the higher target accuracy requires a larger cut-off, a higher interpolation order and more mesh points. However, several exceptions are noticed: taking M = 3 for example, the optimal direct space cut-off 30


Page 30 of 40

Page 31 of 40

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Table 1: The tuned IK-MSME working parameters for System 1. The constrained optimization problem (45)–(46) are solved with M fixed to 1, 2, 3 and 4. The columns, from left to right, are the target accuracy (E ∗ ), the working parameters (M , β, rc , n and Kα ), the actual accuracy (E) and the computational time T . The unit of the accuracy, i.e. E ∗ or E, is kJ/(mol nm). The units of β and rc are nm−1 and nm, respectively. The computational time is recorded in the unit of millisecond. E∗ 100

10−1

10−2

10−3

10−4

10−5

10−6

M 1 2 3 4 1 2 3 4 1 2 3 4 1 2 3 4 1 2 3 4 1 2 3 4 1 2 3 4

β 4.02 4.97 4.68 4.39 5.11 4.71 4.18 4.20 3.99 5.15 3.68 3.48 3.68 4.32 4.39 4.39 4.04 3.87 4.09 3.91 3.41 4.59 4.10 4.10 3.30 4.39 4.75 3.65

rc 0.68 0.55 0.58 0.62 0.61 0.66 0.74 0.74 0.86 0.67 0.93 0.99 1.02 0.87 0.86 0.86 1.00 1.05 0.99 1.04 1.27 0.94 1.06 1.06 1.39 1.05 0.97 1.25

n 4 4 4 4 6 4 4 4 6 6 4 4 6 6 6 6 8 6 6 6 8 8 6 6 8 8 8 6

31

Kα 40 32 30 28 48 42 32 36 50 42 36 32 64 40 40 40 60 48 48 40 64 50 50 48 80 54 54 50


E 0.99 × 100 1.09 × 100 1.11 × 100 1.10 × 100 0.99 × 10−1 0.95 × 10−1 1.03 × 10−1 0.78 × 10−1 1.03 × 10−2 1.07 × 10−2 1.01 × 10−2 0.99 × 10−2 1.03 × 10−3 1.01 × 10−3 0.98 × 10−3 0.98 × 10−3 1.00 × 10−4 1.00 × 10−4 0.77 × 10−4 1.02 × 10−4 1.02 × 10−5 0.82 × 10−5 1.02 × 10−5 1.01 × 10−5 0.99 × 10−6 0.99 × 10−6 0.98 × 10−6 1.02 × 10−6

T 29.7 27.4 33.0 32.9 47.0 39.2 38.6 47.5 60.0 57.0 53.7 54.2 81.9 63.2 73.4 81.9 109.0 84.5 93.9 94.5 143.5 115.2 106.9 113.7 192.3 130.0 154.6 148.0


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

at E ∗ = 10−2 and E ∗ = 10−3 kJ/(mol nm) are 0.93 and 0.86 nm, respectively. The loss of accuracy in the direct space is compensated by a higher interpolation order (4 v.s. 6) and more mesh points (36 v.s. 40) in the reciprocal space accompanied with a larger splitting parameter (3.68 v.s. 4.89 nm−1 ). Thirdly, a clear efficiency improvement is observed when increasing the number of meshes M from 1 to 2, while further increasing M may not accelerate the simulation. In all the cases, we notice that the optimized splitting parameter β is relatively large. For these β, the M > 1 MSME reciprocal error may reach the lower bound, thus further increasing M does not reduce the error. In this situation, the extra cost of using more meshes leads to a decrement in the efficiency. The case of E ∗ = 10−3 kJ/(mol nm) is a perfect example to illustrate the situation. At this accuracy, the M = 2, 3 and 4 MSME use roughly the same splitting parameter (∼ 4.39 nm−1 ), interpolation order (n = 6) and number of mesh points (Kα = 40). We plot the MSME error against the splitting parameter β at n = 6 and Kα = 40 in Fig. 12. As indicated by the arrow in the Figure, the M = 2 MSME error reaches the lower bound at β = 4.39 nm−1 , so the M = 3 and 4 MSME do not further improve the accuracy and thus the efficiency decreases. It is interesting that when the target accuracy decreases to E ∗ = 10−5 kJ/(mol nm), and the M = 2 MSME has to use a higher interpolation order of n = 8, the M = 3 MSME that still uses n = 6 becomes more efficient. In general, more efficiency improvement is expected when the optimized splitting parameter is relatively small, because more accuracy improvement is obtained from using a larger number of staggered meshes. Comparing the MSME with different numbers of meshes, the highest computational efficiency is achieved at M = 2, 3, 3, 2, 2, 3 and 2 for E ∗ = 100 , 10−1 , 10−2 , 10−3 , 10−4 , 10−5 and 10−6 kJ/(mol nm), respectively. Therefore, there is no uniform answer on which M is the best. Both M = 2 and 3 can be optimal, depending on the required accuracy. Furthermore, the optimal M may be different for different applications, computers and software packages, because the performance models are different. The significance of generalizing the SPME to MSME is the extended parameter space with the extra degree of freedom M , in which a

32


Page 32 of 40

Page 33 of 40

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


better computational efficiency may be achieved by tuning the parameters. The advantage of the MSME over SPME in terms of efficiency

1

is 8.4%, 22%, 12%, 30%, 29%, 34%, and

48% for E ∗ = 100 , 10−1 , 10−2 , 10−3 , 10−4 , 10−5 and 10−6 kJ/(mol nm), respectively. There is a clear, but not monotonic, trend that the advantage is more significant for a higher target accuracy. Therefore, it is especially worth trying the MSME for the simulations that requires a high accuracy in the computation of electrostatic interactions.

7

Conclusion and discussion

In this work, the multiple staggered mesh Ewald (MSME) method is proposed to improve the accuracy of the smooth particle mesh Ewald (SPME) method with a modified function B. It takes the average of the SPME reciprocal forces computed on M staggered meshes, and achieves, in a certain range of the splitting parameter, the same accuracy as the SPME that uses M 3 times more mesh points. In the complementary parameter range, it is almost as accurate as the SPME that uses twice order of the B-spline interpolation. The reduction of the necessary FFT mesh points is particularly interesting for the massively parallel computation of the electrostatic interaction, because the three dimensional FFT is a well known bottleneck in the parallel implementation due to the intensive all-to-all data communications among the processors. The accuracy of the MSME is understood by a systematical error estimate. We prove that the multiple staggered meshes reduce the first order part of the error as much as refining the FFT mesh in the SPME, and does not change the second order part, which is roughly the same as doubling the interpolation order in the SPME. Therefore, the M = 2 MSME is always more accurate than the M = 1 MSME (i.e. SPME). Further increasing the number of meshes improves the accuracy, but the reciprocal error will eventually reach the lower bound when the first order error is reduced to lower than the second order error. The error 1

The efficiency is proportional to the inverse of the computational time. Taking M = 3 for example, the advantage of efficiency of is calculated by [1/T (M = 3) − 1/T (M = 1)]/[1/T (M = 1)].

33



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

estimate and the theoretical analysis on different orders of the error are validated both by a uniform and uncorrelated charge system and by a three-point-charge rigid water model. The error estimate developed in this work is significant not only because it explains the numerical phenomena of MSME, but also because it is of key importance in the parameter tuning. The optimal number of meshes M and other working parameters of MSME (i.e. the splitting parameter, the cut-off radius, the order of interpolation and the mesh size) are determined by the proposed parameter tuning algorithm, which minimizes the time-to-solution under the constraint of a prerequisite accuracy. As an example, we tune the parameters for the homogeneous and uncorrelated charge system using MOASP running on one core of an Intel Xeon E5-2692 v2 CPU. We show that the non-trivial MSME is always more efficient than SPME, and the optimal number of meshes is 2 or 3, depending on the required accuracy. A larger number of meshes does not improve the accuracy, because the lower bound of error is reached at the optimal splitting parameter. The optimal parameters are likely to change for different applications, different software packages running on different machines, nevertheless it is worth trying MSME, which may achieve a higher efficiency than the SPME by tuning the parameters in a larger space with M as an extra degree of freedom. Last but not least, the difference between the function B in our implementation of MSME and that in the original SPME should be noticed. We show that the MSME computation on one mesh is more accurate than the original SPME, and is almost as accurate as the PPPM method that is designed to minimize the reciprocal force error in homogeneous and uncorrelated charge systems.

Acknowledgment The authors gratefully acknowledge the financial support from National High Technology Research and Development Program of China under Grant 2015AA01A304. X.G. is supported by the National Science Foundation of China under Grants 91430218 and 61300012.

34


Page 34 of 40

Page 35 of 40

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


H.W. is supported by the National Science Foundation of China under Grants 11501039 and 91530322.

Supporting Information The error estimate for the multiple staggered mesh Ewald method. This information is available free of charge via the Internet at http://pubs.acs.org.

A

Comparison of the function B

The function B of the original SPME (20) is rewritten as

BO-SPME (m) =

Y α

1 ϕˆn (mα ) +

P

lα 6=0

ϕˆn (mα + lα Kα )

.

(51)

The function B of the PPPM method (21) is rewritten as

BPPPM (m) =

Y α

ϕˆn (mα ) +

P

1 . ˆ2n (mα + lα Kα )/ϕˆn (mα ) lα 6=0 ϕ

(52)

The difference between the definitions is thus the summation term in the denominator. We notice that the ϕˆn (mα ) decreases fast as mα increases, and −Kα /2 ≤ mα < Kα /2, then the term ϕˆn (mα +lα Kα ) converges to zero as |lα | increases. The difference between the function B used by the current work (11) and that of the original SPME is dominated by ϕˆn (mα ± Kα ), while the difference between the current work and the PPPM is dominated by ϕˆ2n (mα ± Kα )/ϕˆn (mα ). When mα is far from −Kα /2 or Kα /2, we have |ϕˆ2n (mα ± Kα )/ϕˆn (mα )| ≪ |ϕˆn (mα ± Kα )|. When mα is close to −Kα /2 or Kα /2, the difference in function B is suppressed by F or G (see the energy scheme (13) and the force schemes (16) and (18)), which decays exponentially fast for large m. Therefore, by using the current function B, our computation on a single mesh is closer to the PPPM method than to the original SPME

35



method. We numerically check this conclusion by investigating the reciprocal errors of the M = 1 MSME, the original SPME and the PPPM in Systems 1 and 2, and the results are plotted in Figs. 13 and 14, respectively. In both systems, the M = 1 MSME is in almost prefect agreement with the PPPM method, and is more accurate than the original SPME. 1


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

100 10

-1

10

-2

10

-3

10

-4

System 1

Kα = 32 Kα = 64 Kα = 128 Kα = 256 this work original SPME PPPM

10-5 10

-6

1

2

3

4

5

β [nm-1]

Figure 13: The reciprocal error of the M = 1 MSME (this work, crosses), the original SPME (solid lines) and the PPPM (dotted lines) as a function of the splitting parameter in the System 1. The force scheme is the ik-differentiation. The order of interpolation is n = 4, and the number of mesh points takes Kα = 32, 64, 128 and 256. The difference between the methods lies in the function B.

References (1) Ewald, P. P. Ann. Phys. 1921, 369, 253–287. (2) Perram, J.; Petersen, H.; De Leeuw, S. Mol. Phys. 1988, 65, 875–893. (3) Pollock, E.; Glosli, J. Comput. Phys. Commun. 1996, 95, 93–110. (4) Darden, T.; York, D.; Pedersen, L. J. Chem. Phys. 1993, 98, 10089. (5) Essmann, U.; Perera, L.; Berkowitz, M.; Darden, T.; Lee, H.; Pedersen, L. J. Chem. Phys. 1995, 103, 8577. (6) Hockney, R. W.; Eastwood, J. W. Computer Simulation Using Particles; IOP: London, 1988. 36


Page 36 of 40

Page 37 of 40


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


0

10 10

-1

10

-2

System 2

Kα = 32 Kα = 64 Kα = 128

10-3

Kα = 256

10-4 10

-5

10

-6

this work original SPME PPPM 1

2

3

4

5

-1

β [nm ]

Figure 14: The reciprocal error of the M = 1 MSME (this work, crosses), the original SPME (solid lines) and the PPPM (dotted lines) as a function of the splitting parameter in the System 2. The force scheme is the ik-differentiation. The order of interpolation is n = 4, and the number of mesh points takes Kα = 32, 64, 128 and 256. The difference between the methods lies in the function B. (7) Deserno, M.; Holm, C. J. Chem. Phys. 1998, 109, 7678. (8) Pippig, M. SIAM J. Sci. Comput. 2013, 35, C213–C236. (9) Alejandre, J.; Chapela, G. J. Chem. Phys. 2010, 132, 014701. (10) Isele-Holder, R. E.; Mitchell, W.; Ismail, A. E. J. Chem. Phys. 2012, 137, 174107– 174107. (11) in ’t Veld, P.; Ismail, A.; Grest, G. J. Chem. Phys. 2007, 127, 144711. (12) Isele-Holder, R. E.; Mitchell, W.; Hammond, J. R.; Kohlmeyer, A.; Ismail, A. E. J. Chem. Theory Comput. 2013, 9, 5412–5420. (13) Wennberg, C.; Murtola, T.; Hess, B.; Lindahl, E. J. Chem. Theory Comput. 2013, 9, 3527–3537. (14) Cerutti, D.; Duke, R.; Darden, T.; Lybrand, T. J. Chem. Theory Comput. 2009, 5, 2322–2338. (15) Neelov, A.; Holm, C. J. Chem. Phys. 2010, 132, 234103. 37



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(16) Wang, H.; Zhang, P.; Schütte, C. J. Chem. Theory Comput. 2012, 8, 3243–3256. (17) Abraham, M. J.; Gready, J. E. J. Comput. Chem. 2011, 32, 2031–2040. (18) Wang, H.; Dommert, F.; Holm, C. J. Chem. Phys. 2010, 133, 034117. (19) Kolafa, J.; Perram, J. Mol. Simul. 1992, 9, 351–368. (20) Hummer, G. Chem. Phys. Lett. 1995, 235, 297–302. (21) Petersen, H. J. Chem. Phys. 1995, 103, 3668. (22) Deserno, M.; Holm, C. J. Chem. Phys. 1998, 109, 7694. (23) Stern, H.; Calkins, K. J. Chem. Phys. 2008, 128, 214106. (24) Wang, H.; Schütte, C.; Zhang, P. Phys. Rev. E 2012, 86, 026704. (25) Frenkel, D.; Smit, B. Understanding molecular simulation; Academic Press, Inc. Orlando, FL, USA, 2001. (26) Ballenegger, V.; Cerdà, J.; Holm, C. J. Chem. Theory Comput. 2012, 8, 936–947. (27) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. J. Chem. Phys. 1983, 79, 926–935. (28) Gao, X.; Fang, J.; Wang, H. J. Chem. Phys. 2016, 144, 124113. (29) Berendsen, H. J.; Postma, J. P.; van Gunsteren, W. F.; Hermans, J. Intermolecular forces; Reidel: Dordrecht, 1981; pp 331–342. (30) Berendsen, H.; Grigera, J.; Straatsma, T. J. Phys. Chem. 1987, 91, 6269–6271. (31) Ballenegger, V.; Cerdà, J.; Holm, C. Comput. Phys. Commun. 2011, 182, 1919–1923. (32) Pekurovsky, D. SIAM J. Sci. Comput. 2012, 34, C192–C209.

38


Page 38 of 40

Page 39 of 40

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(33) Richards, D. F.; Glosli, J. N.; Chan, B.; Dorr, M.; Draeger, E. W.; Fattebert, J.-L.; Krauss, W. D.; Spelce, T.; Streitz, F. H.; Surh, M. Beyond homogeneous decomposition: scaling long-range forces on Massively Parallel Systems. In Proceedings of the Conference on High Performance Computing Networking, Storage and Analysis. Portland, Oregon, November 14-20, 2009; ACM: New York, NY, USA, 2009. (34) Jung, J.; Kobayashi, C.; Imamura, T.; Sugita, Y. Comput. Phys. Commun. 2016, 200, 57–65. (35) Gao, X.; Mo, Z.; Fang, J.; Song, H.; Wang, H. Comput. Phys. Commun. 2016,

39



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Graphical TOC Entry

2/3 1/3

40


Page 40 of 40

Multiple Staggered Mesh Ewald: Boosting the Accuracy of the Smooth

Recommend Documents