Coarse-Graining in Interaction Space: A Systematic Approach for

Mar 27, 2008 - On the Investigation of Coarse-Grained Models for Water: Balancing Computational Efficiency and the Retention of Structural Properties...
0 downloads 0 Views 470KB Size
Download PDF
J. Phys. Chem. B 2008, 112, 4711-4724

4711

Coarse-Graining in Interaction Space: A Systematic Approach for Replacing Long-Range Electrostatics with Short-Range Potentials Sergei Izvekov,† Jessica M. J. Swanson,† and Gregory A. Voth* Center for Biophysical Modeling and Simulation and Department of Chemistry, UniVersity of Utah, 315 South 1400 East, Room 2020, Salt Lake City, Utah 84112-0850 ReceiVed: October 25, 2007; In Final Form: December 24, 2007

A short-range effective potential for long-range electrostatic interactions in homogeneously disordered condensed phase systems has been determined with a novel approach to coarse-graining in interaction space. As opposed to coarse-graining the system resolution, this approach “coarsens” the system’s interactions by mapping multiple configurations of an accurate long-range atomistic potential onto a more efficient, shortrange effective potential with a force-matching (FM) method. Developing an empirical potential in this manner is fundamentally different from existing strategies because it utilizes condensed-phase (as opposed to gasphase) atomistic interactions to determine general pair potentials defined on distance meshes (as opposed to fitting predetermined functional forms). The resulting short-range (∼10 Å) effective potential reproduces structural, dynamical, and many thermodynamic properties of liquid water, ions in water, and hydrophobes in water, with unprecedented accuracy. The effective potential is also shown to be transferable to a nonaqueous molten salt system. With continued development, such effective potentials may provide an accurate and highly efficient alternative to Ewald-based long-range electrostatics methods.

I. Introduction The past three decades have marked an era of increasing computational resources and landmark breakthroughs in computer modeling and simulation. Despite these advances, the size and time scales relevant to many important biological and other soft matter phenomena are still orders of magnitude beyond current simulation capabilities. Developing new computational methodologies that allow simulations to explore systems on much larger length and time scales will usher in a new era of collaboration and communication between theory and experiment, making this arguably the most important challenge facing computational chemistry and biophysics today. Here, a new strategy is presented for surmounting this challenge. It relies on the determination of optimal short-range effective potentials for condensed phase systems that can significantly increase simulation efficiency while maintaining structural, dynamical, and thermodynamic accuracy. One of the main computational bottlenecks for simulating complex biomolecular and soft matter systems is the accurate treatment of long-range electrostatic interactions. Accounting for electrostatic interactions in condensed phase systems is challenging because the Coulomb sum is conditionally convergent, and because of the cost of including all pair interactions, such that a brute-force calculation scales as the square of the number of atoms in the system, O(N2). Great effort has gone into developing accurate and efficient electrostatic treatments, resulting in explicit (e.g., Ewald, lattice-sum, and potential truncation),1-4 implicit (i.e., continuum dielectric),5 and mixed explicit-implicitmethods(e.g.,reactionfieldandfastmultipole).6-8 The Ewald method, initially used to account for electrostatics in infinite periodic systems (i.e., the classic Madelung problem), implements a clever mathematical trick that splits the condition* Corresponding author. E-mail: [email protected]. † These authors contributed equally to this work.

ally convergent Coulomb potential into two convergent summations. It is generally accepted as the most accurate, though most expensive, electrostatics treatment for condensed-phase and biomolecular simulations. Lattice sum approaches, such as particle mesh Ewald (PME; see ref 2 and references therein), improve the O(N2) Ewald efficiency considerably; the fastest scales as O(N ln N). They are still, however, considerably slower and more difficult to parallelize than cutoff approaches, in which the potential is simply truncated at a given cutoff length, and which scale linearly with system size O(N). Although efficiency once made truncated potentials the method of choice, there now exists substantial literature demonstrating that they can introduce significant artifacts for polar fluids,9-11 highly charged systems,12,13 and solvated macromolecules.14-19 More sophisticated truncation schemes have been developed employing various functional forms that smoothly shift or switch the electrostatic potential, and sometimes forces, to zero at the cutoff distance.3,9,20,21 Several of the atom-based, force-shifted (FS) forms have been successfully employed in stable simulations of protein and DNA systems on the nanosecond or subnanosecond time scale.20,22,23 Nevertheless, it has been shown that even these ideal smooth cutoff approaches lead to artifacts such as artificially increased structure in liquid water10,11 as well as interfacial and transport artifacts at the solvent-vapor boundary and in membrane simulations.24 Furthermore, accelerated sampling techniques have shown that twin-range and abrupt cutoff techniques lead to long-time scale artifacts that are hidden on the nanosecond time scale,18,25 suggesting that FS cutoffs may result in similar artifacts in time. Over two decades of complications with cutoff approaches have culminated in the widely accepted conclusion that realistic models of condensed-phase systems require the explicit treatment of long-ranged electrostatic interactions via an Ewaldlike method.9,10,16,22 More recently, however, there has been a resurgence in the argument that effective electrostatic interac-

10.1021/jp710339n CCC: $40.75 © 2008 American Chemical Society Published on Web 03/27/2008

4712 J. Phys. Chem. B, Vol. 112, No. 15, 2008

Izvekov et al.

Figure 1. Comparison of various electrostatic potentials highlighting the different damping rates of the FS, FM effective, and DFS potentials. Only the FM effective potential for Rcut ) 0.635 nm (filled circles) falls off more quickly than the strongly damped DFS2 (empty circles). Conversely, all of the FM effective potentials fall off more quickly than the FS potential for Rcut ) 1.0 nm (thin dashed). The energy units are “nm-1”, i.e., U(r ) 1.0 nm) ) 1.0 if U(r) is the Coulomb potential, U(r) ) r-1 (1 nm-1 ) 138.935 kJ/mol).

tions in condensed-phase systems decay considerably faster than r-1,26-28 and several efforts have been directed toward constructing new damped and truncated electrostatic potentials.28-34 The “damped shifted” (DS) potentials, for example,28,31,32,34 were developed in two steps: first, neutralizing the charge within any given cutoff sphere, and second, damping the potential to effectively decrease the range of electrostatic interactions. The first step is equivalent to shifting the potential to zero at the cutoff. The damping function used in the second step could, in principle, take any form. Thus far, however, a complementary error function has been used in mimicry of its use in the Ewald summation. A general form for the DS potential can be expressed as

VDS(r) )

(

erfc(Rr) + P(r - Rcut) r

)

(1)

