Boundary Conditions in Simulations of Aqueous Ionic Solutions: A

J. Phys. Chem. 1995, 99, 1322-1331

1322

Boundary Conditions in Simulations of Aqueous Ionic Solutions: A Systematic Study James E. Roberts and Jurgen Schnitker* Department of Chemistry, University of Michigan, Ann Arbor, Michigan 48109-1055 Received: June 7, 1994; In Final Form: October 4, 1994@

The room temperature properties of aqueous solutions of SPC, TIP4P, and MCY water are examined using either a potential cutoff or Ewald summation. The effects of interaction truncation are compared to the importance of the medium (vacuum or metal) that the entire system is immersed in; for this purpose simulations are also carried out in which the unit cell surface polarization term of De Leeuw, Perram, and Smith (Proc. R. Soc. London 1980, A373, 27) is explicitly accounted for. Both pure water and dilute solutions of C1- and Fe2+ are investigated. It is found that not only solvent-solvent orientational correlations but also ionsolvent energies, ion-induced solvent reorganization energies, and ionic diffusion coefficients can be exceedingly sensitive to the nature of the surroundings. It is asserted that reaction field methods and Ewald summation are the preferred types of boundary conditions, but that the method of parametrization of the potential functions has to be consistent with the eventual simulation conditions.

I. Introduction The evaluation of the slowly decaying Coulombic interactions is a longstanding problem in the computer simulation of condensed phases, both classically and quantum mechanically. Although the literature on the subject is extensive,’ it is only recently that a few papers have explicitly addressed the problem of charged (or uncharged) solutes in molecular solvent^.^-'^ Nevertheless, some dramatic differences between different boundary conditions have already been observed. Thus in a recent study we found a strong dependence of ion-solvent energies on the boundary conditions,2 potentially affecting any of the popular free energy calculations where an uncharged solute is mutated into a charged solute. Other properties that were found to exhibit crucial boundary condition dependencies include the potential of mean force between ferrous and femc ions in water3 and, most importantly, the tertiary structure of protein^.^^^ In this paper, we wish to follow up on some of the scattered observations and systematically investigate the differences between cutoff schemes and Ewald summation in infinitely dilute aqueous ionic solutions, with respect to sensitive measures for energetics, structure, and dynamics. It is important to realize that the two schemes differ fundamentally in more than one way. First, there is the obvious fact that the range of the potential is finite and fairly shortranged in one case and formally infinite in the other. Second, there is another difference between the two approaches that seems to have been widely appreciated only since the publication of a landmark paper by De Leeuw, Perram, and SmithI4 in 1980 (even though the issue is also taken up in the older literat~re).’~-’~ The point is that Ewald summation, if done in the standard way by adding up a reciprocal space sum and a real space sum, is not carried out in vacuum; rather the macroscopic surroundings of the growing sample of boxes are effectively those of a medium of infinite dielectric constant (“conducting” or “tinfoil” boundary conditions). There is no such peculiarity with respect to the surroundings (“vacuum” conditions) in standard cutoff simulations or if so-called direct lattice sums are evaluated. While the effects of different cutoff radii and some other aspects of the boundary conditions in aqueous s ~ l u t i o n ~ and -’~ in pure hydrogen-bonding l i q ~ i d s ‘ ~have - ~ ~occasionally been @

Abstract published in Advance ACS Abstracts, January 15, 1995.

studied, no attempt has been made to explicitly disentangle the effects of the interaction range (finite vs infinite) from those related to the surroundings (vacuum vs metal). The purpose of this paper is thus 2-fold: We would like to (i) obtain per se a more complete picture of the differences between the properties of dilute ionic solutions under cutoff and Ewald conditions and (ii) gain a better physical insight into the nature of these differences by explicitly separating the effects associated with the interaction range from those connected to the nature of the surroundings. Mathematically, the effect of the surrounding medium can be cast in a precise form by comparing direct (=vacuum) lattice sums with Ewald sums, as worked out in the above-mentioned paper by De Leeuw et al.14 It is found that the electrostatic energy for the two cases differs by a “surface term” that is proportional to the square of the unit cell dipole moment. If Ewald summation is used it is-probably for the sake of simplicity-practically always implemented in the traditional way, i.e. in the absence of the surface term. Thus surprisingly little is known concerning the relevance of this term, including its relevance for the most widely studied kind of fluid system, i.e. liquid water and aqueous solutions. In the present study, we consider simulations of aqueous systems under several kinds of boundary conditions where the surface term in the Hamiltonian is either explicitly present or not present. In particular, we examine (1) cubic (or “minimum image”) cutoff conditions, (2)cubic cutoff conditions with a “depolarization correction for surface effects”, (3) spherical cutoff conditions, (4)Ewald summation, and (5) Ewald summation with a “polarization correction for surface effects”. While methods 1, 3, and 5 all correspond to surrounding vacuum conditions, methods 2 and 4 both mimic tinfoil boundary conditions. Comparison of the data along various lines can now provide insight into the nature of the differences between these types of boundary conditions. For example, if the surrounding medium effects were much more important than the effects of periodic replication, method 1 for a short-ranged system and method 5 for an infinitely replicated system should yield exactly the same results, since both correspond to a strictly cubic geometry and surrounding vacuum conditions. Analogously, the tinfoil methods 2 and 4 would yield the same results if the unit cell replication was essentially irrelevant. Another infor-

0022-3654195120991322$09.00/0 0 1995 American Chemical Society

J. Phys. Chem., Vol. 99, No. 4, 1995 1323

Simulations of Aqueous Ionic Solutions mative comparison is the one between simple cubic truncation (1) and simple spherical truncation (3) in order to identify artifacts associated with the anisotropy of the simulation cell. 19,20.25,26 It also turns out that our cutoff/tinfoil boundary conditions (2) are the cubic equivalent of the conventional reaction field m e t h ~ d . ~The ~ - latter ~ ~ can in fact be described as “combining a spherical cutoff with a polarization correction for surface effects”, in obvious deviation from the standard physical reasoning behind the method. The rather drastic effects of reaction fields on the dielectric properties of liquid water have been well documented, but surprisingly the technique has only been implemented in simulations of the pure l i q ~ i d ~ Oand - ~ in ~ simulations with an uncharged solute!b We will see that our renewed emphasis on an alternative view of the reaction field method-as a correction technique for spurious surface effects in a spherical truncation geometry-implies that the method may well be extended to systems with charged species. Naturally one expects some model dependence of the simulation results. It can be shown that the crucial surface term affects the electrostatic potential via the quadrupole moment of the unit cell,15and hence molecular models with larger partial charges and larger multipole moments might conceivably be more sensitive to the boundary conditions. We investigate three different (rigid) water models that are in widespread use, namely, the SPC, TIP4P, and MCY models. After a preliminary examination of the neat liquids, we study solutions of C1- and Fez+, thus considering two of the most common representatives for mono- and divalent ions in water. In section 11, we first summarize the equations relevant for the electrostatics of a unit cell under periodic boundary conditions, with special attention to some issues that specifically arise in simulations of disordered condensed phase systems. In section 111, we describe the molecular dynamics simulation setup. The results of pure water simulations are presented and discussed in section IV, before proceeding to simulations of aqueous ionic solutions in section V. The conclusions are given in section VI.

