Damped reaction field method and the accelerated convergence of the real space Ewald summation
aa r X i v : . [ c ond - m a t . m t r l - s c i ] N ov Damped reaction field method and the accelerated convergence ofthe real space Ewald summation.
Victor H. Elvira and Luis G. MacDowell ∗ Departamento de Qu´ımica F´ısica, Facultad de Ciencias Qu´ımicas,Universidad Complutense de Madrid, 28040, Spain. (Dated: February 18, 2018)
Abstract
In this paper we study a general theoretical framework which allows to approximate the realspace Ewald sum by means of effective force shifted screened potentials, together with a self term.Using this strategy it is possible to generalize the reaction field method, as a means to approximatethe real space Ewald sum. We show that this method exhibits faster convergence of the Coulombenergy than several schemes proposed recently in the literature while enjoying a much more soundand clear electrostatic significance. In terms of the damping parameter of the screened potential,we are able to identify two clearly distinct regimes of convergence. Firstly, a reaction field regimecorresponding to the limit of small screening, where effective pair potentials converge faster thanthe Ewald sum. Secondly, an Ewald regime, where the plain real space Ewald sum converges faster.Tuning the screening parameter for optimal convergence occurs essentially at the crossover. Theimplication is that effective pair potentials are an alternative to the Ewald sum only in those caseswhere optimization of the convergence error is not possible. NTRODUCTION
The increasing availability of large scale computer facilities is allowing to study ever morecomplicated systems, with greater detail as well as longer length and time scales. Despite thisprogress, the accurate evaluation of electrostatic interactions remains the most importantbottleneck in molecular simulations of charged systems.[1, 2]This uncomfortable situation is reflected in the number of different alternatives whichare available in the literature in order to deal with Coulombic interactions.[3–11] Yet, it isclear that the benchmark for both efficiency and accuracy of all studies remain the Ewaldsummation technique.[12–14]In this method, the full electrostatic energy of the system is split into a real space contri-bution, which amounts to a pairwise summation of an effective damped Coulomb potential,and a Fourier contribution, which embodies the long range effects of the Coulomb inter-actions. The latter term features a Fourier transform of the charge distribution, which onthe one hand, brings some conceptual difficulties,[4, 14–16] and on the other, is very timeconsuming to calculate.[17]In the last decade, a number of studies have been devoted to study more or less efficientmethodologies that allow to calculate electrostatic interactions while avoiding the cumber-some Fourier contributions of the Ewald sum.[8, 11, 18–21] Such techniques, named recentlyunder the provocative name of pairwise alternatives to the Ewald sum, have recently receivedconsiderable popularity, but also some degree of controversy as regards efficiency,[20] andaccuracy.[22, 23]The fact is that a pairwise alternative to the Ewald summation has been available eversince the first few simulations of Coulombic systems.[3, 24, 25] Indeed, the Reaction FieldMethod is almost as old as computer simulations of charged systems,[1] yet, it has a cleartheoretical background which more recent approaches lack completely.[11, 18, 20] Despitethis situation, the Reaction Field Method seems to have been largely abandoned in favor ofother techniques, with some exceptions.[26]Recently, Fukuda et al. proposed a heuristic approach to approximate the Coulomb sum.This approach shares advantages of some of the pairwise methods, in the sense that it screensthe Coulomb interactions with a fast decaying function, but has a somewhat more elaborateelectrostatic background.[19] Indeed, it has been recently recognized that this approach may2e considered as a generalization of the reaction field method to screened potentials.[27]In this work we attempt to provide a sound theoretical background for a generalizedreaction field method that achieves fast convergence of the Coulomb sum as with several ofthe most popular pairwise alternatives.[8, 11] The theoretical calculations are supplementedwith numerical results for test systems, which demonstrate the superiority of reaction fieldmethods. Additionally, we perform detailed analysis of the Coulomb sum convergence errorfor either reaction field and Ewald methods. This will allow us to identify the region ofconvergence parameters where each method is advantageous.
THEORETICAL BACKGROUNDEwald Summation
In most simulations of charged systems, the influence of long range electrostatic inter-actions is estimated by assuming the finite sample is surrounded by an infinite number ofreplicas. The energy felt by, say, charge i , may be then estimated as a lattice sum over thereplicas: U i = X n =0 ′ X j q i q j | r ij + n | (1)where n denotes a translation of the unit cell vector, q i is the charge on particle i , and r ij is the position of j relative to i . Furthermore, it is understood that the first sum runs overall possible unit cell translations, while the second sum runs over all charges inside the unitbox. A prime reminds that j must be different from i when n = .This trick only gets rid of the boundary problem, but not of the actual calculation of U i ,since the series is known to be conditionally convergent, with a very slow convergence inthe favorable cases. A particularly inconvenient case is a summation over spherical shells,which is often not convergent.[4, 8]In the classical treatment of Ewald, the conditionally convergent lattice sum of chargesis transformed into two rapidly convergent series, such that: U i = U R i + U F i (2)with U R i = X n =0 ′ X j q i q j e rf c ( α | r ij + n | ) | r ij + n | (3)3nd U F i = 4 π X k = ρ ( k ) e − k α k q i e − i k · r i − q i α √ π (4)where k are vectors in Fourier space and ρ ( k ) is the Fourier transform of the charge density.Additionally, the Ewald sum may contain a surface term, U surf , which accounts for theboundary conditions of the system. Such term arises strictly from the long range electrostaticinteractions with the boundary and surrounding medium, and is therefore essentially aFourier contribution corresponding to the missing k = 0 term of the reciprocal space sum.[4,15, 16] In practice, for bulk systems under metallic boundary conditions the surface termmay be neglected, so that we will henceforth drop this complication.The inverse length, α , plays a key role, dictating the convergence of the series. A largevalue of α leads to a rapidly convergent real space series that can be actually truncatedalready at n = 0, but then the U F contributions converges slowly. Alternatively, a smallvalue of α produces a very fast convergence of U F , but then U R is slowly convergent.An important observation made by Wolf et al. is that for moderately small values of α ,most of the Fourier contribution is actually given by the simple self term, U self = − q i α √ π , (5)so that the expensive reciprocal space summation may be ignored altogether.[8]This observation has allowed for the recent development of efficient methodologies for thecalculation of electrostatic interactions.[8, 11, 18, 19] Accelerating the convergence of the real space sum
Having circumvented the problem of calculating the expensive reciprocal space sum, therestill remains a crucial issue: how fast is the convergence of U R in those cases where thereciprocal space summation may be ignored?For most practical purposes, e rf c ( αr ) decays so fast that the lattice summation requiredfor the evaluation of U R i may be ignored. Rather, a plain spherical cutoff is usually employedfor distances larger than a cutoff radius, R c , with the hope that terms of order e rf c ( αR c ) orsmaller may be neglected. However, it remains desirable to have a cutoff as small as possible.For this purpose, one may try to exploit the techniques of continuum electro–dynamics for4 → FIG. 1. Sketch of the charge neutralization schemes employed in the Wolf and Damped ReactionField methods. In the Wolf method (left), the net charge inside the sphere surrounding the centralion is uniformly spread over its surface. In the Damped Reaction Field Method (right), an addi-tional surface charge due to broken effective dipoles between the central ion and the surroundingcharges is considered. For a perfectly regular crystal lattice, the charge distribution is also uniform,but will usually be non-uniform and therefore produce a net force. the e rf c ( αr ) /r potential in the same spirit as it has been successfully done for the plainCoulomb potential.[6, 8, 24, 28] This would allow to truncate U R i at cutoffs where e rf c ( αr )is actually not negligible.[8, 11, 18, 19]In this regard, another important observation made by Wolf is that Coulomb sums overspherical shells show good convergence whenever the net charge inside the sphere vanishes(Fig. 1).[29] Whence, the poor convergence is due to the the fact that spheres of arbitraryradius centered about an ion i will usually carry a net charge, ∆ q i . The convergence ofthe Coulomb sum over spheres can be much improved by summing over charge neutralizedspheres. This may be achieved in practice by assuming a fictitious uniform charge σ charge i over the surface of the sphere:[8] σ charge i = 14 πR c ∗ X j q j (6)where the asterisk superscript indicates that the summation is restricted to those particles j inside the cutoff sphere (such that r ij < R c ).Whence, for a source charge i with arbitrary electric potential q i φ ( r ), an improved ap-proximation for the energy is: U i ≈ ∗ X j q i q j φ ( r ij ) − q i ∆ q i φ ( R c ) (7)Surprisingly, this charge neutralization may be implemented by means of an effective pairwise5otential:[8] u = q i q j { φ ( r ij ) − φ ( R c ) } (8)where the first term in the right hand side accounts for the bare interaction from the potential q i φ ( r ), while the second term is the energy resulting from the interaction with the fictitiouscharge distribution of Eq. (6). Clearly, such interaction merely shifts the pair potential.Therefore, it does not result in an additional force.However, one expects that the truncation of interactions at R c will not only produce aspurious net charge about ion i , but also a net polarization (Fig. 1).[29] The force stemmingfrom the uniform net charge is zero for reasons of symmetry, but the net polarization willresult in a finite electric field on charge i .[19]This idea may be elaborated quantitatively in terms of the sphere’s polarization and thelaws of electrostatics, which dictate that a polarized dielectric produces a surface chargedensity of magnitude: σ pol i ( e ) = P i · e (9)where P i is the net polarization inside the sphere of radius R c centered on i and e is a unitvector normal to the surface. Since the medium is overall neutral, we assume that the chargesbeyond the sphere act such as to cancel exactly this charge distribution. Accordingly, thefull charge neutralizing distribution is: σ neutral i = − ( σ charge i + σ pol i ( e )) (10)The particle i bears a potential φ ( r ), and thus exerts a force −∇ φ ( r ) on whichever othercharge. Now, consider the infinitesimal force, d F surf i felt by charge i due to the chargeneutralizing distribution on an infinitesimal surface element d S = R sin θ d θ d φ :d F surf i = − q i ∇ φ | R c σ neutral i d S (11)Integrating over all the sphere’s surface, the term stemming from the uniform distribution σ charge i vanishes for reasons of symmetry, but the non uniform σ pol i ( e ) term yields: F surf i = 4 π q i d φ d r (cid:12)(cid:12)(cid:12)(cid:12) R c R c P i (12)After substitution of P i in terms of the explicit charges inside the sphere we obtain: F surf i = q i ∗ X j q j d φ d r (cid:12)(cid:12)(cid:12)(cid:12) R c r ij R c r ij r ij (13)6he total force effectively felt on i then contains the actual interactions between particlesinside the sphere, together with the surface term accounting effectively for interactions withthe remaining charges: F i = − ∗ X j q i q j ∇ φ | r ij + q i ∗ X j q j d φ d r (cid:12)(cid:12)(cid:12)(cid:12) R c r ij R c r ij r ij (14)Clearly, the net force may be cast in terms of an effective pair potential of the form: F ij = − q i q j d φ d r (cid:12)(cid:12)(cid:12)(cid:12) r ij − d φ d r (cid:12)(cid:12)(cid:12)(cid:12) R c r ij R c ! r ij r ij (15)This is the formal result we sought for. It corresponds to a continuous linearly screenedforce which smoothly vanishes at the cutoff. Discussion
In order to understand its significance, it is convenient to recall the result for the effectiveelectrostatic potential of the Reaction Field Method, which assumes a continuous dielectricmedium with dielectric constant ǫ beyond the cutoff:[30] F ij = q i q j (cid:18) r ij − ǫ − ǫ + 1 1 R c r ij R c (cid:19) r ij r ij (16)Comparing this result with Eq. (15), it is clearly seen that our result is recovered for thespecial case where φ ( r ) = 1 /r , as in the Coulomb potential, and the additional assumptionof ǫ = ∞ . Whence, Eq. (15) corresponds to a generalization of the Reaction Field Methodfor arbitrary electric potentials, and a particular choice for the boundary conditions at thesurface of the cutoff sphere . The Reaction Field Method shows that taking into accountexplicitly a finite susceptibility beyond the cutoff requires to assume a surface charge densityof σ pol i = (2 ǫ − / (2 ǫ + 1) P i · e in place of Eq. (9).On the other hand, the assumption embodied in Eq. (9), and other methods,[8, 19, 21]implies that the charges left out beyond the cutoff are fully available to screen the net dipolecreated inside the sphere. This statement may be understood if we consider yet anotherrefinement over the Reaction Field Method, namely, to assume that the charge distributionoutside the cutoff sphere is given as a Boltzmann weighted average dictated by the field of theions inside the sphere. This task may be accomplished in the Debye-H¨uckel approximation,7nd yields a generalization of the Reaction Field Method that accounts for interactions withboth a dielectric continuum and a smooth charge distribution of given concentration:[25, 31] F ij = q i q j (cid:18) r ij − (2 ǫ − κR c ) + ǫ ( κR c ) (2 ǫ + 1)(1 + κR c ) + ǫ ( κR c ) R c r ij R c (cid:19) r ij r ij (17)where κ is the inverse Debye screening length of the free charges. For vanishing concentrationof charges, κ = 0, and this model recovers Eq. (16) exactly. In the opposite limit, κ becomes infinite, and then Eq. (17) becomes equal to the Reaction Field Method withconducting boundary conditions. This illustrates our statement, that the approximationEq. (9) corresponds to an infinite availability of charges for screening of the dipole insidethe cutoff sphere. Unfortunately, the accuracy of Debye-H¨uckel theory is limited to very lowion concentration, so the refinement of Eq. (17) is more a conceptual improvement than anaccurate working equation at typical simulation conditions.The fact is, once the interactions are truncated beyond a cutoff, one can not do withoutan arbitrary approximation as to the charge distribution of the surrounding medium. Thefirst obvious choice is to assume a dielectric response, but then the precise value of thedielectric constant needs to be specified. In simulations of molten salts or ionic fluids, itseems reasonable to choose metallic boundary conditions. On the other hand, in simulationsof polar fluids the choice of ǫ equal to the fluid’s dielectric constant is more natural. If,however, the polar fluid is simulated with explicit charges, as is usually the case, there isthen not an obvious choice.Fortunately, most solvents of interest in studies of charged systems, and particularlywater, have a rather large dielectric constant, so that the ratio (2 ǫ − / (2 ǫ + 1) is veryclose to unity, as in the case of ǫ = ∞ . Furthermore, for finite ǫ , the force of the RFMbecomes discontinuous at the cutoff, and this severely hampers applications in MolecularDynamics. Such inconvenience may be altogether avoided, since it has been shown that theprecise choice of dielectric boundary conditions does not significantly change the outcomeof the simulations, particularly for the case of phase coexistence.[26, 32]For these reasons, we believe the boundary conditions that are implied in Eq. (15) arethe most judicious choice for condensed phases of both i) polar fluids with high dielectricconstant and ii) molten salts that will be studied in the next section. Indeed, they correspondto the accepted choice in a large body of simulations.[8, 19, 21, 26, 32] A word of cautionis required for applications to very low density systems. In such cases, where typically ions8re separated by large distances and the Debye screening length is very large, one cannotexpect that the net charge about a single isolated ion will be neutralized at all within asmall cutoff. In such cases, an Ewald type summation might be the only reliable alternative. Damped force
The advantage of Eq. (15) over the traditional RFM is that it allows to accelerate theconvergence by using a fast decaying electric field φ ( r ) in place of the Coulomb potential.For φ ( r ) = e rf c ( α ) /r , it yields a generalized damped reaction field equation for the force,as: F ij = q i q j (" e rf c ( αr ij ) r ij + 2 α √ π e − α r ij r ij − " e rf c ( αR c ) R c + 2 α √ π e − α R c R c r ij R c ) (18)For the special case where α = 0, we recover the result of Eq. (16) with ǫ = ∞ .It is difficult here to know exactly how Eq. (18) is related to the corresponding expressionsgiven by Wolf et al.[8] The reason is that these authors obtain their forces unconventionallyas an imaginary process where r ij → R c in order to enforce their potential to yield a shiftedforce, which would otherwise not be the case. Whence the authors write: F ij = q i q j ( r ij − R c r ij R c (cid:12)(cid:12)(cid:12)(cid:12) r ij = R c ) (19)We find difficult to interpret what the condition r ij = R c means. If we simply ignore theodd condition, we recover Eq. (15). If, as interpreted by Fennell and Gezelter, we take theequality as written, we then obtain merely a shifted force potential of little electrodynamicsignificance.[11] Formal treatment
In the previous section we have obtained our results using elementary electrodynamics andback of the envelope arguments. Here we derive our results using a more formal treatmentof molecular electrodynamics as formulated by Neumann, Boresch and Steinhauser.[24, 33]Our derivation starts with a basic equation for the electric field that results from thepolarization E ( r ) of a uniform medium due to an external source field, E : E ( r ) = E ( r ) + Z d r ′ T ( r − r ′ ) · P ( r ′ ) (20)9here T ( r ) = ∇∇ φ is the dipole–dipole tensor, and the integration is over the simulationcell under toroidal boundary conditions. This equation was instrumental in the establish-ment of a consistent framework for the calculation of dielectric relaxation phenomena in the1980’s,[24, 34] and has been exploited recently to study the dielectric constant of anisotropicmedia.[35]An important point here is the realization that this equation is also valid for the case where E ( r ) results from atomic charge distributions interacting via a generalized Green’s function φ ( r ), whether it is a Coulomb interaction or some other modified central force potential.The only caution is to keep in mind that the dipole–dipole tensor has to be accordinglymodified as in the definition above. Taking these two considerations into account, we canwrite an equation for the electric field on charge r i that results from the presence of a secondcharge r j interacting with i via the generalized green function φ mod ( r ) = S ( r ) /r , with S ( r )an arbitrary function:[33] E ( r i ) = −∇ φ mod ( r ij ) + Z d rT mod ( r i − r ) · P ( r ) (21)where now T mod ( r i ) = ∇∇ φ mod ( r ). For our purposes, it proves convenient to express T mod in terms of the dipole–dipole tensor corresponding to a Coulomb potential:[33] T mod = ( S − rS ′ + 13 r S ′′ ) T + 13 S ′′ r I (22)where I is a unit matrix and the prime indicates derivation with respect to r .In order to perform the integral of Eq. (21), one needs to take into account that T ( r ) isan odd function, except at the singularity r → , where the tensor is:[24]lim r → T ( r ) = − π δ ( r ) (23)As a result, the integral does not vanish altogether, but rather, yields: E ( r i ) = −∇ φ mod ( r ij ) − π S (0) P ( r i ) + 13 Z d r S ′′ r I · P ( r ) (24)this is a general result for a wide choice of S functions, including whatever polynomial, theexponential or the complementary error function. The integral that remains is difficult tosolve for the general case of a position dependent polarization. However, assuming constantpolarization within the cutoff sphere, the integral can be solved by parts, yielding: E ( r i ) = −∇ φ mod ( r ij ) − π φ ′ mod ( R c ) R c P i (25)10here P i is the uniform polarization inside the sphere about i , that is, P i = q j r ij / πR c .Substitution of P i into the above result and multiplying by the charge q i at i then leadsright away to the general result obtained in the previous section, Eq. (15).The significance of calculating an electric field rather than a potential can be now un-derstood, since a continuous potential at r = R c results ab–initio by integration of theforce: u ( r ) = − Z r ∞ F ( r ) d r (26)considering that F ( r ) vanishes beyond R c , we obtain: u ( r ij ) = q i q j (cid:26) φ mod ( r ij ) − φ mod ( R c ) − φ ′ mod ( R c ) r ij − R c R c (cid:27) (27)Clearly, the resulting potential not only shifts φ mod ( r ) but also produces an additional po-larization force (c.f. Eq. (8)). The equation has the form of a generalized reaction fieldpotential, that may be employed to improve the summation of whatever generalized φ ( r )function. For φ = 1 /r we recover the reaction field model of Neumann,[24] properly mod-ified to produce a continuous force at r = R c as in Ref.6. This is an advantage over othertreatments, where the continuous form is implemented add–hoc on the basis of numericalconvenience (see DL-POLY or Gromacs reference manuals).In order to complete the formulation, we further need to supplement our model with aself term that gauges the total energy, accounting for the difference between the actual andmodified potentials:[33] U self i = lim r → q i ( u ( r ) − φ mod ( r )) (28)Using Eq. (27), we obtain the general expression for the self term of the generalized reactionfield method. U self i = − q i (cid:18) φ mod ( R c ) − φ ′ mod ( R c ) R c (cid:19) (29)Eq. (27) and Eq. (29) are the more important results of this section. For the special casewhere φ ( r ) = 1 /r , we recover exactly the reaction field model of Hummer et al.[6] Here, weshow that the reaction field form is kept generally for whatever Green function φ mod ( r ). Asa result, it may be exploited also to accelerate the convergence of damped coulomb potentialsin Ewald type summations. 11 airwise schemes for the calculation of the real space Ewald summation We finish this section with the pair potential and self terms that result for probably themost significant damped coulomb potential, namely, φ mod ( r ) = e rf c ( αr ) /r . Using Eq. (27)and Eq. (29), we obtain:[36] u ercf ( r ij ) = q i q j ( e rf c ( αr ij ) r ij − e rf c ( αR c ) R c + 12 " e rf c ( αR c ) R c + 2 α √ π e − α R c R c r ij − R c R c ) (30)While for the self term, we have: U self i = − q i (cid:26) e rf c ( αR c ) R c + 12 (cid:18) e rf c ( αR c ) R c + 2 α √ π e − α R c (cid:19)(cid:27) (31)Notice that this self energy corrects for the truncation of the exact real space lattice sum-mation of φ mod ( r ) = e rf c ( αr ) /r , Eq. (3). i.e., it cannot possibly account for the Fourierspace term of Eq. (4). If the purpose is to approximate the full Coulomb sum, Eq. (1), and α is chosen such that the reciprocal space sum may be ignored, one must still account forthe self term of the Fourier sum, Eq. (5).Actually, this point can be illustrated from the formalism afforded by Eq. 21, 27, 28. Tosee this, consider the situation where this approach is applied to approximate the infinitesum Eq. (1) of coulomb potentials, 1 /r , instead of the related e rf c ( r ) /r sum of Eq. (3).Rather than employing Eq. (29) for the self term, we would then need to consider the selfterm as the limiting value of − q i ( u ( r ) − /r ), since now the excess energy is over φ = 1 /r ,rather than over φ mod = e rf c ( r ) /r . The self term would still be as Eq. (31), plus an extraterm lim r → − q i e rf ( r ) /r . One can recognize here exactly the self term of theFourier contribution, Eq. (5). Therefore, the theoretical treatment provides naturally thegeneralized reaction field result of Eq. (30) and Eq. (31), plus the Ewald sum self term asthe best possible approximation to Eq. (1) that can be obtained by performing a truncatedsum of e rf c ( r ) /r terms.The general results Eq. (27), Eq. (29), as well as Eq. (30)–Eq. (31) for the dampedCoulomb potential are exactly as the Zero Charge-Zero Dipole method derived recently byFukuda et al. using a heuristic approach based on the mirror image technique.[19] Ourresults provide an alternative derivation showing clearly the strong connection with theReaction Field technique.[6] As such, we will henceforth refer to Eq. (30)–(31) as DampedReaction Field (DRF) Method. 12he results Eq. (30)–(31) resemble effective pair potentials that have been suggestedrecently.[8, 11, 18, 19, 21] Indeed, they include a shift in the potential that is at the heart ofWolf’s method, but provide also an additional force as in the Fennel and Gezelter methodand related approaches.[11, 19, 21]In fact, all of these methods may be considered as a wide class, that, by use of Eq. 21,27, 28, allows to approximate the real space Ewald summation of Eq. (3). This is achievedusing an effective pairwise potential of the form: u ( r ij ) = q i q j f ( r ij ; α, R c ) r ij ≤ R c r ij > R c (32)together with a self term: U self i = − q i g ( α, R c ) (33)Table I provides a summary of different pairwise methods following the scheme above (Ref.21also discusses the relation amongst different pairwise potential schemes). Recall that for thesake of approximating the full Coulomb sum, Eq. (1), the Fourier contribution Eq. (4) stillneeds to be evaluated. Under favorable cases, that may be achieved by neglecting all of thereciprocal space sum, and accounting only for the remaining self contribution, Eq. (5)).However, the theoretical approach explained here shows that in practice all such meth-ods may be also employed to accelerate the convergence of the real space summation alone,whether one opts to ignore the reciprocal space sum or not. Compared to other method-ologies, however, our results have several advantages: 1) Both the potential and the forceremain continuous at r = R c ; 2) The force and potential are fully consistent with each other,and are obtained in a straight forward manner (c.f. Ref. 8 and 18) 3) The continuity of theexpressions is not merely plugged in for numerical convenience, but results from a clear andwell understood electrodynamic treatment (c.f. Ref. 11). 4) The self term is provided andfollows systematically from the theoretical treatment (c.f. Ref. 11). RESULTS
In the following section, we will test the ability of different real–space lattice summationsas a means to approximate the full electrostatic sum of Eq. (1). Whereas such methods areexpected to provide a good convergence for dense fluids of dipolar molecules,[11] we have13 air potential f ( r ; α, R c ) g ( α, R c ) Limit α → u Ew φ ( r ) - Coulomb u Wolf φ ( r ) − φ ( R c ) φ ( R c ) Shifted Coulomb u WFG φ ( r ) − φ ( R c ) − φ ′ ( R c )( r − R c ) φ ( R c ) − φ ′ ( R c ) R c Force shifted Coulomb u DRF φ ( r ) − φ ( R c ) − φ ′ ( R c ) r − R c R c φ ( R c ) − φ ′ ( R c ) R c Reaction FieldTABLE I. List of pair potentials and self terms used to approximate the real space Ewald sum,Eq. (3), using the general scheme of Eq. (32)-Eq. (33), with φ ( r ) = e rf c ( αr ) /r . The last columnindicates the relation of each model to the Coulomb potential in the limit where α →
0. The entry g ( α, R c ) for u WFG is the self term for the Wolf-Fennell-Gezelter method, not provided in the originalreference.[11] For the case where Eq. (32)-33 are employed to approximate the full Coulomb sum,Eq. (1), an extra self term, 2 α √ π accounting approximately for the Fourier contribution Eq. (4)must be added to g ( α, R c ). chosen to study crystalline and molten ionic salts. This should provide a more stringenttest than merely polar fluids and therefore allow us to obtain conclusions of more generalvalidity.As a test of an ordered ionic solid, we will consider crystalline sodium chloride. Accord-ingly, we will assume an ordered arrangement of + and - unit charges with lattice spacing a , as in the crystal rock salt lattice.[37] For an ionic molten compound, we again considerconfigurations of + and - unit charges thermally sampled from a screened Yukawa fluid withhard sphere diameter σ .The convergence of the average Coulomb energy felt by an ion is monitored by calculating: U R i ( R c ) = ∗ X j u ( r ij ) + U self i (34)with u ( r ) and U self i as given by Eq. (32)-33, and the corresponding choice of f ( r ; α, R c ) and g ( α, R c ) as indicated in Table I for different approximate schemes. Charge neutralization schemes with no damping
In order to make a transparent comparison between DRF and Wolf’s charge neutralizationschemes, let us first consider the case of zero damping, α = 0 so that φ adopts the bare14 IG. 2. Electrostatic potential energy of the NaCl lattice calculated as the direct pairwise sum-mation and by means of Wolf’s and DRF methods. R c / a -15-10-50510 E n e r gy / i on [ q ² / a ] Direct pairwise sumWolf, α = 0DRF, α = 0 Coulomb form 1 /r . In this case, Wolf’s charge neutralization scheme becomes a mere shiftedCoulomb potential, while DRF becomes the Reaction Field Method.Figures 2 and 3 show U i ( R c ) as estimated by bare Coulomb summation, Wolf’s methodand DRF for solid NaCl and the 1:1 molten salt.As expected, in the NaCl lattice (Fig. 2), the direct summation is not convergent, andthe data is so scattered that it is not even possible to guess approximately the exact energy.The introduction of the neutralizing shell proposed by the Wolf’s method leads the energyto converge, and the amplitude of the oscillations decreases remarkably, even though thedamping parameter is zero. However, if the DRF method is used instead, the amplitude ofthe oscillations decreases much more, to the point that they are hardly visible in the scaleof the figure (at least beyond R c ≈ . a ). We have checked that this behavior also occursfor other crystalline structures, such as blende and CsCl.The improved performance of charge neutralization schemes is also seen in the case of amolten ionic compound (Fig. 3). The only difference here is that the direct summationdoes converge, since, due to the lack of long-range order, the charges are effectively shielded.The convergence is, however, very slow. Again, the Wolf’s and DRF methods acceleratesignificantly the convergence, and DRF gives rise to somewhat smaller oscillations thanWolf’s method. 15 IG. 3. Electrostatic potential energy of a molten ionic compound of stoichiometry 1:1 calculatedas the direct pairwise summation and by means of Wolf’s and DRF methods. R c / σ -1-0.8-0.6-0.4 E n e r gy / i on [ q ² / a ] Direct pairwise sumWolf, α = 0DRF, α = 0 Role of damping
The convergence of U i ( R c ) may be much improved by using a damping function φ mod ( r ) =e rf c ( αr ) /r . As noted above, choosing α = 0 transforms φ mod into the bare Coulomb po-tential. For finite α , φ mod becomes a damped Coulomb potential with a decay rate that isgoverned by α .It is important to notice, however, that once α = 0, then immediately the Fourier con-tribution to the Coulomb sum becomes finite, so that even in the limit R c → ∞ , the sum E i ( R c ) is only an approximation for the full Coulomb sum.As noted by Wolf, however, there is a wide choice of finite α where the only significantcontribution to the Fourier sum, Eq. (4), is the self term, Eq. (5). For the NaCl crystal,for example, choosing α = 1 . a − the reciprocal space sum is ≈ · − times smallerthan the self term. The case of NaCl is however a particularly favorable one, and mightnot be always taken for granted. For another simple crystal structure such as CsCl,[38] forexample, the reciprocal space sum is now only ≈ − times smaller than the self term. Astudy of the structure factors reveals that CsCl is a very unfavorable case because here ρ ( k )vanishes only once every two k vectors, while NaCl is probably particularly favorable, since ρ ( k ) vanishes for all but about n+1 every (2 n + 1) k vectors.For a fluid, the relative contribution of the sum cannot generally be determined, but one16ere expects that ρ ( k ) is of finite range, and should therefore decay faster than for the caseof CsCl.Be as it may, one expects there to be a range of sufficiently small α , where, on the onehand, the Fourier contribution is given essentially by the self term, while simultaneously,the convergence of the real sum is much improved. One then hopes that the reciprocal spacesum can be ignored altogether, leaving all of the energy contribution in terms of a relatively fast decaying damped coulomb potential.We test this hypothesis first for a perfect NaCl crystal, then for the 1:1 molten salt. Perfect NaCl crystal
Figures 4 and 5 show the evolution of the real part of the energy of the NaCl crystalwith the cut-off radius for finite α (1.0 a − and 2.0 a − respectively). Results are displayedfor the real space Ewald summation (RSE), Wolf’s method (WM), Wolf-Fennell-Gezelter’smethod (WFG) and DRF.Notice that for finite α , even the real part of the Ewald summation is convergent. How-ever, the oscillations are rather large and exhibit discontinuities, while the convergenceremains relatively slow. The use of Wolf’s method makes the oscillations decrease strongly.The DRF method shows a clear diminution of the oscillations of the energy, not only com-pared with RSE, but also with Wolf’s method, while the performance of WFG is aboutthe same as that afforded by DRF. Further increasing α from 1 . a − to 2 . a − producesa spectacular improvement on the convergence of U i ( R c ). This is obvious right away bycomparing both the y and x scales of figures 4 and 5. i.e., not only the asymptotic valueof U i ( R c ) is reached much faster, but also the amplitude of the oscillations is decreased bymore than an order of magnitude. Obviously, this is at the expense of making the Fouriercontribution larger, and more importantly, increasing the relevance of the reciprocal spacesum. In fact, for α = 2 . a − the reciprocal space sum is now already about 10 − times theself term of Eq. (5), so that it becomes unsafe to approximate the Coulomb sum withouttaking account of the full Fourier contribution.An interesting feature of both WFG and DRF is that, not only is the amplitude of theoscillations smaller, but also the discontinuities that were apparent in RSE and to a smallerextent in WM seem to be considerably smoothed. This smoothing property must be related17 IG. 4. Evolution of the real part of the energy with the cut-off radius. NaCl perfect crystal. α = 1 . a − . R c / a -3.6-3.4-3.2-3 U i R ( R c ) / [ q ² / a ] EwaldWolfWolf-Fennell-GezelterDamped Reaction Field to the presence of a shifted force contribution to the effective pair potential, since only WFGand DRF share this feature among the four methods tested. In fact, WFG and DRF performsimilarly, so that it would seem it is the the shifted force contribution what makes thesemethods perform better than RSE and WM.At this point it is convenient to mention the significance of the self term, Eq. (33).Indeed, at first thought one might consider that a constant contribution U self i merely shiftsthe total energy, but not the underlying dynamics so that it may be completely ignored.However, in order to compare the results of U i for different methods it is essential to accountfor the appropriate self contribution as indicated in Table I. Actually, even for a givenmethod, results with different R c (or α ) can only be compared when U self i is included. Thisis particularly relevant for the case of WFG, which does not include a self term in the originalreference, and could only be compared with the RSE, WM and DRF by including the selfterm of Tab.I. This can also be of great importance in the simulation of open systems, as isthe case of studies performed in the grand canonical ensemble.As a final remark, we note the conclusions drawn from the analysis performed on NaClalso hold for other simple crystal structures such as CsCl and ZnS (not shown).18 IG. 5. Evolution of the real part of the energy with the cut-off radius. NaCl perfect crystal. α = 2 . a − . R c / a -3.51-3.5-3.49-3.48 U i R ( R c ) / [ q ² / a ] EwaldWolfWolf-Fennell-GezelterDamped Reaction Field
EwaldWolfWolf-Fennell-GezelterDamped Reaction Field
FIG. 6. Evolution of the real part of the energy with the cut-off radius. Molten ionic compoundof stoichiometry 1:1. α = 0 . σ − . R c / σ -0.85-0.8-0.75-0.7-0.65-0.6 U i R / [ q ² / σ ] EwaldWolfWolf-Fennel-GezelterDRF
Molten salt
Let us now compare the performance of the pairwise summation schemes for the 1:1molten salt.Depicted in Figures 6 and 7 is the behavior of U i ( R c ) now calculated for a molten ionic19 IG. 7. Evolution of the real part of the energy with the cut-off radius. Molten ionic compoundof stoichiometry 1:1. α = 1 . σ − . R c / σ -0.72-0.71-0.7-0.69 U i R / [ q ² / σ ] EwaldWolfWolf-Fennel-GezelterDRF compound. In a molten system, the oscillations of the energy are smoother than in a crystal,irrespective of the method used. Since the system is disordered, the interaction with newions as the cut-off radius increases does not take place at discrete distances, but rather,continuously.Nevertheless, the oscillations given by RSE are still too high, and the alternative pairwiseschemes clearly improve this situation. In the figures, it is clear that Wolf’s method alreadydiminishes widely the oscillations, and DRF makes this improvement somewhat more sig-nificant.On the contrary, the WFG method now seems to converge slower, and more importantly,does not seem to oscillate but rather approach the asymptotic energy from below.
Calculation of the virial
In the preceding sections we have studied the convergence of the energy. Another impor-tant issue refers to the convergence of the forces. Since, however, the net force exerted onan ion in a perfect crystal is zero, we rather consider the virial per ion, which has similarconvergence issues as the energy. Accordingly, we study the convergence of a virial function,20efined as: W Ri ( R c ) = ∗ X j w ( r ij ) (35)where w ij = −→ f ij · −→ r ij (36)while the forces are given by differentiation of the corresponding pairwise potentials (c.f.Tab.I).Figure 8 shows the evolution of the virial of an ion inside de NaCl lattice with the cut-offradius for a fixed value of α = 1 .
0. Results are shown for RSE, WFG, and DRF. Noticethat we do not employ Wolf’s method here, as it gives a discontinuity in the force at R c .The figure explicitly shows a rather poor convergence of the RSE, which exhibits strongoscillations and discontinuities as a function of R c . At this point, it is worth mentioning thatfor α = 0, the RSE method becomes the bare Coulomb sum and does also not converge at all(not shown). Using the effective pairwise schemes very much improves this situation. Indeed,the amplitude and also the frequency of the oscillations is reduced, while the convergence isalso achieved faster. Apparently, WFG and DRF perform similarly, lending support to theidea that a shifted force is necessary (and perhaps sufficient) to improve the convergence ofthe sum. This is somewhat confusing, given that one can hardly attribute any electrostaticsignificance to the WFG scheme, which merely corresponds to a force shifted potential.Figure 9 displays results for W Ri ( R c ) in the molten salt system, with α = 0 . σ − . Inthis case, the direct RSE summation displays deviations of about the same amplitude thanWFG and DRF, but clearly exhibits stronger discontinuities. Again both WFG and DRFperform similarly. CONVERGENCE AND OPTIMIZATION
In the previous sections we have shown analysis similar to those of Ref.8 and 11, whichindicate that for sufficiently small α , one can devise effective pairwise potentials that allowto approximate the Coulomb sum, Eq. (1), at low computational cost. Furthermore, wehave shown that the damped reaction field method converges at least as good as WM andWFG methods, but has a sounder physical interpretation.However, our theoretical analysis reveals that DRF, as well as WM and WFG, are in21 IG. 8. Evolution of the virial of an ion with the cut-off radius. NaCl perfect crystal. α = 1 . a − . R c / a -5-4-3-2-1 W i R [ q ² / a ] EwaldWolf-Fennell-GezelterDRF
FIG. 9. Evolution of the total virial with the cut-off radius. Molten ionic compound with stoi-chiometry 1:1. α = 0 . σ − . R c / σ -4-3-2-1 W i R [ q ² / σ ] EwaldWolf-Fennell-GezelterDRF fact plausible approximations for the real space Ewald sum, Eq. (3), whether one opts toneglect the reciprocal space sum (Eq. (4)) or not. This is a very relevant issue which hasapparently not been considered previously. It suggest one could perform the full Ewaldsummation, accounting both for the real space and Fourier contributions (whether in thestandard implementation or in the particle mesh approaches 13), but using DRF in order toaccelerate the convergence of the real space term. This would allow to perform the Coulomb22
IG. 10. Deviation of the energy from its convergence value (Rc = ∞ ) as a function of α calculatedfor several cut-off radii. Comparison between Ewald summation and DRF method. ∆ U Fi ( α ; n c = 0)is shown as a black dotted-dashed line in both plots. NaCl perfect crystal. α / a -1 | U i R ( α ; R c ) - U i R ( α ; R c = ∞ ) | / [ q ² / a ] Real Ewald Rc = 1.0Real Ewald Rc = 1.5Real Ewald Rc = 2.0Real Ewald Rc = 2.5Real DRF Rc = 1.0Real DRF Rc = 1.5Real DRF Rc = 2.0Real DRF Rc = 2.5Fourier nc = 0
FIG. 11. Deviation of the energy from its convergence value (Rc = ∞ ) as a function of α calculatedfor several cut-off radii. Comparison between Ewald summation and DRF method. ∆ U Fi ( α ; n c = 0)is shown as a black dotted-dashed line in both plots. Molten ionic compound with stoichiometry1:1. α / σ -1 | U i R ( α ; R c ) - U i R ( α ; R c = ∞ ) | / [ q ² / σ ] Real Ewald Rc = 2.0Real Ewald Rc = 3.0Real Ewald Rc = 4.0Real Ewald Rc = 5.0Real DRF Rc = 2.0Real DRF Rc = 3.0Real DRF Rc = 4.0Real DRF Rc = 5.0Fourier nc = 0 U R i ( α ; R c ) = | U R i ( α ; R c ) − U R i ( α ; R c = ∞ ) | (37)A plot of the real space convergence error is shown in Fig.10 for the case of crystallineNaCl. Full lines indicate the error that results when Eq. (3) is evaluated by sum of plaine rf c ( αr ) /r contributions, i.e., RSE, and dashed lines the case where it is evaluated usingDRF. The figure clearly shows that, for reasonable choices of R c , there is a range of small α , where DRF exhibits much smaller convergence error than RSE. However, as α increases,the difference becomes smaller, and actually exhibits a crossover to a regime where the plainRSE seems to perform better.At any rate, the convergence is always improved for large values of α . However, one cannot increase α arbitrarily, because then the Fourier contribution would become relevant. Inorder to optimize the choice of α for a given fixed R c , the relevant issue is then what theerror of neglecting the Fourier sum is.In order to asses this, we now introduce a measure of the truncation error performed inthe reciprocal space sum of Eq. (4). This error results from the neglect of contributions withreciprocal space vectors k = πL ( n x , n y , n z ) larger than a prescribed cutoff, k c . Similarly toEq. (37), we can therefore introduce a Fourier space truncation error as:∆ U Fi ( α ; n c ) = | U Fi ( α ; n c ) − U Fi ( α ; n c = ∞ ) | (38)where n c defines a cubic cutoff for vectors whose integer components (absolute value) arelarger than n c .Obviously, the pairwise effective potentials as estimated in WM and WFG, as well as inthe zero dipole method of Ref.[19], neglect the Fourier space contribution all together. Thiscorresponds to assuming a Fourier space convergence error for zero cutoff, n c = 0.Fig.10 displays together with ∆ U R i ( α ; R c ) the Fourier space convergence error ∆ U F i ( α ; n c )for the special case where n c = 0.Clearly, it is always possible to choose a value of α sufficiently small that ∆ U F i ( α ; n c = 0) isnegligible. However, the optimal choice of α is one where both ∆ U R i ( α ; R c ) and ∆ U F i ( α ; n c )are of similar order. Obviously, it does not make much sense to achieve an exceptional24onvergence in the real space contribution, if this is at the cost of having a much largerreciprocal space error, because what matters is the overall convergence of the full Coulombsum.Unfortunately, Fig.10 seems to suggest that ∆ U F i ( α ; n c = 0) cuts ∆ U R i ( α ; R c ) at a pointwhere there is no significant advantage of the DRF method over the plain RSE summation.i.e., the point at which ∆ U F i ( α ; n c = 0) is about the same order of magnitude of ∆ U R i ( α ; R c )occurs where DRF and RSE perform similarly. This seems clearly visible, at least for theperfect NaCl crystal and three choices of R c /a = 1 , . α where it exhibits smaller error than the plain RSE is smaller. Infact, in this case ∆ U F i ( α ; n c = 0) cuts the curves ∆ U R i ( α ; R c ) after the crossover where RSEconverges better than DRF, for all but the smallest R c studied.We now study the optimization of the Ewald sum by relaxing the constraint of pairwiseschemes, i.e., by allowing for a finite reciprocal space cutoff. This is still a relevant issue,since the Ewald sum requires to share the computational cost between both the real andFourier contributions.[17] An improvement in the convergence of the real space summationwould allow to shift some of the cost of the reciprocal space sum to the real space term,making the former less expensive.Fig.12-13 display ∆ U R i ( α ; R c ) as estimated from RSE and DRF, compared with ∆ U F i ( α ; n c )for crystalline and molten NaCl, respectively. Notice that the errors are shown in logarithmicscale, with large negative values indicating small errors.From inspection of ∆ U R i ( α ; R c ), it now becomes clear that one can identify two differentregimes. Firstly, a reaction field regime corresponding to small values of the product αR c ,where DRF produces smaller real space errors than the plain Ewald real space sum. Secondly,an Ewald regime , of large αR c , where it is actually the RSE summation which performs betterthan DRF, and presumably, better than whichever pairwise effective potential scheme.Optimization of the Ewald sum requires to choose R c , α and n c such that ∆ U R i ( α ; R c ) ≈ ∆ U F i ( α ; n c ). Thus, a mixed scheme using DRF and reciprocal space sums would be advan-tageous over the Ewald method if the condition above is met at a value of α where theconvergence error of DRF is smaller than that of RSE. Unfortunately, our plots seem toindicate that this is not usually the case. Only for the special choice of n c = 0, i.e., for25 IG. 12. Quadratic deviation of the real and the reciprocal part of the energy from its convergencevalue (infinite Rc or nc) as a function of α calculated for several cut-off radii. Comparison betweenEwald summation and DRF method. NaCl perfect crystal. α / a -1 -25-20-15-10-50 l og { [ U i R ( α ; R c ) - U i R ( α ; R c = ∞ )] ² } Fourier nc = 0Fourier nc = 1Fourier nc = 3Fourier nc = 5Fourier nc = 7Real Ewald Rc = 1.0Real Ewald Rc = 1.5Real Ewald Rc = 2.0Real Ewald Rc = 2.5Real DRF Rc = 1.0Real DRF Rc = 1.5Real DRF Rc = 2.0Real DRF Rc = 2.5 α / a -1 -25-20-15-10-50 l og { [ U i F ( α ; n c ) - U i F ( α ; n c = ∞ )] ² } Fourier nc = 0Fourier nc = 1Fourier nc = 3Fourier nc = 5Fourier nc = 7Real Ewald Rc = 1.0Real Ewald Rc = 1.5Real Ewald Rc = 2.0Real Ewald Rc = 2.5Real DRF Rc = 1.0Real DRF Rc = 1.5Real DRF Rc = 2.0Real DRF Rc = 2.5 complete neglect of the reciprocal space sum, we find that an optimized DRF yields errorssimilar to the optimized RSE.Notice that the conclusions drawn above are restricted to bulk systems. Additionalcare must be taken for inhomogeneous systems, where the charge distribution may exhibitlarge or even diverging correlations of long wavelength. The charge structure factor ρ ( k )in Eq. (4) then becomes very large for small k , and the reciprocal space sum is significanteven for small α . Empirically, this has been observed by Takahashi et al., who studied theliquid-vapor interface of SPCE/E water and noticed the convergence of Wolf’s method couldonly be achieved for cutoffs as large as 20 molecular diameters.[39] Such a poor convergenceindicates very long range, small wave-vector interactions across the interface. These arebest dealt by transferring some of the real space computational burden to the reciprocalspace sum, which precisely is devised to ensure fast convergence of small wave-vector con-tributions. Theoretically, it has been shown that properly taking into account a large scaleinhomogeneity requires to account for the net dipole moment of the system in the directionperpendicular to the interface. Such contribution, which is significant and independent ofthe choice of α , corresponds to a zero wave-vector term of the reciprocal space sum andcannot possibly be accounted for with a real space finite cutoff.[4, 40–42]26 IG. 13. Quadratic deviation of the real and the reciprocal part of the energy from its convergencevalue (infinite Rc or nc) as a function of α calculated for several cut-off radii. Comparison betweenEwald summation and DRF method. Molten ionic compound with stoichiometry 1:1. α / σ -1 -14-12-10-8-6-4-20 l og { [ U i R ( α ; R c ) - U i R ( α ; R c = ∞ )] } Fourier nc = 0Fourier nc = 1Fourier nc = 5Fourier nc = 9Fourier nc = 13Real Ewald Rc = 2.0Real Ewald Rc = 3.0Real Ewald Rc = 4.0Real Ewald Rc = 5.0Real DRF Rc = 2.0Real DRF Rc = 3.0Real DRF Rc = 4.0Real DRF Rc = 5.0 α / σ -1 -14-12-10-8-6-4-20 l og { [ U i F ( α ; n c ) - U i F ( α ; n c = ∞ )] } Fourier nc = 0Fourier nc = 1Fourier nc = 5Fourier nc = 9Fourier nc = 13Real Ewald Rc = 2.0Real Ewald Rc = 3.0Real Ewald Rc = 4.0Real Ewald Rc = 5.0Real DRF Rc = 2.0Real DRF Rc = 3.0Real DRF Rc = 4.0Real DRF Rc = 5.0
CONCLUSIONS
In this work we have considered a number of recent methodologies which allow to approxi-mate the electrostatic energy of charged systems by means of effective pairwise potentials.[8,11, 19]Our theoretical study reveals that these methods may be actually considered as approxi-mations to the real space Ewald sum only, Eq. (3), rather than to the sought Coulomb sum,Eq. (1). This is an important matter, because it means that the pairwise effective poten-tials cannot possibly account for the long range electrostatic response of charged systems.Particularly, such methods cannot account ab-initio for the surface term of the Coulombenergy, which is recovered in the full Ewald sum as the k = 0 limit to the reciprocal spaceterm.[4, 16]For cases where the damping parameter α is small enough, however, all terms besidesthe k = 0 contribution to the reciprocal space sum are negligible and the Ewald self term,Eq. (5) is then a very good approximation for the Fourier space contribution, Eq. (4).In such cases, our numerical study shows that pairwise effective potentials produce in-deed a better convergence of the real space Ewald sum than the Ewald pair potential itself.In practice, this may be achieved by introducing a shifted,[8] or a shifted force,[11] poten-27ial. However, our theoretical analysis shows that the correct electrostatic treatment mayachieve both the potential and force shift by means of a generalization of the reaction fieldmethod,[19] which damps the Coulomb potential via an e rf c ( αr ) function. This dampedreaction field method, used as an approximation to the real space Ewald sum smoothlytransforms from a pure reaction field to an effective Ewald sum, by tuning the dampingparameter α .We have studied the performance of the modified Coulomb sums in a particularly difficultcase of an ionic crystal and its melt. This provides a stringent test of all methods and thusallows to draw more general conclusions. Our results show that the new damped reactionfield method, which has a much more clear electrostatic significance, is also the methodwhich provides better convergence of the energy and the virial for both the crystal and themelt.A close inspection of the convergence errors allows to identify two different regimes.Firstly, a reaction field regime, corresponds to small αR c , where effective pair wise methodssuch as the damped reaction field method converge better than the bare e rf c ( αr ) /r poten-tial. Secondly, an Ewald regime, corresponding to large αR c , where the bare Ewald potentialactually converges better than all effective pair wise methods studied. Our numerical studyshows that the crossover from the reaction field to the Ewald regime occurs precisely forthose values of α where the reciprocal space sum has an error about the same size as the realspace sum. Since this equality is the condition for optimization of the Coulomb sum,[17]the implication is that for an optimized calculation of the electrostatic energy the Ewaldsum method is not improved by effective pair wise methods. Such methods do remain usefuland produce improved convergence of the Ewald sum in those cases where high accuracycalculations, or optimization of the Coulomb sum are not a concern. Alternatively, theymight also be competitive in cases where the availability of parallel computing does notmake competitive the calculation of the reciprocal space sum.This work was supported by Ministerio de Educacion y Ciencia through project FIS2010-22047-C05-05. LGM would like to acknowledge helpful discussions with Sabinne Klapp. ∗ [email protected][1] M. Allen and D. Tildesley, Computer Simulation of Liquids (Clarendon Press, Oxford, 1987).
2] D. Frenkel and B. Smit,
Understanding Molecular Simulation , 2nd ed. (Academic Press, SanDiego, 2002).[3] J. A. Barker and R. O. Watts, Mol. Phys. , 789 (1973).[4] S. W. de Leeuw, J. W. Perram, and E. R. Smith, Proc. R. Soc. Lond. A , 27 (1980).[5] J. Lekner, Physica. A , 826 (1989).[6] D. Hummer, D. M. Soumpasis, and M. Neumann, Mol. Phys. , 769 (1992).[7] T. A. Darden, D. York, and L. Pedersen, J. Chem. Phys. , 10089 (1993).[8] D. Wolf, P. Keblinski, S. Phillpot, and J. Eggebrecht, J. Chem. Phys. , 8254 (1999).[9] S. Tyagi, J. Chem. Phys. , 014101 (2005).[10] Y. Shan, J. L. Klepeis, M. P. Eastwood, R. O. D. r, and D. E. Shaw, J. Chem. Phys. ,054101 (2005).[11] C. J. Fennell and J. D. Gezelter, J. Chem. Phys. , 234104 (2006).[12] K. Esselink, Comp. Phys. Comm. , 375 (1995).[13] A. Y. Toukmaji and J. A. Board, Comp. Phys. Comm. , 73 (1996).[14] G. Hummer, L. R. Pratt, A. E. Garcia, and M. Neumann, AIP Conf. Proc. , 84 (1999).[15] E. R. Smith, Proc. R. Soc. Lond. A , 475 (1981).[16] V. Ballenegger, J. Chem. Phys. , 161102 (2014).[17] J. Kolafa and J. W. Perram, Mol. Sim. , 351 (1992).[18] D. Zahn, B. Schilling, and S. M. Kast, J. Phys. Chem. B , 10725 (2002).[19] I. Fukuda, Y. Yonezawa, and H. Nakamura, J. Chem. Phys. , 164107 (2011).[20] J. S. Hansen, T. B. Schrder, and J. C. Dyre, The Journal of Physical Chemistry B , 5738(2012).[21] I. Fukuda, J. Chem. Phys. , 174107 (2013).[22] F. N. Mendoza, J. Lopez-Lemus, G. A. Chapela, and J. Alejandre, J. Chem. Phys. ,024706 (2008).[23] J. Muscatello and F. Bresme, J. Chem. Phys. , 234111 (2011).[24] M. Neumann, Mol. Phys. , 841 (1983).[25] J. A. Barker, Mol. Phys. , 1057 (1994).[26] J. M. Miguez, D. Gonzalez-Salgado, J. L. Legido, and M. M. Pineiro, The Journal of ChemicalPhysics , 184102 (2010).[27] N. Kamiya, I. Fukuda, and H. Nakamura, Chem. Phys. Lett. , 26 (2013).
28] J. A. Barker and R. O. Watts, Chem. Phys. Lett. , 144 (1969).[29] D. Wolf, Phys. Rev. Lett. , 3315 (1992).[30] M. Neumann, J. Chem. Phys. , 5663 (1985).[31] I. G. Tironi, R. Sperb, P. E. Smith, and W. F. van G unsteren, J. Chem. Phys. , 5451(1995).[32] B. Garzon, S. Lago, and C. Vega, Chem. Phys. Lett. , 366 (1994).[33] S. Boresch and O. Steinhauser, J. Chem. Phys. , 10780 (2001).[34] M. Neumann, O. Steinhauser, and G. S. Pawley, Mol. Phys. , 97 (1984).[35] L. G. MacDowell and C. Vega, J. Phys. Chem. B , 6089 (2010).[36] V. H. Elvira, Evaluaci´on de interacciones de largo alcance en la interfase el´ectrica , Master’sthesis, Universidad Complutense de Madrid (2012).[37] This is an interpenetrated face centered cubic lattice. Plus charges occupy positions(0,0,0),(0,1/2,1/2),(1/2,0,1/2),(1/2,1/2,0), while minus charges occupy analogous positionsthat are displaced by (1/2,0,0).[38] Two interpenetrated cubic primitive lattices, with plus charges occupying (0,0,0) and minuscharges occupying (1/2,1/2,1/2).[39] K. Z. Takahashi, T. Narumi, and K. Yasuoka, J. Chem. Phys. , 174112 (2011).[40] J. Hautman, J. W. Halley, and Y. Rhee, J. Chem. Phys. , 467 (1989).[41] E. Spohr, J. Chem. Phys. , 6342 (1997).[42] I. Yeh and M. L. Berkowitz, J. Chem. Phys. , 3155 (1999)., 3155 (1999).