where P(x) is a polynomial constructed such that VDS(r), and, in the case of the damped force-shifted (DFS) potential, its first derivative approaches zero at Rcut. In the local molecular field theory, introduced by Weeks and co-workers, a DS potential is augmented by a slowly varying component, analogous to the reciprocal summation in the Ewald method.29,30,33 There have also been attempts to use different damping functions, such as the Yukawa potential.35 The two most successful truncated potentials to date are the FS potential, given by eq 1 with R ) 0, and DFS potentials with a properly chosen value of R * 0. As will be shown, the FS potential reproduces the Ewald force distributions more accurately than the strongly DFS potentials (R ) 2.0 nm-1). They also, however, cause structural artifacts and substantially underestimated system densities and energies (even with relatively large cutoffs). The DFS potentials, on the other hand, reproduce structural and certain macroscopic properties more accurately,32 but do so in a system- and property-dependent manner (meaning the ideal damping length seems to be dependent on the condensed-phase environment and the macroscopic property of interest).31 The DFS potentials developed for liquid water,31,34 for example, will be shown herein to be over-damped, resulting in many properties (e.g., configurational energies and densities) that are only slightly better than those produced with the FS potential, and consequently worse than those produced with the Ewald method. Other properties, such

Figure 2. Differences between the Coulomb and various electrostatic potentials 1/r - U(r) (a) and corresponding forces 1/r2 - F(r) (b). The seemingly subtle potential differences in panel a between the fit DFS1 potential (filled circles) and the FM effective potential for Rcut ) 1.0 nm (thin solid) result in more striking force differences in panel b. Likewise, the decreasing “bend” effect at the cutoff distance with increasing cutoff length is more clearly demonstrated in panel b, yielding nearly linear profiles for the FM effective potentials for Rcut ) 1.2 nm (filled triangles) and 1.323 nm (open squares). The energy units are as in Figure 1.

as the diffusion of charged species, are significantly affected by the over-damped DFS potential. The recently proposed additional potential shifting terms32 may alleviate energetic deficiencies, but will have no impact on the electrostatic forces. Thus, most properties, which rely on the forces, can only be improved by increasing the damping length (reducing R in eq 1), which would in turn adversely affect structural properties. These results point to one of the main limitations of the FS and DFS potentials developed to date: the use of predetermined functional forms that have been parametrized to reproduce preselected macroscopic properties.30,31,34 Moreover, it will be shown herein that a pre-defined functional form (as in eq 1) has a limited ability to reproduce long-range electrostatic interactions. The present work seeks to overcome these limitations. Here a novel and systematic “inverse” approach for defining effective potentials is presented and used to determine a shortrange effective potential for long-range electrostatic interactions in disordered homogeneous media such as liquid water. This method is conceptually akin to the recently developed multiscale coarse-graining method (MS-CG)36-40 in that it employs a generalization of the force-matching (FM) algorithm36,41 to determine the effective pair potential that best reproduces the forces from a reference simulation. It is also a unique twist on coarse-graining because only the interactions are “coarsened” while the full atomistic resolution is maintained. Central to the success of the MS-CG method is a statistical mechanical framework that accounts for two- and three-body correlations.38 Developing an effective short-range potential in this manner is fundamentally different from traditional strategies that use the gas phase as a reference state and fit a potential function with some predetermined form to ab initio, cluster, and/or thermodynamic data. Instead, a homogeneous condensed phase is taken as the reference state from which the actual atomistic interactions are distributed over a subset of short-range effective pair potentials defined on distance meshes. Variationally minimizing the functional that averages instantaneous differences between

Coarse-Graining in Interaction Space

J. Phys. Chem. B, Vol. 112, No. 15, 2008 4713

TABLE 1: Coefficients ai from the Least-Squares Fit of FFM(r) to Eq 6 for Rcut ) 1.0 nm (Second Column) Using N ) 7 and for Rcut ) 1.2 nm (Third Column) Using N ) 10 (Units in Atomic Units) i

ai

ai

0 1 2 3 4 5 6 7 8 9 10

-0.165477570871E-03 0.288823451703E-03 -0.122247561247E-03 0.963712701767E-05 0.251954672874E-06 -0.735796273353E-07 0.353601771929E-08 -0.525765995765E-10 0.0 0.0 0.0

-0.547587014180E-04 -0.777061023641E-05 0.649658451885E-04 -0.417496679930E-04 0.926924324623E-05 -0.107542095070E-05 0.710779773104E-07 -0.261040455982E-08 0.439931939572E-10 -0.422656965444E-14 -0.656184357691E-14

the atomistic and effective forces over the atomistic ensemble guarantees that the best possible effective pair potential is found. The same systematic approach might be used to introduce corrections for physical effects that are not well described by the homogeneous reference system, such as the inhomogeneous nature of interfaces. Remarkably, the short-range effective potential determined with this method for a particular model of pure liquid water and at a single temperature is not only charge independent (i.e., consistent for all atomistic electrostatic interactions once scaled by the atomistic charges) but also transferable to other homogeneously disordered (liquid state) systems and/or to other temperatures. The ability of this effective potential to reproduce structural, dynamic, and many thermodynamic properties is demonstrated for liquid water, ions in water, and hydrophobes in water. It is also shown to be transferable to a significantly different system such as molten NaCl. The rest of the manuscript is organized as follows. In section II, the theory and methodology used to determine effective potentials with the FM method are described. In section III, the FM effective potential obtained from atomistic simulations of SPC/E water is described, and the resulting water properties are compared to those for the FS and DFS potentials. The FM effective potential is then tested for its transferability to different models (SPC/E to TIP3P) and temperatures. In addition to structural, dynamic, and thermodynamic properties, the FM effective potential’s reproduction of higher order correlations is also discussed. The transferability analysis is extended to different, more heterogeneous systems, including an ion in water, aqueous ion-ion interactions, and solvated hydrophobes. Finally, the effective electrostatic interaction and the results it provides are tested against a full Ewald MD simulation for molten NaCl. This is done by directly force-matching molten NaCl and comparing these results to those obtained by simply transferring the effective FM interaction from liquid water. In section IV, the main conclusions from this work are presented. II. Theory And Methods A. FM Effective Potentials. The FM method41 has primarily been used to develop effective potentials for coarse-grained molecular systems, those in which the number of atomistic degrees of freedom have been reduced in order to increase simulation efficiency.36-40 In the present work, the same FM algorithm is used to “coarse-grain” a computationally expensive long-range atomistic potential into a more efficient short-range “effective” potential while maintaining the full atomistic resolution. Given a system containing N atoms, the Cartesian coordinates of which are represented by the vector b r N, the FM method

Figure 3. Comparison of DFS1 potential with R ) 1.268 nm-1 (filled circles) obtained through least-squares fit to the FM effective potential for Rcut ) 1.0 nm (thin solid), FS potential with Rcut ) 1.0 nm (thin dashed), and DFS2 potential with R ) 2.0 nm-1 (empty circles). The energy units are as in Figure 1.

determines an optimal effective force field between these atoms by minimizing the quadratic residual,

χ [F B (b r )] ) 2

eff

N

1 3N

〈∑ N

|F B eff bN ) - B Fi(r bN )|2 i (r

i)1

〉

(2)

