J. Phys. Chem. 1996, 100, 2723-2730
2723
Soft Sticky Dipole Potential for Liquid Water: A New Model Yi Liu and Toshiko Ichiye* Departments of Biochemistry and Biophysics and of Chemistry, Washington State UniVersity, Pullman, Washington 99164-4660 ReceiVed: August 8, 1995; In Final Form: December 4, 1995X
A new, efficient potential model for liquid water is presented here. It is based on the hard-sphere sticky dipole potential model for water by Bratko, Blum, and Luzar (J. Chem. Phys. 1985, 83, 6367), referred to as the BBL model. Similar to the BBL model, this new, soft-sphere sticky dipole model has a single interaction site at the molecular center of mass with a spherical repulsive potential, a short-range tetrahedral “sticky” potential, and a point dipolar potential. However, the use of a Lennard-Jones-type soft-sphere, as opposed to the hard-sphere in the BBL model, allows realistic studies of water and aqueous solvation. This is particularly important for the existing parametrizations of biological molecules that use soft-sphere models. The present model gives a liquid water structure comparable to that found by the four-site TIP4P model and also gives an intermolecular energy, a hydrogen bond energy, and a heat capacity of liquid water in good agreement with experimental data and/or results from the TIP3P or the TIP4P model. In addition, Monte Carlo simulations using this model are nearly an order of magnitude faster than those using the TIP3P or TIP4P model. Monte Carlo simulations have also been carried out to study the solvation of a single Na+ or Cl- ion at room temperature, using a hybrid scheme in which the ion-water interaction is modeled by monopole-dipole plus monopole-quadrupole potentials and the water-water interaction is modeled by this new model. The calculated structure of water around the ion and the enthalpy of ionic solvation are in good agreement with those from experiments and from simulations using other water models. Overall, the simplicity, efficiency, and reasonable accuracy of this model make it potentially very useful for studies of aqueous solvation by either computer simulations or integral equation theories.
Introduction Accurate descriptions of liquid water and aqueous solutions are crucial in our understanding of many chemical and biochemical problems. Hence, modeling the behavior of water has been the subject of extensive research. To describe the molecular structure of water, the most widely used effective pair potential models use rigid multiple interaction sites with partial charges, such as the three-site TIP3P,1 SPC,2 and SPC/ E,3 the four-site TIP4P,4 and the five-site ST25 models, whose merits have been well documented.4,6 In particular, the threesite TIP3P, SPC, and SPC/E models are commonly used in biochemical simulations, in part because of their reasonable descriptions of solvation and dielectric properties (dielectric constant 0 ≈ 70).7-10 One apparent drawback of these threesite models is the lack of structure beyond the first coordination shell. On the other hand, the four-site TIP4P model provides structure in excellent agreement with experiment but exhibits a dielectric constant (0 ) 53)11 considerably lower than the experimental value (0 ) 78),12 which is one reason that TIP4P is not as widely used for simulations of solvated biological molecules as the three-site models such as TIP3P, SPC, or SPC/ E. The five-site ST2 model is computationally more expensive than other models cited above. In addition, water models based on these multiple-site potentials with added atomic polarizability have been reported in the literature, such as the recently proposed model by Rick et al.13 Another variation is the various central force models (CFM’s) of water,14-17 in which liquid water is modeled as a 1:2 mixture of oxygen and hydrogen. However, the application of the CFM models to solvation of complex solutes has been very limited, in part because of the * Author to whom correspondence should be addressed. X Abstract published in AdVance ACS Abstracts, January 15, 1996.
0022-3654/96/20100-2723$12.00/0
small time step needed in simulations of this inherently flexible water model. Modeling liquid water explicitly by these potentials in computer simulations is a very powerful technique but has inherent limitations. For example, a major stumbling block for computer simulations of biological macromolecules is that an accurate description of the solvation of biological macromolecules usually requires use of a very large number of water molecules, on the order of 103-106, which in most cases constitutes the majority of the computation and is sometimes forbiddingly expensive. One approach to improving this situation is to devise new water models that are more computationally efficient. For instance, the number of interaction sites can be reduced, since the computational efficiency of a model is approximately proportional to 1/n2, n being the number of sites used for the water monomer of the model. Along these lines are the so-called Stockmayer-type potentials,18-20 consisting of a spherical soft repulsion with an embedded point dipole. The Stockmayer potentials are single-site with an orientation associated with that site. The biggest problem associated with this type of model is the lack of hydrogen bonding, which is a very important determinant of the structure of water. Adding a quadrupole moment substantially slows down the computation, because one must then evaluate the dipole-dipole, dipolequadrupole, and quadrupole-quadrupole interactions. Another one-site model of water is the hard-sphere sticky dipole potential by Bratko, Blum, and Luzar (BBL),21 which uses a tetrahedral octapolar “sticky” potential to describe the hydrogen bonds in the first coordination shell, in addition to the embedded point dipole. The new revised version of the BBL model has the same dipole and octapole terms but includes an exponential term beyond the hard-sphere radius.22 Compared to the dipole plus quadrupole approach, the BBL potential is more efficient © 1996 American Chemical Society
2724 J. Phys. Chem., Vol. 100, No. 7, 1996 because the “sticky” term does not interact with the point dipole explicitly. The BBL model was originally developed for purely analytical calculations23 and has not been used in macromolecular solvation studies because it is a hard-sphere model. Another approach to improving the computational speed is to model the solvent implicitly, rather than calculating the positions of all of the waters explicitly. A common approximation is to use a structureless dielectric continuum24,25 to mimic the dielectric nature of liquid water, which works well for many problems but cannot describe noncontinuum behavior. A different method is to use statistical mechanical integral equation theories.26,27 For the multisite potentials, the reference interaction site model (RISM) theories28-31 are the only reasonably successful site-based integral equation theories. However, many problems are very difficult to tackle because of the large number of equations involved in the formulation of the RISM-type theory for large solutes. Single-site integral equations can be used with the CFM to model pure water32,33 and more recently to model water near planar interfaces.34 Finally, since the BBL model is a reasonable one-site model capable of describing the hydrogen-bonding structure, it naturally becomes a candidate for use in integral equation theories. We have recently developed such a theory that uses the BBL model to describe the three-dimensional structure of water around globular complex solutes,35 albeit there are drawbacks (e.g. hard-sphere) of the BBL potential. From the above, a modified soft-sphere version of the BBL model is potentially very useful because, compared to multisite models, it would be faster in computer simulations and easier to incorporate into integral equation theories. In this paper, we describe such a model, the soft sticky dipole (SSD) potential model of water. It is being developed, although not exclusively, for applications to aqueous solvation of biological systems. The modifications of the BBL model are as follows. First, the hardsphere of the BBL is replaced by a Lennard-Jones (L-J) type potential. Second, since the dielectric properties are of paramount importance for biological systems, the permanent dipole moment is increased from the gas phase value 1.84 D used in the BBL model to 2.35 D used in the TIP3P and SPC/E models, because the larger dipole moment has been attributed as the factor determining their good dielectric properties.7-10 Finally, the sticky potential is similar to that of the original BBL except that a modulating function is used for compatibility with the soft-sphere. Also, an additional potential term is added to the sticky potential to correct the spatial distribution of water in the first shell and the location of the second peak in the O-O radial distribution function that is severely misplaced in the BBL model. In addition, an ion-water potential consisting of a monopole-dipole plus a monopole-quadrupole term is used to describe the interaction between an ion and a SSD water molecule. Here, Monte Carlo (MC) simulations are used to calculate the structure, intermolecular energy, hydrogen bond energy, and heat capacity of liquid water. Also reported are MC results of dilute ionic solutions in which a Na+ or a Cl- is solvated by liquid water at room temperature. The structure of water near the ions and the enthalpy of ionic solvation are calculated and compared to results from experiments. SSD Water Model A. Pure Water Potential. The water molecule of the SSD model is treated as a L-J sphere with an embedded point dipole plus a tetrahedral “sticky” potential, all situated at the molecular center of mass (M) located on the H-O-H bisector, 0.0654 Å from the oxygen toward the hydrogens. The geometry of the
Liu and Ichiye TABLE 1: Parameters for the Ion-Water Potentiala particle type σ (Å) (kcal/mol) µ (D) Qxx (b)b Qyy (b)b Qzz (b)b water (M) Na+ Cl-
3.051 2.579 4.445
0.152 0.118 0.100
2.35
-1.682
1.762
-0.08
a See text for definitions. b Calculated for the TIP3P monomer. Quadrupole moment unit b ) 10-26 esu cm2.
water molecule used in the SSD model is same as that used for TIP3P; the O-H length is 0.9572 Å and the H-O-H angle is 104.52°. The H-M-H angle is 109.47°, the tetrahedral angle. The center of mass is the only interaction site of this model. The explicit coordinates of the three nuclei of a water molecule are not used in the calculation of the intermolecular potential energy but can be indirectly determined from the coordinates of the molecular center and the orientation of the dipole moment and of the molecular plane. The total interaction potential energy between two water molecules i and j is therefore sp νij ) νijLJ(rij) + νijdp(rij,Ωi,Ωj) + νip (rij,Ωi,Ωj)
(2.1)
where rij is the distance between the molecular centers, rij is the separation vector between the two molecular centers, and Ω is the orientation of a water molecule that is determined by both the orientation of the dipole moment vector and the orientation of the molecular plane, all in a coordinate system with an arbitrary origin and an orientation described below. The first term in eq 2.1 is νijLJ(rij), the L-J potential
[( ) ( ) ]
νijLJ(rij) ) 4w
σw rij
12
-
σw rij
6
(2.2)
The parameters w and σw used for the SSD are given in Table 1. The distance parameter, σw, is slightly different from that used for oxygen-oxygen in TIP3P (wherein σw ) 3.1506 Å) because the molecular center, not the oxygen, is used as the interaction site. The second term in eq 2.1 is νijdp(rij,Ωi,Ωj), the point dipole-point dipole potential,
νijdp(rij,Ωi,Ωj) )
µi‚µj rij3
-
3(µi‚rij)(µj‚rij) rij5
(2.3)
where µi and µj are the dipole moment vectors whose magnitude is 2.35 D (like that of a TIP3P or SPC/E monomer). Please note that a factor, 1/4π0, has been (and henceforth will be) omitted for simplicity. The last term in eq 2.1 is νijsp(rij,Ωi,Ωj), the tetrahedral sticky potential,
νijsp(rij,Ωi,Ωj) ) ν°[s(rij)wij(rij,Ωi,Ωj) + s′(rij)wijx(θij)] (2.4) where ν° ) 3.7284 kcal/mol determines the strength of the sticky potential. The tetrahedral coordination in the first shell is satisfied by the function wij(rij,Ωi,Ωj),
wij(rij,Ωi,Ωj) ) sin θij sin 2θij cos 2φij + sin θji sin 2θji cos 2φji (2.5) in which (θij,φij) is the set of spherical polar angles of the position of molecule j in the frame fixed on molecule i and with an orientation such that the Z-axis is parallel to the dipole moment of i and the X-axis is perpendicular to the molecular plane of i. The reader is referred to ref 21 for details. In addition, wijx(θij) in eq 2.4 is an empirical correction term
Soft Sticky Dipole Potential for Liquid Water
J. Phys. Chem., Vol. 100, No. 7, 1996 2725
added to the sticky potential of the original BBL model
wijx(rij,Ωi,Ωj) ) (cos θij - 0.6)2(cos θij + 0.8)2 + (cos θji - 0.6)2(cos θji + 0.8)2 - 2w° (2.6) where w° ) 0.07715. This additional function is necessary because wij(rij,Ωi,Ωj) vanishes at angles θij ) 0° and θij ) 180°, which leads to a 50% probability that a nearest neighbor lies along the dipole axis of the reference monomer. This is an incorrect configuration for the first shell. Indeed, simulation studies of the three-dimensional spatial structure of water have shown that the probability of this aligned structure in the first shell is very small,36,37 especially for the near hydrogen neighbors. Finally, s(rij) and s′(rij) are the modulating functions which interpolate smoothly between 0 and 1 and are given by
s(rij) ) 1 s(rij) )
(rU - rij)2(rU + 2rij - 3rL) (rU - rL)3
s(rij) ) 0
if rij < rL
third term, νmq(r,Ωw), in eq 2.8 is the monopole-quadrupole interaction,
if rL e rij < rU (2.7)
νmq(r,Ωw) ) - q/6r3[Qxx(3 sin2 θ cos2 φ - 1) + Qyy(3 sin2 θ sin2 φ - 1) + Qzz(3 cos2 θ - 1)] (2.12)
if rij g rU
where rL ) 2.75 Å and either rU ) 3.35 Å for s(rij) or rU ) 4.0 Å for s′(rij). The parameters used here are carefully chosen so that the potential, as will be shown below, leads to reasonable structural and energetic results for liquid water. B. Ion-Water Potential. To model solvation of charged solutes, an ion-water potential must be used in supplement to the SSD. This is because the sticky potential used in SSD is formulated to describe the hydrogen bonds in the water-water interaction only, and the point dipole potential cannot mimic hydrogen bonds between an ion (particularly an anion) and the water molecules immediately surrounding it. Also, the ionwater potential must be transferable to all kinds of ions, which essentially eliminates the possibility of using any empirical ionwater sticky potential. To maintain the generality and the consistency of a one-site model, we use a hybrid method in which the ion-water interaction is modeled by ion-dipole plus ion-quadrupole interaction, wherein the dipole and quadrupole are obtained from a center-of-mass point quadrupolar expansion of the charge distribution of the TIP3P monomer. Therefore, the ion-water interaction is given by
νiw(r,Ω) ) νLJ(r) + νmd(r,Ωw) + νmq(r,Ωw)
(2.8)
in which the first term, νLJ(r), is the L-J potential (eq 2.2) between the ion and the center of mass of the water molecule. The indices i and w denote ion and water, respectively. The ion-water cross parameters for energy and distance are determined as follows:
iw ) xiw
(2.9)
and
σiw )
σi + σw 2
(2.10)
The second term, νmd(r,Ωw), in eq 2.8 is the monopole-dipole interaction,
νmd(r,Ωw) )
Figure 1. Schematic representation of the coordinate system used to describe the relative position and orientation of a water molecule with respect to an ion. The origin is located at the center of mass of the water molecule.
qµ cos θ r2
(2.11)
where q is the charge of the ion, µ is the dipole moment of the water molecule, and θ is the angle shown in Figure 1. The
where θ and φ are the angles shown in Figure 1 and the three quadrupole moments, Qxx, Qyy, and Qzz, are calculated by the center-of-mass point quadrupolar expansion of the charge distribution of the TIP3P monomer given the orientation shown in Figure 1. The values of the dipole and quadrupole moments used in this work are given in Table 1. C. Monte Carlo Simulations. We have calculated various structural and energetic properties of the SSD model and two other multisite models, TIP3P1 and TIP4P.4 The simulations here were carried out in the NVT ensemble at the temperature 298 K. A cubic box of N ) 125 water molecules, the experimental density of water (0.03346 molecules/Å3), periodic boundary conditions, and standard Metropolis sampling40 were used. A spherical truncation of the interaction potential was applied at 0.5L, where L ) 15.5165 Å (the edge length of the cube). System size and potential truncation have been found previously to have small effects on the quantities calculated here.4,41,42 The starting configurations were taken from previous simulations and were thoroughly equilibrated again before additional MC runs were used to generate structural and energetic data. We have calculated the average intermolecular energy, the energy pair distribution, the O-O, O-H, and H-H radial distribution functions, and the constant-volume heat capacity, Cv. In the terminology used by Pangali et al.,42 an MC “pass” is defined by N attempted moves, each involving both a translation and a rotation. The number of MC passes used to generate the radial distribution functions, average intermolecular energy, and energy pair distributions was 3200 because these quantities converge rapidly in MC simulations. However, Cv is a fluctuation property and is known to be slowly convergent.42-44 Hence, 120 000 MC passes were used to calculate Cv. These runs were sufficient for the system studied here according to results reported elsewhere.43-46 The acceptance ratio in all MC runs was 40%. To test the hybrid scheme for the ion-water potential, we have carried out MC simulations of two dilute ionic systems, each containing a single Na+ or Cl- fixed at the center of a cubic box containing 215 water molecules with an edge length L ) 18.62 Å. The L-J potential parameters used here for the ions (Table 1) are chosen so that the ion-water dimer potential reproduces the dimer equilibrium distance and the potential well depth as determined by Hartree-Fock calculations47 and the experimental gas phase dimer energy,48 respectively. The choice of these parameters is not critical for our purposes, since the
2726 J. Phys. Chem., Vol. 100, No. 7, 1996
Liu and Ichiye
ionic solvation energies are dominated by the electrostatic interaction. Nevertheless, the Na+ parameters used here are quite close to those used by Pettitt and Rossky,49 and the Clparameters are very close to those used by Dang et al.50 Equations 2.9 and 2.10 were used to determine the ion-water cross potential parameters. The simulation temperature, the density of water, the potential cutoff, and the sampling procedures were the same as above. The effect of potential truncation was not corrected in this study, as in many reported studies of similar systems.50-57 This is because there is no generally accepted way of treating the cutoff effect for a netcharged system without altering the nature of the system. After equilibration, 40 000 MC passes were used for each case to calculate the ion-O and ion-H radial distribution functions and the enthalpy of solvation. The enthalpy of solvation, ∆Hsol, was calculated here as the net-attractive ion-water interaction energy, Eiw, plus the net-repulsive solvent reorganization energy, ∆Eww, i.e.,
∆Hsol ) Eiw + ∆Eww
(2.13)
∆Eww ) Ew - E°w
(2.14)
Figure 2. Comparison of the O-O radial distribution functions for the SSD (solid line) with TIP3P (dotted line), TIP4P (dashed line), and experimental data58 (dotted dash).
where
and Ew and E°w are the energies of solvent in solution and of pure solvent, respectively. We have also examined the influence of the system size on the ionic solution results by carrying out simulations using either ion in 400 water molecules (L ) 22.88 Å). The increase in system size resulted in only small structural changes at distances beyond the first solvation shell. In addition, the system size was found to have little effect on the calculated enthalpy of solvation, which decreased by approximately 5% by going from 215 to 400 water molecules, due almost entirely to the reduced solvent reorganization energy, ∆Eww.
Figure 3. Comparison of the O-H radial distribution functions for the SSD (solid line) with TIP3P (dotted line), TIP4P (dashed line), and experimental data58 (dotted dash).
Results and Discussion In this section, we present results comparing the SSD potential with those of two commonly used water potential functions, namely TIP3P and TIP4P. TIP3P is commonly used in biological studies because it is a simple three-site model for which extensive thermodynamic studies of solvation have been done.1 In fact, TIP3P is used as the default water in widely used molecular modeling software packages such as CHARMM/ QUANTA38 and Discover/Insight.39 Moreover, TIP3P has very similar parameters to SPC and gives an essentially identical structure.4 TIP4P is a four-site model that gives excellent water structure and is very widely used in simulation studies. We have compared the O-O, O-H, and H-H radial distribution functions, the average intermolecular energy, Ei, the energy pair distribution, n(E), and the constant volume heat capacity, Cv. The results are also compared to experimental data and results for other potential models where available. We also present the structural and energetic results of the dilute Na+ and Claqueous solutions. A. Radial Distribution Functions of Pure Water. The O-O (Figure 2), O-H (Figure 3), and H-H (Figure 4) radial distribution functions of SSD, TIP3P, and TIP4P from MC simulations, along with the experimental neutron diffraction results by Soper and Phillips,58 are compared. The overall structure given by SSD is reasonably good and is much improved from that of the original BBL model. As shown in Figure 2, the O-O distribution displays two well-defined peaks, similar to the results from experiment and from simulation using the four-site TIP4P results. A comparison of the coordination
Figure 4. Comparison of the H-H radial distribution functions for the SSD (solid line) with TIP3P (dotted line), TIP4P (dashed line), and experimental data58 (dotted dash).
numbers calculated by integrating the first peak of the O-O curve to 3.3 and 3.5 Å (Table 2) shows good agreement with other models and experiment. The O-H radial distribution function (Figure 3) is also similar to the experimental results and that from simulations using the TIP3P and TIP4P models. The average number of hydrogen-bonded neighbors per monomer calculated by integrating the O-H curve to 2.35 and 2.5 Å (Table 2) shows good agreement with other models. The H-H distribution (Figure 4) is more different from the experimental and the TIP3P and TIP4P results, especially in the third peak, although it should be noted that the precise location and height of the peaks of the H-H curve from neutron diffraction are uncertain, as frequently noted in the literature.4,58-62 The second peak of the O-O radial distribution function bears further discussion, since a peak near 4.5 Å is indicative of tetrahedral ordering. The O-O function for the SSD shown in
Soft Sticky Dipole Potential for Liquid Water
J. Phys. Chem., Vol. 100, No. 7, 1996 2727
TABLE 2: Comparison of Structural and Energetic Results for SSD with TIP3P, TIP4P, SPC, and Experiments propertya numberb
coordination Ei (kcal/mol) E* Hbd (kcal/mol) E h Hbd (kcal/mol) H bonds per monomere H bonds per monomerf Cv (cal/(mol K))
SSD
TIP3P
TIP4P
4.4 (5.2) -9.60 -2.7 -4.84 3.6 (4.0) 4.9 20.0
4.5 (5.5) -9.67 -3.2 -4.62 3.5 (3.9) 4.1 15.6f
4.2 (5.0) -10.02 -2.7 -4.3 3.5 (3.8) 5.1 16.0g
SPC
experiment
(5.1) 4.5 (5.3)c -10.18 -9.92d
(3.9)
3.6 (4.0)c
21.7g
17.9h
a See text for definitions. Unless otherwise specified, data for SSD, TIP3P, and TIP4P are results of this work. Data for SPC are from ref 4. b Calculated by integrating O-O radial distributions to 3.3 Å (or 3.5 Å for numbers in parentheses). c Calculated from data reported in ref 58. d From ref 1. e Calculated by integrating the O-H radial distribution function to 2.35 Å (or 2.5 Å for numbers in parentheses). f Calculated by integrating the energy pair distribution to E* . Hbd g Recalculated from data in ref 4. See text concerning accuracy. h From ref 63, p 172.
Figure 2 exhibits a small, less obvious peak at 3.9 Å, inside the more pronounced peak at 5.0 Å. Apparently, the overall second peak is due to the tetrahedral structure, while the split into the inner and outer peaks seen here for SSD is due to the complex combination of the sticky and dipole interactions. This may be partially corrected by increasing the strength of the sticky potential (ν° in eq 2.4) and/or the dipole moment. However, stronger sticky and dipolar interactions would result in the first O-O peak being too high, sacrificing the more important primary structure. Nevertheless, the structure of the SSD is very much improved from that of the original BBL21 and the newer version of the BBL structure22 in that both of these models overestimate the first O-O and O-H peaks due to the hardcore potential. In addition, another major difference is that while the newer BBL model, and to a lesser extent the original BBL one, exhibits some sort of tetrahedral structure, i.e., small peaks near 3.9 Å and 5.0 Å, the height of these peaks is actually less than 1 and so the first major second peak is at almost 6 Å. This may indicate that the simple packing of nonpolar fluid predominates, or it may be simply a consequence of the hard-sphere potential, as the authors indicate in ref 22. Moreover, the threesite TIP3P model displays no structure beyond the first peak, which has also been shown for the SPC model.4 Overall the structure of SSD water appears comparable to experiment and TIP4P and more structured than TIP3P, which is much less structured than experiment. These results are quite encouraging. B. Average Intermolecular Energy and Energy Pair Distribution of Pure Water. The average intermolecular energy, Ei, is the total potential energy averaged over all intermolecular pairs in the simulation box. From the experimental value of the heat of vaporization at 298 K and 1 atm, the “experimental” value of Ei, -9.92 kcal/mol, can be obtained.4 The results from our calculations (Table 2) are comparable to those obtained both from the other models studied here and from experiment. The SSD gives a much improved intermolecular energy as compared to the original BBL model, which has a too high Ei at about -5.2 kcal/mol. Energy pair distribution (Figure 5), n(E), sometimes referred to as the pair interaction distribution, is the distribution of pair interaction energies experienced by a monomer. This distribution provides useful energetic information that helps define hydrogen bond energies and has been calculated for many water potential models,4,5 although it cannot be obtained experimentally. Figure 5 shows a portion of this distribution calculated from MC simulations for SSD, TIP3P, and TIP4P in such a way that the number of pairs interacting with potential energy
Figure 5. Comparison of the energy pair distribution, n(E), for the SSD (solid line) with TIP3P (dotted line) and TIP4P data (dashed line).
E, n(E), is plotted as a function of E. This plot can also be viewed as a pair potential energy spectrum, which generally shows two regions, a hydrogen-bonding region at largely negative pair energies and a bulk region at near-zero energies, separated by a minimum at E ) E* Hbd, the threshold of hydrogen bond such that a pair of molecules are considered to be hydrogen-bonded if the pair interaction energy is below E*Hbd. The very large population at E ) 0 is simply due to weakly interacting distant molecules (the bulk water). As shown in Figure 5, the shape of the n(E) for the SSD model is somewhat different from that for TIP3P or TIP4P, but the overall trend is quite similar. The wider distribution for SSD compared to those for TIP3P and TIP4P implies slightly stronger hydrogen bonds that are determined largely by the strength of the tetrahedral sticky potential, although the L-J and point dipole potential terms also contribute. The average hydrogen bond strength, E h Hbd, calculated by
∫-∞E * n(E)E dE E h Hbd ) E * ∫-∞ n(E) dE Hbd
Hbd
(4.1)
shows that the hydrogen bonds of SSD are indeed somewhat stronger than that of TIP3P or TIP4P (Table 2). However, reducing the sticky potential strength degrades the structural results. Therefore, modifying the parameters to yield a hydrogen bond strength closer to that of TIP3P or TIP4P does not seem to be warranted. Nevertheless, the difference is only about 0.5 kcal/mol, and similar hydrogen bond behavior has been observed E H*bd n(E) dE gives the for the ST2 potential.5 Finally, ∫-∞ average number of hydrogen-bonded neighbors (Table 2), which is somewhat larger compared to that calculated by integrating the first peak of the O-H radial distribution function for all potentials. C. Heat Capacity of Pure Water. The constant volume heat capacity, Cv, can be calculated from the fluctuations in potential energy from MC simulations by
Cv ) C°v +
〈E2〉 - 〈E〉2 kBT2
(4.2)
where 〈E2〉 - 〈E〉2 is the potential energy fluctuation, kB is the Boltzmann constant, and T is the absolute temperature. C°v ) 1.5NkB is the ideal gas heat capacity contribution. The Cv for SSD calculated from 120 000 MC passes is compared to the experimental value in Table 2. SSD slightly overestimates Cv (20.0 cal/(mol K), as compared to the experimental value, 17.9 cal/(mol K)), but is comparable to the very popular SPC model, which gives Cv ) 21.7 cal/(mol K), and is much better than
2728 J. Phys. Chem., Vol. 100, No. 7, 1996
Liu and Ichiye TABLE 3: Structural and Energetic Results for Aqueous Solvation of Na+ and Cl- Ions propertya
this work other reported workb
experiment
+
ion-O distance (Å) ion-H distance (Å) coordination number Eiw (kcal/mol) ∆Eww (kcal/mol) ∆Hsol (kcal/mol)
Water around a Na Ion 2.30 2.30-2.35 2.9 2.8-3.0 6.3d 5.4-6.0 -221 90 -131 -126f
2.4c 4-6c -106g
-
Figure 6. Radial distribution functions for Na+-O (solid line) and Na+-H (dashed line) calculated from Monte Carlo simulations.
ion-O distance (Å) ion-H distance (Å) coordination number Eiw (kcal/mol) ∆Eww (kcal/mol) ∆Hsol (kcal/mol)
Water around a Cl Ion 3.4 2.7-3.5 2.5 1.7-2.55 9.9e 5.6-8.36 -140 50 -90 -80f
3.10-3.35c 2.22-2.26c 5-11c -82g
a See text for definitions. b Data from simulations using other multisite models reported in the literature, see refs 49-51. c From the X-ray study of Newsome67 and the neutron diffraction study of Neilson and Enderby.68 d Calculated by integrating g(r) to r ) 3.1 Å. e Calculated by integrating g(r) to r ) 4.2 Å. f From simulations using the TIP4P model, see ref 51. g From ref 69.
Figure 7. Radial distribution functions for Cl--O (solid line) and Cl--H (dashed line) calculated from Monte Carlo simulations.
ST2, which gives Cv ) 34.7 cal/(mol K), all reported previously.4-6 For comparison, we have also recalculated Cv for TIP3P and TIP4P (see Table 2) from the results of constantpressure heat capacity Cp reported in a previous work by Jorgensen et al.4 We must point out that these values were calculated by Jorgensen et al. from only 8000 MC passes, 15 times shorter than what was used to calculate Cv for SSD here. In fact, Cv is commonly seen to slowly increase during MC runs until convergence.43 Indeed, the Cv for SSD calculated from 50 000 MC passes (instead of 120 000) is 18.7 cal/(mol K). D. Solvation of Na+ and Cl- Ions. The Na+-water (Figure 6) and Cl--water (Figure 7) distribution functions are similar to those reported in the literature.50-57 The structure of water near the Na+ is predominately controlled by the strong cationdipole interaction, which results in a sharp Na+-O peak at r ) 2.3 Å and a Na+-H peak at r ) 2.9 Å. An opposite pattern was found for water near Cl-; the first Cl--O peak (r ) 3.4 Å) is outside the first Cl--H peak (r ) 2.5 Å). The structure of water near an anion is of particular interest because of the formation of anion-water hydrogen bonds. It has been noted that the local structure of solvent near anions in aqueous solutions is very similar to that found for small anion-water clusters;50,64,65 i.e., one hydrogen of a water molecule in the first solvation shell is close to the anion, and the other hydrogen serves as a hydrogen bond donor to another water molecule, which determines the location of the second shell. This geometry is often termed a linear hydrogen bond structure, as opposed to a bifurcated hydrogen bond structure wherein the two hydrogens are equally bonded to the anion. The bifurcated geometry is a direct consequence of the monopole-dipole potential. The introduction of a quadrupole in our ion-water potential helps correct the tendency to bifurcate. A noticeable difference between the present work and those reported in the
Figure 8. P(cos θ) as a function of cos θ for the first shell (r ) 2-4 Å, solid line), the second shell (r ) 4-6 Å, dashed line), and the third shell (r ) 6-8 Å, dotted dash) of water around the Cl-, calculated from Monte Carlo simulations. See Figure 1 for the definition of angle θ.
literature (such as refs 51 and 57) is in the shape of Cl--water radial distribution functions (Figure 7). The second Cl--O peak is somewhat shifted in some cases,51,57 and the first Cl--H minimum shown in Figure 7 is not deep enough. Although including higher order moments would probably improve the situation, it must be weighted against the loss in computational speed. Also, the shift in the second Cl--O peak was observed in a polarizable water model,57 and this SSD model has no electronic polarizability. Moreover, it is not clear whether including polarizability would give this shift, since another polarizable model gives a second Cl--O peak close to that of the SSD.50 Nevertheless, the results shown here suggest that the primary solvation shell and solvation energy (see Table 3) are well represented by our current approach. Figure 8 shows the calculated orientational probability of water with respect to the ion, P(cos θ), in three consecutive shells (r ) 2-4, 4-6, and 6-8 Å), where P(cos θ) is the probability that a water has an orientation angle, θ, as defined in Figure 1. In the first shell, the most probable orientation is about cos θ ) 0.6, which is in qualitatively good agreement with results from simulations using the TIP3P and SPC potentials.66 Also, the orientation dependence decreases and approaches random at distances near 8 Å. We summarize the structural and thermodynamic results of ionic solutions in Table 3, in which the positions of the first maxima of both ion-O and ion-H radial distribution functions,
Soft Sticky Dipole Potential for Liquid Water the coordination numbers, the ion-water interaction energy, Eiw, the solvent reorganization energy, ∆Eww, and the enthalpy of solvation, ∆Hsol, are given for comparisons to results from experiments and from simulations using other water potentials. The reader is referred to papers by Chandrasekhar et al.51 by Pettitt and Rossky,49 and by Dang et al.50 for details. Although the calculated enthalpy of solvation for either ion is more exothermic than the experimental value, similar trends have been observed for other existing water models as well.50,51 Further, as has been pointed out by Chandrasekhar et al.,51 the uncertainties of the experimental ∆Hsol may be as much as (10 kcal/ mol. Thus, the overall good agreement with available experimental data suggests that the quadrupolar expansion based ionwater potential used here is a reasonable approach. E. Model Efficiency. It should be noted that the SSD model is aimed at applications to solvation of biological macromolecules through either simulations or integral equation theories, and therefore, simplicity and efficiency are as important as accuracy. Being a one-site model, SSD can be readily used in an integral equation theory.35 Moreover, SSD is also observed to be extremely efficient. For example, a 3200-pass MC run used to calculate the mean intermolecular energy and energy pair distribution requires CPU times of 0.44 h for SSD, 2.9 h for TIP3P, and 3.2 h for TIP4P, all on an IBM RISC6000 workstation. Thus, for these calculations, SSD is approximately 6.6 and 7.3 times faster than TIP3P and TIP4P, respectively. The fast speed of SSD is attributed to its simplicity, because, to describe the interaction between two water molecules, one needs to calculate only one distance (M-M) and two angles (θij,φij) for SSD but nine distances for TIP3P or SPC and ten for TIP4P. The simplicity, reasonable accuracy and the nearly 1 order of magnitude gain in computational speed clearly suggest that SSD should be useful in modeling solvent effect in many aqueous systems, especially water around biological macromolecules, because solute-solvent-type simulations are generally very slowly convergent. We must point out that the computational efficiency of the SSD model has not been tested by molecular dynamics simulations, which requires incorporating a quaternion dynamics algorithm.70 However, for models such as TIP3P, quaternion dynamics is at least as fast as normal MD and the implementations are easier,71 since SHAKE72 is no longer needed. Therefore, we estimate that SSD will also be faster than other multisite models in MD, just as in MC simulations. In calculations of the radial distribution functions in which the exact locations of the O and H are needed explicitly, SSD is still faster, with the CPU time ratio being 1:2.7:2.9 for SSD:TIP3P:TIP4P. Moreover, these calculations are usually necessary only when comparisons to multisite models are being made. Conclusions In this paper we have presented a new one-site effective pair potential model for water intended for modeling the solvent effect in biological systems. The present model, based on modifications of the hard-sphere sticky dipole potential of Bratko, Blum, and Luzar, contains a Lennard-Jones-type repulsive soft-sphere, plus a tetrahedral “sticky” potential and an embedded point dipole. Although a one-site model, it gives water structure comparable to that of the four-site TIP4P and better than that of the three-site TIP3P and the similar SPC (data not shown here for the latter). Also, this model gives an intermolecular energy, a hydrogen bond energy, and a heat capacity in good agreement with results from experiments and other models. Therefore, it represents an improvement over the BBL in both structure and intermolecular energy. In
J. Phys. Chem., Vol. 100, No. 7, 1996 2729 addition, MC simulations using this model are nearly 1 order of magnitude faster than those using TIP3P or TIP4P. The efficiency is mainly attributed to its simplicity. Moreover, this potential model can be readily used in numerical integral equation theories. Finally, a hybrid potential in which ionwater interaction is modeled by ion-dipole plus ion-quadrupole terms has been shown to give reasonable descriptions of the structure and energetics of dilute ionic solutions involving a Na+ or a Cl- ion. Overall, these results are very encouraging, since the simplicity, efficiency, and very reasonable accuracy may be utilized in many systems, particularly solvation of macromolecular solutes, such as a protein. In forthcoming papers, we will report results concerning the static dielectric properties of this model and results of applications to an integral equation theory for aqueous solvation. Acknowledgment. We are grateful to the National Science Foundation for support of this work through Grant MCB9118085 and MCB-9506796. This research was sponsored in part by the Phillips Laboratory, Air Force Material Command, USAF, through the use of the MHPCC under cooperative agreement number F29601-93-2-0001. The view and conclusions contained in this document are those of the authors and should not be interpreted as necessarily representing the official policies or endorsements, either expressed or implied, of Phillips Laboratory or the U.S. Government. Y.L. thanks Satheesan Babu, Brian Beck, Jin Kee Hyun, John Koerner, Paul Swartz, and Robert Yelle for helpful discussions. References and Notes (1) Jorgensen, W. L. J. Chem. Phys. 1981, 103, 335-340. (2) Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; Hermans, J. In Intermolecular Forces; Pullman, B., Ed.; Reidel: Dordrecht, Holland, 1981; pp 331-341. (3) Berendsen, H. J. C.; Grigera, J. R.; Straatsma, T. P. J. Phys. Chem. 1987, 91, 6269-6971. (4) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D. J. Chem. Phys. 1983, 79, 926-935. (5) Stillinger, F. H.; Rahman, A. J. Chem. Phys. 1974, 60, 15451557. (6) Zhu, S.-B.; Singh, S.; Robinson, G. W. In Modern Nonlinear Optics, Part 3; Evans, M., Kielich, S., Eds.; John Wiley & Sons: New York, 1994; Vol. LXXXX. (7) Watanabe, K.; Klein, M. L. Chem. Phys. 1989, 131, 157-167. (8) Dang, L. X.; Pettitt, B. M. J. Phys. Chem. 1990, 94, 4303-4308. (9) Reddy, M. R.; Berkowitz, M. Chem. Phys. Lett. 1989, 155, 173176. (10) Alper, H. E.; Levy, R. M. J. Chem. Phys. 1989, 91, 1242-1251. (11) Neumann, M. J. Chem. Phys. 1986, 85, 1567-1580. (12) Hasted, J. B. In Water. A comprehensiVe treatise; Franks, F., Ed.; Plenum: New York, 1972; Vol. 1; pp 255-309. (13) Rick, S. W.; Stuart, S. J.; Berne, B. J. J. Chem. Phys. 1994, 101, 6141-6156. (14) Lemberg, H. L.; Stillinger, F. H. J. Chem. Phys. 1975, 62, 16771690. (15) Lemberg, H. L.; Stillinger, F. H. Mol. Phys. 1976, 32, 353-362. (16) Stillinger, F. H.; Rahman, A. J. Chem. Phys. 1978, 68, 666-670. (17) Duh, D.-M.; Perera, D. N.; Haymet, A. D. J. J. Chem. Phys. 1995, 102, 3736-3746. (18) Kusalik, P. G. Mol. Phys. 1989, 67, 67-80. (19) Barnes, P. In Progress in Liquid Physics; Croxton, C. A., Ed.; Wiley-Interscience: New York, 1978. (20) Barnes, P.; Finney, J. L.; Nicholas, J. D.; Quinn, J. E. Nature 1979, 282, 459-464. (21) Bratko, D.; Blum, L.; Luzar, A. J. Chem. Phys. 1985, 83, 63676370. (22) Blum, L.; Vericat, F.; Bratko, D. J. Chem. Phys. 1995, 102, 14611462. (23) Blum, L.; Vericat, F. Mol. Phys. 1995, 86, 809-817. (24) Gilson, M. K.; Sharp, K. A.; Honig, B. J. Comput. Chem. 1988, 9, 327-335. (25) Gilson, M. K.; Honig, B. Proteins 1988, 4, 7-18. (26) Hansen, J.-P.; McDonald, I. R. Theory of Simple Liquids, 2nd ed.; Academic Press: London, 1986. (27) Gray, C. G.; Gubbins, K. E. Theory of Molecular Fluids; Clarendon Press: Oxford, 1984; Vol. 1.
2730 J. Phys. Chem., Vol. 100, No. 7, 1996 (28) Chandler, D.; Andersen, H. C. J. Chem. Phys. 1972, 57, 19301937. (29) Chandler, D.; Singh, Y.; Richardson, D. M. J. Chem. Phys. 1984, 81, 1975-1982. (30) Hirata, F.; Rossky, P. J. Chem. Phys. Lett. 1981, 83, 329-334. (31) Hirata, F.; Rossky, P. J.; Pettitt, B. M. J. Chem. Phys. 1983, 78, 4133-4144. (32) Ichiye, T.; Haymet, A. D. J. J. Chem. Phys. 1988, 89, 4315-4324. (33) Ichiye, T.; Haymet, A. D. J. J. Chem. Phys. 1990, 93, 8954-8962. (34) Booth, M. J.; Duh, D.-M.; Haymet, A. D. J. J. Chem. Phys. 1994, 101, 7925-7933. (35) Liu, Y.; Ichiye, T. Chem. Phys. Lett. 1994, 231, 380-386. (36) Svishchev, I. M.; Kusalik, P. G. J. Chem. Phys. 1993, 99, 30493058. (37) Liu, Y.; Ichiye, T. Unpublished work, 1994. (38) CHARMM is molecular modeling software distributed by Harvard University. QUANTA is a molecular modeling software trademark of Molecular Simulations Inc. (39) Discover and Insight are molecular modeling software trademarks of Biosym Technologies Inc. (40) Metropolis, N.; Rosenbluth, A. W.; Metropolis, M. N.; Teller, A. H.; Teller, E. J. Chem. Phys. 1953, 21, 1087-1092. (41) Owicki, J. C.; Scheraga, H. A. J. Am. Chem. Soc. 1977, 99, 74037412. (42) Pangali, C.; Rao, M.; Berne, B. J. Mol. Phys. 1980, 40, 661-680. (43) Mezei, M.; Swaminathan, S.; Beveridge, D. L. J. Chem. Phys. 1979, 71, 3366-3373. (44) Jorgensen, W. L. Chem. Phys. Lett. 1982, 92, 405-410. (45) Pangali, C. S.; Rao, M.; Berne, B. J. In Computer Modeling of Matter; Lykos, P., Ed.; American Chemical Society: Washington, DC, 1978. (46) Jorgensen, W. L. J. Am. Chem. Soc. 1979, 101, 2011-2016. (47) Kistenmacher, M.; Popkie, M.; Clementi, E. J. Chem. Phys. 1973, 59, 5842-5848. (48) Arshadi, M.; Yamdagni, R.; Kebarle, P. J. Phys. Chem. 1970, 74, 1475-1482. (49) Pettitt, B. M.; Rossky, P. J. Chem. Phys. 1986, 84, 5836-5844. (50) Dang, L. X.; Rice, J. E.; Caldwell, J.; Kollman, P. A. J. Am. Chem. Soc. 1991, 113, 2481-2486. (51) Chandrasekhar, J.; Spellmeyer, D. C.; Jorgensen, W. L. J. Am. Chem. Soc. 1984, 106, 903-910. For experimental data, see for example:
Liu and Ichiye Enderby, J. E.; Cummings, S.; Herdman, G. J.; Neilson, G. W.; Salmon, P. S.; Skipper, N. J. Phys. Chem. 1987, 91, 5851-5858. Powell, D. H.; Barnes, A. C.; Enderby, J. E.; Neilson, G. W.; Salmon, P. S. Faraday Discuss. Chem. Soc. 1988, 85, 137-146. (52) Brooks, C. L., III. J. Phys. Chem. 1986, 90, 6680-6684. (53) Impey, R. W.; Madden, P. A.; McDonald, I. R. J. Phys. Chem. 1983, 87, 5071-5083. (54) Mezei, M.; Beveridge, D. J. Chem. Phys. 1981, 74, 6902-6910. (55) Clementi, E.; Barsotti, R. Chem. Phys. Lett. 1978, 59, 21-25. (56) Heinzinger, K. Pure Appl. Chem. 1985, 57, 1031-1042. (57) Sprik, M.; Klein, M. L.; Watanabe, K. J. Phys. Chem. 1990, 94, 6483-6488. (58) Soper, A. K.; Phillips, M. G. Chem. Phys. 1986, 107, 47-60. (59) Narten, A. H.; Thiessen, W. E.; Blum, L. Science 1982, 217, 10331034. (60) Narten, A. H. J. Chem. Phys. 1972, 56, 5681-5687. (61) Palinkas, G.; Kalman, E.; Kovacs, P. Mol. Phys. 1977, 34, 525537. (62) Soper, A. K.; Silver, R. N. Phys. ReV. Lett. 1982, 49, 471-474. (63) Eisenberg, D.; Kauzmann, W. The Structure and Properties of Water; Oxford University Press: New York, 1969. (64) Rao, M.; Berne, B. J. J. Phys. Chem. 1981, 85, 1498-1505. (65) Israelachvili, J. Intermolecular and Surface Forces, 2nd ed.; Academic Press: San Diego, CA, 1992. (66) Hyun, J.-K.; Babu, C. S.; Ichiye, T. J. Phys. Chem. 1995, 99, 51875195. (67) Newsome, J. R.; Neilson, G. W.; Enderby, J. E. J. Phys. C 1980, 13, L923-L926. (68) Neilson, G. W.; Enderby, J. E. Annu. Rep. Prog. Chem., Sect. C 1979, 76, 185-220. (69) Friedman, H. L.; Krischnan, C. V. In Water: A ComprehensiVe Treatise; Franks, F., Ed.; Plenum: New York, 1973; Vol. 6. (70) Allen, M. P.; Tildesley, D. J. Computer Simulations of Liquids; Oxford University Press: New York, 1987. (71) Kuharski, R. A. Private communication. (72) Ryckaert, J. P.; Ciccotti, G.; Berendsen, H. J. C. J. Comput. Phys. 1977, 23, 327-341.
JP952324T