11. The Electrostatics of a Replicated System In the following, we first give the fundamental result of De Leeuw et a1.14 for the true electrostatic energy of an infinitely replicated unit cell of point charges. Because of the conditional convergence of Coulombic lattice sums, we anticipate that this energy does not have a unique value. Subsequently, we indicate how to resolve some complications that arise in simulations because of use of the minimum image convention and because of truncation of the interactions; it will be pointed out that these complications can effectively be ignored. De Leeuw et al. calculated the energy per unit cell Evacuum of a perfect crystal of N electroneutral unit cells (with N c-) that has spherical geometry and is surrounded by vacuum. Thus the crystal is infinite, but still has a well-defined center, and the energy is obtained by summing over all pair interactions inside the crystal (without invoking any sort of minimum image boundary condition). According to De Leeuw et al., the “direct lattice sum energy” Evacuum can be written as the sum of two contributions: l4

-

Equivalently, we can write

which obviously implies

The energy EEwald is a unique reference energy, reflecting the fact that it corresponds to carrying out the lattice summation with a suitable convergence factor. The lattice sum energy under vacuum conditions, Evacuum, and correspondingly the socalled “surface term”, Le. either Epolarization or ,!&epolarization, vary with the geometrical situation, in particular the specific choice of unit cell. As shown by De Leeuw et al.14-and as also implicit in Felderhof s derivation of fluctuation theorems for dielectk@--the surface term is given by

(3) Here ri and qi are the positions and charges, respectively, of all the particles in the unit cell of volume v (the usual factor of 1/(4n€0) has been omitted). We see that the surface term is proportional to the square of the unit cell dipole moment. The original derivation of De Leeuw et al. is for a cubic Bravais lattice of unit cells, but more complicated expressions for the surface term-for arbitrary, noncubic lattices and arbitrary, nonspherical shapes of the growing crystal of unit cells-were subsequently given by Smith3’ and Olives.38 We emphasize that it is the dependence of the surface term, Epolarization or Edepolarization, on the precise specification of the system that reflects the famous “dependence on the order of summation” of the value of a lattice sum. The designation as a “surface term” relates to the fact that any charges in the interior of the unit cell will cancel out, and thus the unit cell dipole moment has to be attributed to an imbalance or asymmetry of the unit cell surface charge distribution. Since EEwdd is the energy of a crystal under screening of all surface effects, it physically corresponds to the energy of a crystal surrounded by a hypothetical perfectly conducting medium that neutralizes all surface charges (tinfoil boundary conditions). Mathematically, however, this is not-as might initially be expected on the basis of eq lb-due to the presence of a term that has somehow been “added on” after evaluating the two sums in real space and reciprocal space in the normal manner. The tinfoil boundary conditions are inherent to the standard procedure of calculating Ewald sums, and probably it is precisely because of this circumstance that there has been such widespread disregard of the physical issue. Once the association of direct lattice sum energies with vacuum conditions, and of Ewald energies with tinfoil conditions, has been made, eqs l a and l b can be used to go from one set of conditions to the other. Thus vacuum conditions can be recovered by adding a “polarization correction”, eq 3, to the Ewald energy, and tinfoil conditions can be obtained not only by Ewald summation but also by correcting a direct lattice sum energy with a “depolarization term”, thus removing the effects associated with the imbalance or asymmetry of the surface charge distribution. From eq 3 it is also clear that the total electrostatic energy under tinfoil conditions is the same as under surrounding vacuum conditions if the unit cell dipole moment vanishes. An apparent problem with the mathematically sophisticated treatment of De Leeuw et a l . I 4 is that their specific direct lattice sum Hamiltonian, while theoretically suitable for a Monte Carlo simulation, is not really relevant for molecular dynamics simulations. In the case of the latter, the dynamic implementation of periodic boundary conditions requires that the mimimum image condition is applied to each and every pair interaction.

1324 J. Phys. Chem., Vol. 99, No. 4, 1995 (To see this, consider a particle j that is leaving on one side of the unit cell and its image j‘ that is entering on the opposite side. Obviously, at the moment of the surface crossing both images have to be at the same distance from the reference particle i if a discontinuous jump of the energy is to be avoided.) Thus, we do not have a distinct “center” of the replicated system as implied by the formalism of De Leeuw et al. Furthermore, interactions are normally truncated, either after the nearest neighbor or within a sphere of a given radius, succinctly described as “toroidal” instead of periodic boundary conditions. De Leeuw et al. indeed suggest evaluating the surface-related part of the electrostatic potential from the mimimum image positions r, of all charges in the unit cell, relative to the reference point r:

Roberts and Schnitker posed to site-by-site-basi~,’~>~~-~~ since the correct quadrupole moments in Eq 4 are those of slightly distorted (or “corrugated’) unit cells. Third, the apparent violation of the electroneutrality condition has to be addressed that is posed by the absence of counterions in an ionic solution under infinite dilution conditions. The entailing problems can be circumvented by assuming the presence of a compensating background charge and rigorous enforcement of the minimum image condition.* Finally, in a rigid solvent model the intramolecular contribution to the surface potential (4) can effectively be neglected, thus simplifying practical calculations.2 While simulations using Ewald summation are always carried out under tinfoil boundary conditions, whether knowingly or unknowingly, this is not necessarily true for simulations under cutoff conditions. In the reaction field method, each molecular dipole pi is at the center of a cutoff sphere of radius r, which in turn is surrounded by a dielectric continuum of permittivity E W . This gives rise to an additional term in the Hamiltonian 0f30a

Using the electroneutrality condition, it can easily be shown that summation over all charges in the system, each of which sees the potential 4, in fact leads to the energy expression 3.2,14 Evidently, the surface potential depends on the trace of the quadrupole moment tensor of the complete charge distribution of the (minimum image) unit cell. However, the potential 4 is not the only one that is compatible with the energy expression 3, as we have shown in our previous paper.2 The potential that literally corresponds to De Leeuw et al.’s direct lattice sum Hamiltonian is indeed slightly more complicated, with both a dipolar and a quadrupolar contribution.2 This follows from the often ignored work of Redlack and Grindlay on the electrostatic potential in finite systems of replicated and electroneutral unit ~ e 1 l s . lIt~ can also be shown,2 however, that the expression 4 is indeed rigorously true for any finite system of replicated unit cells if minimum image conditions are invoked and if the system has cubic or spherical symmetry. Of course these are exactly the conditions we are interested in, and in summary we find that the surface potential (4)-even though inconsistent with the direct lattice sum Hamiltonian that is the point of departure of De Leeuw et al.’s derivation-and also the surface energy expression (3) are both applicable to simulation cells under minimum image conditions with truncation of the interactions. Analogous to what was said before about the total system energies, Eq 4 allows the recovery of the electrostatic potential under surrounding vacuum (or tinfoil conditions) from a simulation under the opposite kind of conditions by polarization (or depolarization) correction.’* Most interestingly, we find that potential-dependent single-particle properties like solvation energies2 or positions of energy bandsI6 vary with the unit cell quadrupole moment and thus can show a dependence on the boundary conditions even i f the unit cell dipole moment vanishes (Le. if the total system energy is independent of the boundary conditions). Standard source codes for doing Ewald sums do not allow for inclusion of the correction terms 3 and 4 and thus for a choice between surrounding vacuum and tinfoil boundary conditions. In fact, a number of minor subtleties have to be resolved if eqs 3 and 4 are to be implemented in a molecular dynamics simulation of an ionic solution at infinite dilution. These are addressed in detail in our previous paper,2 and thus we only give a brief summary. First, in reflection of what was outlined above, it is important to take the distances (r, - r) in the surface potential (4)as minimum image distances. Second, it is important to do the minimum image shifting on a molecule-by-molecule-as op-

1ERF-1 1 E , = - ; T r c 3 -&:Mi

(5)

E,+-

2 with Mi being the dipole moment of the spherical cavity. If we replace rc3= (3/47c)Vsphere by (3/4x)Vcube, and if we use the fact that the cavity dipole moment for strict minimum image truncation is a constant,