This residual compares the total force acting on each atom i r N ), and the according to the original (reference) force field, B Fi( b eff N r ), and averages this effective short-range force field, B F i (b difference over the configurations sampled with the original force field in the canonical ensemble at temperature T. Minimizing eq 2 with respect to the coefficients that describe B F eff rN) i (b produces an effective force field that is the best possible approximation (in a variational least-squares sense42) to the original more complicated force field. It should be readily apparent that when B F eff r N ) is described by functions that span i (b the full range of the reference force field, the FM procedure will recover the reference force field exactly, and the residual in eq 2 will be zero. For example, if the reference potential was the Ewald lattice summation, and B F eff r N ) was described i (b by pairwise radial functions over all atoms up to the Ewald real space cutoff and k-space vectors thereafter, then the exact Ewald potential would be recovered. The goal of this work, however, is to find a more efficient r N ), may be short-range effective potential. In general, B F eff i (b many-bodied in character and therefore be described by manybody functions. In the present work it is assumed that the effective electrostatic potential is a pairwise-additive, radially symmetric sum over a subset of atoms in the system, yielding a force on atom i given by NC

r NC) ) B F eff i (b

fij(rij)n bij ∑ j*i

(3)

where NC is the number of atoms within a certain cutoff distance. The pairwise force fij(r bij)n bij can be further assumed to be of the form f Rβ(rij) acting between atoms i and j of type R and β, a function of their separation |r bij | ) |b ri - b rj | and acting along rij/rij. Thus, the FM method provides the the unit vector b nij ) b optimal approximation to the original force field within this class of potentials. If the effective forces, f Rβ(r), are defined by functions that depend linearly on their parameters, e.g., numerical basis functions such as discrete delta, linear spline, or cubic spline

4714 J. Phys. Chem. B, Vol. 112, No. 15, 2008

Izvekov et al.

Figure 4. Time evolution of ∆F defined in eq 7 (upper row) and [∆FI]conf max defined in eq 8 (lower row) for FS potential (left), DFS2 potential (middle), and FM effective potential (right). The cyan results are for Rcut ) 1.0 nm, and the black are for Rcut ) 1.2. The white lines show the running time averages of each quantity.

Figure 5. The same data as in Figure 4 for the ions in simulations of four Cl- ions solvated in 512 SPC water molecules.

functions, then eq 2 can be reduced to an overdetermined set of linear equations, and a unique solution can be determined. In the present work, f Rβ(r) is described by cubic splines that connect k points on a discretized mesh {rk} up to point rkmax ) Rcut, thus preserving continuity of the forces and their first two derivatives across the mesh points.36,43 Other basis functions having a finite interaction length scale (i.e., cutoff) are certainly possible,40 but the spline basis has proven to be a highly accurate way to represent the numerical force data. The spline representation also facilitates ensuring that the forces are continuous at the cutoff distance, a condition that is crucial for energy conservation and simulation stability. This is accomplished by forcing the f Rβ(r) functions and their derivatives to be zero at Rcut with additional constraints,

f Rβ(Rcut) ) 0 and

d2f Rβ(Rcut)/dr2 ) 0

(4)

in the FM equations. Ensuring that dfRβ(Rcut)/dr ) 0 is already a requirement of the spline representation.

Generally, there are too many configurations sampled in an ensemble (i.e., too many equations) to solve the variationally derived, overdetermined set of linear FM equations41 with conventional techniques (e.g., single-value decomposition or QR algorithms44). In the present work, a previously described block averaging approach36,37 was used to solve the FM equations in sets (blocks) of configurations. Averaging over block solutions results in an excellent approximation to the exact solution as long as the fluctuations in the residual curvature and those in force field parameters are uncorrelated.39,40 Correlations may occur in complex systems when different blocks correspond to different parts of conformational phase space, but this correlation can be monitored and minimized by swapping or randomizing configurations between blocks. For the homogeneous condensed phase systems studied in the present work, swapping was not necessary, as indicated by roughly equivalent solutions for all the blocks in the ensemble. In summary, the FM procedure is implemented as follows. First, the reference system is selected, and atomistic trajectories are generated. Second, the desired reference forces (e.g., the Ewald electrostatic forces) are collected from each sampled configuration. Third, the FM method is used to fit a chosen

Coarse-Graining in Interaction Space

J. Phys. Chem. B, Vol. 112, No. 15, 2008 4715

TABLE 2: Statistical Analysis of Force Deviations for Data Shown in Figures 4 and 5 〈∆Q〉 (〈δ∆Q〉)e (a.u.)

model, Rcut (nm) potential

〈∆F〉 (〈δ∆F〉)a

[∆F]tmax b

c 〈[∆FI]conf max 〉

SPC/E, 1.0 FS DFS2 FM

0.031 (0.001) 0.032 (0.002) 0.016 (0.002)

0.039 0.044 0.028

0.130 0.157 0.094

2.000 2.944 1.611

1.84 (2.32) 1.94 (2.44) 1.94 (2.44)

SPC/E, 1.2 FS DFS2 FM

0.021 (0.001) 0.031 (0.002) 0.010 (0.002)

0.028 0.043 0.019

0.089 0.152 0.061

1.060 2.767 0.842

2.14 (2.70) 2.28 (2.87) 2.28 (2.87)

TIP3P, 1.0 FS DFS2 FM

0.034 (0.002) 0.036 (0.003) 0.018 (0.003)

0.044 0.050 0.032

0.144 0.177 0.107

2.591 3.460 1.776

1.81 (2.29) 1.91 (2.41) 1.91 (2.41)

TIP3P, 1.2 FS DFS2 FM

0.024 (0.001) 0.035 (0.003) 0.012 (0.002)

0.032 0.049 0.022

0.100 0.171 0.069

1.964 3.349 1.389

2.12 (2.68) 2.26 (2.84) 2.26 (2.84)

Ion, SPC, 1.0 FS DFS2 FM

0.088 (0.027) 0.110 (0.034) 0.074 (0.024)

0.272 0.359 0.257

0.189 0.239 0.163

9.394 13.793 6.224

1.61 (0.69) 1.26 (0.96) 1.26 (0.96)

Ion, SPC, 1.2 FS DFS2 FM

0.061 (0.019) 0.106 (0.033) 0.049 (0.016)

0.208 0.349 0.189

0.131 0.231 0.108

5.257 13.632 5.454

1.61 (0.69) 1.26 (0.96) 1.26 (0.96)

t d [[∆FI]conf max ]max

a Average (and standard deviation) of the relative deviation of the Euclidean norm per configuration. b Maximum relative deviation of the Euclidean norm per ensemble. c Time average of the maximum force deviation per configuration. d Ensemble maximum force deviation. e Ensemble average (and standard deviation) of the uncompensated charge within the cutoff sphere surrounding atom with maximum force deviation.

functional form (e.g., pairwise radial cubic spline meshes) to the reference forces. Finally, the FM effective potential for each unique pair of possible interactions, hereafter denoted as UFM(r), can be obtained by integrating the force profile, FFM(r), and then shifting the result to zero at Rcut. As described above, the FM methodology can be used to, in effect, “coarse-grain” the long-range electrostatic interactions into shorter-range interactions. While this methodology alone could prove to be both powerful and insightful, its impact might be rather limited if every condensed phase system of interest required a new and unique application of the FM algorithm in order to develop a system-specific short-range effective potential. Here, we have gone a step further to explore whether or not there exists an atom type (i.e., charge) independent and transferable short-range representation of electrostatics, i.e., one that can be determined once for a single pair interaction (e.g., O-H) in a single system (e.g., water) and then transferred without new force-matching to other pair interactions and other liquid-state systems. Finding such an effective potential might seem impossible, but evidence is presented below that this concept may not be as unlikely as one might assume. At a minimum, a reasonably accurate and transferable short-range effective potential for condensed phase electrostatics could serve as a “zeroth-order” description to which corrections for physical effects not well described by the zeroth-order approximation can be included. Indeed, the f Rβ(r) solutions (cf. eq 3) are found to be independent of charge for the systems studied in the present work, meaning that the FM effective electrostatic force between unit charges, FFM(r) ) f Rβ(r)/qRqβ, is equiValent for all atom types R and β. Moreover, this effective potential (and force) results in remarkable accuracy when transferred to other aqueous systems and to a molten ionic liquid, though corrections are expected to be necessary for inhomogeneous environments such as interfaces and macromolecules.

