An Extension of Wolf's Method for the Treatment of Electrostatic

Jan 22, 2015 - An Extension of Wolf's Method for the Treatment of Electrostatic Interactions: Application to Liquid Water and Aqueous Solutions. Georg...
0 downloads 12 Views 1MB Size
Download PDF
Article pubs.acs.org/JPCB

An Extension of Wolf’s Method for the Treatment of Electrostatic Interactions: Application to Liquid Water and Aqueous Solutions George S. Fanourgakis*,† Computation-based Science and Technology Research Center (CaSToRC), The Cyprus Institute, P.O. Box 27456, 1645 Nicosia, Cyprus ABSTRACT: We report the results of molecular dynamics simulations, which indicate that, by carefully splitting the electrostatic interactions in short- and long-range contributions and employing the charge-neutralization method of Wolf, accurate predictions of various properties of liquid water and aqueous solutions can be achieved without the need for the Ewald summation. In order to assess the accuracy of the proposed approach, several molecular dynamics simulations under different thermodynamic conditions are performed, employing various rigid, flexible, pairwise additive, and many-body polarizable water models. The predictions of the new approach are compared to the benchmark results obtained with the Ewald summation. It is found that while in the new approach there are no adjustable parameters, such as a damping parameter, the obtained results are more accurate than the results of similar approaches that are based on the Wolf method, while at the same time less or no additional computational effort is required. It is also concluded that the error of the results is smaller or at least comparable to the statistical error of a typical molecular dynamics simulation.

1. INTRODUCTION Electrostatic and dispersion van der Waals (vdw) interactions are among the most common types of interactions between molecules. These type of interactions are found in almost any empirical force field that is currently employed for the description of a molecular system of chemical or biological interest. Both of them are usually modeled as pairwise additive ones, and in order to compute the total potential energy of a system, the energy between all pairs of atoms is summed up. However, while the van der Waals energy decays fast with the interatomic separation (as r−6 with r being the distance between two atoms) the electrostatic one decays much slower (as r−1). Therefore, in a system of infinite size in which, for practical reasons, the interactions of each atom should be restricted with atoms up to a certain distance (cutoff), the two types of interactions exhibit different behavior and thus require completely different treatment. During a molecular simulation, and for a reasonable choice of the cutoff, it is usually accurate enough to truncate the vdw interactions at the cutoff and to apply an analytically computed correction to the energy and pressure.1 It is not possible, however, to apply the same technique to the electrostatic interactions, as it can be shown that a similar correction in this case has an infinite value. Also, completely ignoring interactions beyond the cutoff during a molecular simulation usually leads to dramatic fluctuations of the obtained energy and therefore to inaccurate or nonphysical results. The Ewald summation2 is considered today the most accurate method for the calculation of the electrostatic energy in a periodic system. The main idea is to decompose the electrostatic potential into a short- and a long-range component. The first one is evaluated in real space, while the latter in Fourier (reciprocal) © XXXX American Chemical Society

space. The separation of the initial potential into the two components is usually controlled by a parameter α, which, although arbitrary, governs the relative rate of convergence of the two components. By choosing an appropriate α, both components converge rapidly in their respective spaces. Approaches such as the particle-mesh Ewald (PME)3 and the particle−particle particle−mesh (P3M)4 have been introduced to reduce the computational cost of the reciprocal space calculation. An excellent review of several of the most modern exact and approximate approaches for the treatment of the electrostatic interactions can be found in ref 5. Approximate methods, such as the spherical-cutoff methods, have been extensively employed for the treatment of the electrostatic interactions in macromolecular systems and have been implemented in major packages such as CHARMM6 since the early days of molecular simulations. Several of these approaches have been discussed and evaluated.7 Wolf et al.,8 based on the observation that in a condensed ionic system the net Coulomb potential is rather short ranged,9 had some interesting insights, while an alternative approach for the calculation of the electrostatic energy was proposed. A key point is that when the direct summation of the spherically truncated r−1 pair potential is used for the energy estimation, problems are encountered since the system summed up is practically never neutral. In order to overcome this difficulty, Wolf et al. proposed a simple and elegant approach to neutralize the system. It was shown that charge neutralization can be achieved by shifting the original Received: October 22, 2014 Revised: January 6, 2015

A

DOI: 10.1021/jp510612w J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B Coulomb potential by 1/Rc (where Rc is the cutoff), making in this way all pair interactions at the spherical cutoff equal to zero. Application of this approach gave accurate results for several ordered and disordered ionic systems. Faster convergence was achieved by splitting the interaction in a short- and long-range contribution and applying a damping parameter, α, in the same fashion as in the Ewald summation. However, in contrast to the latter method, in the Wolf approach the long-range interactions are considered small with no contribution to the total electrostatic energy. Also, the parameter α, which in the Ewald summation controls the rate of convergence, controls the accuracy of the results in the Wolf approach and, therefore, should be carefully chosen. In a recent theoretical and computational work10 it was found that reasonable results can be obtained by choosing α = 2/Rc. The optimum value of α, however, is usually determined by comparing the results for different α’s with the results of the Ewald method. Despite the previous drawback, the fact that no calculations in the reciprocal space are required makes the Wolf method significantly faster than the Ewald summation. Over the past few years, several extensions and improvements in the original formulation have been proposed.11−14 Zahn et al.11 first noticed that in the original work of Wolf et al. the suggested expression for the forces did not correspond to the potential energy function. In order to remove this inconsistency, they suggested a different expression for the energy. Alternative expressions were also given later by Fennell et al.12 and will be discussed later in this work. The separation of the Coulomb potential into a short- and long-range contribution has been previously described and discussed in the framework of the local molecular field (LMF) theory.15 In the present work, an alternative partition of the Coulomb potential in a short- and long-range contribution is proposed. In the spirit of the Wolf approach, it is required that the short-term contribution of the potential vanishes at the cutoff distance in order to satisfy the charge-neutralization condition. However, the splitting function is chosen such as, at the same time, higher derivatives of both short- and long-term of the potential are continuous at any interparticle separation. The accuracy of the present scheme is examined by performing several molecular dynamics simulations for liquid water and aqueous solutions at various thermodynamic conditions using several different force fields for the description of the intermolecular interactions. The paper is organized as follows: in section 2 the approach of Wolf for the treatment of the electrostatic interactions is briefly described along with some of its more recent extensions. In the same section, details regarding the suggested new approach and its application to simple pairwise-additive and many-body polarizable force fields similar to that examined here are also presented. In section 3 the details of the molecular simulations aimed at comparing the new with the previous approaches are described. The results of the molecular simulations are presented in section 4, while final conclusions are drawn in the last section 5.

V Coul =

1 2

qiqj

N Ni(R c)

∑∑

rij

i = 1 j ′= 1

(1)

where rij is the distance between the charges i and j, while j′ indicates that j ≠ i. For simplicity of notation the Coulomb’s constant, ke = 1/4πε0, was set to one. It was also assumed that each atom i interacts only with the Ni(Rc) atoms that are inside a sphere of radius Rc centered on i, and that Rc is in any case less or equal to half of the box length, Rc ≤ L/2. Despite the simplicity of the previous expression, it is wellknown, that direct estimation of the Coulomb energy from eq 1 is difficult even for a very large spherical radius, as huge fluctuations in VCoul are encountered even for very small changes in the value of Rc. On the other hand, it is also known that the effective Coulomb potential of the ions in a condensed system is short ranged (see for example ref 9). Wolf et al. pointed out that when the net charge inside a sphere of radius Rc is close to zero, the Coulomb energy VCoul of atom i (which is in the center of the i sphere) is close to the correct one. Based on this observation, for an approximate estimation of the Madelung energy, Wolf et al.8 proposed to subtract a charge-neutralization term, Vneutr , from i the energy VCoul of each atom. The total energy of the system is i written then as V ≈ V Coul − V neutr

(2)

where V neutr =

1 2

N

