Effect of the treatment of long-range forces on the dynamics of ions
in aqueous solutions
Lalith Perera, Ulrich Essmann, and Max L. Berkowitz
Department of Chemistry, University of North Carolina, Chapel Hill, North Carolina 27599
͑Received 25 July 1994; accepted 15 September 1994͒
The goal of the present work is to study the dependence of the limiting ionic mobility of such anions
as ﬂuoride, chloride, and bromide in water on the way the long-range forces are treated in the
computer simulations. With this in mind we have performed molecular dynamics computer
simulations where the long-range electrostatic forces were treated using: ͑a͒ simple truncation
procedure, ͑b͒ energy switching procedure, ͑c͒ reaction ﬁeld method, and ͑d͒ Ewald summation
technique. Our analysis shows that the switching procedure with the short-range switching function
introduces artifacts into the simulations. These artifacts are responsible for the faster decay and
oscillations in the velocity autocorrelation function of the ions and therefore for the lower value of
the diffusion coefﬁcients. © 1995 American Institute of Physics.
I. INTRODUCTION
We are witnessing now a rapid development of computational
hardware that will allow us to perform molecular
dynamics computer simulations on systems with large number
of particles and over long simulation times. To utilize the
potential of the hardware, we have to develop reliable software
and therefore resolve many challenging problems and
answer many questions, old and new. Among the questions
that need to be answered we ﬁnd, for example, the question
of how the treatment of long-ranged Coulomb forces inﬂuences
the structure and dynamics of the system when we
simulate an aqueous solution in the bulk. Many recent computer
simulations clearly demonstrated that the values of the
structural and thermodynamic properties of aqueous solutions
strongly depend on the treatment of the long-range
electrostatic forces.1–6
The usually adopted methods to treat long-ranged forces
are to use some kind of cutoff procedure or the Ewald summation
technique. In addition to these two methods, new
ones, based on multipole expansions are now developed,7
but
we are not going to address them in this work. The cutoff of
the electrostatic interactions can be achieved by ͑i͒ applying
a simple truncation procedure, ͑ii͒ multiplying the electrostatic
potential or force by some kind of a switching function
or a shifting function ͑when the width of a switching function
is equal to zero, the switching procedure is equivalent to
the truncation procedure͒, ͑iii͒ using the reaction ﬁeld
method. All of these cutoff procedures and the Ewald summation
method introduce approximations into the treatment
of the real systems and have their advantages and disadvantages.
The use of switching or shifting functions modiﬁes
Coulomb interaction and forces it to zero after a certain
distance.8
In the reaction ﬁeld method the particles beyond
the sphere with the certain cutoff radius Rc are represented
by a dielectric continuum.9
The advantage of the cutoff
schemes is that they reduce computational demands on the
system. In the Ewald summation method the periodic boundary
conditions are fully utilized and the lattice sum is
calculated.10
Sometimes this technique may produce artiﬁcial
results due to the enhancement of the periodicity, while the
real system is not periodic.
In many simulations the results depend on the way the
electrostatic forces are treated. An interesting example,
which illustrates dramatic inﬂuence of such a treatment on
the structural and thermodynamic properties of aqueous solutions,
is the calculation of the potential of mean force between
two chloride anions in water. Earlier computer simulations
predicted the ClϪ
–ClϪ
pairing in aqueous solutions.11
Later studies showed that the pairing disappeared when the
Ewald summation procedure was used.3
A similar effect was
observed for the ferrous-ferric system.12
While most of the previous work was devoted to the
study of how the treatment of long-range forces inﬂuences
structural and thermodynamics properties, much less is
known of how these treatments inﬂuence the dynamical
properties of the system. If one starts the study of thermodynamic
properties of aqueous solutions with the potential of
mean force, the proper place to start the study of the dynamical
properties of aqueous solutions is to investigate the inﬂuence
of long-range force treatment on the behavior of the
ionic mobility, particularly the limiting ionic mobility. This
was done recently by Rossky and collaborators13
who
showed that for the bromide ion the velocity autocorrelation
function and the value of the diffusion coefﬁcient of this ion
in water strongly depend on the treatment of long-range
forces. They also found that the velocity autocorrelation
function of bromide looked similar to the one of the Brownian
particle when the Ewald summation technique was used.
When the switching procedure was applied to the system the
velocity autocorrelation function of the ion had an oscillatory
character. Such oscillations in the velocity autocorrelation
function were also observed previously for chloride ion in
the simulations where the switching of the potential was
used.14,15
In addition, the value of the diffusion coefﬁcient
obtained in Ref. 14 was around half the value measured in
the experiment. The above-mentioned studies13–15
prompted
us to ask the following questions: ͑a͒ how does the value of
the chloride diffusion coefﬁcient change when the Ewald
summation is used to treat the long-ranged electrostatic
force? ͑b͒ how does the qualitative character of the chloride
450 J. Chem. Phys. 102 (1), 1 January 1995 0021-9606/95/102(1)/450/7/$6.00 © 1995 American Institute of Physics
velocity autocorrelation function change with the change in
the long-range force treatment? While in the present work
we will attempt to answer the above-asked questions related
to chloride, the main goal of the present study is to understand
how the treatment of the long-range forces inﬂuences
the limiting ionic mobility of the anions.
II. TREATMENT OF THE LONG-RANGE FORCES
To treat the long-range forces we use four procedures in
this work. They are: ͑i͒ the simple truncation procedure, ͑ii͒
a switching function procedure, ͑iii͒ the reaction ﬁeld
method, and ﬁnally ͑iv͒ the Ewald sum technique. The truncation
procedure is the simplest one: in it we consider two
molecules as interacting if the distance between their centers
of mass is less than a cutoff radius Rc . If the distance is
larger than Rc the interaction is neglected. That means that
we introduce an artiﬁcial discontinuity in our potential using
this procedure. In the simulations with the switching procedure
the discontinuity in the potential is avoided through the
multiplication of the interaction with the switching function
that starts reducing the potential at some distance RL and
reduces it continuously to zero at some cutoff distance Rc .
Therefore, the interaction between two particles V1,2 is written
as
V1,2ϭ͚i,j
Vel͑ri,j͒S͑R͒, ͑1͒
where the sum runs over all sites i on particle 1 and all sites
j on particle 2. R is the distance between two reference
points on particle 1 and particle 2. In our simulations we
used the following switching function:16
S͑R͒ϭͭ1 if RрRL
1Ϫ͑RϪRL͒2
͑3RcϪRLϪ2R͒/͑RcϪRL͒3
0 if RуRc
. ͑2͒
Other forms of the switching function exist, but we use this
particular one to be able to compare our present results with
the previous simulations where the ionic mobility of the
chloride ion was investigated.14,15
In our present simulations
R is the distance between oxygen atoms of water molecules
when the interaction between two water molecules is considered.
It is the distance between the oxygen of water and the
ion, when the interaction between the ion and a water molecule
is considered.
As was recently shown by Steinbach and Brooks17
͑and
as we shall see below͒ the switching procedure may result in
the creation of large sudden changes in the force. To avoid
these changes, the shifting procedure is used to turn off the
potential continuously over long distance until it becomes
zero at Rc . For example, some of the shifting functions used
in the simulations of water and aqueous solutions have the
general form18
S͑R͒ϭ͓1Ϫ2͑R/Rc͒n
ϩ͑R/Rc͒2n
͔␪͑RcϪR͒. ͑3͒
In Eq. ͑3͒ ␪ is the Heaviside step function and n is some
integer. Usually this shifting function is applied to every
site–site Coulomb interaction, thus modifying it over the entire
range of distances. For nϭ1 and nϭ2 the resulting
modiﬁed electrostatic interactions were called MEI4 and
MEI5, respectively, by Brooks et al.8
The shifting functions
in these cases were called F1 and F2 by Prevost et al.18
Steinbach and Brooks17
recently introduced the force shift
functions, that monotonically turn off the force at the distance
Rc . With the application of such a function the force
acting on site i ͑with charge qi͒ due to its Coulomb interaction
with site j ͑charge qj͒ at a distance R, has the following
form:
F͑R͒ϭ
qiqj
R2 ͓1Ϫ͑R/Rc͒␤ϩ1
͔
R
R
␪͑RcϪR͒, ͑4͒
where ␤ is some integer.
In the reaction ﬁeld method that we use in this work we
check the distance between centers of mass of two molecules.
If it is smaller than Rc we assume that molecules
interact and the force of interaction between sites i and j is
given by the following expression:3
F͑rij͒ϭ
qiqj
rij
2 ͓1Ϫ͑rij /Rc͒3
͔
rij
rij
␪͑RcϪR͒. ͑5͒
Equation ͑5͒ shows that the use of a reaction ﬁeld method is
equivalent to an application of a force shifting function of
Eq. ͑4͒ with ␤ϭ2.
Finally in the Ewald summation method we take full
advantage of the periodic boundary conditions by periodically
repeating our simulation box. To achieve convergence
of the Coulomb sum we build up the periodic system in
roughly spherical layers. Outside the sphere we have to
specify the nature of the dielectric continuum surrounding
the sphere, i.e., specify the dielectric constant ⑀r of the continuum.
There is a rather simple connection between the energy
of the very large sphere of boxes surrounded by an ideal
conductor V(⑀rϭϱ͒ and the energy of the system embedded
into vacuum V(⑀rϭ1͒. The connection is given by the following
expression:19
V͑⑀rϭϱ͒ϭV͑⑀rϭ1͒Ϫ
2␲
3␯ ͚ͯi
qiriͯ
2
, ͑6͒
where ␯ is the volume of the unit cell. When the large sphere
is in the vacuum, it has a dipolar layer on its surface. The last
term in Eq. ͑5͒ is responsible for removing this layer, when
the sphere is embedded in perfect conductor ͑tin foil boundaries͒.
We run our Ewald simulations with these tin foil
boundary conditions.
III. MOLECULAR DYNAMICS SIMULATIONS
In the present study, we have selected three systems involving
ﬂuoride, or chloride, or bromide ion solvated in 500
water molecules as well as a system involving one bromide
ion solvated in 215 water molecules. Such a selection of
different systems allows us to observe how the properties
change in a molecular dynamics ͑MD͒ simulation due to the
variation of the cutoff distance and the size of the ion. In the
simulations with the ﬂuoride and bromide ions the SPC water
model20
was chosen, while for the chloride ion, the
TIP4P water model21
was used. The choice of different water
451Perera, Essmann, and Berkowitz: Dynamic of ions in aqueous solutions
J. Chem. Phys., Vol. 102, No. 1, 1 January 1995
models was primarily due to the fact that we can compare
our current results with some of the previously obtained results
from solvation studies of the above-mentioned ions under
various simulation schemes.13–15
Also, we stress here
that it is not our intention in the present study to emphasize
the model performance relative to the reality, but to critically
evaluate methods that are involved in the treatment of longrange
interactions. All three ions were represented by a point
charge having a Lennard-Jones ͑LJ͒ center on it. The
Lennard-Jones parameters of the ﬂuoride–SPC water22
and
the chloride–TIP4P water23
interactions were the ones reported
previously in the literature. The LJ parameters of the
bromide–SPC water interaction were obtained after combining
the parameters reported previously24
for bromide with
the SPC water parameters. In all these cases we can write the
interaction potential between the ion and a water molecule as
VwIϭ
AwI
RoI
12Ϫ
CwI
RoI
6 ϩ͚j
qIqj
rIj
͑7͒
and the interaction potential between two water molecules as
Vwwϭ
Aww
Roo
12 Ϫ
Cww
Roo
6 ϩ͚ij
qiqj
rij
. ͑8͒
In Eqs. ͑7͒ and ͑8͒ RoI is the distance between the oxygen of
the water molecule and the ion center. Roo is the distance
between oxygens of two water molecules. All the potential
parameters along with the charges on the ions used in the
current study are presented in Table I.
All our simulations were carried out in the NVT ensemble
and the density was always selected to be 1 g/cc at
temperature 300 K. For simulations with any one of the
above-mentioned ions solvated in 500 water molecules, the
required box length was 24.66 Å, whereas in the simulations
with 215 water molecules the value was 18.62 Å.
Numerical evaluations of the trajectories were done using
the leap-frog algorithm and the internal geometries of
water molecules were constrained through the Shake procedure.
To achieve the desired temperature, all the systems
studied were coupled to a thermal bath25
with a coupling
constant of 0.1 ps. All the trajectory calculations were performed
with a time step of 2 fs. In all simulations the equilibration
period of 40 ps was followed by a period of 120 ps
for the data collection.
For each ion solvated in 500 water molecules, we have
performed our molecular dynamics simulations using ͑1͒ the
simple truncation procedure, ͑2͒ the reaction ﬁeld method,
͑3͒ a switching of the potential ͑which alters the force in the
switching region͒, and ͑4͒ the Ewald summation procedure
to treat long-range interactions. In the case of 216 particles in
the primary box, we have selected to study only the case of
bromide, since as it turned out in this case the four different
methods produced the largest deviations in the large size
system. This system with 216 particles was also studied using
the same four methods to treat long-range forces as the
ones mentioned above for the systems with 501 particles. In
the Ewald scheme, the convergence parameter ␣ had the
value of 0.372 ÅϪ1
in order to conﬁne the real space part of
the sum to the primary simulation box. In the evaluation of
the reciprocal space part, the maximum number of the components
of the k vectors had the value of 7 and k2
Ͻ51.
Reaction ﬁeld simulations were carried out with ⑀rfϭϱ and
the cutoff distances Rc were 12 and 9.2 Å in the simulations
involving 500 and 215 water molecules, respectively. The
same cutoff values were used in the simulations with the
switching as well as the ones with the truncation. When the
switching method was used, RL of Eq. ͑2͒ had a value of
0.95 Rc yielding a switching region of 0.6 Å for the cutoff
radius of 12 Å.
IV. RESULTS AND DISCUSSION
The main goal of this work is to determine how the
treatment of long-range forces inﬂuences the dynamics of
ions in aqueous solutions, more speciﬁcally, the limiting
ionic mobility. With this in mind, we have calculated the
velocity autocorrelation functions of the ions. The results of
the calculations are presented in Figs. 1͑a͒, 1͑b͒, and 1͑c͒
where the velocity autocorrelation functions for ﬂuoride,
chloride, and bromide ions are shown. In each case, the correlation
function was calculated from MD simulations with
500 water molecules, and using four different ways to treat
long-range forces. The results from the present simulations
for chloride, using the switching function, and for bromide,
using both switching function and the Ewald sum method,
are in agreement with the previously published results.13–15
A
general feature common to Figs. 1͑a͒–1͑c͒ is that the velocity
autocorrelation function obtained from the simulation
with the switching is displaying more oscillations compared
to the correlation functions obtained from the MD with the
other three methods. The difference is weak for ﬂuoride,
stronger for chloride, and very pronounced for the bromide
ion. Also in cases of chloride and bromide the initial decay
of the velocity autocorrelation function is enhanced when the
switching function method is used in the simulations. The
correlation functions obtained from the simulations with the
Ewald method and with the reaction ﬁeld are, within the
statistical accuracy, practically the same. Surprisingly, the
correlation functions from simulations with truncation look
very similar to the correlation functions obtained from the
simulations with the Ewald sum and with the reaction ﬁeld.
To illustrate the inﬂuence of the cutoff value Rc on the
dynamic properties of the system, we display in Fig. 2 the
velocity autocorrelation functions for bromide ion obtained
from the simulations with 215 water molecules. From Figs. 2
and 1͑c͒ we observe that the largest change in the character
TABLE I. Potential parameters used in the present study.
Water model Aw-w kcal/mol Å12
Cw-wkcal/mol Å6
qH
SPC 629 400.0 625.5 0.41
TIP 4P 600 000.0 610.0 0.52
Ion Aw-I kcal/mol Å12
Cw-I kcal/mol Å6
qI
Fluoride 343 100.0 426.2 Ϫ1.0
Chloride 3 950 000.0 1461.2 Ϫ1.0
Bromide 9 118 300.0 1911.8 Ϫ1.0
452 Perera, Essmann, and Berkowitz: Dynamic of ions in aqueous solutions
J. Chem. Phys., Vol. 102, No. 1, 1 January 1995
of the correlation function occurs in the simulations where
the switching function is used.
The enhanced oscillations and the fast initial decay in
the velocity autocorrelation function obtained from the simulations
with the switching functions are expected to result in
a lower value of the diffusion constant, compared to the values
of this constant obtained from the simulations where the
other three methods were used. Indeed, when we calculated
the values of the diffusion coefﬁcients for the ions using the
mean square displacement formula, that is what we observed.
The mean square displacement curves, calculated for
the whole 120 ps of the simulations are displayed in Fig. 3
for the chloride ion. Similar curves are obtained for the other
ions. The values of the diffusion coefﬁcients obtained from
the ﬁt of mean square displacement curves to straight lines in
the interval of 1–3 ps are presented in Table II͑a͒. The
change in the value of the diffusion coefﬁcient cannot be
attributed to the change in the coordination number of the
ion, when different methods of long-range force treatments
were used. The latter ones nearly do not change, as Table
II͑b͒ shows.
To see the effect of the ion size on the character of the
ionic motion in water, we display the velocity autocorrelation
functions for the ﬂuoride, chloride, and bromide ions in Fig.
4. The correlation functions in Fig. 4 are obtained from the
simulations with the Ewald summation. While the correlation
function for ﬂuoride displays strong oscillations, the oscillations
are reduced for chloride and completely disappear from
bromide. This different character of the ionic motion is rather
simple to explain, if we neglect the effect of the ionic mass
FIG. 1. Velocity autocorrelation functions of ͑a͒ ﬂuoride, ͑b͒ chloride, and
͑c͒ bromide ions from the simulations with 500 water molecules.
FIG. 2. Velocity autocorrelation functions of the bromide ion from the simulations
with 215 water molecules.
FIG. 3. Mean square displacements of the chloride ion from the simulations
with 500 water molecules.
453Perera, Essmann, and Berkowitz: Dynamic of ions in aqueous solutions
J. Chem. Phys., Vol. 102, No. 1, 1 January 1995
on its motion. Due to the small size of the ion the electric
ﬁeld emanating from ﬂuoride is strong. This keeps the water
molecules around the ion for a long residence time, thus
creating a cage around the ion. The ion therefore rattles in
this cage, and the signature of this motion is observed as the
oscillations in the velocity autocorrelation function. The cage
effect is weaker for chloride and disappears for bromide.
Figure 4 also shows that we cannot describe the motion of
ﬂuoride and chloride ions as Brownian, but the trend is
clearly visible.
Since we know now that the simulations with the switching
function result in an enhancement of the oscillatory motion,
the main question is: why does this enhancement happen?
A clue to the answer is obtained from Figs. 5͑a͒, 5͑b͒,
and 6. In Figs. 5͑a͒ and 5͑b͒ we show the plots for the
bromide–water oxygen ͓Fig. 5͑a͔͒ and the bromide–water
hydrogen ͓Fig. 5͑b͔͒ radial distribution functions. ͑In what
follows we show the distribution functions for bromide only,
since we get the same qualitative information from the consideration
of the other ion–water distribution functions.͒ As
we can see from Fig. 5͑a͒ the four methods produce statistically
similar ion–water oxygen distribution functions in general,
but when we use the potential switching method in the
simulation, a spurious peak at around 11 Å ͑0.4 Å away from
RL͒ can be observed. Correspondingly, we observe the enhancement
of hydrogen density at around 10.3 Å in Fig. 5͑b͒.
In Fig. 6 we display the average cosine of the angle between
the water dipole and the vector connecting the water oxygen
with the ion center. As we can see from Fig. 6 a spurious
peak is also observed for this distribution at the distance just
before the switching function is turned on. Together, Figs. 5
and 6 indicate that the effect of the introduction of a switching
function into the ion–water interaction potential is to
create a buildup of water just before the distance RL from the
ion, where the switching of the potential begins. Also, at this
distance the water molecules in the simulations with the
switching function display a preferential orientation of their
dipoles towards the ion, which is easily inferred from the
buildup of hydrogen and directly from Fig. 6. Thus we see,
that when we use a switching function in our simulations, we
create a shell of water around the ion close to the distance
RL . This shell outlines a large sphere of water around the
TABLE II. ͑a͒ Diffusion coefﬁcients ͑in 10Ϫ5
cm2
/s͒. ͑b͒ Coordination num-
bers.
Ion Switching
Reaction
ﬁeld
Ewald
sum Truncation Experimenta
͑a͒
Fluoride 1.56 2.06 1.72 1.80 1.46
Chloride 0.78 1.57 1.57 1.91 2.03
Bromide 1.37 2.43 2.34 2.50 2.08
͑b͒
Fluoride 6.10 6.16 6.20 6.12
Chloride 7.28 7.42 7.34 7.36
Bromide 7.87 7.68 7.60 7.76
a
Reference 27.
FIG. 4. Velocity autocorrelation functions of ﬂuoride, chloride, and bromide
ions from the simulations with 500 water molecules using the Ewald summation
technique.
FIG. 5. ͑a͒ The bromide–oxygen radial distribution functions from simulations
with 500 water molecules. The insert is the enlargement of the last 3 Å
segment. ͑b͒ The hydrogen–Bronide radial distribution functions obtained
from the same set of simulations.
454 Perera, Essmann, and Berkowitz: Dynamic of ions in aqueous solutions
J. Chem. Phys., Vol. 102, No. 1, 1 January 1995
ion, and the ion is performing a rattlelike motion in this large
sphere, which is why we observe enhancement of the oscillations
in the velocity autocorrelation function.
Why then do we observe the enhancement of the oxygen
and particularly hydrogen density in the region where the
action of the switching function starts? The answer can be
obtained from the consideration of Fig. 7. Figure 7 displays
how the truncation, reaction ﬁeld, and switching procedures
modify the Coulomb interaction between two similar
charges. As we can see from Fig. 7, the reaction ﬁeld modiﬁes
the interaction in a smooth manner and turns it off at
RϭRc ; the truncation procedure turns off the interaction
abruptly at RϭRc . As we observe from Fig. 7, a large force
is artiﬁcially created in the region of the switch, which is
totally spurious. This large force can be either attractive or
repulsive ͑depending on what atom it acts͒, but in any case it
leads to a depletion of water in the switching region, which
is responsible for the artifacts of the simulation.
Why then is the switching function method often used in
the simulations? It seems that it was introduced into MD to
avoid the warm up of the system, due to the energy discontinuity
when the simple truncation method is employed
alone. Although the switching function helps in stabilizing
the value of the system’s temperature, it can create more
serious problems. Our results show that for the study of limiting
ionic mobility the usage of the Ewald summation technique
and of the reaction ﬁeld technique ͑which is equivalent
to the usage of the shifting of the force technique, proposed
recently by Steinbach and Brooks͒17
produced basically the
same results. Previous simulation by Prevost et al.18
also
demonstrated that the application of the shifting function ͓the
one described in Eq. ͑3͒ with nϭ1͔ to a Coulomb interaction
potential produced the same results as the Ewald summation
for pure water. The simple truncation method when applied
to simulations with just one ion produced results which are
very similar to the results from the simulations where Ewald
and reaction ﬁeld methods were used. However we do not
recommend use of this simple truncation technique in simulations
with more than one ion, since this cutoff scheme produces
artiﬁcial correlations between ions.26
Our present study is in no way a comprehensive study of
how the treatment of long-range forces inﬂuences the structure,
thermodynamics, and dynamics of aqueous solutions.
But on the basis of it we would recommend use of the Ewald
summation or the reaction ﬁeld techniques for the study of
aqueous ionic solutions in computer simulations. At least
these two methods produce consistent results and can be
given a simple physical interpretation.
With respect to the questions asked in Sec. I and related
to the change in the behavior of the velocity autocorrelation
function and the value of the diffusion coefﬁcient of the
chloride ion, we observed that: ͑a͒ the velocity autocorrelation
function of chloride is still oscillatory, even when Ewald
summation or reaction ﬁeld techniques are used to treat longrange
interactions, ͑b͒ the value of the diffusion coefﬁcient
for chloride calculated using the Ewald summation technique
is very close to the value obtained from the simulation using
the reaction ﬁeld method. Although these values are still below
the experimental one by about 20%, considering the uncertainty
in the calculated value ͑estimated to be around
10%–15%͒ the agreement between the calculation and experiment
can be now considered as reasonable.
FIG. 6. Distributions of the cosine of the angle between the line joining the
bromide ion and the water oxygen and the dipole vector as a function of
distance from the ion.
FIG. 7. The force between two charges as a function of R*(ϭR/Rc under
different cutoff schemes. The switching of the force begins at 0.95 Rc .
455Perera, Essmann, and Berkowitz: Dynamic of ions in aqueous solutions
J. Chem. Phys., Vol. 102, No. 1, 1 January 1995
ACKNOWLEDGMENTS
This work was supported by grants from the Ofﬁce of
Naval Research. Part of the simulations were performed on
the Cray YMP at the North Carolina Supercomputer Center.
We are grateful to Professor P. Rossky and Dr. B. Brooks for
making preprints of Refs. 13 and 17 available to us prior to
their publication.
1
R. J. Loncharich and B. R. Brooks, Proteins: Structure, Function Genet. 6,
32 ͑1989͒.
2
P. E. Smith and B. M. Pettitt, J. Chem. Phys. 95, 8430 ͑1991͒.
3
G. Hummer, D. K. Soumpasis, and M. Neumann, Mol. Phys. 77, 769
͑1992͒; 81, 1155 ͑1993͒.
4
H. Alper and R. Levy, J. Chem. Phys. 91, 1242 ͑1989͒.
5
D. M. York, T. A. Darden, and L. G. Pedersen, J. Chem. Phys. 99, 8345
͑1993͒.
6
H. Schreiber and O. Steinhauser, Chem. Phys. 168, 75 ͑1992͒; J. Mol.
Biol. 228, 909 ͑1992͒; Biochemistry 31, 5856 ͑1992͒.
7
L. Greengard and V. Rokhlin, J. Comp. Phys. 73, 325 ͑1987͒.
8
C. L. Brooks III, B. M. Pettitt, and M. Karplus, J. Chem. Phys. 83, 5897
͑1985͒.
9
J. A. Barker and R. O. Watts, Mol. Phys. 26, 789 ͑1973͒.
10
S. W. de Leeuw, J. W. Perram, and E. R. Smith, Proc. R. Soc. London, Ser.
A 373, 27 ͑1980͒.
11
L. X. Dang and B. M. Pettitt, J. Am. Chem. Soc. 109, 5531 ͑1987͒; J.
Phys. Chem. 94, 4303 ͑1990͒.
12
J. S. Bader and D. Chandler, J. Phys. Chem. 96, 6423 ͑1992͒.
13
G. S. Del Buono, T. S. Cohen, and P. J. Rossky, J. Mol. Liq. 60, 221
͑1994͒.
14
M. R. Reddy and M. Berkowitz, J. Chem. Phys. 88, 7104 ͑1988͒.
15
M. Berkowitz and W. Wan, J. Chem. Phys. 86, 376 ͑1987͒.
16
O. Steinhauser, Mol. Phys. 45, 335 ͑1982͒.
17
P. J. Steinbach and B. R. Brooks, J. Comp. Chem. 15, 667 ͑1994͒.
18
M. Prevost, D. van Belle, G. Lippens, and S. Wodak, Mol. Phys. 71, 587
͑1990͒.
19
M. P. Allen and D. J. Tildesley, Computer Simulation of Liquids ͑Oxford
University Press, Oxford, 1987͒.
20
H. J. C. Berendsen, J. P. M. Postma, W. F. van Gunsteren, and J. Hermans,
in Intramolecular Forces, edited by B. Pullman ͑Reidel, Dordrecht, 1981͒.
21
W. L. Jorgensen, J. W. Madura, R. W. Impey, and M. L. Klein, J. Chem.
Phys. 79, 926 ͑1983͒.
22
T. P. Staatsma and H. J. C. Berendsen, J. Chem. Phys. 89, 5876 ͑1988͒.
23
J. Chandrasekhar, D. C. Spellmeyer, and W. L. Jorgensen, J. Am. Chem.
Soc. 106, 903 ͑1984͒.
24
G. Palinkas, W. O. Riede, and K. Heinzinger, Z. Naturforsch. Teil A 32,
1137 ͑1977͒.
25
H. J. C. Berendsen, J. P. M. Postma, W. F. van Gunsteren, A. DiNola, and
J. R. Haak, J. Chem. Phys. 81, 3684 ͑1984͒.
26
W. F. van Gunsteren ͑private communication͒.
27
P. W. Atkins, Physical Chemistry ͑Freeman, New York, 1986͒.
456 Perera, Essmann, and Berkowitz: Dynamic of ions in aqueous solutions
J. Chem. Phys., Vol. 102, No. 1, 1 January 1995