B. Determining an Effective Potential for Electrostatic Interactions. In order to determine the optimal short-ranged effective pair potential for long-ranged electrostatic interactions in condensed-phase environments, we began by force-matching liquid water at ambient conditions. Although ab initio simulations could and have provided the atomistic information,41,45 using an empirical potential permitted a clean decomposition of Lennard-Jones and electrostatic nonbonded interactions. Only the latter were force-matched. 1. Reference Simulation. In the reference simulation, bulk water was modeled with the SPC/E potential using a rigid molecular geometry. To ensure that our results were independent of periodic artifacts, two system sizes were simulated; the first had 215 water molecules in a cubic box with a length of 1.86 nm, and the second held 1720 water molecules and had a box length of 3.72 nm. The FM potentials presented below were identical for the two box sizes. In each case, the system was simulated for 2 ns in the constant NVT ensemble at T ) 300 K using the Nose´-Hoover thermostat with a relaxation time of 0.2 ps. The electrostatic interactions were calculated with the PME summation with a convergence tolerance of 10-5. Configurations were recorded every picosecond for a total of 2000 reference configurations. For each recorded configuration, the forces were then recalculated with the Ewald summation using a smaller tolerance of 10-6, which corresponds to the following Ewald parameters: R ) 4.0514 nm-1 and kmax ) 8 for the smaller system and kmax ) 15 for the larger system, where kmax is a maximal number of reciprocal vectors in each positive direction. The two different electrostatics treatments, which were chosen for efficiency in the simulation and accuracy in the force calculations, produced nearly identical forces with relative errors well below the tolerance used in the FM procedure (10-5). Therefore, the sampled configurational ensemble was representative of the reference atomistic potential.

4716 J. Phys. Chem. B, Vol. 112, No. 15, 2008

Figure 6. (a) Comparison of the oxygen-oxygen RDFs for SPC/E water (T ) 300 K) simulated with various electrostatic potentials. Shown are the reference Ewald (cyan solid line), the FS (Rcut ) 1.0 nm; dashed line), and the FM effective (Rcut ) 1.0 nm; black solid line) potentials. (b) Same results shown in panel a, highlighting the deviations with a higher resolution.

Izvekov et al.

Figure 8. Same RDFs shown in Figure 6 for the TIP3P water model.

The FM effective forces could not be determined for r < rcore, where rcore ) 0.148 nm was the smallest O-H intermolecular separation observed in the reference simulation. One might expect the FFM(r) to be similar to the Coulomb force for core separations. Indeed, a least-squares fit of FFM(r) to γ/rη over the interval 0.148 < r < 0.265 nm resulted in γ ) 1.0262 and η ) 2.0236. Thus, the FM force can be extrapolated onto the core region as

FFM(r < rcore) )

1 1 + FFM(rcore) r2 r 2core

(5a)

γ rη

(5b)

or alternatively

FFM(r < rcore) )

using a core value of rcore ) 0.148 nm for all cutoffs and the above values of γ and η. Although both extrapolations resulted in identical simulated properties for SPC/E water, the rest of the simulations reported here used eq 5a. III. Results Figure 7. (a) Comparison of the RDFs shown in Figure 6 with those obtained with the FS (open circles) and FM effective (filled circles) potentials with Rcut ) 0.9 nm. (b) Comparison of the RDFs shown in Figure 6 with those obtained with the strongly damped DFS2 potential with Rcut ) 1.0 nm (dashed line).

2. FM Procedure. After obtaining the atomistic ensemble and corresponding Ewald forces, the FM procedure was used to determine the O-O, O-H, and H-H effective electrostatic interactions. All calculations used a distance mesh defined by 0.005 nm bins and a block size of 4. The three resulting effective potentials were identical within tolerance of the FM procedure, indicating that FFM(r) ) f Rβ(r)/qRqβ is indeed charge-independent. Five different cutoffs were used to calculate FFM(r), Rcut ) 0.635, 0.9, 1.0, 1.2, and 1.323 nm (the fractional units are simply the result of using atomic units in the original calculations). Atomic units are also used in the formulas such that the Coulomb potential and force adopt the form 1/r and 1/r2, respectively. However, the distance in the figures shown later are in nanometers for convenience and clarity.

A. FM Effective Potential. 1. Nature of the FM EffectiVe Potentials. The FM effective potentials, UFM(r), for the five different cutoff lengths are shown in Figure 1. The results from the large and small reference simulations were identical within the FM tolerance, proving that our results were independent of the box size and periodic artifacts. Clearly, all UFM(r) fall to zero more rapidly than r-1 because of the imposed boundary conditions (eq 4). Interestingly, the effective potentials seem to evolve with increasing cutoff length, as shown by the differences between the Coulomb potential and the FM potentials (1/r - UFM(r)) and respective data for forces (1/r2 - FFM(r)) in Figure 2. For all cutoff distances, the force difference is almost linear for small separations but then bends as it approaches the cutoff. This bend becomes less pronounced as Rcut increases, resulting in nearly linear differences for the FM effective potentials with Rcut > 1.2 nm. The decay rate of the FM potentials accelerates with distance. For example, fitting FFM(r) for Rcut ) 1.323 nm to γ/rη over the range 0.5 < r < 0.8 nm resulted in η ) 2.7. The same fit

Coarse-Graining in Interaction Space

J. Phys. Chem. B, Vol. 112, No. 15, 2008 4717

TABLE 3: Liquid-Phase Properties of the SPC/E and TIP3P Models from Simulations Run with Various Electrostatic Treatments at T ) 300 K and P ) 1 bara Fb (kg/m3)

Dc (10-9 m2/s)

Ud (kJ/mol)

RPe (10-4 1/K)

κTf (10-5 1/bar)

SPC/E: Ewald FS, 1.0 FS, 0.9 FM, 1.0 FM, 0.9 FM, 0.635 DFS1, 1.0g DFS2, 1.0h

998.4 (11.2) 994.4 (11.2) 991.8 (11.3) 997.6 (11.2) 997.7 (11.2) 996.3 (11.5) 997.3 (11.2) 995.1 (10.9)

2.50 (0.01) 2.85 (0.01) 2.65 (0.03) 2.50 (0.01) 2.55 (0.03) 2.68 (0.05) 2.48 (0.01) 2.70 (0.01)

-46.60 (0.30) -45.01 (0.30) -44.81 (0.30) -46.65 (0.30) -46.74 (0.30) -46.84 (0.30) -46.10 (0.30) -45.50 (0.30)