∑ qiΔqi i=1

1 Rc

(3)

The net charge, Δqi, which is located at the surface of a sphere centered on atom i is given by the expression Ni(R c)

Δqi =

∑

qj

(4)

j=1

By substitution of eqs 4 and 3 in eq 2 we arrive at SP SP V SP = V short + V self

(5)

where SP V short

1 = 2

N Ni(R c)

∑∑ i = 1 j ′= 1

⎛1 1⎞ qiqj⎜⎜ − ⎟⎟ Rc ⎠ ⎝ rij

(6)

and SP V self =−

1 2Rc

N

∑ qi2 i=1

(7)

From eq 6 it becomes clear that the charge neutralization at the system surface is equivalent to the shifting of all Coulomb pairs of interactions u(rij) = qiqj/rij so that u(Rc) = 0. For this reason, the previous expressions for energy are denoted by SP (Shifted Potential). For rij > Rc the pair interactions u(rij) = 0. The last term of eq 5, which represents the self-energy of the system, VSP self, remains constant and, therefore, does not affect the derivatives of energy. It was recently shown16 that the previous expression converges to the Ewald energy although the convergence is very slow. One of the problems with the previous expression of the energy, VSP short, is that the associated force (or potential derivative) on ion i due to its interaction with ion j is discontinuous for a separation rij = Rc, since Fij(Rc) = qiqj/R2c ≠ 0 while Fij(rij > Rc) = 0.

2. THEORY 2.1. Shifted and Charge-Neutralized Coulomb Potential. Lets consider first a system of N charges qi inside a cubic box of linear dimension L for which periodic boundary conditions (PBC) have been imposed in the three directions. The total Coulomb energy of the system can be written as B

DOI: 10.1021/jp510612w J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B This may introduce problems during a molecular dynamics simulation such as poor energy conservation or even more severe artifacts. For this reason, Fennell et al.12 employed the shifted force function, the general form of which is given by ⎛ du(r ) ⎞ ⎟ (r − R c ) u shift(r ) = u(r ) − u(Rc) − ⎜ ⎝ dr ⎠ r = R c

Vlong

SF V short

N Ni(R c)

∑∑ i = 1 j ′= 1

(8)

⎤ ⎡1 ⎛ 1 ⎞ 1 qiqj⎢ − + ⎜ 2 ⎟(rij − Rc)⎥ Rc ⎢⎣ rij ⎝ Rc ⎠ ⎦⎥

DSP Vshort =

1 Rc

N

∑ qi2

V=

1 2

N

∑ ∑ qiqj 1 i = 1 j ′= 1

rij

=

1 2

N

⎡

N

∑ ⎢⎢ ∑ qiqj i=1

⎣ j ′= 1

i=1 j=1

1 2

N Ni(R c)

∑∑ i = 1 j ′= 1

DSF Vshort =

1 2

N Ni(R c)

∑∑ i = 1 j ′= 1

rij

(12)

s(rij)

i = 1 j ′= 1

1 − s(rij) 1 Vself = − lim 2 rij → 0 rij

(17)

⎡ erfc(αrij) erfc(αRc) qiqj⎢ − rij Rc ⎢⎣

in which both potential and forces are continuous at the cutoff radius. Notice again that in Fennell’s work VDSF self = 0. However, in contrast to the force-shifted potential described before, it was found in this study that, for the damping parameters α chosen in the present simulations, better accuracy in energy is obtained by using eq 19. Therefore, in the reported results below, the selfterm has been considered during the simulations. In a recent study,19 the self-term was included in DSF as well, while in a study20 of liquid water with the SPC/E model in which the selfterm was not considered, in agreement with our observations, no satisfactory values of energy were obtained, even for very long values of Rc (up to 20 Å). 2.3. Alternative Forms of Splitting Functions. It is trivial to rewrite the expressions for the shifted potential (SP) in eq 6 and shifted force (SF) in eq 9 in the form of eq 11, by using the following splitting functions respectively

where qiqj

(16)

⎛ erfc(aRc) ⎞ N α α DSF = −⎜ + + Vself exp( −α 2Rc2)⎟ ∑ qi2 π π Rc ⎝ ⎠ i=1

s(rij)

(11)

V = Vshort + Vself + Vlong

N Ni(R c)

⎛ erfc(αrij) erfc(αRc) ⎞ ⎟⎟ qiqj⎜⎜ − rij Rc ⎝ ⎠

(18)

In the previous expression s(rij) is a splitting function that separates the Coulomb interaction in two parts. Moreover, by assuming that the s(rij) function is such that the [1 − s(rij)]/rij does not diverge as rij → 0, the above expression can be written as

∑∑

(15)

(19)

1 − s(rij) ⎤ ⎥ + ∑ qiqj ⎥⎦ rij j ′= 1

1 2

rij

2 2 ⎤ ⎛ erfc(αR ) 2α exp( −α Rc ) ⎞ c ⎥ ⎟ r R ( ) +⎜ + − ij c 2 Rc π ⎥⎦ ⎝ Rc ⎠

N

Vshort =

1 − s(rij)

Regarding the term, Wolf et al. showed that it is equivalent to the reciprocal-space energy in the Ewald sum. For a wide range of values of the parameter α, VDSP long is very small DSP compared to VDSP short + Vself and remains almost constant for any atomic arrangement. For this reason, in the original development,8 and also in the present work, Vlong will be neglected. For several crystals examined, it was found that the Madelung constant obtained with the DSP method converges much faster with the spherical cutoff compared to the SP. However, as in the case of SP, DSP suffers from a discontinuity of the force at rij = Rc. Starting from DSP and using the general form of the shifted force function (eq 8), one arrives at the damped shifted force expression (DSF)12

The previous expressions not only maintain the charge neutrality required in the Wolf method but also meet the requirement for continuous forces since F(rij ≥ Rc) = 0. For this reason, the approach is denoted by SF (Shifted Force). Notice, however, that traditionally the self-term of the force-shifted potential is set to zero. Later, this approach will be examined further and the previous recommendation will be followed as it was also found in this work that, for the systems under study, the obtained energies with VSF self = 0 are closer to the energies obtained with the Ewald summation. In a recent study, Hansen et al.17 employed this approach in simulations of liquid water and a model of a molten salt. They arrived at the conclusion that the accuracy of the results obtained with SF is comparable to that of DSF (see next subsection) with a carefully chosen damping parameter α. 2.2. Damping the Pair Interaction. In order to improve the convergence properties of the previous approach and to decrease the oscillations in energy with increasing Rc, the Coulomb potential was damped according to the scheme: N

∑ ∑ qiqj

VDSP long

(10)

i=1

N

N ⎛ erfc(αRc) α ⎞ DSP = −⎜ + Vself ⎟ ∑ qi2 π ⎠ i=1 ⎝ 2Rc

(9) SF V self =−

N

Due to the close connection to the Ewald method, Wolf chose as the splitting function the complementary error function s(r) = erfc(αr), with α as the damping parameter, although he states that, in principle, any splitting function could be chosen.18 Notice that 1 − s(r) = 1 − erfc(αr) = erf(αr), with erf(αr) as the error function. When the neutralization condition is applied, one arrives at the following expression for the energy

Application of the previous expression to the Coulomb pair interaction leads to the following expressions for the short-range and self-energy respectively: 1 = 2

1 = 2

rij

(13) N

∑ qi2 i=1

⎛ rij ⎞ sSP(rij) = 1 − ⎜ ⎟ ⎝ Rc ⎠

(14)

and C

(20) DOI: 10.1021/jp510612w J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B and ⎛ rij ⎞ ⎛ rij ⎞2 sSF(rij) = 1 − 2⎜ ⎟ + ⎜ ⎟ ⎝ Rc ⎠ ⎝ Rc ⎠

(21)