it can be seen that Eq 5 is an equivalent of the previous surface energy #!?depolarization(i.e. the negative of Eq 3), namely, for a cubic truncation geometry and for a dielectric constant of the surrounding medium of EW = 00. Conversely, the surface term of De Leeuw et al. can be said to provide a “cubic reaction field”.30b This is also evident from the extra term in the pair potential in standard reaction field simulations; that is, there is a quadratic dependence on the pair distance,30ajust as is present in Eq 4 for the surface potential. We now recall that the formalism of Redlack and Grindlay15 and of De Leeuw et al.14 is completely general for any electroneutral system. In particular, calculation of the surface terms 3 and 4 for charged species (such as ions in solution) does not pose any fundamental difficulties, and even the physical interpretation (elimination of spurious surface effects) remains the same. Hence, we conclude that at least formally the reaction field method can well be applied to systems with charged species, contrary to what is often assumed.

111. Molecular Dynamics Simulations Room temperature simulations of bulk water and of dilute aqueous solutions of C1- and Fe2+ were carried out with three different water models and four alternative kinds of boundary conditions. We note that the emphasis is not on realistic modeling but on an illustration of the model dependencies of the differences between systems under different boundary conditions. Only at the very end of this paper are we concerned with the suitability of the various simulation setups for actual modeling purposes. The simulated system consisted of 200 rigid water molecules (or 200 molecules plus 1 ion) in a simple cube, using the threesite SPC40 model or either one of the four-site TIP4P‘“ and MCY42 models. To keep the comparison between the models

J. Phys. Chem., Vol. 99, No. 4, 1995 1325