5.2 5.0 4.9 5.1 5.0 4.3 5.1 4.2

4.6 4.7 4.9 4.6 4.7 4.9 4.7 4.4

TIP3P: Ewald FS, 1.0 FS, 0.9 FM, 1.0 FM, 0.9 FM, 0.635

1014.4 (11.5) 1010.6 (11.9) 1007.7 (12.3) 1015.3 (11.5) 1015.2 (11.9) 1016.3 (12.0)

5.30 (0.01) 5.51 (0.01) 5.55 (0.04) 5.30 (0.01) 5.38 (0.04) 5.42 (0.07)

-41.14 (0.28) -39.77 (0.28) -39.40 (0.28) -41.21 (0.28) -41.23 (0.28) -41.32 (0.28)

8.8 9.0 9.0 8.8 8.8 7.6

4.8 5.0 5.5 4.8 5.0 5.1

model: potential, Rcut (nm)

a Standard deviations in parenthesis. b Density. c Self-diffusion coefficient. d Average configuration energy. e Thermal expansion coefficient. Isothermal compressibility. g The DFS1 potential (R ) 1.268 nm-1) that was fit to the FM effective potential. h The DFS2 potential (R ) 2.0 nm-1).

f

over larger distances, 0.8 < r < Rcut and 0.9 < r < Rcut, yielded η ) 5.3 and η ) 6.2, respectively. The last fit supports the notion that electrostatic forces effectively decay as r-6, as observed in crystal-like polarizable media.27,28 A more accurate representation of FFM(r) can be obtained using the polynomial expansion

FFM(r g rcore) )

1 r2

∑ i)0

airi

(6)

As previously described, FFM(r) for r e rcore can be represented as eq 5a or 5b. The coefficients ai that best fit FFM(r) for Rcut ) 1.0 nm using N ) 7 and for Rcut ) 1.2 using N ) 10 are given in Table 1. These values of N were found to be the minimal polynomial powers necessary to reproduce the original numerical FM results. 2. Comparisons to Other Potentials. In order to clarify how UFM(r) is related to the DS and DFS potentials, the FM potential for Rcut ) 1.0 nm was fit to eq 1 by varying the parameters R and β ) P(x). The resulting damping parameter R was 1.268 nm-1, regardless of a variable or fixed value of β ) 0. Figure 3 compares UFM(r) to this damped potential, which will be referred to as DFS1. Note that the damped potential and force can be shifted to avoid the small jump at the cutoff (shown more clearly in Figure 2) in order to improve the energy conservation without affecting the macroscopic properties in the simulation. Interestingly, most authors31,32,34 have used a faster decaying DFS potential for liquid water, selecting R ) 2.0 nm-1, which will be referred to as the DFS2 potential, to better match Ewald simulation densities. As shown in Figures 1-3, the DFS2 potential is more damped than all of the effective potentials except for UFM(r) for Rcut ) 0.635 nm. Although the DFS2 potential leads to surprisingly good structural properties for liquid water (as shown below), it shows little improvement over the FS potential in density and energetic properties. In contrast, the DFS1 potential results in unsatisfactory structural properties for liquid water, but improved density and energetic properties. 3. Force DeViations. To gain further insight into the accuracy of these short-range potentials, the force deviations for the FS, DFS2, and FM potentials were analyzed for Rcut ) 1.0 and 1.2

x∑ x∑

||F Bcutoff -B FEwald ||2 I I

∆F )

N

+

nm. The first set of quantities analyzed were the ratios of the Euclidean norm of the difference between the two sets of forces relative to the Euclidian norm of the Ewald forces,

I∈conf.

(7) ||F BEwald ||2 I

I∈conf.

where the summation is over all atoms in a single atomistic configuration. The temporal distribution of this quantity was characterized by its time average 〈∆F〉, standard deviation 〈δ∆F〉, and maximum [∆F]tmax. The second quantity analyzed was the maximum relative force deviation acting on a single atom in a configuration,

[∆FI]conf max ) maxI∈conf.

||F Bcutoff -B FEwald || I I ||F BEwald || I

(8)

in addition to its time average 〈[∆FI]conf max 〉 and time maximum t ] . The final quantity analyzed was the total (i.e., [[∆FI]conf max max uncompensated) charge ∆Q within a sphere of radius Rcut centered on the atom with the maximum force deviation measured in eq 8 along with this quantity’s time average 〈∆Q〉 and standard deviation 〈δ∆Q〉. All statistics were calculated from 6 ns simulations of 512 SPC/E and TIP3P water molecules using the NVT ensemble. Configurations were sampled every 0.1 ps for a total of 6 × 107. A similar analysis was carried out for the four Cl- ions dissolved in 512 SPC water molecules from a 16 ns simulation. As further described below, the SPC water model was chosen as an additional test of transferability. In Figure 4 ∆F and [∆FI]conf max are plotted as functions of time for the SPC/E model using the FS, DFS2, and FM potentials with Rcut ) 1.0 nm in cyan and Rcut ) 1.2 nm in black. The white lines are the running total averages of each respective quantity. In Figure 5 the same quantities from Figure 4 are shown for the Cl- ions. The statistical analysis of the distributions ∆F, [∆FI]conf max , and ∆Q is compiled in Table 2.

4718 J. Phys. Chem. B, Vol. 112, No. 15, 2008

Figure 9. Same SPC/E water RDFs shown in Figure 6 for T ) 275 K (a) and T ) 340 K (b).

As previously appreciated,32 the FS potential gives better agreement with the Ewald method than the DFS potentials. For example, the DFS2 potential with Rcut ) 1.0 nm produced slightly larger deviations and noticeably broader distributions (i.e., 〈∆F〉 and 〈δ∆F〉), as shown in Figure 4 and summarized in Table 2. Increasing the cutoff to 1.2 nm causes a very moderate decrease in force deviations for the DFS2 potential (∼3% for 〈∆F〉), but a significant decrease for the FS and FM potentials (∼30% for 〈∆F〉). Finally, Figure 5 and Table 2 show that the same trends hold for the Cl- ions, although the magnitude of deviations is almost an order of magnitude larger. An interesting pattern emerges from the 〈∆Q〉 and 〈δ∆Q〉 values in Table 2. The FM and DFS2 results consistently coincide, while the FS results differ. This suggests that the atoms with the largest force deviations are the same for the DFS2 and FM potentials, but different for the FS potential. The average uncompensated charge 〈∆Q〉 is actually smaller for all of the short-range electrostatic potentials in water, but larger for the ions. The increased structure and polarization of the water molecules surrounding the ions likely leads to the larger uncompensated charge. B. Accuracy in Various Liquid-State Environments. The FM effective potential was first tested on SPC/E bulk water (i.e., the same model that was used in its determination from eq 2). It was then used to simulate TIP3P water to test the transferability of the effective potential on a different water model. Finally, it was used to simulate SPC/E and TIP3P water at multiple temperatures to test temperature transferability. If not otherwise stated, each system was composed of 512 H2O molecules and was simulated for 2 ns using periodic boundary conditions and a constant NPT ensemble implemented with the Melchionna modifications of the Hoover algorithm,46 which employs a Nose´-Hoover thermostat and a barostat. 1. SPC/E Water. The structural properties of the SPC/E water simulations (T ) 300 K) are described by the oxygen-oxygen radial distribution functions (RDFs), gOO, in Figures 6 and 7. As shown in Figure 6, the FM (Rcut ) 1.0 nm) and Ewald results are in almost perfect agreement. The only difference is a small dip in the FM gOO that spans 0.01 nm just beyond the cutoff. In contrast, the FS (Rcut ) 1.0 nm) gOO shows much larger deviations at all distances. The largest FS deviation is in the vicinity of the cutoff where the oxygen density is significantly structured, resulting in a hump preceding and a dip following