Similar expressions can be written for the DSP (eq 16) and DSF (eq 18). Undoubtedly, there are more splitting functions s(rij) that maintain the charge neutrality condition as all previous expressions and provide continuous forces for rij = Rc as the SF and DSF. A necessary condition is that any s(rij) should satisfy, s(0) = 1 and s(Rc) = 0. The question to be asked, therefore, is whether there is another function s(rij) to separate the Coulomb interaction in a short and long term such as to minimize the contribution of the long term which is going to be omitted from the calculation. A possible choice is that the s(rij) should be constructed such that higher derivatives of the pair potential ushort(rij) = qiqjs(rij)/rij decay to zero for rij → Rc. However, it should be also required that the long term of the pair potential ulong(rij) = qiqj(1−s(rij))/rij varies smoothly. A function that satisfies these requirements is ⎛ rij ⎞ ⎛ rij ⎞3 ⎛ rij ⎞4 sSP2(rij) = 1 − 2⎜ ⎟ + 2⎜ ⎟ − ⎜ ⎟ ⎝ Rc ⎠ ⎝ Rc ⎠ ⎝ Rc ⎠

Figure 1. Illustration of the short-range (ushort) and long-range (ulong) contributions to the Coulomb pair energy 1/rij as computed using the sSP2(rij), sSP3(rij), and the erfc( π rij/Rc) splitting functions.

(22)

2.4. Polarizable Force Fields. For the sake of completeness, the mathematical expressions needed for the energy evaluation of a molecular system containing charges qi and polarizable sites λi are presented below (the polarizabilities are denoted here with the symbol λ instead of the usual α to avoid confusion with the damping parameter α used before). The expressions given below are similar to that of the Ewald summation (see for example ref 23) but without the reciprocal-space part. It will be assumed that there are no electrostatic interactions between the atoms of the same molecule, which is usually the case for small molecules such as water. The electrical energy may be separated to the electrostatic (Coulomb) and the induction energy:

uSP2 short(rij)

This function provides continuous potential up to its second derivatives at rij = Rc (SP2 indicates shifted potential and up to its second derivative). At the same time, the first derivative of uSP2 long(rij) is zero at rij = 0. This function was first suggested by Markland et al.21 to separate the electrostatic interactions in short and long contributions in the course of a path integral molecular dynamics simulation. However, in that case a different distance σ < Rc instead of Rc was used in eq 22 to separate the two contributions, while, in contrast to the present work, the longterm contribution was approximately computed. It was shown21 that the sSP2(rij) splitting function is similar to the error function erfc( π rij/Rc). Following the same procedure, another splitting function sSP3(rij) was constructed22 that provides continuous values for the pair potential uSP3 short(rij) up to its third derivative (denoted as SP3). The first and second derivatives of the uSP3 long(rij) are zero at rij = 0. ⎛ rij ⎞ 7 ⎛ rij ⎞ 21 ⎛ rij ⎞ 5 ⎛ rij ⎞ ⎜ ⎟+ ⎜ ⎟ − 7⎜ ⎟ + ⎜ ⎟ 4 ⎝ Rc ⎠ 4 ⎝ Rc ⎠ 2 ⎝ Rc ⎠ ⎝ Rc ⎠ 5

sSP3(rij) = 1 −

6

V = V Coul + V ind

where V Coul =

7

i=1

7 8Rc

i=1

(28)

i

∑ Tij̃ qj}

j , m′j

j , mj

(29)

and α

φi μ = φi μ ,real + φi μ ,self = {∑ Tijαμj , α } + {−∑ Tij̃ μj , α } j , m′j

j , mj

(30) (24)

In the previous expressions it was assumed that, in the real terms of the potential ϕi, there are no intramolecular interactions between the atoms of the molecule mj where atom j belongs (denoted in the summation by mj′). In contrast, for the self-terms only interactions between different atoms of the same molecule mj are considered. Tij and Tαij are electrostatic tensors defined as

N

∑ qi2

∑ μi ,α Eiq,α

