The optimal P3M algorithm for computing electrostatic energies in periodic systems
TThe optimal P3M algorithm for computing electrostatic energiesin periodic systems
V. Ballenegger
Institut UTINAM, Universit´e de Franche-Comt´e, UMR 6213,16, route de Gray, 25030 Besan¸con cedex
France . J. J. Cerda and O. Lenz
Frankfurt Inst. for Advanced Studies,J.W. Goethe - Universit¨at, Frankfurt, Germany
Ch. Holm
Frankfurt Inst. for Advanced Studies,J.W. Goethe - Universit¨at, Frankfurt, Germany andMax-Planck-Institut f¨ur Polymerforschung, Mainz, Germany (Dated: October 26, 2018)
Abstract
We optimize Hockney and Eastwood’s Particle-Particle Particle-Mesh (P3M) algorithm toachieve maximal accuracy in the electrostatic energies (instead of forces) in 3D periodic chargedsystems. To this end we construct an optimal influence function that minimizes the RMS errorsin the energies. As a by-product we derive a new real-space cut-off correction term, give a trans-parent derivation of the systematic errors in terms of Madelung energies, and provide an accurateanalytical estimate for the RMS error of the energies. This error estimate is a useful indicatorof the accuracy of the computed energies, and allows an easy and precise determination of theoptimal values of the various parameters in the algorithm (Ewald splitting parameter, mesh sizeand charge assignment order). a r X i v : . [ phy s i c s . c o m p - ph ] A ug . INTRODUCTION Long range interactions are ubiquitously present in our daily life. The calculation ofthese interactions is, however, not an easy task to perform. One needs indeed to resort tospecialized algorithms to overcome the quadratic scaling with the number of particles, assoon as the simulated system includes more than a few hundred particles, see for example thereview of Arnold and Holm . In Molecular Dynamics simulations, one is mainly interestedin the accuracy of the force computation, since they govern the dynamics of the system. Incontrast, in Monte Carlo (MC) simulations, the concern is to compute accurate energies.If the potential is of long range ( e.g. a Coulomb potential or dipolar interaction), andone has chosen to use periodic boundary conditions, the computation of both observablesis quite time consuming if one uses the traditional Ewald sum. Since the seminal workof Hockney and Eastwood it has been common to resort to a faster way of calculatingthe reciprocal space sum in the Ewald method with the help of Fast-Fourier-Transforms(FFTs). These algorithms are called mesh-based Ewald sums, and various variants exist .They all scale as N log N with the number of charged particles N , and the algorithms arenowadays routinely used in simulations of bio-systems, charged soft matter, plasmas, andmany more areas. The most accurate variant is still the original method of Hockney andEastwood, which they called particle-particle-particle-mesh (P3M), and into which variousother improvements like the analytical differentiation used in other variants of the mesh-based Ewald sum can be built in. In addition, an accurate error estimate for P3M exists,so that one can tune the algorithm to a preset accuracy, thus maximizing the computationalefficiency before doing any simulations .While in the standard P3M algorithm , the lattice Green function, called the “influencefunction”, is optimized to give the best possible accuracy in the forces, the electrostaticenergy is usually calculated with the same force-optimized influence function. However,there are certainly situations where one needs a high precision of the energies, for instancein Monte Carlo simulations, and the natural question arises whether one can optimize theinfluence function to enhance the accuracy of the P3M energies. The main goal of this paperis to derive the energy-optimized influence function, and to derive an analytical estimate forthe error in the P3M energies. This error estimate is a valuable indicator of the accuracyof the calculations and allows a straightforward and precise determination of the optimal2alues of the various parameters in the algorithm (Ewald splitting parameter, mesh size,charge assignment order).The present derivation of the optimal influence function, and the associated error esti-mate, is concise and entirely self-contained. The present paper can thus also serve as apedagogical introduction to the main ideas and mathematics of the P3M algorithm.The paper is organized as follows. In Sec. II, we briefly review the ideas of the standardEwald method and provide the most important formulae. In Sec. III, we derive directand reciprocal space correction terms which compensate, on average, the effects of cut-offerrors in the standard Ewald method. We interpret the formulae in terms of the direct andreciprocal space components of the Madelung energies of the ions. In Sec. IV, the calculationof the reciprocal energy according to the P3M algorithm ( i.e. with a fast Fourier transformand an optimized influence function) is presented. The mathematical analysis of the errorsintroduced by the discretization on a grid is performed in Sec. V. This analysis is used inSec. VI to derive the energy-optimized influence function and the associated RMS errorestimate. The derivation shows that the P3M energies must be shifted to compensate forsystematic cut-off and aliasing errors in the Madelung energies of the ions. Finally, ouranalytical results are tested numerically in Sec. VII. II. THE EWALD SUM
We consider a system of N particles with charges q i at positions r i in an overall neutraland (for simplicity) cubic simulation box of length L and volume V = L . If periodicboundary conditions are applied, the total electrostatic energy of the box is given by E = 12 (cid:88) n ∈ Z N (cid:88) (cid:48) i,j =1 q i q j v ( r ij + n L ) , (2.1)where v ( r ) = 1 / | r | is the Coulomb potential, r ij = r i − r j , and n is a vector with integercomponents that indexes the periodic images. The prime indicates that the (divergent)summand for i = j has to be omitted when n = 0.Because of the slow decay of the Coulomb interaction, the sum in (2.1) is only condition-ally convergent: its value is not well defined unless one specifies the precise way in whichthe cluster of simulation boxes is supposed to fill R . Often, one chooses a spherical orderof summation, which is equivalent to the limit of a large, spherically bounded, regular grid3f replicas of the simulation box, embedded in vacuum. The simulation box can then bepictured as the central LEGO brick in a huge ball made up of such bricks. If this “lego ball”is surrounded by a homogeneous medium with dielectic constant (cid:15) (cid:48) ( (cid:15) (cid:48) = 1 if it’s vacuum)and if the simulation box has a net dipole moment M = (cid:80) i q i r i , the particles in the ballwill feel a depolarizing field created by charges that appear on the surface of the uniformilypolarized ball. It can be shown that the work done against this depolarizing field whencharging up the system is E ( d ) = 2 π M (1 + 2 (cid:15) (cid:48) ) L (2.2)in the case of a spherical order of summation (for other summation orders, see the articlesof Smith and Ballenegger and Hansen ). The energy E ( d ) is contained, even if not easilyseen, in the total electrostatic energy (2.1) (at least when (cid:15) (cid:48) = 1 since such a vacuumboundary condition was assumed in writing (2.1)). Obviously, the energy E ( d ) vanishes ifwe employ metallic boundary conditions defined by (cid:15) (cid:48) = ∞ .The fact that E ( d ) depends on the order of summation, and hence on the shape of macro-scopic sample under consideration, is a consequence of the conditional convergence of thesum (2.1). Due to the energy cost E ( d ) , the fluctuations of the total dipole moment of thesimulation box (and hence of the considered macroscopic sample) depend on the dielectricconstant (cid:15) (cid:48) and on the shape of the sample. The energy E ( d ) is crucial to ensure, for ex-ample, that the dielectric constant (cid:15) of the simulated system obtained from the Kirkwoodformula , which relates (cid:15) to the fluctuations of the total dipole moment, is independent ofthe choices made for the sample shape and for the dielectric boundary condition .Ewald’s method to compute the energy (2.1) is based on a decomposition of the Coulombpotential, v ( r ) = ψ ( r ) + φ ( r ), such that ψ ( r ) contains the short-distance behavior of theinteraction, while φ ( r ) contains the long-distance part of the interaction and is regular atthe origin. The traditional way to perform this splitting is to define φ ( r ) = erf( αr ) /r, r = | r | , (2.3)and ψ ( r ) = v ( r ) − φ ( r ) = erfc( αr ) /r (2.4)With this choice, ψ ( r ) corresponds to the interaction energy between a unit charge at adistance r from another unit charge that is screened by a neutralizing Gaussian charge dis-4ribution whose width is controlled by the Ewald length α − . Following this decompositionof the potential, the electrostatic energy can be written in the well-known Ewald form : E = E ( r ) + E ( k ) + E ( d ) (2.5)where the real-space energy E ( r ) contains the contributions from short-range interac-tions ψ ( r ), i.e. E ( r ) = 12 (cid:88) n ∈ Z N (cid:88) (cid:48) i,j =1 q i q j ψ ( r ij + n L ) , (2.6)and the reciprocal space energy E ( k ) contains contributions from long-range interactions φ ( r )(apart from the contributions that are responsible for the conditional convergence which areincluded in the term E ( d ) in (2.5)). The fact that the surface term (or “dipole term”) E ( d ) is independent of the Ewald parameter α shows that this contribution is not specific tothe Ewald method, but more generally reflects the problems inherent to the conditionalconvergence of the n sum in Eq. (2.1). Contrary to E ( r ) which can be computed easily inreal space thanks to the rapid decay of the ψ interaction, E ( k ) is best computed in Fourierspace, where it can be expressed as E ( k ) = E ( ks ) − E ( s ) (2.7)where E ( ks ) = 12 L (cid:88) k ∈ K k (cid:54) =0 | (cid:101) ρ ( k ) | (cid:101) φ ( k ) (2.8) E ( s ) = Q α √ π (2.9)with Q = N (cid:88) i =1 q i . (2.10)In (2.8), (cid:101) φ ( k ) is the Fourier transform of the reciprocal interaction (2.3), (cid:101) φ ( k ) = (cid:90) e − i k · r φ ( r )d r = 4 πk exp( − k / α ) , (2.11)and (cid:101) ρ ( k ) is the Fourier transformed charge density (cid:101) ρ ( k ) = N (cid:88) i =1 q i e − i k · r i . (2.12)5he sum in (2.8) is over wave vectors in the discrete set K = { π n /L : n ∈ Z } . The term k = 0 is excluded in the sum because of the overall charge neutrality. The self-energy term E ( s ) compensates for the self-energies (the reciprocal interaction of each particle with itself q i φ ( r = ) = q i α/ √ π ) that are included in E ( ks ) .The energy (2.1) converges only for systems that are globally neutral. For systems witha net charge, the sum can be made convergent by adding a homogeneously distributedbackground charge which restores neutrality. In that case, an additional contribution E ( n ) = − π α L (cid:16) N (cid:88) i =1 q i (cid:17) (2.13)must be added to (2.5) to account for the interaction energies of the charges with theneutralizing background.The reciprocal energy E ( ks ) , defined by the Ewald formula (2.7), is the starting point ofmesh-based Ewald sums, which are methods to compute efficiently that energy in many-particle systems. Notice that (2.7) can also be written in an alternative form in terms of apair potential and a Madelung self-energy, see Appendix A. The inverse length α tunes therelative weight of the real space E ( r ) and the reciprocal space E ( k ) contributions to the energy,but the final result is independent of α . In practice, E ( r ) and E ( k ) can be computed usingcut-offs, because the sum over n in (2.6) and the sum over k in (2.8) converge exponentiallyfast. Typically, one chooses α large enough to employ the minimum image convention inEq. (2.6).At given real and reciprocal space cut-offs r cut and k cut , there exists actually an optimal α such that the accuracy of the approximated Ewald sum is as high as possible. This optimalvalue can be determined with the help of the estimates for the cut-off errors derived byKolafa and Perram , by demanding that the real and reciprocal space contributions to theerror are equal. Kolafa and Perram’s root-mean-square error estimates are∆ E ( r ) (cid:39) Q (cid:114) r c L e − α r ( αr cut ) (2.14)and ∆ E ( k ) (cid:39) Q α e − ( πk cut /αL ) π k / . (2.15)These error estimates make explicit the exponential dependence of the error on the real andreciprocal space cut-offs. 6ormula (2.15) is actually valid only when a correction term (given by Eq. (3.8) below)is added, to compensate the systematic error that affects the reciprocal energies when thesum over wave vectors in (2.8) is truncated. The origin of this correction term is explainedin detail in Sec. III. A similar term must also be introduced in the P3M algorithm whenone computes the electrostatic energy. Similarly, the direct-space energy (2.6) also containsa systematic error when the pair-wise interaction is truncated at the cut-off distance r cut .The derivation in the next section will also provide a correction term for this effect.To summarize, the final Ewald formula for the total electrostatic energy reads E = E ( r ) [eq . (2.6)]+ E ( ks ) [eq . (2.8)] − E ( s ) [eq . (2.9)]+ E ( d ) [eq . (2.2)]+ E ( n ) [eq . (2.13)] (2.16)Furthermore, when the sums in E ( r ) and E ( k ) are evaluated numerically using cut-offs, anadditional correction term E cut , defined in Eq. (3.13) below, must be added to the truncatedenergy, as shown in the next section. III. CORRECTION TERM FOR TRUNCATED EWALD SUMS
If we consider electroneutral systems where the charged particles are located at random,we expect the electrostatic energy to vanish on average, because there is an equal probabilityto find a positive or negative charge at any relative distance r . However, when periodicboundary conditions (PBC) are applied, the average energy of random systems does notvanish, because each charge interacts with its own periodic images (and with the uniformneutralizing background provided by the other charges).Since this interaction energy E img of an ion with its periodic images and with the neu-tralizing background does not depend on the position of the ion in the simulation box, itplays the role of a “self-energy”. We will refer to E img as the Madelung (self-)energy of anion, to avoid confusion with the self-energy q φ (0) already defined in the Ewald method asthe reciprocal interaction of a particle with itself.We denote by angular brackets the average over the positions of the N charged particles:7 · · ·(cid:105) = 1 V N (cid:90) V N · · · d r ... d r N . (3.1) A. Madelung energy
The Madelung energy of an ion takes the form E img = q ζ , where ζ is a purely numericalfactor in units of 1 /L that depends only on the size and shape of the simulation box.Let us calculate the average electrostatic energy of random charged systems in PBC,to find the value of ζ and derive a correction term for cut-off errors in truncated Ewaldsums (some results derived here will be used in Sec. VI A). On the one hand, the averageCoulomb energy of the random systems is by definition Q ζ/
2, while on the other hand, itcan be calculated as the sum of a direct space contribution (cid:10) E ( r ) (cid:11) and a reciprocal spacecontribution (cid:10) E ( k ) (cid:11) . The average reciprocal energy is, using (2.9) and (2.8), (cid:10) E ( k ) (cid:11) = 12 L (cid:88) i,j q i q j (cid:88) k ∈ K k (cid:54) =0 (cid:10) e − i k · ( r i − r j ) (cid:11) (cid:101) φ ( k ) − Q α √ π . (3.2)Since (cid:104) exp( − i k · r j ) (cid:105) = δ k , , all terms with j (cid:54) = i vanish (this is due to the fact that theEwald pair potential averages to zero, see Appendix A). By contrast, “self” terms ( i = j )remain and lead to (cid:10) E ( k ) (cid:11) = Q (cid:16) L (cid:88) k ∈ K k (cid:54) =0 (cid:101) φ ( k ) − α √ π (cid:17) = Q ζ ( k ) , (3.3)where the second equality defines ζ ( k ) . The average real-space energy of a single ion ofcharge q i in periodic random systems is (cid:68) E ( r ) i (cid:69) = q i (cid:16) (cid:88) n ∈ Z n (cid:54) =0 ψ ( n L ) − L (cid:90) R ψ ( r ) d r (cid:17) (3.4)where the first term is the sum of the direct interactions of the ion with all its periodicimages, while the second term corresponds to its interaction with the uniform backgroundcharge density − q i /L provided by the other particles in the system. Since (cid:90) R ψ ( r ) d r = 4 π (cid:90) ∞ r ψ ( αr ) d r = πα , (3.5)8e can write the average total real-space energy as (cid:10) E ( r ) (cid:11) = Q (cid:16) (cid:88) n ∈ Z n (cid:54) =0 ψ ( n L ) − πα L (cid:17) = Q ζ ( r ) , (3.6)which defines ζ ( r ) . The second term in ζ ( r ) is, not surprisingly, identical to the energy E ( n ) defined in (2.13). Notice that the above result for (cid:10) E ( r ) (cid:11) may also be obtained by splitting(2.6) into self ( i = j ) and interaction terms, and using for the latter (cid:80) j (cid:54) = i q j = − q i whichfollows from the electro-neutrality condition. The expression of the factor ζ = ζ ( r ) + ζ ( k ) istherefore ζ = (cid:16) (cid:88) n ∈ Z n (cid:54) =0 ψ ( n L ) − πα L (cid:17) + (cid:16) L (cid:88) k ∈ K k (cid:54) =0 (cid:101) φ ( k ) − α √ π (cid:17) . (3.7)Eq. (3.7) can be computed for a number of different box geometries . For a cubicsimulation box of size L , it yields ζ (cid:39) − . /L. The above calculation shows that, when a charged system is simulated using PBC, theelectrostatic energy (2.1) contains an additional contribution Q ζ/
2. The existence of thisMadelung self-energy can be made more apparent in the Ewald formula for E , as shown inAppendix A. B. Madelung cut-off error correction terms
The Ewald sums (2.6) and (2.8) are necessarily truncated when evaluated in a simula-tion. These truncations introduce systematic cut-off errors in the total energy, because theMadelung self-energies of the ions are then not fully accounted for. This systematic error istypically of the same order of magnitude, or even larger, than the fluctuating error, due tothe use of cut-offs, in the Ewald pair interaction energy . Note, that no similar systematicerror affects the electrostatic forces, because the Madelung energy does not depend on theposition of the ion.Fortunately, it is easy to suppress the systematic bias in the computed energies. Wesimply have to add the cut-off correction E ( k )cut = Q L (cid:88) k ∈ K k>k cut (cid:101) φ ( k ) (3.8)9o the computed k -space energies, which Kolafa and Perram termed the diagonal correc-tion . The value of E ( k )cut does not depend on the configuration and may thus be computedin advance using a sufficiently large second cut-off k (cid:48) cut > k > k cut . Using definition (3.3)of ζ ( k ) , we can rewrite (3.8) as E ( k )cut = Q (cid:16) ζ ( k ) − ζ ( k )cut (cid:17) (3.9)where ζ ( k )cut = 1 L (cid:88) k ∈ K k (cid:54) =0 , k The idea of particle-mesh algorithms is to speed up the calculation of the reciprocalenergy E ( ks ) with the help of a Fast-Fourier-Transform (FFT). To use a FFT, the chargedensity must be assigned to points on a regular grid. There are several ways of discretizingthe charge density on a grid, and to get the electrostatic energy from the Fourier transformedgrid. We will use the P3M method of Hockney and Eastwood (but with the standard Ewaldreciprocal interaction (2.3)), because this method surpasses in efficiency the other variantsof mesh based Ewald sums (PME, SPME) .For simplicity, we assume the number of grid points M to be identical in all three direc-tions. Let h = L/M be the spacing between two adjacent grid points. We denote by M theset of all grid points: M = { m h : m ∈ Z and 0 ≤ m x,y,z < M } .The mesh based calculation of the reciprocal energy is made in the following steps: A. Assign charges to grid points The charge density ρ M ( r ) at a grid point r is computed via the equation ρ M ( r ) = (cid:90) U ( r − r (cid:48) ) ρ ( r (cid:48) )d r (cid:48) , r ∈ M , (4.1)where U ( r ) = h − W ( r ) with W the charge assignment function (the factor h − ensuresmerely that ρ M ( r ) has the dimensions of a density). A charge assignement function is clas-sified according to its order P , i.e. between how many grid points per coordinate directioneach charge is distributed. Typically, one chooses a cardinal B-spline for W , which is apiece-wise polynomial function of weight one. The order P gives the number of sections inthe function. In P3M, we only need the Fourier transform of the cardinal B-splines, whichare (cid:102) W ( P ) ( k ) = h (cid:18) sin( k x h/ k x h/ k y h/ k y h/ k z h/ k z h/ (cid:19) P . (4.2)Notice that ρ M ( r ) = h − (cid:80) i q i W ( r − r i ), apart at the boundaries where the periodicityhas to be properly taken into account. 11 . Fourier transform the charge grid Compute the finite Fourier transform of the mesh-based charge density (using the FFTalgorithm) (cid:101) ρ M ( k ) = h (cid:88) r ∈ M ρ M ( r ) e − i k · r = FFT { ρ M } , k ∈ (cid:101) M . (4.3)Here k is a wave vector in the reciprocal mesh (cid:101) M = { π n /L : n ∈ Z , | n x,y,z | < M/ } .We stress that (cid:101) ρ M ( k ) differs from (cid:101) ρ ( k ) for k ∈ (cid:101) M , because sampling of the charge densityon a grid introduces errors (see Sec. V). C. Solve Poisson equation (in Fourier space) The mesh-based electrostatic potential Φ M is given by the Poisson equation, which re-duces to a simple multiplication in k-space: (cid:101) Φ M ( k ) = (cid:101) ρ M ( k ) (cid:101) φ ( k ) , k ∈ (cid:101) M , (4.4)with (cid:101) φ ( k ) the Fourier transformed reciprocal interaction (2.11). However, instead of using (cid:101) φ ( k ) in the above equation, it is better to introduce an “influence” function (cid:101) G ( k ). Wereplace therefore Eq. (4.4) by (cid:101) Φ M ( k ) = (cid:101) ρ M ( k ) (cid:101) G ( k ) , k ∈ (cid:101) M . (4.5)where (cid:101) G ( k ) is determined by the condition that it leads to the smallest possible errors inthe computed energies (on average for uncorrelated random charge distributions). (cid:101) G ( k ) willbe determined later (see Eq. (6.21)); it can be computed once and for all at the beginningof a simulation since it depends only on the mesh size and the charge assignment function. (cid:101) G ( k ) plays basically the same role as the reciprocal interaction (cid:101) φ ( k ), except that it is tunedto minimize a well defined error functional in (cid:101) ρ M ( k ). We stress that (cid:101) G ( k ) is defined onlyfor k ∈ (cid:101) M (we dropped the subscript M on the influence function to alleviate the notation).The idea of optimizing (cid:101) G ( k ), which is a key-point of the P3M algorithm, ensures that themesh based calculation of the reciprocal energy gives the best possible results . Get total reciprocal electrostatic energy Expression (2.8) is approximated on the mesh by E ( ks )P3M = 12 L (cid:88) k ∈ e M k (cid:54) =0 | (cid:101) ρ M ( k ) | (cid:101) G ( k ) . (4.6)The total reciprocal energy follows from subtracting the self-energies from the above quan-tity: E ( k )P3M = E ( ks )P3M − E ( s ) . E. Electrostatic energy of individual charges (optional) If the reciprocal energy of each individual particle is needed (and not only their sum asin step 4), the potential mesh must be transformed back to real space via an inverse FFT, i.e. Φ M ( r m ) = 1 L (cid:88) k ∈ e M (cid:101) φ M ( k ) e i k · r m = FFT − { (cid:101) φ M } . (4.7)The mesh-based potential is then mapped back to the particle positions using the samecharge assignment function:Φ( r ) := (cid:88) r m ∈ M p W ( r − r m )Φ M ( r m ) . (4.8)In this equation, M p = { m h : m ∈ Z } is the mesh extended by periodicity to all space,and Φ M ( r ) is assumed to be periodic (with period L ). The interpretation of Eq. (4.8) isthe following: due to the discretization each particle is replaced by several “sub-particles”which are located at the surrounding mesh points and carry the fraction W ( r − r m ) of thecharge of the original particle. The potential at the position of the original particle is givenby the sum of the charge fraction times the potential at each mesh points. The reciprocalelectrostatic energy of the i th particle is then q i Φ( r i ) / 2, and the total reciprocal energy(including self-energies) is the sum E ( ks )P3M = 12 (cid:90) V ρ ( r )Φ( r )d r = 12 (cid:88) i q i Φ( r i ) . (4.9)This formula gives the same result for the total energy as Eq. (4.6). A mathematical proofof the equivalence is given in Appendix C. 13 . ANALYSIS OF DISCRETIZATION ERRORS If the fast Fourier transform has the benefit of speed, it has the drawback of introducingerrors in the k -space spectrum of the charge density: (cid:101) ρ M ( k ) differs from the true Fouriertransform (2.12)(times a trivial factor (cid:101) U ( k )) because of the discretization on a finite grid.The difference is two-fold. Firstly, (cid:101) ρ ( k ) is defined for any vector in the full k-space K ,whereas (cid:101) ρ M ( k ) is defined only for k ∈ (cid:101) M , i.e. in the first Brillouin zone. This is a firstnatural consequence of discretization: if the grid spacing is h , it necessarily introduces acut-off | k x,y,z | < π/h in k -space. Secondly, the act of sampling the charge density at gridpoints, which is mathematically embodied in Eq. (4.3) by the presence of a discrete Fouriertransform instead of a continuous FT, introduces aliasing errors. While a continuous FTwould simply transform the convolution Eq. (4.1) intoFT { ρ M } ( k ) = (cid:101) U ( k ) (cid:101) ρ ( k ) , k ∈ K , (5.1)the finite Fourier transform results in (see proof in Appendix B) (cid:101) ρ M ( k ) = FFT { ρ M } ( k ) = (cid:88) m ∈ Z (cid:101) U ( k + m k g ) (cid:101) ρ ( k + m k g ) , k ∈ (cid:101) M , (5.2)where k g = 2 π/h . The sum over m shows that spurious contributions from high frequenciesof the full spectrum (cid:101) U ( k ) (cid:101) ρ ( k ) are introduced into the first Brillouin zone (cid:101) M . These unwantedcopies of the other Brillouin zones into the first one are known as aliasing errors .To avoid aliasing errors, the spectrum needs to be entirely contained within the firstBrillouin zone. Since (cid:101) ρ ( k ) may contain arbitrary high frequencies, this can only be achievedby choosing U ( k ) to be a low-pass filter satisfying (cid:101) U ( k ) = 0 for k ∈ K \ (cid:101) M . But thecharge assignment function would then have a compact support in k -space, and hence aninfinite support in r -space. This is not acceptable, as it would require the grid to have aninfinite extension. The need to keep the charge assignment function local in r -space meansthat (cid:101) U ( k ) cannot be a perfect low pass filter. Aliasing errors are therefore unavoidable,and the impact of these errors must be minimized, by choosing a good compromise for thecharge assignment function and optimizing the influence function. The influence functioncan indeed compensate partially for the aliasing errors, because the spectrum of (cid:101) U ( k ) isknown exactly at all frequencies.The error in reciprocal energy, for a given configuration ρ ( r ) of the charges, is defined by14he difference∆ E ( k ) = E ( k )P3M − E ( k ) (5.3)where E ( k ) is the exact reciprocal energy (see (2.8) and (2.9)). The above analysis of dis-cretization errors results in the explicit formula for this error∆ E ( k ) = 12 L (cid:88) k ∈ e M k (cid:54) =0 | (cid:101) ρ M ( k ) | (cid:101) G ( k ) − L (cid:88) k ∈ K k (cid:54) =0 | (cid:101) ρ ( k ) | (cid:101) φ ( k ) , (5.4)where (cid:101) ρ M ( k ) is given by (5.2). The error ∆ E ( k ) is due to the finite resolution h offered by themesh. The finiteness of h introduces the cut-off π/h in k -space ( k ∈ (cid:101) M ) and causes aliasingerrors ( (cid:101) ρ M ( k ) (cid:54) = (cid:101) ρ ( k ) (cid:101) U ( k )) that cannot be entirely eliminated by the charge assignmentfunction. VI. OPTIMIZATION OF THE P3M ALGORITHM We derive in this section the influence function (cid:101) G ( k ) that minimizes the error (5.4) onaverage for uncorrelated systems, and give a formula for the associated RMS errors. Theaverage over random systems is denoted by angular brackets, as in Sec. III.Notice that the assumption of the absence of correlations is never satisfied in practice(even for uniform systems because negative charges tend to cluster around positive chargesand vice-versa). The error estimate proves however to predict quite accurately the error inreal systems with correlations, notably in liquids where the pair distribution function g ( r )decays rapidly to one. A. Shift in the energies to avoid systematic errors The P3M energies (4.6) contain in general systematic errors, i.e. (cid:10) ∆ E ( k ) (cid:11) (cid:54) = 0, becausethe Madelung energies of the ions obtained in the mesh calculation contain cut-off andaliasing errors. The average error K = (cid:10) ∆ E ( k ) (cid:11) = (cid:68) E ( k )P3M (cid:69) − (cid:10) E ( k ) (cid:11) (6.1)is a constant that must be subtracted from the P3M energies, to ensure that the energiesare right on average. The corrected P3M energies are thus obtained by applying a constant15hift to the original P3M energies: E ( k )P3M , corr = E ( k )P3M − K, (6.2)where the constant K depends on the various P3M parameters like mesh size, charge as-signment order (CAO) and Ewald splitting parameter.Let us determine analytically the constant (6.1). Writing it as K = (cid:68) E ( k )P3M (cid:69) − (cid:10) E ( k ) (cid:11) ,we can use the result (3.3) for (cid:10) E ( k ) (cid:11) : it is nothing but Q ζ ( k ) / i.e. the k -space Madelungenergies of the ions. The other term (cid:68) E ( k )P3M (cid:69) can be calculated in the same way as (3.2).Using (4.6), (5.2) and (2.9), we find (cid:68) E ( k )P3M (cid:69) = Q (cid:16) L (cid:88) k ∈ e M k (cid:54) =0 (cid:101) G ( k ) (cid:88) m ∈ Z (cid:101) U ( k + k g m ) − α √ π (cid:17) = Q ζ ( k )P3M , (6.3)which defines ζ ( k )P3M . The result (6.3) can be interpreted as the average k -space Madelungenergies of the ions as obtained from the mesh calculation , i.e. including cut-off and aliasingerrors. The explicit expression of the correction constant (6.1) is thus K = Q (cid:16) ζ ( k )P3M − ζ ( k ) (cid:17) = Q L (cid:16) (cid:88) k ∈ e M k (cid:54) =0 (cid:101) G ( k ) (cid:88) m ∈ Z (cid:101) U ( k + k g m ) − (cid:88) k ∈ K k (cid:54) =0 (cid:101) φ ( k ) (cid:17) (6.4)In the last sum in (6.4), the terms with | k x,y,z | > π/h are equivalent to the k -space cut-off correction defined in (3.8). These terms compensate for the fact that the Madelungenergies of the ions are underestimated in the mesh calculation because of the cut-off π/h introduced by the finite size of the mesh. The remaining terms in (6.4) compensate, onaverage, the aliasing errors that affect the Madelung energies of the ions obtained from themesh calculation.Notice that the two correction terms E ( r )cut and − K can be combined together in the simpleexpression E cutP3M = E ( r )cut − K = Q (cid:16) ζ − ζ ( k )P3M − ζ ( r )cut (cid:17) (6.5)where ζ is defined by (3.7) and ζ ( r )cut is given in (3.12). We stress that E cutP3M has the samestructure as the correction term (3.13) for truncated Ewald sums. The difference lies in thereplacement of ζ ( k )cut by the quantity ζ ( k )P3M defined in (6.3), which accounts for both the cut-offand aliasing errors that affect the reciprocal energies computed on the mesh.16n summary, the final formula for computing the total electrostatic energy with the P3Malgorithm is E ≈ E ( r ) [eq . (2.6)]+ E ( ks )P3M [eq . (4.6)] − E ( s ) [eq . (2.9)]+ E ( d ) [eq . (2.2)]+ E ( n ) [eq . (2.13)]+ E cutP3M [eq . (6.5)] (6.6)The correction term E cutP3M is necessary to compensate on average for systematic errors in themesh calculation. It can be computed once for all before the start of a simulation, since itdepends only on the size of the simulation box, the size of a mesh cell, the charge assignmentfunction and the influence function. B. RMS error estimate for energy The result (5.4) is an exact measure of the error in the P3M energies for a given configu-ration ρ ( r ) of the particles. Let us average this expression over all possible positions of theparticles to get a useful overall measure of the accuracy of the algorithm. The RMS errorof the corrected P3M energies is, by definition,(∆ E ( k )RMS ) := (cid:68) ( E ( k )P3M , corr − E ( k ) ) (cid:69) = (cid:10) (∆ E ( k ) − K ) (cid:11) (6.7)where we used (5.3) and (6.2). We can isolate in ∆ E ( k ) “interaction” terms ( i (cid:54) = j ) from selfterms ( i = j ):(∆ E ( k )RMS ) = (cid:28)(cid:16) ∆ E ( k )int + ∆ E ( k )self − K (cid:17) (cid:29) . (6.8)We recall from Sec.III that the interaction terms vanish on average for random systems: (cid:68) ∆ E ( k )int (cid:69) = 0. The correlation (cid:68) ∆ E ( k )self · ∆ E ( k )int (cid:69) = 0 (6.9)vanishes as well for the same reason (this is due to the fact that the average Ewald interactionenergy between a fixed particle i and a particle j (cid:54) = i is zero, see Appendix A). Eq. (6.8)17educes therefore to(∆ E ( k )RMS ) = (cid:68) (∆ E ( k )int ) (cid:69) + (cid:68) (∆ E ( k )self − K ) (cid:69) , (6.10)where the first term accounts for fluctuating errors in the interactions energies, and thesecond term accounts for fluctuating errors in the corrected Madelung self-energies of theions. Since the latter term may be written as (cid:68) (∆ E ( k )self − K ) (cid:69) = (cid:68) (∆ E ( k )self ) (cid:69) − K , (6.11)we remark that the shift − K derived in the previous section, in addition to removing thesystematic bias in the k -space energies, also reduces the fluctuating errors of the k -spaceself-energies by an amount − K .In the substraction ∆ E ( k )self − K , it can be seen, from (5.4) in which only i = j terms arekept and (6.4), that all terms containing (cid:101) φ ( k ) cancel out, so we have (cid:68) (∆ E ( k )self − K ) (cid:69) = (cid:42)(cid:16) L (cid:88) k ∈ e M k (cid:54) =0 (cid:101) G ( k ) (cid:88) i q i (cid:104) (cid:88) m (cid:88) m (cid:101) U ( k m ) (cid:101) U ( k m ) e ik g ( m − m ) · r i − (cid:88) m (cid:101) U ( k m ) (cid:105)(cid:17) (cid:43) (6.12)where we used the symmetry (cid:101) U ( − k ) = (cid:101) U ( k ) and introduced the shorthand notation k m = k + k g m . When the square is expanded, the summation over particles (cid:80) i becomes a doublesummation (cid:80) i,i (cid:48) . All terms with i (cid:48) (cid:54) = i vanish, because (cid:104) exp( ik g ( m − m ) · r i ) (cid:105) = δ m , m ,leaving identical sums over m which cancel each other. The remaining terms i (cid:48) = i evaluateto (cid:68) (∆ E ( k )self − K ) (cid:69) = 14 L ( (cid:88) i q i ) (cid:88) k ∈ e M k (cid:54) =0 (cid:88) k (cid:48) ∈ e M k (cid:48) (cid:54) =0 (cid:101) G ( k ) (cid:101) G ( k (cid:48) ) ×× (cid:110) (cid:88) m (cid:88) m (cid:88) m (cid:101) U ( k m ) (cid:101) U ( k m ) (cid:101) U ( k (cid:48) m ) (cid:101) U ( k (cid:48) m − m + m ) − (cid:88) m (cid:101) U ( k m ) · (cid:88) m (cid:48) (cid:101) U ( k (cid:48) m (cid:48) ) (cid:111) = 14 L ( (cid:88) i q i ) H (6.13)with H = 1 L (cid:88) k ∈ e M k (cid:54) =0 (cid:88) k (cid:48) ∈ e M k (cid:48) (cid:54) =0 (cid:101) G ( k ) (cid:101) G ( k (cid:48) ) (cid:110) (cid:88) m (cid:88) m (cid:54) = m (cid:88) m (cid:101) U ( k m ) (cid:101) U ( k m ) (cid:101) U ( k (cid:48) m ) (cid:101) U ( k (cid:48) m − m + m ) (cid:111) . (cid:80) i q i with thevalencies of the ions. The prefactor is somewhat complicated since it involves a doublesummation over wave vectors and a triple summation over alias indices m , m , m , but H can be evaluated reasonably fast. The numerical calculation of H can be acceleratedby taking profit of the symmetries (the sum over k can be restricted to only half an octant ofthe reciprocal mesh), and by skipping inner loops in the triple summation over alias indicesif the product of the charge fractions is almost zero.We calculate now the fluctuations of the errors in the interaction energies, i.e. the firstterm of Eq. (6.10). That term reads, using (6.8), (5.4), (5.2) and (2.12) and keeping onlyinteraction terms: (cid:68) (∆ E ( k )int ) (cid:69) = (cid:42)(cid:16) L (cid:88) k ∈ e M k (cid:54) =0 (cid:88) i,ji (cid:54) = j q i q j (cid:88) m ∈ Z e i k m · r i ×× (cid:104) (cid:101) G ( k ) (cid:88) m (cid:48) ∈ Z e − i k m (cid:48) · r j (cid:101) U ( k m ) (cid:101) U ( k m (cid:48) ) − e − i k m · r j (cid:101) φ ( k m ) (cid:105)(cid:17) (cid:43) . (6.15)The calculation of this average is straightforward, though somewhat tedious. We find thatit reduces to (cid:68) (∆ E ( k )int ) (cid:69) (cid:39) Q L H , (6.16)where H = 2 L (cid:88) k ∈ e M k (cid:54) =0 (cid:34) (cid:101) G ( k ) (cid:16) (cid:88) m (cid:101) U ( k m ) (cid:17) − (cid:101) G ( k ) (cid:88) m (cid:101) U ( k m ) (cid:101) φ ( k m ) + (cid:88) m (cid:101) φ ( k m ) (cid:35) . (6.17)The factor 2 in H originates from the fact that each pair of particles appears twice inthe sum over i and j (cid:54) = i in (6.15). Expression (6.17) is the analog for the energy of theparameter Q introduced by Hockney and Eastwood to measure the accuracy of the P3Mforces . Notice that (6.17) is given in real space by H = 2 V cell (cid:90) V cell d r (cid:90) L d r [ φ P3M ( r ; r ) − φ ( r − r )] (6.18)where φ P3M ( r ; r ) is the reciprocal potential at r created by a unit charge located at r , asobtained from the P3M algorithm. (This potential is given in Fourier space by combining19C4) with (5.2) in which we set ρ ( r ) = δ ( r − r ).) H is hence twice the squared deviationbetween the potential φ P3M obtained from the mesh calculation and the exact reciprocalpotential φ , summed over all relative positions r within the simulation box, and averagedover all possible positions of charge r in a mesh cell ( V cell = h ).Inserting the above results in (6.10), our final expression for the RMS error of the (cor-rected) P3M energies is∆ E ( k )RMS = (cid:112) Q H + ( (cid:80) i q i ) H L / (6.19)where H and H are defined in Eqs. (6.17) and (6.14). This error depends on theinfluence function (cid:101) G ( k ). The optimal influence function (the one that minimizes the error)will be determined in the next section. The above error estimate, together with the optimalinfluence function (6.21) and the constant shift (6.4) which must be applied to the P3Menergies, constitute the main results of this paper.The RMS error (6.19) displays two different scalings with the valencies of the ions:( (cid:80) i q i ) for errors coming from pair interactions (such a scaling also governs errors in P3Mforces ) and (cid:80) i q i for errors in Madelung self-energies. Because of these different scalings,the errors from pair interactions are expected to dominate in systems with many chargedparticles ( Q (cid:29) (cid:80) i q i ). Notice that H is, roughly speaking, proportional to ( (cid:80) k (cid:101) G ( k )) ,while H scales like (cid:80) k (cid:101) G ( k ). The errors in the Madelung self-energies increase there-fore more rapidly than the errors in the pair interaction energies when the Ewald splittingparameter α (and hence (cid:101) G ( k )) is increased, or when the size of the mesh is increased. Theimportance of the two source of errors (fluctuations in pair interaction energies versus fluc-tuations in Madelung self-energies) will be compared in Sec. VII for a test system with Q = 100. C. Optimal influence function We can now determine the optimal influence function (cid:101) G ( k ), by imposing the conditionthat it minimizes the RMS error (6.19). Since the errors coming from pair P3M interactionsare expected to dominate the self-interaction errors (except in systems with few particles),we optimize the influence function only with respect to the pair interactions. Setting δH δ (cid:101) G ( k ) = 0 , (6.20)20ives immediately (cid:101) G ( k ) = (cid:88) m ∈ Z (cid:101) U ( k m ) (cid:101) φ ( k m ) (cid:16) (cid:88) m ∈ Z (cid:101) U ( k m ) (cid:17) (6.21)where we recall that the Fourier-transformed reciprocal interaction (cid:101) φ ( k ) is given by (2.11).An optimization of the influence function with respect to the full RMS error could beperformed, but would require solving a linear system of M equations to compute (cid:101) G ( k ).The numerical results shown in Sec. VII will confirm that such a full optimization is notnecessary in typical systems.Since (cid:101) φ ( k ) decays exponentially fast, the optimal influence function is given in goodapproximation by (cid:101) G ( k ) (cid:39) (cid:101) φ ( k ) (cid:101) U ( k ) (cid:16)(cid:80) m ∈ Z (cid:101) U ( k m ) (cid:17) . (6.22) (cid:101) G ( k ) differs thus from (cid:101) φ ( k ) by a factor which is always less than one. This damping of theinteraction compensates as well as possible for the aliasing errors introduced by the use ofa fast Fourier transform. If (cid:101) U ( k ) were a perfect low-pass filter ( (cid:101) U ( k m ) = 0 if m (cid:54) = 0), noaliasing error would occur and the influence function would reduce to (cid:101) G ( k ) = (cid:101) φ ( k ) / (cid:101) U ( k ).This is indeed the result expected from (5.4) when aliasing errors are absent. The trueoptimal influence function (6.21) differs from this simple expression by contributions fromthe high-frequency spectrum of (cid:101) U ( k ) and reciprocal interaction (2.11).Hockney and Eastwood obtained the following optimal influence function by minimizingthe errors in the forces instead of the energy : (cid:101) G (forces) ( k ) = (cid:80) m ( k · k m ) (cid:101) U ( k m ) (cid:101) φ ( k m ) k (cid:16)(cid:80) m (cid:101) U ( k m ) (cid:17) (6.23)Obviously, this function is also given in very good approximation by (6.22). This explainswhy influence functions (6.21), (6.22) and (6.23) all give very similar results when computingenergies and forces.Inserting (6.21) into (6.17), we find that the minimal value of H is H (cid:12)(cid:12) min = 1 L (cid:88) k ∈ e M k (cid:54) =0 (cid:88) m ∈ Z (cid:101) φ ( k m ) − (cid:32) (cid:80) m (cid:101) U ( k m ) (cid:101) φ ( k m ) (cid:80) m (cid:101) U ( k m ) (cid:33) . (6.24)21his is the expression of H to be used in the RMS error estimate (6.19) when the P3Malgorithm is optimized to yield the smallest possible errors in the pair interaction energies.We recall that the errors in the P3M energies originate from aliasing effects (due to thesampling on a grid) and truncation errors (due to the fact that the reciprocal mesh containsonly a finite number of wave vectors). The truncation error can only be reduced by choosinga larger mesh or by using a reciprocal interaction with a faster decay in k -space, whereasthe aliasing errors may be reduced by increasing the order of the charge assignment function(up to the maximum order allowed by the size of the mesh). The intrinsic truncation errorof a given mesh and reciprocal interaction can be obtained by assuming (cid:101) U ( k ) in (6.24) tobe a perfect low-pass filter: H , cut − off (cid:12)(cid:12) min = 1 L (cid:88) k ∈ e M k (cid:54) =0 (cid:34) (cid:88) m ∈ Z (cid:101) φ ( k m ) − (cid:101) φ ( k ) (cid:35) . (6.25)By inserting this formula in (6.19), we get an estimate of the intrinsic RMS cut-off errorin k -space, caused by the finite number of wave vectors in the reciprocal mesh. The RMSerror associated with (6.25) depends only on the size of the mesh and on the choice of thereciprocal interaction, i.e. Ewald parameter α if the standard form (2.3) is used. VII. NUMERICAL CHECK OF ACCURACY In this section, we test the analytical results (optimal influence function, energyshift E cutP3M , RMS error estimate) derived in the previous section. We do this by comparingthe P3M energies with the exact energies calculated in a specific random system. In thefollowing, all dimensions are given in terms of the arbitrary length unit L and charge unit C . In particular, energies and energy errors are given in units of L / C . We choose thesame test system as the one defined in Appendix D of Deserno and Holm : 100 particlesrandomly distributed within a cubic box of length L = 10 L , half of them carry a positive,the other half a negative unit charge. The statistical average (cid:104)· · ·(cid:105) is calculated by aver-aging over at least 100 different configurations of this test system (these configurations aredetermined by using the same random number generator as in Deserno and Holm ). Wellconverged Ewald sums (in metallic boundary conditions) were used to compute the exactenergies of the test systems. The first three systems have energies − . − . − . Q ζ/ (cid:39) − . M =4 , , , r cut = 4 . 95, and different orders of the charge assignmentfunction (from 1 to 7). Our calculations show that using the energy-optimized influencefunction (6.21), instead of the force-optimized influence function (6.23), leaves the energiesalmost unchanged. (A slight improvement in accuracy appears only when the aliasing errorare at their maximum, namely for a charge assignment order of 1 and large values of α .)This behavior could have been expected, since both influence functions are almost equivalentto the simple formula (6.22).We compare in Fig. 1 the measured systematic error (cid:104) E P3M − E exact (cid:105) of the uncorrectedP3M energies [Eq. (6.6) without term E cutP3M ], for CAO’s ranging from 1 to 7, to the expectedbias − E cutP3M . The agreement is perfect for all CAO’s and for all values of Ewald’s splittingparameter α . The energy shift E cutP3M in Eq. (6.6) removes therefore entirely the systematicerror, as it should.Fig. 1 illustrates that the systematic errors in the uncorrected energies, which are dueto cut-off and aliasing errors in the Madelung self-energies of the ions, have two differentcontributions of opposite sign. At small values of α , the r -space cut-off error dominates andleads to an overestimation of the energy because the negative interaction energy of an ionwith the neutralizing background charge provided by the other particles is not fully takeninto account. The cut-off correction (3.11) derived in Sec. III does compensate very well forthis effect. At large values of α , k -space cut-off and aliasing errors dominate, and lead to anunderestimation of the Madelung self-energies (expression (6.4) is indeed always negative).Since the systematic error (cid:104) E ( k )P3M − E ( k )exact (cid:105) in the reciprocal energies arise solely from selfterms (the Ewald interaction between a pair of particles is zero on average), this error canalternatively, and more efficiently, be measured by computing the P3M energy of a systemmade up of a single ion in the box, averaging that energy over different positions of theparticle relative to the mesh. To restore electro-neutrality, the interaction energy (2.13)with the (implicit) neutralizing background must of course be taken into account beforecomparing the result with the exact Madelung self-energy q ζ/ α -10-50 < E P M - E e x ac t > MeasuredPredicted CAO 1 CAO 2 CAO 3CAO 7 FIG. 1: Comparison between the measured systematic error of the uncorrected P3M energies(crosses) and the theoretical prediction − E cutP3M (solid lines) as a function of Ewald parameter α .The average is performed over 1000 test systems consisting in 100 charges located at random in abox of size L = 10. The mesh has size M = 8 and the real-space cut-off is r cut = 4 . easily be transposed to any cubic system with an arbitrary number of ions since the energyshift E cutP3M scales merely as Q /L .Having validated the energy shift (6.5), we test now the accuracy of the RMS errorestimate (6.19). We show in Fig. 2 the theoretical predictions for the RMS error of thecorrected P3M energies for different mesh sizes (thick solid lines), at fixed CAO 2. Thedominant error at small values of α comes from the truncation in the real-space calculation,while k -space cut-off and aliasing errors dominate at large values of α . The plot showsalso separately the contribution ∆ E ( k )self , which accounts for fluctuating errors in the k -spaceMadelung self-energies, and the contribution ∆ E ( k )int which accounts for fluctuating errors inthe P3M pair interaction energies. Near the optimal value of α , the error ∆ E ( k )int dominatesslightly ∆ E ( k )self by half an order of magnitude. This validates the use of the optimal influencefunction (6.21), which was designed to minimize errors in the pair P3M interaction energiesonly. Notice that ∆ E ( k )self overcomes ∆ E ( k )int at large values of α , in agreement with the scalingwith α discussed in Sec. VI. The errors in Madelung self-energies must therefore by included24o predict correctly the full RMS error curve in our test system with 100 charged particles,but they are expected to become negligible when the number of ions is increased above afew hundred. a R M S e rr o r D E RMS D E (k)self D E (k)int m e s h m e s h m e s h FIG. 2: Theoretical predictions for the RMS error of the corrected P3M energies for CAO 2and three different mesh sizes [thick solid lines], for the same system and real-space cut-off as inFig. 1. The two contributions which make up the total k -space error are also shown independently:RMS error in the pair P3M interaction energies (Eq. (6.16), thin solid line) and RMS error in theMadelung self-energies (Eq. (6.13), dashed line). The predicted RMS errors agree very well with the measured RMS errors, as shown inFig. 3. The small deviations at low values of α are due to a loss of accuracy of Kolafaand Perram’s r -space error estimate (2.14), and to the fact that this error estimate doesnot take into account the improvement in accuracy brought by the new cut-off correctionterm (3.11). In the regime where the dominant error comes from the k -space calculation,the agreement with our RMS error estimate is excellent, especially at high values of thecharge assignment order. The errors in the k -space calculation are caused by truncationand aliasing effects. The aliasing errors can be reduced by increasing the charge assignmentorder, but the accuracy cannot go below the minimum k -space cut-off error (6.25) (dashedcurve in Fig. 3), which is intrinsic to the mesh size and choice of reciprocal interaction.25 α -3 -2 -1 ∆ E R M S Measured RMS errorPredicted RMS errorMinimum cutoff error CAO 1CAO 2CAO 3CAO 4CAO 5CAO 6CAO 7 FIG. 3: Comparison between the measured (crosses) and predicted (solid lines) RMS errors of the(corrected) P3M energies, for the same test system, mesh size and real-space cut-off as in Fig. 1.The minimal error due to direct and reciprocal space cut-offs is shown as a dashed line. The pronounced minimum in the RMS error curves stresses the importance of using theoptimal value of α when performing simulations with the P3M algorithm (or with the othervariants of mesh based Ewald sums). Our accurate RMS error estimate for the P3M energiescan be used to quickly find the optimal set of parameters (mesh size, charge assignmentorder, Ewald splitting parameter) that lead to the desired accuracy with a minimum ofcomputational effort . Whatever the chosen parameters, it can serve also as a valuableindicator of the accuracy of the P3M energies. VIII. CONCLUSIONS In this article, we discussed in detail which ingredients are necessary to utilize the P3Malgorithm to compute accurate Coulomb energies of point charge distributions. The usageof a nearly linear scaling method ( ≈ N log N ) like P3M is almost compulsory for systemscontaining more than a few thousand charges.In particular, we derived the cut-off corrections for the standard Ewald sum transparentlyand interpreted the systematic errors in terms of Madelung energies. This route lead us to26n additional real-space cut-off correction term that has so far not been discussed in theliterature. Building on these results, we have deduced the k -space cut-off correction term inthe case of the P3M algorithm, where additional aliasing errors play a role. Furthermore wederived the exact form of the influence function that minimizes the RMS errors in the ener-gies, and showed that this function is not much different from the force-optimized influencefunction, which a posteriori justifies why in most P3M implementations the usage of theforce-optimized influence function does not lead to inaccurate results. Based on the energyoptimized influence function we derive an accurate RMS error estimate for the energy, andperformed numerical tests on sample configurations that demonstrate the validity of ourerror estimates and the necessity to include our correction terms. We also demonstratedthat the electrostatic energy of an individual particle in the system can be obtained in theP3M method, but at the expense of an additional inverse fast Fourier transform.With the help of the newly derived error estimates we can easily tune the desired accu-racy of the P3M algorithm and find suitable parameter combinations before running anysimulation.The P3M algorithm can be generalized along our discussed lines to compute other longrange interactions. Of particular interest are dipolar energies, forces and torques, and theassociated error estimates for these quantities. This will be the content of a forthcomingpublication. Our P3M generalization for the energies will be included in a future version ofthe molecular simulation package Espresso , that is freely available under the GNU generalpublic license. The website provides up-to-date information. Acknowledgements Funds for this research were provided by the Volkswagen Stiftung under grant I/80433and by the DFG within grant Ho-1108/13-1. APPENDIX A: EWALD PAIR POTENTIAL AND MADELUNG SELF-ENERGY The Ewald formula for the electrostatic energy E of a periodic charged system can bewritten in a form that underlines the fact that E includes the Madelung self-energies Q ζ/ Q = (cid:80) i q i and ζ is defined in (3.7)]. We recall from Sec. II that the Ewald27ormula for E reads, if the system is globally neutral and if we employ metallic boundaryconditions, E = 12 (cid:88) i,j (cid:88) (cid:48) n ∈ Z q i q j ψ ( r ij + n L ) + 12 L (cid:88) i,j q i q j (cid:88) k ∈ K k (cid:54) =0 e − i k · ( r i − r j ) (cid:101) φ ( k ) − Q α √ π . (A1)The “self-energy terms” in E , i.e. term E ( s ) and terms i = j , are Q (cid:16) (cid:88) n (cid:54) =0 ψ ( n L ) + 1 L (cid:88) k (cid:54) =0 (cid:101) φ ( k ) − α √ π (cid:17) = Q (cid:16) ζ + πα L (cid:17) . (A2)We can write therefore E = 12 (cid:88) i (cid:54) = j q i q j (cid:16) (cid:88) n (cid:54) =0 ψ ( r ij + n L ) + 1 L (cid:88) k (cid:54) =0 e − i k · r ij (cid:101) φ ( k ) (cid:17) + Q (cid:16) ζ + πα L (cid:17) = 12 (cid:88) i (cid:54) = j q i q j V Ewald ( r ij ) + Q ζ (A3)where we defined the Ewald pair interaction V Ewald ( r ) = (cid:88) n (cid:54) =0 ψ ( r + n L ) + 1 L (cid:88) k (cid:54) =0 e − i k · r (cid:101) φ ( k ) − πα L . (A4)Notice that in writing (A3), we used (cid:80) i (cid:80) j (cid:54) = i q i q j π/ ( α L ) = − Q π/ ( α L ) which followsfrom electro-neutrality. Thanks to the inclusion of this constant in the definition of V Ewald ( r ),the Ewald pair potential does not depend on the parameter α [ ∂/∂α V Ewald ( r ) = 0] and itsaverage over the simulation box is zero : (cid:104) V Ewald ( r ) (cid:105) = 1 L (cid:90) L d r V Ewald ( r ) = 0 . (A5)The latter property is simply a consequence of (cid:104) exp( i k · r ij ) (cid:105) = δ k , and Eq. (3.5).In conclusion, expression (A3) shows explicitly that the electrostatic energy of a periodiccharged system includes the Madelung self-energies Q ζ/ . The fact that theEwald interaction between a pair of particles averages to zero when one particle exploresthe whole simulation box is also noteworthy aspect of Ewald potential . APPENDIX B: PROOF OF EQ. (5.2) Eq. (5.2) is a consequence of the Sampling theorem [refs] and is straightforward to demon-strate. The sum in (4.3) is rewritten as an integral (cid:101) ρ M ( k ) = h (cid:90) d r (cid:48) (cid:90) V d r W ( r ) U ( r − r (cid:48) ) ρ ( r (cid:48) ) e − i k · r (B1)28here we used (4.1) and introduced an infinite mesh of Dirac delta functions W ( r ) = (cid:88) r m ∈ M p δ ( r − r m ) = 1 h (cid:88) m ∈ Z e − ik g m · r . (B2)(We recall that k g = 2 π/h ). Using the above representation of W ( r ) and introducing in (B1)the Fourier series representation of the periodic charge density, ρ ( r (cid:48) ) = 1 L (cid:88) k (cid:48) ∈ K (cid:101) ρ ( k (cid:48) ) exp( i k (cid:48) · r (cid:48) ) , (B3)we recover the result (5.2) after straightforward simplifications. APPENDIX C: PROOF OF EQUIVALENCE BETWEEN EQS. (4.9) AND (4.6) Eq. (4.9) is equivalent to E ( ks )P3M = 12 V (cid:88) k ∈ K k (cid:54) =0 (cid:101) ρ ∗ ( k ) (cid:101) Φ( k ) , (C1)where (cid:101) Φ( k ) is the full Fourier transform ( k ∈ K ) of the back-interpolated potentialmesh (4.8): (cid:101) Φ( k ) = (cid:90) V Φ( r )d r = h (cid:90) V d r e − i k · r (cid:90) d r (cid:48) W ( r (cid:48) ) U ( r − r (cid:48) )Φ M ( r (cid:48) ) . (C2)We replace in this equation Φ M ( r (cid:48) ) and W ( r (cid:48) ) by their expressions (4.7) and (B2), andperform the integration over r (cid:48) : (cid:101) Φ( k ) = 1 L (cid:88) k (cid:48) ∈ e M (cid:101) Φ M ( k (cid:48) ) (cid:88) m (cid:90) V d r e − i k · r (cid:101) U ( k (cid:48) + k g m ) e i ( k (cid:48) + k g m ) · r . (C3)The integration over r introduces a Kronecker symbol δ k , k (cid:48) + k g m . We get therefore the simpleresult (cid:101) Φ( k ) = (cid:101) Φ M ( k ) (cid:101) U ( k ) (C4)where the function (cid:101) Φ M ( k ), which is defined originally only for k ∈ (cid:101) M , is now understoodto be extended periodically to all K space. Notice that the inverse FFT does not introducealiasing errors: the sum over m merely renders (cid:101) Φ M ( k ) periodic. In accordance with (4.3)29nd (4.5), we extend also (cid:101) ρ M ( k ) and (cid:101) G ( k ) periodically, with period 2 π/h . Using the aboveresult and (4.5), the reciprocal energy (C1) can be expressed as E ( ks )P3M = 12 L (cid:88) k ∈ K (cid:101) ρ ∗ ( k ) (cid:101) U ( k ) (cid:101) ρ M ( k ) (cid:101) G ( k ) (C5)= 12 L (cid:88) k ∈ e M (cid:88) m ∈ Z (cid:101) ρ ∗ ( k + k g m ) (cid:101) U ( k + k g m ) (cid:101) ρ M ( k ) (cid:101) G ( k ) . (C6)This may be compared with Eq. (4.6), i.e. E ( ks )P3M = 12 L (cid:88) k ∈ e M k (cid:54) =0 (cid:101) ρ ∗ M ( k ) (cid:101) ρ M ( k ) (cid:101) G ( k ) . (C7)Recalling (5.2) and the fact that (cid:101) U ( k ) is real, we see that both expressions are equivalent. A. Arnold and C. Holm,in Advanced Computer Simulation Approaches for Soft Matter SciencesII , eds. C. Holm and K. Kremer (Springer, Berlin, 2005) R.W. Hockney and J.W. Eastwood, Computer Simulation Using Particles (IOP, Bristol, 1988) M. Deserno and C. Holm, J. Chem. Phys. (1998): 7678 U. Essmann, L. Perera et al., J. Chem. Phys. (1995): 8577 M. Deserno and C. Holm, J. Chem. Phys. (1998): 7694 S.W. de Leeuw, J.W. Perram and E.R. Smith, Proc. R. Soc. London, Ser. A (1980):57 J.-M. Caillol, J. Chem. Phys. (1994): 6080 E.R. Smith, Mol. Phys. (1988): 1089 V. Ballenegger and J.-P. Hansen, J. Chem. Phys. (2005): Art. 114711 J. Kirkwood, J. Chem. Phys. (1939): 911 A. Alastuey and V. Ballenegger, Physica A (2000): 268 V. Ballenegger and J.-P. Hansen, Mol. Phys. (2004): 599 P. Ewald, Ann. Phys. (Leipzig) (1921): 253 G. Hummer, L.R. Pratt and A.E. Garc´ıa, J. Phys. Chem. (1995):14188 M.P. Allen and D.J. Tildesley, Computer Simulation of Liquids Oxford Science Publications.(Clarendon Press, Oxford, 1987) J. Kolafa and J.W. Perram, Mol. Sim. (1992): 351 S.G. Brush, H.L. Sahlin and E. Teller, J. Chem. Phys. (1966): 2102 B.R. Nijboer and T.W. Ruijgrok, J. Stat. Phys. (1988): 361 T. Darden, D. Pearlman and L.G. Pedersen, J. Chem. Phys. (1998): 10921 Z. Wang and C. Holm, J. Chem. Phys. (2001): 6351 H.J. Limbach, A. Arnold, B.A. Mann, and C. Holm, Comp. Phys. Comm. , (2006): 704 G. Hummer, L.R. Pratt and A.E. Garc´ıa, J. Phys. Chem. A (1998):7885(1998):7885