Izvekov et al.

Figure 10. Same TIP3P water RDFs shown in Figure 8 for T ) 275 K (a) and T ) 340 K (b).

Figure 11. Water density as a function of temperature for the SPC/E (a) and TIP3P (b) models obtained using the Ewald method (cyan line) and the FS (dashed line) and FM effective (black solid line) potentials with Rcut ) 1.0 nm.

Rcut. Figure 7a shows that, when a smaller cutoff is used for the FS potential (Rcut ) 0.9 nm), the structural deviations increase significantly. The equivalent FM potential, however, has little effect on the structural deviations (i.e., the post-cutoff dip is only slightly deeper). Surprisingly, the moderately damped DFS1 potential (Rcut ) 1.0 nm) that was fit to the FM potential shows only slight improvement over the FS potential (results not shown), while the strongly damped DFS2 potential produces an RDF that is almost identical to that of the FM potential (cf. Figure 7b). This was unexpected, given the apparent visual similarity between the DFS1 and the FM potentials in Figure 3. Some of the thermodynamic and dynamical properties for SPC/E water at T ) 300 K and P ) 1 bar are summarized in Table 3. The FM effective potential for Rcut ) 1.0 nm outperforms every other tested potential for every measured property, and even that for Rcut ) 0.9 nm generally outperforms the FS potential for Rcut ) 1.0 nm. The density discrepancy between the Ewald and FM potentials with Rcut ) 1.0 nm is only 0.08%, while those for the FS and the DFS2 potentials are ∼5-fold larger (∼0.4%,). Similarly, the configurational energy

Coarse-Graining in Interaction Space

J. Phys. Chem. B, Vol. 112, No. 15, 2008 4719

TABLE 4: Liquid-Phase Properties of the SPC/E and TIP3P Models from Simulations Run with Various Electrostatic Treatments at T ) 275 K and T ) 340 K, and P ) 1 bara model: potentialb, Rcut (nm) SPC/E: Ewald FS, 1.0 FM, 1.0 TIP3P: Ewald FS, 1.0 FM, 1.0

Fc (kg/m3)

Dd (10-9 m2/s)

Ue (kJ/mol)

RPf (10-4 1/K)

κT g (10-5 1/bar)

1008.9 (10.6) 971.8 (11.9) 1006.1 (10.8) 966.7 (12.4) 1007.5 (10.6) 973.5 (11.7)

1.39 (0.02) 5.25 (0.01) 1.70 (0.02) 5.45 (0.01) 1.40 (0.02) 5.22 (0.01)

-48.21 (0.28) -44.24 (0.33) -46.56 (0.28) -42.64 (0.33) -48.22 (0.28) -44.29 (0.33)

4.7 7.8 4.2 8.9 4.6 8.0

4.0 5.7 4.3 6.4 4.1 5.5

1034.1 (11.6) 977.4 (13.3) 1031.3 (11.4) 969.7 (13.3) 1036.0 (11.6) 978.5 (13.3)

3.72 (0.02) 8.55 (0.01) 3.71 (0.02) 9.14 (0.01) 3.72(0.02) 8.88 (0.01)

-42.47 (0.26) -39.10 (0.31) -41.08 (0.26) -37.67 (0.31) -42.50 (0.26) -39.15 (0.26)

6.5 10.3 7.8 10.1 7.0 10.1

4.5 6.7 4.4 7.2 4.5 6.7

a Standard deviations in parenthesis. b Two numbers refer to T ) 275 K and T ) 340 K. c Density. d Self-diffusion coefficient. e Average configuration energy. f Thermal expansion coefficient. g Isothermal compressibility.

for the FM effective potential deviates less from the Ewald results than those of the DFS2 and FS potentials. The agreement in the self-diffusion constants from the Ewald and FM simulations is excellent and, again, better than all of the other potentials. Similarly, in a constant NVT simulation, where density discrepancies are no longer an issue, the FM effective potential (Rcut ) 1.0 nm) reproduces the Ewald self-diffusion rates within 0.01 × 10-9 m2/s, while FS and DFS2 still differ by 0.1-0.2 × 10-9 m2/s. In contrast to the structural results, the moderately damped DFS1 potential resulted in better values for the density, diffusion, and average internal energy than the strongly damped DFS2 potential. The dielectric properties of systems simulated with the various potentials were also explored by measuring the static dielectric constant, a property that converges very slowly,47 from 10 ns simulations run in the constant NVT ensemble. As described in the Appendix, a correction to the dielectric constant was calculated for the Ewald simulation, but no corrections were appropriate for the FS, DFS, or FM effective potentials. For the purposes of calculating the dielectric constant, atomic partial charges were treated in the normal way, despite the nonCoulombic (i.e., effective) treatment of the electrostatic interactions. The Ewald (including the proper correction), FS, and FM potentials resulted in r ) 72.6, 71.1, and 71.5, respectively. 2. Model Transferability (TIP3P Water). The ability of the FM potential to reproduce the structure of bulk water transfers well to the TIP3P water potential, which has different partial charges qR/β on the hydrogen and oxygen sites. Figure 8 shows the oxygen-oxygen RDFs, gOO, for the Ewald, FM, and FS potentials (the latter two use Rcut ) 1.0 nm). The FS potential produced a better structure at small separations than it did for the SPC/E model, but still resulted in substantial deviations from the second minimum to beyond the cutoff value. The thermodynamic and self-diffusion properties for the TIP3P model are also well reproduced by the FM effective potentials. The density discrepancies for the short-range potentials were almost identical to those for SPC/E water. The FM effective potentials again show excellent agreement in self-diffusion rates; the FS and DFS2 potentials do not. Similar to the SPC/E model, constant NVT simulations for TIP3P water (results not shown) resulted in Ewald and FM effective potential self-diffusion coefficients that are within 0.01 × 10-9 m2/s, while FS and DFS2 potentials deviate by ∼0.2 × 10-9 and 0.3 × 10-9 m2/s, respectively. Finally, the trend in dielectric constants observed in the SPC/E