Simulations of Aqueous Ionic Solutions simple, the density was always set to the experimental value of 0.997 g cmW3,thus ignoring the fact that a different value (in particular in the case of MCY water) may actually be preferred if an optimum overall reproduction of experimental properties is called for. For the simulations of C1- in water, the short-ranged ionwater interaction was taken as Lennard-Jones-like, using the parameters of Chandrasekhar et al. for the interaction with TIP4P water.43 For simplicity, we ad hoc used the same potential type and the same parameters also in conjunction with the partial charges of the SPC and MCY models. The shortranged interaction of Fe2+ with SPC water was described by the potential function of Curtiss et aL4 neglecting the fact that the parameters were actually optimized in simulations with a flexible version of the SPC model. For the interaction of Fe2+ with TIP4P and MCY water, we again combined the distributed charges of these water models with the original short-range interaction of Fe2+ with SPC water. In all simulations of ions in water, we use the same unit cell length as in the pure water case, implying a negligible partial molar volume of the solute. Five main ways for the evaluation of the Coulombic interactions were studied: a minimum image (cubic) cutoff, a minimum image cutoff plus depolarization correction, a spherical cutoff of half the box length, Ewald summation, and Ewald summation plus polarization correction. All water-water and water-ion interactions, including polarization and depolarization terms, were evaluated in the molecular minimum image convention as defined by the oxygen positions. In order to have energy conservation in the cubic cutoff simulations, the Coulombic interactions and the correction terms were smoothly tapered off over a range of 0.2 8, with the product of three separate spline that were functions of the pair distances IAxl, IAyl, and I A Z I . ~A~ corresponding spherical tapering function was employed in the simulations using a spherical cutoff. Ewald summation was carried out with a convergence parameter of a = 5 . 6 L that allows the truncation of the real space sum within the original minimum image cube. The reciprocal space part was evaluated using fast single sums,45with k,, = 5 and k2 < 27 for a total of 587 terms. Standard rigid body molecular dynamics simulations were carried out, with a time step of 3 fs and using the Verlet algorithm in conjunction with the SHAKE method.46 After preparing sufficiently equilibrated configurations, most runs were carried out as NVT runs of 10000 steps duration, resampling the velocities from a Boltzmann distribution (T = 25 "C) every 100 steps.46 If dynamic quantities were of interest, longer NVE runs of at least 50000 steps duration were performed. Diffusion coefficients were obtained from the mean square displacements, and reorientational correlation times were calculated with suitable analytic continuations of the relevant time correlation functions to infinity. In particular among the shorter NVT runs there is a small spread of the average temperatures, due to the finite amount of velocity sampling. Since solvent reorganization energies are exceedingly sensitive to very small temperature differences between the systems that are to be compared (for our system size z15 kJ mol-' for just 1 K of difference in temperature), an explicit temperature correction was invoked for the calculation of solvent reorganization energies. Thus the energies of the solvent plus ion runs and of the pure water reference runs were always corrected to exactly 25 OC by making use of approximate heat capacity values for SPC, TIP4P, and MCY water, as determined from the potential energy fluctuation in NVE runs of the neat l i q ~ i d s . ~ 'The , ~ ~ solvent reorganization

TABLE 1: Dipole Moments p (in D) and Components ax, and Q, of the Traceless Quadrupole Moment Tensor (in esu cm2,x = H-H Vector, z = Bisector, Relative to Center-of-Mass) of Water Models

fF26

P

SPC TIP4P MCY

QXX

2.27 2.18 2.19

2.12 2.20 3.14

QW

QZ

-1.82 -2.09 -2.78

-0.29 -0.11 -0.36

TABLE 2: Static Properties of Pure Water at T = 25 "C and e = 0.997 g c ~ n with - ~ the SPC, TIP4P, or MCY Model, Using Either a Cubic Cutoff or Ewald Summation SPC SPC TIP4P TIP4P MCY MCY

cubic cutoff Ewald cubic cutoff Ewald cubic cutoff Ewald

-44.31 -41.55 -43.55 -41.17 -37.20 -35.41

f 0.04 -870 f 180 f 0.02 +510 f 80 f 0.25 -1100 f 140 f 0.27 +190 f 210 f 0.18 +7160 f 70 f 0.15 $8290 f 140

0.13 f 0.03 5.06 f 2.49 0.16 f 0.02 1.91 f 0.97 0.15 f 0.02 1.45 f 0.72

energy was then calculated in the usual way from the difference in the total system energies with and without the presence of the solute. Error bars were obtained from subaverages over five consecutive segments of a given simulation. Because of our special temperature correction procedure, we did not calculate error bars for solvent reorganization and hydration energies. We can estimate, however, that the uncertainties of the latter are smaller than f 1 0 kJ mol-' if both runs are of the NVE type and smaller than f 5 0 kJ mol-' if one or both are of the NVT type.

IV. Results for Pure Water As a preliminary to our investigation of aqueous solutions, we examine how the properties of pure SPC, TIP4P, and MCY water are influenced by the specific variations of the boundary conditions that are the subject of this paper. We consecutively address the thermodynamics, the orientational structure, and the dynamics, wherever possible also comparing to data and observations from older, primarily reaction field simulation studies. According to Eqs 2 and 3, the surface effects depend on the dipolar and quadrupolar moments of the unit cell. These in tum will somehow reflect the size of the multipole moments of the individual molecules. Table 1 lists the dipole moments and the (traceless) quadrupole moment tensors of the three water models investigated here.49 While the dipole moments of the three models are-not surprisingly-very similar, it can be seen that there is still some clear variation among the models, with MCY water having a larger quadrupole moment than either SPC or TIP4P water. Our initial comparison for a given property and a given model is always between a minimum image cutoff (Le. a cubic cutoff that includes all molecules of the unit cell) and Ewald summation. We can think of the cubic cutoff as the limiting case of a lattice sum that has been truncated after the biggest term. Thus the two cases mentioned do not differ in anything other than the extent of replication (finite vs infinite, respectively) and the nature of the surrounding medium (vacuum vs tinfoil conditions, respectively). Thermodynamics. Table 2 gives values of some basic static properties at 25 "C, for either one of the three water models and for the two basic boundary conditions "cubic cutoff' and "Ewald". Table 3 gives more detailed information for the arbitrarily selected case of SPC water, comparing a total of five different boundary conditions and making it possible to explicitly infer the role of the surface effects.

1326 J. Phys. Chem., Vol. 99, No. 4, 1995

Roberts and Schnitker

TABLE 3: Static Properties of SPC Water at T = 25 "C and e = 0.997 g cm-3 under Various Types of Boundary Conditions CvlJ K-' E d k J mol-' Phar mol-' RK cubic cutoff -44.31 i0.04 -870 i 180 77.1 0.13 f 0.03 cubic cutoff/ -44.13 f 0.06 -730 f 170 81.6 0.82 f 0.29 tinfoil spherical cutoff -42.29 f 0.02 f 3 5 0 k 30 70.4 0.15 i 0.02 Ewald -41.55 f 0.02 f 5 1 0 i 80 75.9 5.06 i 2.49 Ewaldvacuum -41.71 f 0.03 f 7 6 0 f 70 78.7 0.16 f 0.02 Analysis of the data shows that even though potential energy, pressure, and constant volume heat capacity all vary with the boundary conditions, they are only marginally affected by changing the surface conditions. Rather, there is a significant difference between cubic and spherical cutoff conditions (see, for example, the pressure in Table 3), supplementing observations from a few older s t ~ d i e s . ~ ~We - * also ~ ~ *notice ~ that there are two opposing effects in the sense that the long-range interactions with particles in the comers of the cubic unit cell are on average attractive, while infinite replication beyond the nearest neighbors gives rise to a net repulsion. Use of vacuum conditions in conjunction with a small spherical cutoff thus implies significant error cancellation, and it is apparent how this popular approach can fortuitously yield basic thermodynamic data that are not too different from those of the infinitely replicated system under tinfoil boundary conditions. Orientational Structure. While the pair correlation functions (not shown) are rather insensitive to the boundary conditions, higher order correlations are strongly influenced, as noticed before for other water models.z0,21This is particularly evident in the orientational correlation functions that are relevant for the dielectric properties and familiar from reaction field simulation s t ~ d i e s . ~ ~We , ~ ,note ~ ~ that - ~ ~the five kinds of boundary conditions studied here are complementary to these previous studies in that they precisely exclude the standard reaction field case with a spherical truncation geometry. As an overall measure for the dielectric properties, first Tables 2 and 3 give the Kirkwood g-factors gK, Le. the mean square fluctuations of the unit cell dipole moment. Dielectric constants could be calculated from these data using appropriate fluctuation formulas, but we refrain from such calculations because our simulation runs are not of the extremely long duration that in this particular context is considered to be necessary for a meaningful analysis.31 An interesting observation in Table 3 is the fact that going from Ewald to Ewaldvacuum conditions lowers g K virtually to the cubic cutoff value, while going from cutoff to cutoff/tinfoil conditions leads to an increase of the value of gK, but still does not recover the full value under Ewald conditions. This can be rationalized if we recall that the cubic cutoff/tinfoil simulation is for technical reasons carried out with a tapering of the depolarization term, cf. eq 4,by the product of three cubic spline functions, leading to a decrease of the depolarization energy, cf. eq 3, by about 10%. Using the language of the reaction field approach, this is equivalent to a lowering of the infinite dielectric constant of the surrounding m e d i ~ m . The ~ ~ ~value ,~ 6~ is in our case obtained from

EW+1

I ,

-3

-4

I

ttir

-5 L

2

cubic cutoff / tinfoil

4

6

8

10

----

12

Y 14

1

r/A Figure 1. Orientational correlation functions gmllO(r)for SPC water at 25 "C under various types of boundary conditions, as annotated.

cutoff/tinfoil case, and there is hence a certain "asymmetry" of polarization and depolarization correction. Disregarding the latter discrepancy, it seems from comparison of the g-factors in Tables 2 and 3 that only the surface effects are really of any importance, with gK being fairly large in the two cases of tinfoil boundary conditions and small in the three cases of surrounding vacuum conditions. However, correlation function analysis shows that the situation is considerably more complicated, with the microscopic orientational correlations exhibiting exactly the opposite kind of dependencies. Figures 1 shows, for SPC water and for all five kinds of boundary conditions, the orientational correlation function goo

l0(r),

g,,"O(r)

=- 3 f i

(COS

eii)

i.e. essentially the average cosine between any two water dipoles i and j . Negative oscillations of this function correspond to a parallel alignment of the dipoles. Due to basic hydrogenbonding, nearest neighbor correlations ( r % 2.8 A) are always expected to be of this type, irrespective of the boundary conditions. Whether there is a tendency for an antiparallel alignment (positive oscillations of goo1lo) at larger distances depends on the specific boundary conditions, as discussed in several pla~es~O-'~ and not repeated here. We rather notice in Figure 1 that the influence of the surrounding medium on the microscopic orientational order is not only very subtle-with the functions for the two cubic cutoff conditions also being hardly distinguishable and the ones for the two Ewald conditions being fairly close-but actually less of a factor than both the extent of replication (the Ewald and cutoff functions are clearly different) and the mode of truncation (the spherical and cubic cutoff functions are clearly different, in agreement with early observations for ST2 waterz0). By radial integration of the function goo1l0(r), we obtain the r-dependent Kirkwood factor gK(r), shown in Figure 2 (which in the Ewald case is in agreement with what has been reported in ref 32 for SPC water',:

L' 2

With a box length of L = 18.17 A and a tapering distance of Q = 0.2 A, we have ERF x 45. This only applies to the cubic

We now see five pattems that are all distinctly different. Contact with the seemingly contradictory data for the asymptotic g-factor from above is made through the long-distance behavior of this

Simulations of Aqueous Ionic Solutions

. I Phys. . Chem., Vol. 99, No. 4, 1995 1327

TABLE 5: Dynamic Properties of TIP4P Water at T = 25 "C and = 0.997 g ~ m - Using ~ , Either a Cubic Cutoff, Ewald Summation, or a Spherical Cutoff dipole vector H-H vector tl/ps tdps 71/72 slips rdps ~1172 cubic cutoff 10.8 3.3 3.3 4.1 2.4 1.8 sphericalcutoff 4.1 1.2 3.3 4.0 1.6 2.5 Ewald" 2.9 1.1 2.6 3.4 1.4 2.4 T = 21 "C.

-4

t 2

4

6

8

10

12

14

r/A Figure 2. Integrated orientational correlation functions gK(r) for SPC water at 25 "C under various types of boundary conditions, as annotated.

TABLE 4: Dynamic Properties of SPC Water at T = 25 "C - ~ Various Types of Boundary and = 0.997 g ~ m under Conditions olio-5 dipole vector H-H vector cmz s-' tl/ps tdps 71/72 rllps tdps tlltz cubic cutoff 3.2f0.210.1 cubic cutoffltinfoil 3.3 f 0.3 9.9 sphericalcutoff 3.5 f 0.5 3.7 Ewald 4.2 k 0.1 2.4 Ewaldvacuum 3.9 k 0.2 2.7

3.0 3.0

3.3 3.3

3.0

2.1 2.0

1.3

2.8 2.7 2.9

2.9

1.3

2.2

2.6 2.4

1.1 1.1

2.3 2.1

0.87 0.94

3.4

1.6 1.5

correlation function. As recognized by S t e i n h a ~ s e rit, ~is~the ~ weighting factor of 3 in the integral of eq 9 that makes the function gK( I ) exceedingly sensitive to even the tiniest changes of the orientational correlations at longer distances. Finally, we emphasize the very pronounced difference between the spherical truncation geometry in standard reaction field simulations and the cubic reaction field conditions studied here. The former have been shown, for a number of water models, to give orientational functions gmllo(r)and r-dependent Kirkwood factors that are similar to the ones obtained under Ewald condition^,^^-^^ in particular without the big, sine-like oscillation in gK(r) that we find under the respective cubic conditions. Thus we have another example-in this case for tinfoil boundary conditions in spherical truncation geometry-where the error cancellation is such that the properties of the infinitely replicated system happen to be rather closely reproduced. Dynamics. Table 4 tabulates data for the translational and rotational diffusion of SPC water under the five kinds of boundary conditions. Reorientational times are given for two molecule-fixed axes, the dipole vector and the H-H connecting vector, and for the autocorrelation functions of both the first (tl) and second (t2) order Legendre polynomials. For ideal Brownian motion we expect ~ 1 1 ~=23. It can be seen that there is only a slight increase in the translational mobility if the surrounding medium is changed from vacuum to tinfoil conditions. The rotational mobility also has a tendency to increase, but not each and every correlation time calculated follows this trend. A generally increased mobility under tinfoil conditions would be consistent with isolated findings in some other s t u d i e ~ . ~ ~ ~ ~ ~ ~ More interesting is a prominent discrepancy which makes a distinction between cubic cutoff and Ewald conditions-regardless of the nature of the surrounding medium-and for which the

spherical cutoff case, once again, resembles much more the Ewald system. Thus we see under either one of the cubic cutoff conditions that the reorientational motion is (i) rather drastically slowed down, (ii) much more anisotropic, as clear from the fact that always t(dipo1e) >> t(H-H), and (iii) considerably less ideal, as clear from the fact that the ratios tlIt2 for the motion of both dipole and H-H vectors deviate more from the reference value of 3. These observations concerning the reorientational motion are not only characteristic for SPC water, as can be seen from the corresponding data for TIP4P water in Table 5. It is natural to attribute these substantial effects to differences in the long-range correlations with molecules in the corners of the cubic simulation cell. In Figure 1 it can in fact be seen that under cubic cutoff conditions the dipole-dipole correlations for distances r L 12 8, are on average those of clearly parallely oriented molecules. This feature is much less pronounced under spherical cutoff conditions and virtually absent under Ewald conditions. Hence, once again the behavior seen with a "suitable" kind of cutoff conditions happens to resemble the properties of the infinitely replicated system, at least to some extent.

V. Results for Aqueous Ionic Solutions Building on our results for pure SPC, TIP4P, and MCY water, we now investigate if the presence of either a C1- ion or an Fez+ ion in each of these systems leads to any new dependencies on the specification of the boundary conditions. We note that the cubic cutoff/tinfoil conditions amount to an extension of the reaction field approach-usually reserved for the simulation of pure molecular liquids-to aqueous ionic solutions, with an explicit ion-solvent reaction field term in cubic geometry. Thermodynamics. Tables 6 and 7 summarize a number of static and dynamic results for C1- and Fe2+, respectively, in either one of the three water models at 25 "C under cubic cutoff or Ewald conditions. Tables 8 and 9 show data for the entire spectrum of boundary conditions, but only for C1- in TIP4P water and for Fe2+in SPC water, respectively. As far as applicable, the energy values are consistent-within error bars-with values given in the l i t e r a t ~ r e (with ~ ~ , ~one ~ exceptions1). The ion-water potential energy is found to become systematically less attractive while going from cubic cutoff to Ewald conditions, with changes of +63, +56, and +32 kJ mol-' for C1- in SPC, TIP4P, and MCY water, respectively, and +697, $723, and +781 kJ mol-' for Fez+ in the same three water models, respectively. This is in agreement with observations in the Br- in water system.13 We also find that the differences are basically undiminished if the ion-solvent energies of a given simulation-regardless of how the latter was carried out-are calculated with both methods and subsequently compared.2 The effect is entirely due to the difference in the nature of the surrounding medium and can therefore easily be corrected for by polarization (or depolarization) correction. This is clear from the data in Tables 8 and 9, where we see, for example, that the energy between Fez+ and SPC water is - 1975 kJ mol-' under Ewald conditions and almost exactly the same (- 1978

1328 J. Phys. Chem., Vol. 99, No. 4, 1995

Roberts and Schnitker

TABLE 6: Simulations of C1- in SPC, TIP4P, or MCY Water at 25 "C, Using Either a Cubic Cutoff or Ewald Summation C1- in SPC C1- in SPC C1- in TIP4P C1- in TIP4P C1- in MCY C1- in MCY

cubic cutoff Ewald cubic cutoff Ewald cubic cutoff Ewald

-606 f 11 -543 f 8

+302 $250

-304 -293

0.15 i 0.01 6.49 i 1.57

-605 f 2 -549 i 2

$341 +257

-264 -292

0.18 i 0.01 3.20 i 0.49

-672 f 9 -640 f 4

+315 f247

-357 -393

0.14 i 0.01 2.32 f 0.83

TABLE 7: Simulations of Fez+ in SPC, TIP4P, or MCY Water at 25 "C, Using Either a Cubic Cutoff or Ewald Summation Fez+in SPC Fez+in SPC Fez+in TIP4P Fe2+in TIP4P Fez+in MCY Fez+in MCY

cubic cutoff Ewald cubic cutoff Ewald cubic cutoff Ewald

-2672 f. 1 -1975 f 3

$998 $810

- 1674 -1165

0.15 f.0.01 2.83 f.0.63

-2515 f 15 - 1792 4z 6

$982 +754

- 1533 -1038

0.15 f 0.03 3.29 i0.55

-2640 f 12 - 1859 f 30

+975 $804

-1665 - 1055

0.16 i 0.01 2.42 f 0.75

TABLE 8: Simulations of C1- in TIP4P Water at 25 "C under Various Types of Boundary Conditions &-watfld mol-' E,,,,lkT mol-' EhydlkT mol-' gK cubic cutoff -605 i 2 $341 -264 0.18 f 0.01 cubic cutoff/tinfoil -547 f 3 +272 -275 1.11 f 0.36 $326 -264 0.18 f 0.01 spherical cutoff -590 f 5 -549 i 2 1-257 -292 3.20 f 0.49 Ewald Ewaldvacuum -609 f 6 +296 -313 0.17 i 0.01

D,,,/10-5 cm2SKI

0.5 f 0.1 1.0 f 0.1 0.6 i 0.1 2.2 i 0.8 0.6 i 0.1

TABLE 9: Simulations of Fe2+ in SPC Water at 25 "C under Various Types of Boundary Conditions cubic cutoff cubic cutoff/tinfoil Ewald Ewaldvacuum

-2672 i 1 -1978 f 9

+998 +819

-1674 -1159

0.15 f 0.01 0.71 f 0.11

0.12 i 0.03 0.52 i 0.10

-1975 f 3 -2695 f 15

+810 +loo8

-1165 -1687

2.83 f 0.63 0.17 f 0.02

0.73 f 0.19 0.02 f 0.01

kJ mol-') under cubic cutoff/tinfoil conditions. Also, polarization and depolarization correction are in this case completely "symmetrical". Besides this major effect, there is also a slight difference between cubic and spherical cutoff conditions, as clear from the C1- data in Table 8. This difference is consistent with a certain dependence of the thermodynamics of hydration on the size of the cutoff radius that has been reported for several system^.^-^.'^ The fact remains that this dependence is still fairly minor relative to the surface effects, at least for the cutoff radii that are typical for standard simulations. From Tables 6 and 7 it can be seen that a switch from cubic cutoff to Ewald conditions also affects the solvent reorganization energy in a systematic way. Under Ewald conditions it is consistently smaller, namely, by -52, -82, and -68 kJ mol-' for C1- and by -188, -228, and -171 kJ mol-' for Fe2+ in SPC, TIP4P, and MCY water, respectively. (Apparently, the TIP4P model is most sensitive to this particular change of boundary conditions.) From the data for the Fe2+ ion in Table 9 and-although somewhat less clearcut-also from the data for C1- in Table 8, one concludes that this variation is again predominantly a surface effect and can potentially be accounted for by polarization or depolarization correction. If we now consider the experimentally relevant hydration energy (=difference of the total energies of the systems with and without solute), we notice that the two contributing terms (ion-solvent energy and reorganization energy) show surface effects of opposite sign. In the case of C1-, they are also of the same magnitude, with the result that the hydration energy is relatively insensitive to the boundary conditions, changing by $11, -28, and -36 kJ mol-' if going from cubic cutoff to Ewald conditions, for each of the three water models, respec-

tively. In the case of Fe2+, on the other hand, the surface effects associated with the ion-water potential energy are clearly dominant so that there is a manifest change also of the hydration energy (by +509, +495, and +610 kJ mol-' for the three water models, respectively). Structure. Paralleling the pure solvent behavior, we find the ion-solvent pair correlation functions (not shown) to be rather insensentive to the specification of the boundary conditions. For C1- in TIP4P water, the ion-oxygen g(r) has its first peak always at r w 3.2 A, with a peak height of ~ 3 . and 8 a first shell coordination number of ~ 7 . 2 .There is a very small variation in the height of the second peak, which is found at r z 4.9 A with a height of ~ 1 . 2 - 1 . 3and a second shell (running) coordination number of ~ 2 9 - 3 0 . These observations are consistent with other simulation results for halide ions in ~ a t e r . ~We . ' ~do not see anything reminiscent of the pronounced change in the height of the second peak with changing cutoff radius that was found in an integral equation study of the C1in water system.6 The ion-oxygen correlations for the doubly charged Fe2+ ion are equally uniform. In SPC water, the corresponding g(r) has its first peak always at r x 2.1 A, with a peak height of ~ 1 2 . and 5 a first shell coordination number of exactly 6, and a second peak at r x 4.3 A, with a height of ~ 1 . 9 - 2 . 0and a second shell coordination number of ~ 2 0 - 2 2 . The relatively most pronounced difference we can identify is that the second and third shells of the ion-hydrogen g(r) are somewhat more clearly separated under surrounding vacuum conditions than under tinfoil boundary conditions. Even this is only a weak effect. The solute-solvent angular distribution functions, as represented by the angles of molecule-fixed vectors relative to the

J. Phys. Chem., Vol. 99, No. 4,1995 1329

Simulations of Aqueous Ionic Solutions ion-oxygen vector, are equally invariable, with molecules that are on average bond-oriented relative to C1- and reverse-dipoleoriented relative to Fez+. Differences between the various boundary conditions only become evident if the solvent-solvent correlations very close to the ion are examined. For example, the orientational correlation functions goo1l0(r)for all molecules that are within 4 8, of the ion (approximately the first hydration shell) show some variation with the type of boundary conditions (not shown) that in all likelihood can simply be attributed to the varying water structure as such of these systems. The variation of the structure of the “bulk water” component of the ionic solutions is in fact entirely analogous to what was described before for the neat liquids. For example, if we only include molecules at distances I 2 8 8, from the ion, we find (not shown) essentially the same orientational correlation functions gw”O(r) and gK(T) as in the previous Figures 1 and 2. Also, the variation of the Kirkwood g-factor, averaged over all molecules and given in Tables 6-9, is similar to what was seen before for pure water. Dynamics. A truly dramatic dependence on the type of boundary conditions is found for the ionic diffusion coefficient Di,. The values given in Table 8 for C1- in TIP4P water and in Table 9 for Fez+ in SPC water span a range of a factor of %5 in either case. This is an effect that is much more pronounced than the one for the translational diffusion of the pure solvent (see Tables 2 and 3). A similar strong speedup of the ionic diffusion while proceeding from cutoff to Ewald conditions was also seen in simulations of Br- in SPC water under both sets of boundary ~0nditions.l~ For unknown reasons, all our calculated diffusion coefficients for Fez+ in SPC water (Table 9) are smaller than what has been reported for this system with a 7.5 8, cutoff, regardless whether simulated with or without inclusion of solvent f l e ~ i b i l i t y . ~ ~ Interestingly enough, cubic and spherical cutoff conditions for C1- in TIP4P water yield almost the same result, contrary to what is seen for some other properties. From the data in Tables 8 it rather becomes clear that it is basically the nature of the surrounding medium that determines the ionic diffusion behavior. D,,,is equally small under cubic cutoff, spherical cutoff, and Ewaldvacuum conditions, and it is large under Ewald conditions. An approximately intermediate value is found for cubic cutoff/tinfoil conditions, probably to be rationalized by the finite dielectric constant of the surrounding medium (see section IV). In Figures 3 and 4, the ionic velocity autocorrelation functions and the corresponding power spectra under cubic cutoff conditions are opposed to those under Ewald conditions, for C1- in TIP4P water and for Fez+ in SPC water, respectively. In all cases we find an oscillatory, non-Brownian behavior of the time correlation functions (insets of figures), consistent with previous reports for small ions in various water models, including C1under c ~ t o f f ,C1~ ~under ,~~ and Fez+ under cutoff condition^.^^ In no case does the decay of the velocity correlations look as Brownian as reported for some heavier ions, like Br- under Ewald13 and Cs+ under cutoff conditions.s3b The power spectra (main panels of figures) show that the accelerated diffusion under Ewald conditions goes along with a clear shift of the vibrational spectrum to lower frequencies, for both anion and cation. However, even under Ewald conditions the line shapes still have to be described as a superposition of several Lorentzians. VI. Conclusions

Prior to the work of De Leeuw et al.,14 there was very little awareness that Ewald summation as a convenient tool for the

-0.4

0.2

0.1

0.0

tips

I 0

25

75

50

100

0 f ps-’ Figure 3. Power spectra (main panel) as obtained from velocity autocorrelation functions (inset) for C1- in TP4P water at 25 “C, using either a cubic cutoff (solid lines) or Ewald summation (dashed lines).

‘--‘‘,

, 0

25

50

75

100

0 1ps-1

Figure 4. Power spectra (main panel) as obtained from velocity autocorrelation functions (inset) for Fez+in SPC water at 25 “C, using either a cubic cutoff (solid lines) or Ewald summation (dashed lines). enforcement of lattice sum convergence does not only involve a mathematical trick but also implies a fundamental change of the physical situation, namely, a switch from nonconducting to conducting boundary Conditions. It seems that in more recent applications there will at most be a disclaimer that the technique is “used in absence of the surface term of De Leeuw et al.”. As can be seen from the summary of our results in Chart 1, there are in fact a few properties that depend quite significantly on the presence of this term in the Hamiltonian. What is the nature of the unit cell surroundings that is physically “desirable” in a simulation? In crystals, particularly pyroelectric or ferroelectric crystals, surface polarization effects can truly be expected and hence vacuum conditions are to be preferred. Paradoxically, these conditions are not realized in the standard solid-state physics technique for calculating lattice sums, Le. Ewald summation. In any disordered system of macroscopic size, however, there is no distinct surface structure and accordingly polarization effects should vanish. Thus surface polarization in simulated liquids is only present as a consequence of our having to deal

1330 J. Phys. Chem., Vol. 99, No. 4, 1995 CHART 1: Summary of Boundary Condition Dependencies of Aqueous Ionic Solution Properties Primarily dependent on nature of surrounding medium (vacuum or tinfoil): asymptotic Kirkwood factorg, ion-solventpotential energy solvent reorganization energy ionic mobility Primarily dependent on unit cell replication (toroidal or infinitely periodic) and mode of interaction truncation (cubic or spherical): pressure of solvent orientational order of solvent,g,”O (r) translationalmobility of solvent reorientational mobility of solvent Sensitive to all features of the boundary conditions: distance-dependentKirkwood factor g, (r) Minor sensitivity to boundary conditions: potential energy of solvent heat capacity of solvent solvent-solventpair correlation functions ion-solventpair correlation functions

with unit cells of very small size. Tinfoil boundary conditions successfully correct for this small system artifact. Of course, the spurious dipolar and quadrupolar moments of a liquid state unit cell automatically become less and less important as system size and cutoff are increased. Eventually the nature of the surrounding medium becomes irrelevant as we approach the thermodynamic limit. However, this smooth connection between cutoff and Ewald conditions is not easily established even with what is currently considered to be “large” system sizes and “large” cutoff radii. For example, in a recent simulation of an aqueous solution that included 2200 water molecules, results obtained with a spherical cutoff of 20 8, and with Ewald summation were still considerably different. l3 We thus arrive at the conclusion that Ewald summation-while usually considered to be the correct method for the description of crystals-is on the contrary expected to be particularly appropriate for the simulation of liquids in small unit cells, due to the elimination of the surface term. In the same vein, it can be argued that the reaction field method is also an appropriate method for the simulation of liquids. Now, do the systems modeled in this study under tinfoil conditions really give the best description of the actual systems? Unfortunately, a comparison of simulated and experimental data is not straightforward, both because of general shortcomings of the models-not necessarily related to choosing the “right” calculational method for the long-range interactions-and because of uncertainty in the experimental information. In addition, there is the compounding problem of a certain bias that is tacitly introduced through the boundary conditions that are inherent in the method for the construction of the potential functions. A good example is the thermodynamics of ionic hydration. Our simulated C1- ion always exhibits a coordination number that is presumably too large,56thus making a direct comparison of calculated and measured hydration energies questionable. Just as problematic is the solvation in the Fe2+ system since for multiply charged ions it is particularly unclear if they can be adequately modeled with simple pair potential^.^^ Finally, experimental estimates for the hydration energy can be quite inaccurate. If we nevertheless compare with standard hydration enthalpies from the literatureSRof -381 W mol-’ for C1- and -1946 W

Roberts and Schnitker mol-’ for Fe2+, we find that these values are consistently underestimated in all simulations (Tables 6-9). In particular for the case of Fe2+, however, it is striking that the calculated hydration energies under vacuum conditions are always much closer to the experimental figure than those under tinfoil conditions (Tables 7 and 9). However, before we interpret these findings as providing support for the conventional cutoff method, we need to recall that ion-solvent potential energies are (normally) also used in the parametrization of the potential functions. These energies originate either in quantum mechanical calculations of the energy surface of the dimer or-less commonly-in energetic information extracted from series of simulation runs of the complete system (typically with a short spherical cutoff for computational efficiency) that scan the space of potential function parameters. In either case, the potential function is effectively parametrized for vacuum conditions. Thus any relative closeness of experimental hydration energies and simulated values under cutoff conditions is easily rationalized. It can also be suspected that the experimental data could equally well be reproduced in simulations under Ewald conditions if we used potential functions parametrized in explicit anticipation of the subsequent simulation conditions. Nevertheless, even with the current potential functions, our systems exhibit some properties that are best reproduced using Ewald summation. Of course, the often studied dielectric properties have to be mentioned first, the central point being that simulations carried out with Ewald summation are devoid of spurious long-range orientational correlations (see Figure 1). Another property that is only obtained “correctly” with the Ewald technique is the diffusion coefficient of the C1- ion (cf. our simulated value of 2.2 x cm2 s-l with the experimental59 figure of 2.064 x cm2 s-l). Hence at least some observations are directly in line with what has been asserted above (and by other investigators), namely, that Ewald summation generally corresponds to the physically most desirable situation. Interestingly enough, the analysis carried out here implies that we can argue for the superiority of Ewald summation in a rather nontraditional way, namely, by specifically attributing it to the tinfoil boundary conditions that are inherent to the technique. As ample computational resources are becoming available to a larger and larger community, the interest in the Ewald technique in fact only appears to be growing, with a steady flow of suggestions for improvements60 and generalization@ of the method. Still, it is a computationally expensive technique, making the reaction field method appealing as a cheap altemative that can also provide tinfoil boundary conditions. An interesting choice to be made in reaction field simulation studies is the value of the continuum dielectric constant ERF in the H a m i l t ~ n i a n . ~ The ~ * ~traditional ~ - ~ ~ view is that the reaction field provides a long-range correction to the truncated Coulombic interactions, and therefore-in order to avoid discontinuities on the length scale of the cutoff radius-the input dielectric constant EW should be the same as the calculated value E of the molecular model itself (Le. ERF should ideally be determined in a self-consistent way). We have seen, however, that the reaction field serves the equally important purpose of eliminating the spurious surface effects in a small simulation cell, thus ideally calling for ERF = =. It is clear that both goals-proper long-range correction and elimination of all surface effects-cannot be simultaneously achieved, and either choice (ERF = E or ERF = =) implies some compromise with respect to the physically desirable situation. A most interesting byproduct of our surface-potential-based

J. Phys. Chem., Vol. 99, No. 4, 1995 1331

Simulations of Aqueous Ionic Solutions interpretation of the reaction field expression is the conclusion that the approach is in fact not limited to the study of pure liquids and-contrary to some belief-may well be extended to systems with charged species, i.e. ionic solutions. Naturally, for a good overall reproduction of experimental properties ion-solvent potential functions should be used that are explicitly compatible with this specification of the boundary conditions. By far the most ambiguous issue is the suitability of the simplest and still most popular of all methods, i.e. a simple spherical truncation with surrounding vacuum conditions. On the one hand, there can be little argument that this technique leads to incorrect long-range orientational correlations, and thus all properties that are sensitive to the orientational structure, like the dielectric constant30and possibly even phenomena like the stability of protein helices in aqueous s0lution,4~~ cannot be realistically modeled. Nevertheless, it is also true that systems simulated with a short spherical cutoff can sometimes-by virtue of error cancellation-qualitatively mimic the properties of Ewald systems (see Tables 3-9, as occasionally emphasized by other investigators.21.62Of course, spherical truncation also continues to be an appealing technique for the simple reason of computational simplicity. Finally, we need not forget that this technique-as discussed above-is under normal circumstances the one that is most compatible with the standard methods for the parametrization of potential functions. Thus there are at least a few reasons why the use of small spherical cutoffs can occasionally still be considered as a viable altemative. In general, however, cutoffs should be used in conjunction with reaction field terms for both solvent-solvent and ionsolvent interactions and together with potential functions that have been parametrized in anticipation of this choice of the boundary conditions. Acknowledgment is made to the Donors of The Petroleum Research Fund, administered by the American Chemical Society, for partial support of this research (Grant No. 25299-G6). The computations were performed with resources provided by an NSF Instrumentation Grant (No. CHE-8917309) to the Department of Chemistry of the University of Michigan. References and Notes (1) For an excellent review of the older literature see: The Problem of Long Range Forces in the Computer Simulation of Condensed Media; Ceperley, D., Ed.; NRCC Proc. 9; Lawrence Berkeley Laboratory, 1980. (2) Roberts, J. E.; Schnitker, J. J . Chem. Phys. 1994, 101, 5024. (3) Bader, J. S.; Chandler, D. J . Phys. Chem. 1992, 96, 6423. (4) (a) Schreiber, H.; Steinhauser, 0. Biochemistry 1992, 31, 5856; (b) J . Mol. Biol. 1992, 228, 909; (c) Chem. Phys. 1992, 168, 75. (5) York, D. M.; Darden, T. A.; Pedersen, L. G. J . Chem. Phys. 1993, 99, 8345. (6) Brooks, C. L., 111 J . Chem. Phys. 1987, 86, 5156. (7) Madura, J. D.; Pettitt, B. M. Chem. Phys. Letr. 1988, 150, 105. (8) Straatsma, T. P.; Berendsen, H. J. C. J . Chem. Phys. 1988, 89, 5876. (9) Kim, K. S. Chem. Phys. Lett. 1989, 156, 261. (10) Smith, P. E.; Pettitt, B. M. J . Chem. Phys. 1991, 95, 8430. (11) Lee, F. S.; Warshel, A. J . Chem. Phys. 1992, 97, 3100. (12) Hummer, G.; Soumpasis, D. M.; Neumann, M. Mol. Phys. 1992, 77, 769. (13) Del Buono, G. S.; Cohen, T. S . ; Rossky, P. J. Mol. Liq. 1994, 60, 221. (14) De Leeuw, S. W.; Perram, J. W.; Smith, E. R. Proc. R . SOC.London 1980, A373, 27. (15) Redlack, A.; Grindlay, J. J . Phys. Chem. Solids 1975, 36, 73.

(16) Euwema, R. N.; Surratt, G. T. J . Phys. Chem. Solids 1975,36, 67. (17) Massidda, V. Physica 1978, 9 5 4 317. (18) Stewart, S. N. J . Comput. Phys. 1978, 29, 127. (19) Klein, M. L.; McDonald, I. R. J . Chem. Phys. 1979, 71, 298. (20) Pangali, C.; Rao, M.; Beme, B. J. Mol. Phys. 1980, 40, 661. (21) An&ea, T. A.; Swope, W. C.; Andersen, H . C. J . Chem. Phys. 1983, 79, 4576. (22) Brooks, C. L., III; Pettitt, B. M.; Karplus, M. J . Chem. Phys. 1985, 83, 5897. (23) Prevost, M.; Van Belle, D.; Lippens, G.; Wodak, S. Mol. Phys. 1990, 71, 587. (24) Alper, H. E.; Bassolino, D.; Stouch, T. R. J . Chem. Phys. 1993, 98, 9798. (25) Adams, D. J. Chem. Phys. Lett. 1979, 62, 329. (26) Evans, M. W. Comput. Phys. Commun. 1990, 59, 495. (27) Barker, J. A,; Watts, R. 0. Mol. Phys. 1973, 26, 789. (28) Adams, D. J.; Adams, E. M.; Hills, G. J. Mol. Phys. 1979, 38, 387. (29) De Leeuw, S. W.; Perram, J. W.; Smith, E. R. Proc. R. SOC.London 1983, A388, 177. (30) (a) Steinhauser, 0. Mol. Phys. 1982, 45, 335; (b) Ber. BunsenGes. Phys. Chem. 1983, 87, 128; (c) Chem. Phys. 1983, 79, 465. (31) (a) Neumann, M. J. Chem. Phys. 1985, 82, 5663; (b) 1986, 85, 1567. (32) Anderson, J.; Ullo, J. J.; Yip, S. J . Chem. Phys. 1987, 87, 1726. (33) Alper, H. E.; Levy,R. M. J . Chem. Phys. 1989.91, 1242. Belhadj, M.; Alper, H. E.; Levy, R. M. Chem. Phys. Lett. 1991, 179, 13. (34) Zhu, S.-B.; Wong, C. F. J . Chem. Phys. 1993, 98, 8892. (35) Smith, P. E.; Van Gunsteren, W. F. J . Chem. Phys. 1994, 100, 3169. (36) Felderhof, B. U. Physica 1980, 101A, 275. (37) Smith, E. R. Proc. R. SOC.London 1981, A375, 475. (38) Olives, J. J . Phys. Lett. 1985, 46, L-1143. (39) Tasaki, K.; McDonald, S.; Brady, J. W. J . Comput. Chem. 1993, 14, 278. (40) Berendsen, H. J. C.; Postma, J. P. M.; Van Gunsteren, W. F.; Hermans, J. In Intermolecular Forces; Pullman, B., Ed.; Reidel: Dordrecht, 1981. (41) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. J . Chem. Phys. 1983, 79, 926. (42) Matsuoka, 0.;Clementi, E.; Yoshimine, M. J . Chem. Phys. 1976, 64, 1351. (43) Chandrasekhar, J.; Spellmeyer, D. C.; Jorgensen, W. L. J . Am. Chem. SOC.1984, 106, 903. (44) Curtiss, L. A,; Halley, J. W.; Hautman, J.; Rahman, A. J . Chem. Phys. 1987, 86, 2319. (45) Sangster, M. J.; Dixon, M. Adv. Phys. 1976, 25, 247. (46) Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids; Clarendon: Oxford, 1987. (47) Cheung, P. S. Y. Mol. Phys. 1977, 33, 519. (48) We note that we found many heat capacity values in the literature to be notoriously wrong, presumably because of confusion about the degrees of freedom and the use of the correct fluctuation formula. (49) Watanabe, K.; Klein, M. L. Chem. Phys. 1989, 131, 157. (50) Retardation of the translational water motion under Ewald conditions has been reported, but was later shown to be erroneous; see: Teleman, 0.;Wallqvist, A. Int. J . Quantum Chem. 1992, 41, 381. (51) Sprik, M.; Klein, M. L.; Watanabe, K. J . Phys. Chem. 1990, 94, 6483, give an ion-water potential energy for C1- in TIP4P water that is very close to our cutoff values even though the simulations were carried out using Ewald simulation. A difference between simulation and analysis conditions seems to be the most likely explanation. (52) Gukdia, E.; Padr6, J. A. Chem. Phys. 1990, 144, 353. (53) (a) Berkowitz, M.; Wan, W. J . Chem. Phys. 1987, 86, 376. (b) Reddy, M. R.; Berkowitz, M. J . Chem. Phys. 1988, 88, 7104. (5i) Impey, R. W.; Madden, P. A.; McDonald, I. R. J . Phys. Chem. 1983, 87, 5071. (55) Neuven. H. L.: Adelman. S. A. J . Chem. Phvs. 1984. 81.4564. (56) Eiddrby, J. Ann. Rev. Phys. Chem. 1983, 34; 155. (57) Probst, M. M. Chem. Phys. Lett. 1987, 137, 229. (58) Cotton, F. A.; Wilkinson, G. Advanced Inorganic Chemistry, 5th ed.; Wiley: New York, 1988. (59) Hawlicka, E. Z. Naturforsch. 1986, 4Ia, 939. (60) Darden, T.; York, D.; Pedersen, L. J . Chem. Phys. 1993, 98, 10089. (61) Hautman, J.; Klein, M. L. Mol. Phys. 1992, 75, 379. (62) Linse, P.; Andersen, H. C. J . Chem. Phys. 1986, 85, 3027. JF'941374F