ϕiq = ϕiq ,real + ϕiq ,self = {∑ Tijqj} + {−T̃(0)qi −

and SP3 V self =−

1 2

(27)

i

The total electric potential ϕi on atom i is also separated to the electrostatic ϕqi and induced ϕμi parts as ϕi = ϕqi + ϕμi , where

N

∑ qi2

∑ qiϕiq

V ind = −

The short- and long-term contributions to the pair potential obtained using sSP2(rij), sSP3(rij), and erfc( π rij/Rc) are illustrated in Figure 1. It can be seen that while both sSP2(rij) and sSP3(rij) functions were constructed in order to satisfy continuous derivatives of the potential, the corresponding forms of the short and long contributions of the potential still closely resemble the one obtained using the complementary error function. For the SP2 and SP3 functions the corresponding selfterms as obtained from eq 14 are written as 1 Rc

1 2

and

(23)

SP2 V self =−

(26)

(25) D

DOI: 10.1021/jp510612w J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B Tij = Tijα

s(rij) rij

according to which an effective density replaces the point charges/dipoles. For all previous water models, in order to determine the accuracy that each method predicts the liquid density, constant temperature and pressure (NPT) simulations were performed using the Berendsen thermostat and barostat29 in the range T = 260−370 K and P = 1 atm for a system of N = 256 water molecules in a cubic box. Periodic boundary conditions (PBC) were applied using the minimum image convention. The van der Waals and short-range electrostatic interactions were truncated at a distance Rc = 9 Å. For the former, standard long-range corrections1 were included in the energy and the virial tensor. As regards the Ewald summation, the implementation for polarizable force fields described in ref 23 was employed. For the treatment of the real-space electrostatic interactions the Ewald parameter was set to 0.38 Å−1, while, in the reciprocal space, the largest k-vector considered was k = 8. A total of 1108 k-vectors were included in the calculation. In order to satisfy the internal constrains of the water molecules during the simulations with the rigid models, the Rattle algorithm30 was employed for integrating the equations of motion. The time step was 2 fs, while the total integration time was 10 and 8 ns in the simulations with the TIP4P/2005 and POL3 models, respectively. For the flexible TTM3-F model, a multiple time step integration algorithm was used.31 In this case, every δt = 0.25 fs, the fast intramolecular component of the force was used for the integration of the equations of motion, while the full force was computed every three steps. The total integration time was 7.5 ns. The velocityverlet algorithm32 was employed in this case. In order to examine the dependence of the obtained results from the real-space cutoff, Rc, a few additional NPT simulations were performed with the previous cutoff-based approaches and the TIP4P/2005 model using N = 500 water molecules and Rc = 12 Å. In addition to the density, the dielectric constants ε0 at T = 300 K were also obtained from 10 ns long simulations at the constant temperature and volume (NVT) ensemble. The density was the one determined before for each water model and method for the treatment of the electrostatic interactions. In these simulations the temperature was kept constant using a canonical velocityrescaling thermostat.33 The total dipole moment M of the system was used for the calculation of ε0 according to 4π (⟨M2⟩ − ⟨M⟩·⟨M⟩) ε0 = 1 + 3kBTV (35)

,

= ∇α Tij ,

Tijαβ = ∇α ∇β Tij

(31)

The T̃ ij term appearing in the self-terms is defined as T̃ ij = 1/rij − Tij, while the higher terms in the expansion are obtained similarly to Tij in eqs 31. Finally, the term T̃ (0) is defined as T̃ (0) = limrij → 0 T̃ ij and the splitting function s(rij) should be such that the value of T̃ (0) is finite. The Einstein summation convention for Greek indices was followed in eqs 30 and 28. The α-component of the electric field, Ei,α, on atom i is also μ written as Ei,α = Eqi,α + Ei,α with Eiq, α = −

∂φiq ∂ri , α

α

= {∑ Tijαqj} + {−∑ Tij̃ qj} j , m′j

j , mj

(32)

and Eiμ, α = −

∂φi μ ∂ri , α

αβ = {∑ Tijαβμj , β } + {−∑ Tij̃ μj , β } j , m′j

j , mj

(33)

Finally, the dipole moment μi,α needed for the calculation of Vind in eq 28 is μi , α = λiEi , α

(34)

For a system of D polarizable sites, eqs 32−34 form a 3D × 3D system of linear equations. During the course of a molecular dynamics simulation the system is usually solved using iterative methods. It should be also mentioned that, in the Ewald summation as well as in the present study, while there are no electrostatic interactions between atoms that belong to the same molecule, there is still an intramolecular contribution in the selfterm, which affects both the energy and forces. For the latter it was found that the contribution is very small. Notice also that, apart from a different splitting function, the same approach described above was recently employed19 for the study of the liquid−vapor interface of water using a flexible version of the SPC/E model.

3. METHODOLOGY In order to assess the accuracy of the previously described approaches, several molecular dynamics simulations of liquid water as well as of aqueous NaCl solutions have been performed. The obtained results have been compared to the results of the benchmark Ewald technique. For liquid water the following water models have been employed for the description of the intermolecular interactions: (i) The TIP4P/2005 model24 which has been found to predict very accurately several properties of liquid water and ice compared to other rigid, nonpolarizable water models;25 (ii) the rigid and polarizable POL3 model,26 in which oxygen and hydrogen atoms are polarizable and carry partial charges; and (iii) the TTM3-F model,27 which is an ab initio based flexible, polarizable model. The hydrogen atoms and an extra site (M-site), which is located at the bisector of the HOH bend angle, carry the partial charges. The M-site is also polarizable, and it is used for the purpose of describing manybody polarization effects. The TTM3-F model employs the Thole approach for the treatment of electrostatic interactions,28

The self-diffusion coefficient D was also determined using the equation D=

1 3

∫0

∞

dt ⟨v(t ) ·v(0)⟩

(36)

where ⟨v(t)·v(0)⟩ is the velocity autocorrelation function. In the present study its length was 2 ps, and it was computed as the average of the autocorrelation functions obtained from 1000 constant volume and energy (NVE) trajectories of 4 ps length each. All NVE trajectories initiated with momenta obtained from the Maxwell distribution. During the simulations with the TIP4P/2005 water model and the Ewald summation technique, 1000 atomic configurations were stored on the disk and used afterward for a first evaluation of the accuracy of the various schemes. More specifically, using as a reference the benchmark results of the Ewald method, the average relative errors of the energy, forces, and the trace of the virial tensor were computed for these configurations. The error in forces for each configuration was defined as14 E

DOI: 10.1021/jp510612w J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B Ferr =

1 N

N

∑ i=1

|f iew − f imethod | f iew

(37)

ew method where f ew = |fmethod | is the magnitude of the force i = |fi | and f i i acting on atom i as obtained by the Ewald summation and one of the approximate methods examined in this study, respectively. Based on the results of this analysis at T = 300 K, it was also possible to determine a damping parameter α for the DSF method which is expected to be close to the optimum one and it was used later during the MD simulations. As it will be shown later, the results obtained with the TIP4P/ 2005 model and the SP3 function were in better agreement with the Ewald method than the other approximate schemes; for this reason, in the studies of liquid water with the TTM3-F and POL3 as well as in the NaCl solution, the simulations were performed only with the Ewald and the SP3 function. In all simulations of liquid water with the polarizable models, N = 256 water molecules and Rc = 9 Å were used. For the study of the NaCl solution, the SPC/E model34 was employed for the description of water−water interactions, while the Lennard−Jones parameters used for the description of dispersion interactions between ion−water and ion−ion were taken from Smith and Dang.35 It should be noted that the same potential model was recently employed by Aragonez et al.36 for the study of solubility of NaCl in water. In order to examine the accuracy of the present calculations, the same system parameters were used. Therefore, during the NPT simulations at T = 298 K, the simulation box comprises N = 270 water molecules and NNaCl = 1−25 molecules of salt. When the accurate Ewald summation technique was employed, the results for the density of the solution as a function of the salt concentration were found to be in excellent agreement with the previously reported ones. All results reported here have been obtained using our own MD computer codes. For the simulations with polarizable water models, the parallelization scheme described in ref 37 was employed. In a few cases the dl-poly38 suite of codes was employed for the validation of the present calculations.

Figure 2. Comparison of the Madelung energies against Rc obtained using the DSP, SP2, and SP3 methods. The Rc is expressed in units of the lattice parameter, a.

employed for the treatment of the electrostatic interactions during the NPT simulations performed in the temperature range T = 260−370 K. Based on the analysis of 1000 configurations obtained during the 10 ns trajectory at T = 300 K for a system of N = 500 water molecules, the relative errors in energy, force norm, and virial tensor were computed using all previously described approaches. The results obtained with Rc = 9 and 12 Å are summarized in Table 1. By comparing the results of the SP and SF forms, it is seen that shifting the force of the potential leads to more accurate estimations of the forces and virial. However, at the same time, the accuracy in the energy estimation is reduced. In the same table, it can also be seen that the relative Table 1. Relative Errors for the Energy, Force Norm, and Virial (Last Three Columns Respectively) with Respect to the Ewald Technique for Different Methods (First Column) and Damping Parameters, α (second column)a

4. RESULTS AND DISCUSSION 4.1. Madelung Constant of NaCl Lattice. As a first test to examine the accuracy and efficiency of the new approach, the Madelung constant for a perfect NaCl crystal has been calculated. The obtained values of energy as a function of the real-space cutoff distance, Rc, as computed with DSP, SP2 and SP3 are shown in Figure 2. For the DSP the damping parameters α = 0.8/ a and α = 1.5/a (a is the lattice constant) previously used by Wolf8 were chosen. It can be seen that the DSP with α = 0.8/a shows oscillatory behavior up to 2.0−2.5a while after 3.0a it is practically converged. The same behavior is exhibited by the SP2 and SP3 although the oscillations are about 1 order of magnitude smaller compared to the DSP and the obtained results converge after 2.5a. Finally, using DSP with α = 1.5/a, faster convergence of energy is achieved even for a very short cutoff radius, Rc = 1.5a. We may conclude that SP2 and SP3 converge relatively fast to the Madelung constant; however, DSP with a carefully chosen damping parameter may converge even faster. 4.2. Pairwise Additive TIP4P/2005 Model. Molecular dynamics simulations with the pairwise additive TIP4P/2005 water model are computationally much less demanding compared to that with many-body polarizable models. Therefore, a comprehensive analysis of all approaches has been performed for this model. At the beginning, the Ewald summation was

method DSP

DSF

α[Å−1] 0.21a 0.23b 0.26c 0.26a,c 0.22b

SP SF SP2 SP3 DSP

DSF SP SF SP2 SP3

0.16a 0.17b 0.21c 0.21a,c 0.17b

energy Rc = 9.0 Å 0.16 0.17 0.23 0.37 0.91 0.87 2.99 0.21 0.06 Rc = 12.0 Å 0.04 0.05 0.11 0.16 0.41 0.51 1.70 0.08 0.02

force

virial

1.70 1.59 1.81 1.93 1.30 17.01 3.08 1.20 1.20

2.08 1.17 0.67 0.60 0.76 31.44 4.35 0.49 0.18

1.08 1.00 1.26 1.27 0.82 12.61 1.75 0.69 0.69

1.70 1.11 0.27 0.31 0.54 28.00 1.55 0.25 0.10

a

The superscripts a,b,c indicate the optimal damping parameter for the energy, force norm, and virial, respectively.

F

DOI: 10.1021/jp510612w J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

Rc = 12 Å are compared to the Ewald method, while afterward, for the polarizable models and aqueous solutions, results obtained with the Ewald summation and the SP3 with Rc = 9 Å only are presented and compared. Molecular dynamics simulation in the NPT ensemble for the temperature range T = 260−370 K were performed using the SF, DSF, SP2, and SP3 approaches as well as the Ewald summation. For the DSF method the damping parameters α = 0.25 Å−1 for Rc = 9 Å and α = 0.20 Å−1 for Rc = 12 Å were chosen. According to the results in Table 1, these damping parameters are closer to the parameters that minimize the errors in energy and virial, while also providing reasonable accuracy in forces. The results for the liquid density and energy as a function of temperature are illustrated in Figures 4 and 5 respectively, along with the

errors are significantly reduced in the damped SP and SF forms in which the optimal damping parameter is chosen. However, the optimal values of α are different for the energy, force norm, and virial in the DSP, while in the DSF the optimal α for the force norm is different from the one of energy and virial. The same conclusion, namely that SP provides more accurate energies than SF, holds here, since when the optimal damping parameters are chosen, the DSP predicts more accurate energies and the DSF predicts more accurate forces and virial. The SP2 form, although it does not use any adjustable parameters, gives as accurate results as the DSP and DSF with optimal damping parameters. Finally, the SP3 provides results that, when compared to those obtained with the Ewald method, are more accurate than the other approximate schemes examined in this work. Results of a similar analysis for the SF, DSF, and SP3 are also illustrated in Figure 3 for a real-space cutoff Rc between 9 and 12

Figure 4. Temperature dependence of the density ρ(T) (in g/cm3), for the TIP4P/2005 water model as predicted with the Ewald summation technique, and the SF, DSF, SP2, and SP3 methods with Rc = 9 Å. The results of SF, DSF, and SP3 methods with Rc = 12 Å are also shown.

Figure 3. Relative error in energy, forces, and virial of the TIP4P/2005 model as a function of the cutoff, for the SF, SP3, and DSF with various damping parameters α (in Å−1) with respect to the Ewald technique.

benchmark results obtained with the Ewald method. The corresponding numerical values are tabulated in Tables 2 and 3. The estimated errors for these quantities vary as a function of

Å. For the DSF three different damping parameters, α, have been used. All methods provide more accurate results for longer cutoffs. It can be easily seen that the optimum value of the damping parameter α used in DSF is different for each of the quantities examined, while it varies significantly with the cutoff. SP3 appears to provide more accurate estimates for all three quantities and for all values of the cutoff. In particular, the error in energy and virial of the SP3 with Rc = 9 Å is comparable to that of DSF with Rc = 12 Å. The results presented so far provide a first estimate for the accuracy of each approach; the real issue is, however, whether, in the course of a molecular simulation, the predictions of previous schemes are accurate enough so they can be used as an alternative to the Ewald summation. According to the results of several previous studies (see for example refs 7 and 12) a real-space cutoff of 12 Å or more is suggested. In particular, Steinbach et al.7 state that the cutoff is more important than the shifting function employed. It is therefore interesting to investigate if the remarkable accuracy that SP3 exhibits so far can lead to accurate results during molecular simulations using the relative short cutoff distance, Rc = 9 Å. For this reason, for the TIP4P/2005 model, the predictions of all previously described approaches obtained from molecular dynamics simulations with Rc = 9 Å and

Figure 5. Same as Figure 4 for the average potential energy (in kcal/ mol). G

DOI: 10.1021/jp510612w J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

Table 2. Temperature Dependence of the Density ρ(T) (in g/cm3), for the TIP4P/2005 Water Model As Predicted with the Ewald Summation Technique, and the SF, DSF, SP2, and SP3 Methods with Rc = 9 Åa

a

temp

Ewald

SF

DSF

SP2

SP3

SF

DSF

SP3

260 270 280 290 300 310 340 370

0.9997 1.0018 1.0018 1.0001 0.9979 0.9944 0.9794 0.9593

0.9943 0.9956 0.9969 0.9957 0.9928 0.9892 0.9733 0.9514

0.9984 0.9998 0.9992 0.9975 0.9954 0.9914 0.9755 0.9541

0.9995 1.0006 1.0006 0.9990 0.9970 0.9926 0.9775 0.9568

1.0015 1.0032 1.0022 1.0012 0.9989 0.9952 0.9801 0.9599

0.9995 1.0005 0.9997 0.9980 0.9951 0.9914 0.9755 0.9543

0.9983 1.0002 1.0006 0.9986 0.9963 0.9924 0.9769 0.9561

0.9998 1.0014 1.0016 1.0001 0.9974 0.9938 0.9788 0.9586

The results of SF, DSF, and SP3 method with Rc = 12 Å are also shown, in the last three columns of the table.

Table 3. Same as in Table 2 for the Average Potential Energy (in kcal/mol) temp

Ewald

SF

DSF

SP2

SP3

SF

DSF

SP3

260 270 280 290 300 310 340 370

−12.025 −11.857 −11.694 −11.537 −11.382 −11.231 −10.789 −10.362

−11.640 −11.468 −11.293 −11.125 −10.964 −10.804 −10.347 −9.909

−12.059 −11.891 −11.726 −11.569 −11.413 −11.258 −10.814 −10.382

−12.039 −11.871 −11.709 −11.548 −11.393 −11.240 −10.798 −10.370

−12.012 −11.846 −11.688 −11.529 −11.375 −11.225 −10.786 −10.360

−11.744 −11.576 −11.418 −11.263 −11.109 −10.960 −10.519 −10.094

−12.026 −11.858 −11.696 −11.537 −11.384 −11.232 −10.789 −10.361

−12.017 −11.852 −11.687 −11.532 −11.376 −11.225 −10.783 −10.356

temperature; they are 0.002 g/cm3 for the density and 0.01 kcal/ mol for the energy at T = 260 K, while at T = 370 K the corresponding errors are 0.0005 g/cm3 and 0.002 kcal/mol for the density and energy, respectively. Regarding the liquid density, it can be seen in Figure 4 that, in agreement with the Ewald technique, all methods predict a maximum value during the simulation at T = 270 K. However, the SF with Rc = 9 Å underestimates the density by a large amount (0.006−0.008 g/cm3), while with Rc = 12 Å the results are significantly improved and they differ by 0.001−0.003 g/cm3 to that obtained with the Ewald summation. The DSF and SP2 with Rc = 9 Å, also, systematically underestimate the liquid density by a smaller amount: between 0.001−0.004 g/cm3 for the DSF and 0.001−0.002 g/cm3 for the SP2. Similarly to the SF, DSF provides significantly more accurate results by increasing the cutoff to 12 Å. Among all methods, the SP3 shows the best agreement with the Ewald results, as the density differs by an amount usually less than 0.001 g/cm3. However, while this error is similar to the statistical errors of the previous MD simulations, the densities predicted with the SP3 with Rc = 9 Å are always larger than those found with the Ewald method, suggesting a systematic difference. By increasing the cutoff to Rc = 12 Å the agreement between the SP3 and the Ewald method becomes even better, and no systematic differences in the obtained values for the density were observed. According to the results presented in Table 3 and illustrated in Figure 5, similar conclusions with the ones found for the density can be drawn for the average potential energy of water as a function of temperature. The SF shows poor agreement with the Ewald method, as it overestimates the energy by ∼0.4 kcal/mol for Rc = 9 Å and ∼0.25 kcal/mol for Rc = 12 Å. In contrast, the results of all other methods differ by less than 0.04 kcal/mol compared to the results obtained with the Ewald summation, even with Rc = 9 Å. SP2 and SP3 yield slightly better estimates for the energy than DSF for Rc = 9 Å while for Rc = 12 Å the energies of SP3 and DSF are practically the same as those of the Ewald summation.

It is worth mentioning that the accuracy of the results obtained from the MD simulations for the energy and density using all previous approximate schemes are in agreement with the conclusions drawn before, during analysis of the 1000 configurations. The previously presented approximate schemes were further evaluated by comparing the radial distribution functions (RDFs) with those obtained with the Ewald method. All RDFs were computed during a constant temperature and volume (NVT) simulation at T = 300 K. The density of the system was in all cases the one determined before with the Ewald method (ρ = 0.9979 g/cm3). Figure 6 (upper panel) shows the RDF between oxygen atoms, gEwald OO , obtained with the Ewald method. For a convenient

Figure 6. (a) Oxygen−oxygen RDF obtained from a NVT simulation at T = 300 K with the TIP4P/2005 water model and the Ewald method for the treatment of electrostatics. (b) Difference between the oxygen− oxygen RDFs obtained with the SF, DSF, SP2, SP3 forms and the Ewald method. H

DOI: 10.1021/jp510612w J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

4.3. Polarizable TTM3-F and POL3 Models. The results presented so far for the pairwise additive water model indicate that not only the energies but also the virial and the forces are accurately predicted with the SP3 splitting function. It is not surprising, therefore, that for the many-body polarizable POL3 and TTM3-F force fields, for which the Vind is based on the derivative of the electrostatic potential (eq 32), the agreement between the Ewald technique and the SP3 is good as well. In Figure 7 the temperature dependence of the energy and density

visualization, in the lower panel of Figure 6 only the differences ΔgOO between the gOO obtained with the SF, DSF, SP2, and SP3 forms and the Ewald method are illustrated. It was found that comparing the RDFs was a less sensitive test for the evaluation of all approaches. Only the RDF obtained with the SF form appears to differ slightly from that of the Ewald method both at small and at large interatomic separations. However, these differences are still small and insignificant. Therefore, it may be concluded that all SF, DSF, SP2, and SP3 approaches predict gOO in excellent agreement with the Ewald method. Similar conclusions to that of the oxygen−oxygen RDF were drawn for the oxygen−hydrogen and hydrogen−hydrogen RDFs (not shown here). As a final test for the accuracy of all previous approaches, the values of the dielectric constant and the diffusion coefficient were computed at T = 300 K and their values are given in Table 4. Table 4. Dielectric Constants and Self-Diffusion Coefficients for the TIP4P/2005 Liquid Water at T = 300 K for a System with N = 256 Water Molecules and Rc = 9 Å and a System with N = 500 Water Molecules and Rc = 12 Å as Obtained with the Ewald Summation, SF, DSF, and SP3 method Ewald SF DSF SP3 Ewald SF DSF SP3

ε0 Rc = 9.0 Å 57 55 57 58 Rc = 12.0 Å 59 59 55 59

D[Å2 ps−1] 0.225 0.225 0.231 0.225

Figure 7. Temperature dependence of the density ρ (in g/cm3) and average potential energy (in kcal/mol) for the TTM3/F water model as predicted with the Ewald summation technique and the SP3.

0.232 0.241 0.233 0.223

for the TTM3-F water model is illustrated, while the corresponding results for the POL3 model are shown in Figure 8. Numerical values for the previous quantities are presented for

They have been computed at the density that each method predicts. The error in the dielectric constant was estimated by the size of fluctuations in the running average during the last 4 ns of simulation to be about 3. For the diffusion coefficients the errors are 0.004 and 0.003 Å2 ps−1 for N = 256 and 512 water molecules, respectively. Even though the densities that the simulations were performed at are different from each other, the obtained results for both the dielectric constant and self-diffusion coefficient are in all cases the same within the statistical error. Also, the cutoff does not seem to affect the accuracy of the results. These findings are in agreement with previous conclusions drawn17 using the SPC/E model and the SF method. The previous results clearly show that in the course of a molecular dynamics simulation of liquid water with a pairwise additive force field, treatment of electrostatic interactions with the SP3 form provides accurate estimates for all properties examined. In particular, the predictions of SP3 for the density of liquid water are significantly closer to those of the Ewald summation, compared to the other approaches examined here (SF, DSF, and SP2). It is worth mentioning that several of the properties examined before are accurately predicted by DSF and SF as well. At the same time, the functional form of SP3 is simpler than the one of DSP and DSF, while it does not require optimization of any parameter. This is a very important point, as the present and several previous studies have shown that in DSF the damping parameter should be separately optimized for each property of interest. For these reasons, in the rest of this study, only the results of SP3 will be compared to the ones with the Ewald method.

Figure 8. Same as Figure 7 for the POL3 model.

the TTM3-F and POL3 in the Tables 5 and 6, respectively. It can be seen that for the TTM3-F model the agreement between the Ewald and the SP3 in the potential energy is within 0.01 kcal/ mol. Regarding the density, the SP3 overestimates its value by ∼0.002 g/cm3, showing a similar trend to what was found earlier for the TIP4P/2005 model. It is therefore expected that a slightly larger value of the real-space cutoff Rc will further improve the accuracy of the results obtained with the SP3. The self-diffusion coefficient and the dielectric constant at T = 300 K were I

DOI: 10.1021/jp510612w J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B Table 5. Temperature Dependence of the Average Potential Energy (in kcal/mol) and Density (in g/cm3) for the TTM3-F Water Model As Predicted with the Ewald Summation Technique and the SP3 density

Table 7. Comparison of the Results Obtained with the Ewald Method and the SP3 for the Density and the Average Potential Energy of an Aqueous Solution with 270 Waters and NNaCl Molecules of Salt

energy

density

energy

temp

Ewald

SP3

Ewald

SP3

NNaCl

Ewald

SP3

Ewald

SP3

260 270 280 290 300 310 330 350 370

1.0010 1.0023 1.0025 1.0028 1.0024 1.0001 0.9946 0.9862 0.9759

1.0035 1.0044 1.0057 1.0051 1.0037 1.0021 0.9956 0.9874 0.9771

−10.479 −10.307 −10.130 −9.955 −9.782 −9.610 −9.278 −8.951 −8.627

−10.469 −10.294 −10.116 −9.945 −9.778 −9.608 −9.277 −8.956 −8.632

1 2 3 4 5 7 10 12 15 17 20 22 25

1.0076 1.0158 1.0239 1.0318 1.0390 1.0538 1.0747 1.0875 1.1061 1.1206 1.1342 1.1444 1.1590

1.0083 1.0164 1.0243 1.0323 1.0401 1.0545 1.0756 1.0887 1.1071 1.1189 1.1352 1.1455 1.1607

−11.77 −12.37 −12.96 −13.53 −14.10 −15.22 −16.84 −17.87 −19.38 −20.34 −21.75 −22.66 −23.96

−11.77 −12.37 −12.95 −13.53 −14.10 −15.21 −16.83 −17.87 −19.37 −20.33 −21.74 −22.63 −23.95

Table 6. Same as Table 5 for the POL3 Model density

energy

temp

Ewald

SP3

Ewald

SP3

260 270 280 290 300 310 340 370

1.0308 1.0243 1.0172 1.0080 0.9999 0.9910 0.9605 0.9255

1.0318 1.0250 1.0172 1.0092 1.0005 0.9914 0.9611 0.9265

−10.497 −10.318 −10.142 −9.971 −9.802 −9.636 −9.143 −8.652

−10.491 −10.317 −10.142 −9.971 −9.803 −9.636 −9.143 −8.655

estimated at D = 0.220 ± 0.005 Å2 ps−1 and ε0 = 95 respectively, with the SP3 in excellent agreement with the Ewald method (D = 0.222 ± 0.005 Å2 ps−1 and ε0 = 96 respectively). The POL3 model, despite the fact that it incorporates three polarizable sites instead of one, shows an even better agreement between the results obtained with the Ewald method and SP3. In this case, the obtained results differ by less than 0.01 kcal/mol for the energy and 0.001 g/cm3 for the density. The diffusion coefficients at T = 300 K were calculated as D = 0.268 ± 0.005 Å2 ps−1 (Ewald) and D = 0.265 ± 0.005 Å2 ps−1 (SP3), while the dielectric constants are ε0 = 121 (Ewald) and ε0 = 123 (SP3). 4.4. Aqueous Solutions. As a final test of the present approach, mixtures consisting of 270 waters and 1−25 NaCl molecules were examined. The results of the simulations with the Ewald method and the SP3 are summarized in Table 7 and are illustrated in Figure 9. It is worth noting that the density of the solution as a function of the number of the salt molecules obtained with the Ewald method practically coincides with the published results obtained with Monte Carlo techniques.36 Direct comparison of the Ewald and SP3 leads to exactly the same conclusions as in the previous studies for liquid water. The density of the solution is systematically slightly overestimated by the SP3. However, the difference is usually smaller than 0.001 g/ cm3 which, for practical purposes, is accurate enough. Based on the previous results for liquid water with the TIP4P/2005 model, it is expected that a slightly longer cutoff will further improve the agreement between the Ewald and SP3. Regarding the average potential energy, the SP3 and Ewald methods provide almost the same values for all concentrations of the solution examined. 4.5. Timings. In order to evaluate the efficiency of the SP3 with respect to the Ewald summation, additional simulations with the three water models have been performed and the CPU time required for the evaluation of the electrostatic energy and

Figure 9. Density ρ (in g/cm3), and average potential energy (in kcal/ mol) for a NaCl aqueous solution as a function of the number of NaCl molecules in 270 water molecules predicted using the Ewald summation technique and the SP3.

derivatives was recorded. In all benchmarks a system of N = 256 water molecules with a 9 Å real space cutoff was employed for 1000 integration steps in the NVE ensemble, while all other technical details of the simulations were the same as those before. For the TIP4P/2005 model it was found that SP3 requires for the evaluation of the electrostatic interactions about half the time that Ewald summation needs. About 35% of the gain comes from the reciprocal-space calculations; the remainder comes from the evaluation of the pair interactions between charges, which is faster using the SP3 function with respect to the error function. For the polarizable TTM3-F and POL3, that incorporate one and three polarizable sites respectively, calculation of the induced dipoles as well as the induction energy and its derivatives is also required. In our implementation of the Ewald summation, reciprocal space represents 45%−55% of the total time needed for computations involving induced dipoles. These computational gains combined with the faster calculations in real space make the calculation of induced dipoles with the SP3 function about 2.5−3 times faster than the Ewald technique. Overall for the two polarizable models, calculations of the electrostatic and J

DOI: 10.1021/jp510612w J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

model that incorporates point charges and three polarizable sites and the TTM3-F model that incorporates smeared charges (according to the treatment Thole suggested) and one polarizable site. Moreover, even for a system comprised of more than one component as an aqueous solution, the accuracy of the obtained results was satisfactory. In a recent study,40 few interesting theoretical developments based on the DSF function were presented for the efficient calculation of the many-body polarization interactions. The application of these approaches in several molecular systems gave very accurate results demonstrating in this way that cutoff-based approaches can accurately describe polarization. As the systems and the methodology followed in that study were different from the present one, a direct comparison of the two approaches is yet not possible. However, it would be interesting in future to apply some of the ideas presented in that work with the SP3 function as well. It should also be mentioned that truncation methods have been employed in the past for the study of biological systems and water interfaces (for a quick review see ref 5 and references therein). It was found that while for homogeneous bulk systems these methods are able to make accurate predictions, for inhomogeneous systems inaccurate or incorrect results that strongly depend on the truncation method and choice of parameters may be obtained. An example is the water interface in which the dielectric constant significantly varies across the interface. In this case an extremely long cutoff was necessary,19,41 in order to obtain satisfactory agreement with the results of the Ewald method. In that respect, for the previous systems, the present approach that also takes into account only the local environment of each atom is expected to have similar drawbacks as in the case of the other methods. In summary, the proposed SP3 function retains the simplicity of the Wolf method for the treatment of electrostatic interactions while, for the systems examined, providing results comparable to those of the Ewald method without the need for adjustable parameters. For further validation of the present approach, work for a wider range of systems is currently in progress.

induced interactions are more than 2.5 times faster when the SP3 function instead of the Ewald summation is used. It should be emphasized that the previous numbers are very sensitive not only to the method used for the calculations in the reciprocal space (e.g., PME is faster to the Ewald technique employed here) but also to several other technical details. For example, by employing an approximate error function (see for example ref 39), or using look up tables (a technique that is used in dl-poly and several other major packages), the computational cost in real space between the Ewald summation and the SP3 becomes comparable.

5. CONCLUSIONS It is already known that, in a condensed system, accurate estimations of the electrostatic energy can be made by shifting the potential and maintaining charge neutrality. In this work we investigated the effect of shifting not only the potential and its corresponding first derivative but also higher potential derivatives which now are vanishing when the interatomic separation approaches the spherical cutoff. The accuracy of two splitting functions SP2 and SP3, which maintain charge neutrality and provide continuous second and third derivatives respectively, was investigated by performing molecular dynamics simulations of liquid water and aqueous solutions of sodium chloride. We have employed pairwise additive and many-body polarizable force fields for the description of the intermolecular interactions. For liquid water we have performed simulations for a wide-range of temperatures, while for the salt solutions for a wide range of salt concentrations. In all cases, it was found that employing the SP3 function to split the electrostatic interactions in short- and long-range contributions leads to results which are in reasonable agreement with the Ewald summation. Based on the results with the pairwise additive TIP4P/2005 water model, similar conclusions with previous studies12,14,17 were drawn. More specifically, it was found that in the DSP and DSF methods the damping parameter α should be optimized independently for each property in order to reproduce accurate results during molecular simulations (see Table 1). Despite the fact that the SP2 and SP3 splitting functions do not incorporate any adjustable parameters, they provide higher accuracy than the other approaches, especially the SF (which also does not have any adjustable parameters). A direct comparison of the results obtained with the SP2 and SP3 functions showed that the latter are more accurate. At the same time, the evaluation of an additional term in the SP3 polynomial compared to the SP2 does not increase significantly the computational cost. Therefore, among the two proposed functions the SP3 is recommended. In almost all simulations performed in the present study, the spherical cutoff distance was set to Rc = 9 Å. It was found that, in the case the SP3 function was employed, the obtained values of potential energy and density were slightly higher than those of the Ewald method. By increasing this value to Rc = 12 Å, the results obtained with the two methods are practically indistinguishable for the case of the TIP4P/2005 model. It is emphasized again that a value of Rc in the range 9−12 Å is used in almost all studies where the contribution of the reciprocal space is taken into account. The quality of the results obtained with the SP3 function was not affected by the complexity of the potential model employed for the description of the interatomic interactions. Therefore, the excellent agreement between the SP3 and Ewald summation found for a simple pairwise additive model (TIP4P/2005) was also found for the many-body polarizable models: the POL3

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Phone: +30 2810-545067. Fax: +30 2810-545166. Present Address †

Environmental Chemical Processes Laboratory, Department of Chemistry, University of Crete, P.O. Box 2208, 70013, Heraklion, Crete, Greece. Notes

The authors declare no competing financial interest.

■

REFERENCES

(1) Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids; Oxford University Press: New York, 1989. (2) Ewald, P. P. Die Berechnung Optischer und Elektrostatischer Gitterpotentiale. Ann. Phys. 1921, 369, 253−287. (3) Darden, T.; York, D.; Pedersen, L. Particle Mesh Ewald: An N Log(N) Method for Ewald Sums in Large Systems. J. Chem. Phys. 1993, 98, 10089−10092. (4) Eastwood, J. W.; Hockney, R. W.; Lawrence, D. N. P3M3DP − The Three-Dimensional Periodic Particle-Particle/Particle-Mesh Program. Comput. Phys. Commun. 1980, 19, 215−261. K

DOI: 10.1021/jp510612w J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B (5) Cisneros, G. A.; Karttunen, M.; Ren, P.; Sagui, C. Classical Electrostatics for Biomolecular Simulations. Chem. Rev. 2014, 114, 779− 814. (6) Brooks, B. R.; Bruccoleri, R. E.; Olafson, B. D.; States, D. J.; Swaminathan, S.; Karplus, M. CHARMM: A Program for Macromolecular Energy, Minimization, and Dynamics Calculations. J. Comput. Chem. 1983, 4, 187−217. (7) Steinbach, P. J.; Brooks, B. R. New Spherical-Cutoff Methods for Long-Range Forces in Macromolecular Simulation. J. Comput. Chem. 1994, 15, 667−683. (8) Wolf, D.; Keblinski, P.; Phillpot, S. R.; Eggebrecht, J. Exact Method for the Simulation of Coulombic Systems by Spherically Truncated, Pairwise r−1 Summation. J. Chem. Phys. 1999, 110, 8254−8282. (9) Wolf, D. Reconstruction of NaCl Surfaces from a Dipolar Solution to the Madelung Problem. Phys. Rev. Lett. 1992, 68, 3315−3318. (10) Demontis, P.; Spanu, S.; Suffritti, G. B. Application of the Wolf Method for the Evaluation of Coulombic Interactions to Complex Condensed Matter Systems: Aluminosilicates and Water. J. Chem. Phys. 2001, 114, 7980−7988. (11) Zahn, D.; Schilling, B.; Kast, S. M. Enhancement of the Wolf Damped Coulomb Potential: Static, Dynamic, and Dielectric Properties of Liquid Water from Molecular Simulation. J. Phys. Chem. B 2002, 106, 10725−10732. (12) Fennell, C. J.; Gezelter, J. D. Is the Ewald Summation Still Necessary? Pairwise Alternatives to the Accepted Standard for LongRange Electrostatics. J. Chem. Phys. 2006, 124, 234104. (13) Yonezawa, Y. A Long-Range Electrostatic Potential Based on the Wolf Method Charge-Neutral Condition. J. Chem. Phys. 2012, 136, 244103. (14) Yonezawa, Y. Electrostatic Properties of Water Models Evaluated by a Long-Range Potential Based Solely on the Wolf Charge-Neutral Condition. Chem. Phys. Lett. 2013, 556, 308−314. (15) Chen, Y.-G.; Weeks, J. D. Local Molecular Field Theory for Effective Attractions Between like Charged Objects in Systems with Strong Coulomb Interactions. Proc. Natl. Acad. Sci. U.S.A. 2006, 103, 7560−7565. (16) Angoshtari, A.; Yavari, A. Convergence Analysis of the Wolf Method for Coulombic Interactions. Phys. Lett. A 2011, 375, 1281− 1285. (17) Hansen, J. S.; Schrøder, T. B.; Dyre, J. C. Simplistic Coulomb Forces in Molecular Dynamics: Comparing the Wolf and Shifted-Force Approximations. J. Phys. Chem. B 2012, 116, 5738−5743. (18) Heyes, D. M. Electrostatic Potentials and Fields in Infinite Point Charge Lattices. J. Chem. Phys. 1981, 74, 1924−1929. (19) Mendoza, F. N.; Lòpez-Lemus, J.; Chapela, G. A.; Alejandre, J. The Wolf Method Applied to the Liquid-Vapor Interface of Water. J. Chem. Phys. 2008, 129, 024706. (20) Takahashi, K. Z. Truncation Effects of Shift Function Methods in Bulk Water Systems. Entropy 2013, 15, 3249−3264. (21) Markland, T. E.; Manolopoulos, D. E. A Refined Ring Polymer Contraction Scheme for Systems with Electrostatic Interactions. Chem. Phys. Lett. 2008, 464, 256−261. (22) Fanourgakis, G. S.; Markland, T. E.; Manolopoulos, D. E. A Fast Path Integral Method for Polarizable Force Fields. J. Chem. Phys. 2009, 131, 094102. (23) Nymand, T. M.; Linse, P. Ewald Summation and Reaction Field Methods for Potentials with Atomic Charges, Dipoles, and Polarizabilities. J. Chem. Phys. 2000, 112, 6152−6160. (24) Abascal, J. L. F.; Vega, C. A General Purpose Model for the Condensed Phases of Water: TIP4P/2005. J. Chem. Phys. 2005, 123, 234505. (25) Vega, C.; Abascal, J. L. F. Simulating Water with Rigid NonPolarizable Models: A General Perspective. Phys. Chem. Chem. Phys. 2011, 13, 19663−19688. (26) Caldwell, J. W.; Kollman, P. A. Structure and Properties of Neat Liquids Using Nonadditive Molecular Dynamics: Water, Methanol, and N-Methylacetamide. J. Phys. Chem. 1995, 99, 6208−6219. (27) Fanourgakis, G. S.; Xantheas, S. S. Development of Transferable Interaction Potentials for Water. V. Extension of the Flexible,

Polarizable, Thole-Type Model Potential (TTM3-F, v. 3.0) to Describe the Vibrational Spectra of Water Clusters and Liquid Water. J. Chem. Phys. 2008, 128, 074506. (28) Thole, B. T. Molecular Polarizabilities Calculated with a Modified Dipole Interaction. Chem. Phys. 1981, 59, 341−350. (29) Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; DiNola, A.; Haak, J. R. Molecular Dynamics with Coupling to an External Bath. J. Chem. Phys. 1984, 81, 3684−3690. (30) Andersen, H. C. Rattle: A “velocity” Version of the Shake Algorithm for Molecular Dynamics Calculations. J. Comput. Phys. 1983, 52, 24−34. (31) Martyna, G. J.; Tuckerman, M. E.; Tobias, D. J.; Klein, M. L. Explicit Reversible Integrators for Extended Systems Dynamics. Mol. Phys. 1996, 87, 1117−1157. (32) Swope, W. C.; Andersen, H. C.; Berens, P. H.; Wilson, K. R. A Computer Simulation Method for the Calculation of Equilibrium Constants for the Formation of Physical Clusters of Molecules: Application to Small Water Clusters. J. Chem. Phys. 1982, 76, 637−649. (33) Bussi, G.; Donadio, D.; Parrinello, M. Canonical Sampling Through Velocity Rescaling. J. Chem. Phys. 2007, 126, 014101. (34) Berendsen, H. J. C.; Grigera, J. R.; Straatsma, T. P. The Missing Term in Effective Pair Potentials. J. Phys. Chem. 1987, 91, 6269−6271. (35) Smith, D. E.; Dang, L. X. Computer Simulations of NaCl Association in Polarizable Water. J. Chem. Phys. 1994, 100, 3757−3766. (36) Aragones, J. L.; Sanz, E.; Vega, C. Solubility of NaCl in Water by Molecular Simulation Revisited. J. Chem. Phys. 2012, 136, 244508. (37) Fanourgakis, G. S.; Tipparaju, V.; Nieplocha, J.; Xantheas, S. S. An Efficient Parallelization Scheme for Molecular Dynamics Simulations with Many-Body, Flexible, Polarizable Empirical Potentials: Application to Water. Theor. Chem. Acc. 2007, 117, 73−84. (38) Todorov, I. T.; Smith, W.; Trachenko, K.; Dove, M. T. DL_POLY_3: New Dimensions in Molecular Dynamics Simulations Via Massive Parallelism. J. Mater. Chem. 2006, 16, 1911−1918. (39) Press, W. H.; Flannery, B. P.; Teukolsky, S. A.; Vetterling, W. T. Numerical Recipes in Fortran: The Art of Scientific Computing, 2nd ed.; Cambridge University Press: 1992. (40) McLaughlin, K.; Cioce, C. R.; Pham, T.; Belof, J. L.; Space, B. Efficient Calculation of Many-Body Induced Electrostatics in Molecular Systems. J. Chem. Phys. 2013, 139, 184112. (41) Takahashi, K. Z.; Narumi, T.; Yasuoka, K. Cutoff Radius Effect of the Isotropic Periodic Sum and Wolf Method in Liquid-Vapor Interfaces of Water. J. Chem. Phys. 2011, 134, 174112.

L

DOI: 10.1021/jp510612w J. Phys. Chem. B XXXX, XXX, XXX−XXX