simulations was also observed for the TIP3P model (r ) 91.5, 91.0, and 91.2 for the Ewald, FS, and FM potentials, respectively). 3. Temperature Transferability. The FM potential also demonstrated excellent temperature transferability for the SPC/E and TIP3P water models, as indicated by the RDFs for T ) 275 K and T ) 340 K (cf. Figures 9 and 10) and thermodynamic properties (Table 4). This transferability is further validated by the density maximum between T ) 150 and T ) 250 K (Figure 11). Screening effects in condensed phase systems are generally thought to strongly depend on the thermodynamic conditions, so the FM effective potential’s temperature transferability is a remarkable feature. Moreover, this transferability is quite encouraging in that the “missing entropy” effects, which result in poor temperature transferability in resolution coarse-grained models, are less of an issue when all atom resolution is maintained. Thus, the present “coarse-graining in interaction space” approach may offer distinct advantages over resolution coarse-graining in the preservation of atomistic resolution and the associated entropy within the length scale of the effective potential cutoff. 4. Higher Order Correlations. The pair correlations that define an RDF are not sufficient for describing the detailed structure of a condensed-phase system. In fact, models that are developed to perfectly reproduce an RDF often fail to reproduce higher order correlations (e.g., three-body effects).48 A poignant example of this was recently presented for the specific case of short-range electrostatic interactions.49 This work demonstrated that a short-range electrostatics potential developed using the inverse Monte Carlo (IMC) method to reproduce the RDF of TIP3P water resulted in skewed orientational structure, as indicated by the distributions of the angle between the centers of mass of three neighboring water molecules and the orientational order parameter

q)1-

3

3

4

∑∑ 8 j)1 k)j+1

(

cos θOjOiOk +

)

1 3

2

(9)

where θOjOiOk is the angle connecting the oxygen of the ith reference molecule with the oxygens of its four closest neighbors (Oj‚‚‚Oi‚‚‚Ok).50 The orientational order parameter, q, is a good measure of the local tetrahedral structure; a value of 1 corresponds to a perfect tetrahedral network. In Figure 12 the angle distributions of θOjOiOk, θHjOiHk, and θHjHiHk (measured

4720 J. Phys. Chem. B, Vol. 112, No. 15, 2008

Izvekov et al.

Figure 12. Distributions of angles θOjOiOk, θHjOiHk, θHjHiHk, and θCM3 (a,b) and of the order parameter q in eq 9 (c,d) for Ewald, DFS2, and FM effective potentials (cyan line) and FS potential (black line) from simulations using SPC/E (a,c) and TIP3P (b,d) water models at T ) 300 K and P ) 1 bar. The short-ranged potentials were those with Rcut ) 1.0 nm. Deviations in the FS potential are highlighted by dotted circles in panel a and arrows in panel c. Figure 14. Comparison of Na+-O (a) and Cl--O (b) RDFs obtained from simulations run with the Ewald method (cyan solid line) and the FS (dashed line), FM effective, and DFS2 (black solid line) potentials with Rcut ) 1.0 nm. Insertions highlight the deviations around the cutoff.

Figure 13. Distributions of the hydrogen bond angle θHB from the same simulations represented in Figure 12.

between the four closest neighbors of the ith molecule) are presented for both the SPC/E and TIP3P simulations using the Ewald, FS, FM effective, and DFS2 potentials. The fourth curve in Figure 12 (labeled as “CM”) shows the distribution of the angle formed between centers of mass of the reference molecule and its two closest neighbors. The latter provides a direct comparison with ref 49. Figure 13 shows the distribution of the hydrogen bond angles θHB formed between the oxygens of any two hydrogen-bonded molecules and the donor O-H bond. The hydrogen bond was limited to neighbors with O‚‚‚O and O‚‚‚H separations that were less than 0.35 and 0.25 nm, respectively. All distributions were obtained from 10 ns simulations in a constant NVT ensemble. The FS, DFS2, and FM effective potentials all do an excellent job of reproducing the full Ewald angular distributions. For both SPC/E and TIP3P, the distribution curves from Ewald, DFS2, and FM potentials were indistinguishable on the scale of graphs. However, the FS potential produced a slightly skewed structure (as highlighted in Figure 12 by dotted circles and arrows) for the SPC/E model. No difference was observed for the TIP3P model. The FS potential also caused a decrease in the maximum hydrogen bond angle (θHB ) 0°) indicating softer hydrogen bonds for both water models, though the effect was more pronounced for SPC/E. The increased sensitivity of the SPC/E model compared to the TIP3P model is a likely due to the fact that SPC/E has larger partial charges, which would lead to larger

deviations in the hydrogen-bonding structure due to a shifted Coulomb potential at the cutoff distance. The average values presented in Table 5 show that the FM effective potential is the most accurate for the SPC/E model. The FM effective and DFS2 results are virtually identical for the TIP3P model (with minor differences in the θHB and CM distributions). The CM angle and order parameter distributions for the TIP3P model (Figure 12) show far better agreement with the Ewald results than those previously presented for the short-range IMC potential.49 The ability of the FM effective potential to capture the orientational structure emphasizes the power of the FM methodology and the importance of a statistical mechanical framework that builds higher order correlations (three-body in this case) into the resulting potential energy function.38 5. SolVated Ions. In order to assess transferability of the FM effective potential to less homogeneous systems, it was compared to the Ewald, FS, and DFS2 potentials in simulations of solvated ions. Two systems, one with a single sodium ion solvated in 512 TIP3P water molecules and the other with four chloride ions solvated in 512 rigid SPC water molecules, were simulated for 20 ns using periodic boundary conditions and the constant NPT ensemble defined by a modified Hoover algorithm46 and ambient pressure and temperature. The chloride ion, modeled with the OPLS-AA parameters, was solvated in SPC water as an additional test of model transferability. The oxygen solvation structures (RDFs) around each ion are depicted in Figure 14. The structure on the first and second solvation shells is almost identical for all three potentials, but the third solvation shell and the vicinity of the cutoff are more accurately captured by the FM effective potential. Similar to bulk water, the gNa+O and gCl-O from the Ewald and FM effective potentials differ only in a small region just after the cutoff. Also similar to bulk water, the DFS2 potential produced structural results that are comparable to those of the FM effective potential. The ion-oxygen RDFs along with some thermodynamic and dynamical properties are summarized in Table 6. The configurational energies and the densities are best for the FM effective potential, and only

Coarse-Graining in Interaction Space

J. Phys. Chem. B, Vol. 112, No. 15, 2008 4721

TABLE 5: Average Value (and Standard Deviation) of the Angle (in deg), Hydrogen Bond Angle (in deg), and Order Parameter Distributions model: angle

Ewald

FS

DFS2

FM

SPC/E: OOO HOH HHH CM HB q

104.31 (25.85) 103.26 (25.35) 99.72 (31.83) 108.76 (22.36) 18.51 (10.50) 0.632 (0.190)

104.14 (26.13) 103.11 (28.56) 99.59 (31.93) 108.67 (22.61) 18.79 (10.64) 0.624 (0.191)

104.28 (25.91) 103.23 (28.39) 99.70 (31.86) 108.74 (22.39) 18.54 (10.51) 0.612 (0.190)

104.31 (25.85) 103.26 (25.35) 99.72 (31.83) 108.76 (22.36) 18.51 (10.50) 0.632 (0.190)

TIP3P: OOO HOH HHH CM HB q

102.86 (28.35) 102.01 (30.13) 99.04 (32.27) 107.98 (25.60) 22.05 (11.96) 0.555 (0.191)

102.77 (28.44) 101.93 (30.21) 98.96 (32.30) 107.93 (25.70) 22.22 (12.01) 0.551 (0.191)

102.86 (28.36) 102.01 (30.14) 99.04 (32.27) 108.00 (25.60) 22.08 (11.97) 0.555 (0.191)

102.86 (28.35) 102.01 (30.13) 99.04 (32.27) 107.98 (25.60) 22.05 (11.96) 0.555 (0.191)

TABLE 6: Properties of Na+ in TIP3P Water (First Row) and Four Cl- in SPC Water (Second Row) from the Ewald Method and the FS, DFS2, and FM Effective Potentials with Rcut ) 1.0 nma property max (gion-O)b min(gion-O)c 〈U〉/N (kJ/mol)d Fw (kg/m3)e 〈V〉/N (10-2 nm-3)f 〈δV〉/N (10-4 nm-3)g D (10-9 m2/s)h

Ewald

FS

DFS2

FM

7.31 3.69 0.145 0.516 -41.82 (0.28) -44.24 (0.29) 1013.5 (12.1) 968.2 (11.2) 2.946 3.066 3.51 3.56 1.07 (0.02) 1.17 (0.01)

7.24 3.65 0.152 0.533 -40.15 (0.28) -41.66 (0.29) 1008.6 (12.3) 961.0 (11.2) 2.960 3.089 3.60 3.66 1.04 (0.02) 1.36 (0.01)

7.30 3.67 0.146 0.516 -40.40 (0.28) -41.90 (0.29) 1011.2 (12.1) 964.3 (11.2) 2.953 3.078 3.54 3.58 4.28 (0.02) 3.00 (0.01)

7.31 3.69 0.145 0.516 -41.65 (0.28) -43.33 (0.29) 1014.6 (12.1) 966.0 (11.2) 2.943 3.073 3.52 3.58 1.08 (0.02) 1.12 (0.01)

a

e

Standard deviations in parenthesis. b First maximum in ion-O RDF. c First minimum in ion-O RDF. d Average configuration energy per particle. Water density. f Average volume per particle. g Standard deviation in the volume per particle. h Ion diffusion coefficient.

slightly better for the DFS2 than the FS potential. The most striking difference is that the ion self-diffusion rate is significantly overestimated (3-4-fold) by the DFS2 potential, but well described by the FM potential. The FM potential’s ability to describe diffusion is rather surprising, given the fact that the ion-H2O distances can fall into the core region of the FM potential, which was approximated by bare Coulomb potential (eq 5a). Using the optimal FFM(r) for this region, which could be different from the Coulomb potential, should only improve ion dynamics and diffusion. Next, the association of two like ions (Cl--Cl-) in rigid SPC water was simulated. Potentials of mean force (PMFs), denoted by ∆G(r) in Figure 15, were obtained with force constrained ensemble simulations51 using a step size of 0.01 nm for the ionion separation. Each point was simulated for 100 ps after 30 ps of equilibration. As shown in Figure 15, the PMFs obtained with the FS, DFS2, and FM potentials are all very similar. Close inspection reveals that the FM effective potential result agrees better with the Ewald PMF at shorter separations (e.g., first decay and first (contact pair) minimum), but becomes closer to the FS potential near the desolvation (first) barrier. However, the differences between these PMFs are approximately as large as the error bars in each calculation (shown for the Ewald potential in Figure 15). 6. SolVated Hydrophobes. Finally, the FM potential for Rcut ) 1.0 nm was tested on an aqueous solution of methane molecules. The interactions in the system were modeled with the OPLS-AA force field52 and the flexible TIP3P water potential. The latter was chosen in order to assess whether the

FM effective potential transfers to a flexible water model. The simulations employed periodic boundary conditions and the constant NPT ensemble using the modified Hoover algorithm46 at ambient conditions. First, the solvation structure of four methane molecules dissolved in 512 water molecules was studied. The most visible structural differences, shown in Figure 16, were observed in the gOHMe and gHHMe distribution functions. The RDFs from methane atoms to water oxygens appeared to be insensitive to the way the electrostatics were calculated. As opposed to polar and charged systems, the water-methane RDFs do not exhibit perturbed structure in the region of cutoff for either the FS or the FM potentials. Rather, the first solvation shell peaks were most affected and, as shown in Figure 16, were slightly underestimated by the FS potential. 7. Nonaqueous Systems: Molten Salts. To assess the importance of the aqueous environment on the nature of the effective electrostatic interactions, the FM effective potential was explored in a molten salt composed of liquid NaCl ions simulated at 1200 K. There is no reason to assume a priori that the present effective electrostatic potential, originally force-matched to a particular model of pure liquid water at a particular thermodynamic state point, will accurately represent electrostatic interactions in disordered nonaqueous systems, especially one as different as a molten salt. This test is therefore an important one. Furthermore, the entire FM procedure can also be repeated for the molten salt system and the resulting effective shortranged electrostatic interaction specific to that system compared to the one transferred from the water system (i.e., both directly compared in its numerical form but also in its predictions of

4722 J. Phys. Chem. B, Vol. 112, No. 15, 2008

Izvekov et al.

Figure 15. Potential of mean force between two Cl- ions in SPC water simulations using the Ewald method (thick solid line) and the FS (dashed line) and FM effective (thin solid line) potentials with Rcut ) 1.0 nm. The errors (shown for the Ewald curve) were similar for all three potentials. Figure 17. Comparison of FM effective potentials from liquid water, U wFM(r), (solid line) with that fit directly to molten NaCl salt, U slt FM(r), (dashed line) obtained with Rcut ) 1.0 nm (a) and their absolute w slt w U wFM(r) - U slt FM(r) and relative (U FM(r) - U FM(r))/U FM(r) differences (b). Units of kJ/mol are used to show the slight deviations and can be converted to the units used in Figures 1 and 2 with the conversion factor 1 nm-1 ) 138.935 kJ/mol.

Figure 16. RDFs for methane in flexible TIP3P water at T ) 300 K and P ) 1 bar obtained with the Ewald method (cyan line) and the FS (dashed line) and FM effective (black solid line) potentials with Rcut ) 1.0 nm. Dotted circles highlight the differences.

the NaCl liquid properties relative to the predictions from the transferred interaction and to the full Ewald case). The molten NaCl salt was represented by 256 ionic pairs placed in a periodic cube. The ion-ion interactions were described by the Born-Mayer-Huggins potential with TosiFumi parameters.53,54 The system properties were collected from a 2 ns simulation run in the constant NPT ensemble at 1200 K and ambient (1 bar) pressure. The same FM procedure that was presented for liquid water was again used to directly force-match the molten salt system. Reference configurations were sampled from a 1 ns simulation run in the constant NVT ensemble at equilibrium density. Structures were sampled every 0.5 ps, resulting in 2000 reference configurations. The molten salt FM slt (r), was only obtained for effective electrostatic force, F FM Rcut ) 1.0 nm. slt w (r) and the original water U FM (r) The resulting salt U FM effective potentials are compared in Figure 17a. The two potentials are hardly distinguishable in the scale of the plot. The absolute and relative differences between the two are shown in Figure 17b. Although the absolute difference gets as large as 8 kJ/mol, the relative differences are