6548
J. Phys. Chem. 1984, 88, 6548-6556
A Molecular Dynamics Method for Calculating the Solubility of Gases in Liquids and the Hydrophobic Hydration of Inert-Gas Atoms in Aqueous Solution William C. Swopet and Hans C. Andersen* Department of Chemistry, Stanford University, Stanford, California 94305 (Received: June 1 1 , 1984; In Final Form: September 21, 1984)
We have developed a coupling parameter method for calculating the chemical potential of a solute in a solution using molecular dynamics calculations and applied it to the study of inert-gas atoms in water. The calculations used a modified version of the revised central force model for the water-water interaction and a Lennard-Jones interaction between the solute atom and the oxygen atom of a water molecule. For some values of the Lennard-Jones E and u parameters, the calculated solubilities are similar in magnitude and temperature dependence to the measured results for the noble-gas atoms in water. The implications of these results for the computer simulation study of hydrophobic interactions are discussed.
I. Introduction Nonpolar molecules and nonpolar groups on polar molecules are often called ”hydrophobic” (Le. they fear water), and the word hydrophobic is used to describe two different aspects of the properties of aqueous solutions of such m o l e c ~ l e s . ~ - ~ “Hydrophobic hydration” refers to the way in which a single nonpolar solute interacts with and changes the structure of the water in its vicinity. Its macroscopic manifestation is the unusual temperature dependence of the solubilities of atoms and nonpolar molecules in water. Its microscopic basis is usually assumed to be a tendency for the solute to promote the formation of additional “structuie” in the water, and the nature of this structure formation has been the object of much discussion. “Hydrophobic interaction” refers to the supposed tendency for nonpolar solutes and groups dissolved in water to come close together. It is usually assumed to be a major causal factor for the immiscibility of water and nonpolar substances over a wide range of concentrations, the thermodynamic properties of aqueous solutions, the conformations of proteins and nucleic acids in solution, association equilibria in aqueous solution, micelle formation by amphiphilic molecules, and the stability of lipid bilayers and biological membranes. The basis for hydrophobic interaction is usually assumed to be the unusual way in which nonpolar molecules and groups are hydrated in water, i.e. hydrophobic hydration. Computer simulation methods, such as the molecular dynamics and Monte Carlo methods, have become the primary theoretical tool for the study of water and aqueous and many of the applications of these methods to solutions have been concerned with understanding hydrophobic The other theoretical method that has led to major insights into hydrophobic effects is the theory of Pratt and Char~dler,*~-~l some of whose predictions have been strikingly confirmed by subsequent computer simulations. The results of computer simulation calculations are dependent on the assumed forms of the intermolecular and intramolecular potentials. (Indeed their major virtue is that this is all they depend on, other than the validity of classical statistical mechanics.) Thus, it is important to test and validate the potentials used in a calculation in order to know how much to trust the results as valid for real materials. The need for this has been underscored by Morse and Rice,32 who found that some of the proposed intermolecular potentials for water are not consistent with the structure of proton ordered ice Ih. They also tested a number of intermolecular potentials for their consistency with the structures of two other forms of ice. In particular, the virtues and deficiencies of the MCY potential,33the ST2 p~tential?~ and the revised central were discussed in some detail. force This paper is concerned with whether a simple model potential for the water-water interaction and a simple Lennard-Jones ‘Present address: IBM Instruments Inc., Orchard Park, P.O. Box 332, Danbury, CT 06810.
0022-3654/84/2088-6548$01.50/0
atom-water potential are sufficient to yield simulations that contain the essential features of the hydrophobic effects. Concern for the validity of the potentials is stimulated not only by the work of Morse and Rice but also by the results of Rapoport and Scheraga,26who observed no tendency for solute atoms to aggregate in their simulations of a solution, and by the theory of Pratt and Chandler,31which predicts that the potential of mean force for two atomic solutes at infinite dilution depends sensitively on the relative magnitudes of the solute-solute and solute-water attractive interactions.
(1) Franks, F. In “Water: A Comprehensive Treatise”; Franks, F., Ed.; Plenum Press: New York, 1973; Vol. 2, p 1 . (2) Franks, F.; Reid, D. S. In “Water: A Comprehensive Treatise”; Franks, F., Ed.; Plenum Press: New York, 1973; Vol. 2, p 323. (3) Ben-Naim, A. “Water and Aqueous Solutions: Introduction to a Molecular Theory”; Plenum Press: New York, 1974. (4) For reviews of computer simulation studies of water and aqueous solutions, see ref 5-7. (5) Stillinger, F. H. Adv. Chem. Phys. 1975, 31, 1 . (6) Wood, D. W. In ”Water: A ComprehensiveTreatise”; Franks, F., Ed.; Plenum Press: New York, 1979; Vol. 6, p 279. (7) Stillinger, F. H. Science 1980, 209, 451. ( 8 ) Dashevsky, G . ;Sarkisov, G. N. Mol. Phys. 1974, 27, 1271. (9) Owicki, J. C.; Scheraga, H. A. J . Am. Chem. SOC.1977, 99, 7413. (10) Swaminathan, S.;Harrison, S. W.; Beveridge, D. L. J . Am. Chem. SOC.1978, 100, 5705. (11) Geiger, A,; Rahman, A.; Stillinger, F. H. J . Chem. Phys. 1979, 70, 263. (12) Okazaki, S.; Nakanishi, K.; Touhara, H.; Adachi, Y. J. Chem. Phys. 1979, 71, 2421. (13) Pangali, C.; Rao, M.; Berne, B. J. J . Chem. Phys. 1979, 71, 2975. (14) Pangali, C.; Rao, M.; Berne, B. J. J . Chem. Phys. 1979, 71, 2982. (15) Rossky, P. J.; Karplus, M. J . Am. Chem. SOC.1979, 101, 1913. (16) Swaminathan, S.; Beveridge, D. L. J . Am. Chem. SOC.1979, 101, 5832. (17) Alagona, G.; Tani, A . J . Chem. Phys. 1980, 72, 580. (18) Okazaki, S.;Nakanishi, K.; Touhara, H.; Adachi, Y . J . Chem. Phys. 1980, 72, 4253. (19) Bolis, G . ;Clementi, E. Chem. Phys. Lett. 1981, 82, 147. (20) Nakanishi, K.: Okazaki. S.: Ikari. K.: Touhara. H. Chem. Phvs. Lett. 1981, 84, 428. (21) Okazaki, S.; Nakanishi, J.; Touhara, H.; Watanabe, N. J . Chem. Phys. 1981, 74, 5863. (22) Bok, G.;Corongiu, G.;Clementi, E. Chem. Phys. Lett. 1982,86,299. (23) Alagona, G.; Tani, A. Chem. Phys. Lett. 1982,87, 337. (24) Jorgensen, W. L. J. Chem. Phys. 1982, 77, 5757. (25) Chandrasekhar,J.; Jorgensen, W. L. J. Chem. Phys. 1982, 77,5080. (26) Rapoport, D. C.; Scheraga, H. A. J. Phys. Chem. 1982, 86, 873. (27) Jorgensen, W. L.; Madura, J. D. J . Am. Chem. SOC.1983,105, 1407. (28) Rosenberg, R.; Mikkilineni, R.; Berne, B. J. J . Am. Chem. SOC.1982, 104, 7647. (29) Pratt, L. R.; Chandler, D. J. Chem. Phys. 1977, 67, 3683. (30) Pratt, L. R.; Chandler, D. J . Chem. Phys. 1980, 73, 3430. (31) Pratt, L. R.; Chandler, D. J. Chem. Phys. 1980, 73, 3434. (32) Morse, M. D.; Rice, S.A. J. Chem. Phys. 1982, 76, 650. (33) Matsuoka, 0.;Clementi, E.; Yoshimine, M. J . Chem. Phys. 1976,64, 1351. (34) Stillinger, F. H.; Rahman, A . J . Chem. Phys. 1974, 60, 1545. (35) Stillinger, F. H.; Rahman, A . J. Chem. Phys. 1978, 68, 666.
0 1984 American Chemical Society
Inert-Gas Atoms in Water
The Journal of Physical Chemistry, Vol. 88, No. 26, 1984 6549
One of the simplest quantitative manifestations of hydrophobic effects is in the solubility of inert-gas elements and small organic molecules in aqueous solution.’ These solubilities are much smaller and have a more negative temperature derivative and a more positive second temperature derivative than do the solubilities of the same solutes in a nonpolar solvent such as cyclohexane. One way to test a set of water-water and atom-water potentials is to see whether they can predict the observed features of the solubilities of inert-gas atoms in water. The solubility is more difficult to calculate than the energy, pressure, and pair correlation functions normally obtained from molecular dynamics and Monte Carlo calculations. This is because the solubility is related to the chemical potential, which has entropic as well as energetic contributions. Typically, entropic quantities like free energies, chemical potentials, and the entropy itself must be obtained by a thermodynamic integration or coupling parameter integration. Evaluation of such quantities by averaging some dynamical property over a single molecular dynamics trajectory or Monte Carlo chain is in principle possible but subject to severe practical d i f f i c ~ l t i e s . ~These ~ difficulties are especially severe for calculating partial molar quantities for solutes, where the desired quantity is calculated as the small difference between two large number^.^ We have developed a coupling parameter method for calculating the chemical potential of a solute in a solution and applied it to the case of an atomic solute in water. We use a modification of the revised central force model potentiaP of Stillinger and Rahman as the water-water intermolecular potential and a Lennard-Jones potential for the atom-water interaction. The calculations are for a range of values of the Lennard-Jones E and u parameters and for a range of temperatures. We find that for some, but not all, values o f t and u, this simple model system has the solubility properties of real inert gases in water and thus exhibits the macroscopic manifestations of hydrophobic hydration. 11. Theory of the Solubility Calculation The basis for the present calculation of the solubility is a formula for the chemical potential of a solute based on the use of a coupling parameter. We imagine adding one more solute atom to a solution but that the interaction between the added atom and the other atoms and molecules in the solution is multiplied by a coupling parameter, t, which can take on any value from zero to unity. We average the interaction between the added solute and the rest of the solution over an ensemble whose Hamiltonian contains the coupling parameter. Let U ( t )be this average. (Note that U ( t ) is the average of the full potential of interaction of the solute with the solution for a system in which that interaction is only partially effective in determining the dynamics and equilibrium distributions.) Then the chemical potential, ph, of the atom in the solution is given by
where pAs is the number density of solute atoms in the solution and A is the usual thermal deBroglie wavelength of the solute. This result is derived for the special case of a one-component atomic fluid by Hill;’ but the same method can be used to derive the result for a solution of atoms in a molecular fluid. The derivation uses a canonical ensemble, but in the thermodynamic limit any of the customary ensembles can be used to evaluate U(f). If the solution is very dilute, then U ( t ) can be calculated by considering a system containing just one solute atom and a large number of solvent molecules. A common definition of the solubility s of a gas in a liquid is the mole fraction of the solute in a solution that is in equilibrium with a gas phase containing the same solute with a partial pressure of 1 atm. Let us assume that such a gas is ideal. Then the chemical potential of the solute in the gas phase can easily be calculated to be ~~
where pAg is the density of A in the gas. When the gas phase and the solution are in equilibrium, the chemical potential of A must be the same in both phases. Equating the chemical potentials, we find
When the solution is very dilute, the mole fraction of solute in the solution is accurately given by pAs/pwo,where the quantity in the denominator is the number density of pure solvent. Let p A 2 be the density of solute in the gas phase when its partial pressure is 1 atm. It follows from this last equation and the definition of the solubility s that
This is the basic theoretical equation used in the present work for calculation of solubilities. To evaluate U ( t ) ,we assume some forms for the water-water potential and the solute-water potential and perform constant temperature-constant pressure molecular dynamics calculation^^^ for a system containing one solute atom and several water molecules. In the calculations, the solute-water potential is multiplied by a factor of when it is used to calculate the forces of interaction between the solute and the water. At each time step, the total energy of interaction of the solute with all the water molecules (without the factor of E ) is evaluated and the results is averaged over time. The time average, for a long enough trajectory, is then assumed to be equal to the ensemble average, and so the time average is just U ( t ) . A practical problem associated with the implementation of this procedure arises when the solute-water interaction is assumed to diverge more rapidly than r-3 for small separations r of the solute and the water molecule. In the motion generated by the dynamics calculation for small values of E, the solute can come very close to water molecules since the forces that would prevent this have been multiplied by 5. Such close encounters give large instantaneous values of the solute-water interaction energy. This energy then fluctuates greatly during the trajectory, and the resulting time average has a large statistical error. A strategy that overcomes this problem to a large extent is to modify the potential by eliminating the divergence in the potential at short distance. This is done by assuming a form for the potential that is very large but finite at small distances but that is equal to the desired potential for all distances that can be physically realized in the fully coupled ( E = 1) system. The calculations therefore give, in principle, the solubility of a modified solute whose interaction potential with water is different from the potential originally postulated. Nevertheless, the “modified” solute and the desired “original” solute have the same solubility. (Indeed, all the thermodynamic properties of the two solutes are identical, as long as they are in thermodynamic states that do not sample configurations in which the two potentials are different.) Since the solute-water potential actually used in the simulation is finite, U ( t ) does not diverge but stays finite at approaches 0. This reduces the statistical error in its evaluation. (For 6 = 0, Ucan actually be evaluated analytically if the potential is of a simple enough form.) 111. Water-Water and Solute-Water Interactions In this section we define the intermolecular and intramolecular potentials used in the molecular dynamics calculations. We regard a water molecule as a set of three vibrating mass points. The intermolecular potential we used is a modification of the revised central force model of Stillinger and R a h n ~ a n . ~ ~ For the revised central force model, there is an oxygen-oxygen potential, an oxygen-hydrogen potential, and a hydrogen-hydrogen potential. Each of the three potentials is a function of the scalar
~~
(36) Owicki, J. C.; Scherage, H. A. J . Am. Chem. SOC.1977, 99, 7403. (37) Hill, T. L. “Statistical Mechanics”; McGraw-Hill: New York, 1956.
(38) Andersen, H. C. J. Chem. Phys. 1980, 72, 2384.
6550 The Journal of Physical Chemistry, Vol. 88, No. 26, 1984 distance between the atoms. The energy of interaction between two water molecules is a sum of interactions between the atoms on each molecule. The intramolecular vibrational potential energy of a molecule is also a sum of the interactions between each pair of atoms in the molecule. The behavior of the OH and HH potentials at short distances was adjusted to give a local minimum in the intramolecular energy at the observed geometry of the water molecule. The three atomatom potentials at larger distances were adjusted to fit known properties of the water dimer and other properties of water. The potential we used was the revised central force model described above, with three modifications. First, the atom-atom potentials for the intramolecular energy were simple parabolic potentials with minima at the observed intramolecular distances and curvatures that matched those of the revised central force potentials. (This was done for simplicity in the present calculation and in previous ~ o r k . ~ ~Second, , ~ ~ )in order to eliminate the long-ranged part of the potential, the sum of atom-atom intermolecular potentials was multiplied by a switching function designed to leave the potential unmodified for oxygen-oxygen distances less than 5.5 A, to make the potential zero for distances larger than 6 A, and to make the potential go smoothly to zero in that range of distances. (The characteristics of central force potentials switched off in this way are discussed in detail by Andrea et aL40 They found that such a modification of the potential has a small effect on the thermodynamic properties and on the radial distribution functions. This type of potential has the advantage of allowing energy conservative dynamics, unlike the usual scheme of setting the force to zero for large intermolecular distances.) Using the potentials described above, constant temperatureconstant volume molecular dynamics calculation^^^ on a system of 512 molecules were performed for a temperature of 29.5 OC and a density of 1 g/cm3. The average pressure of the fluid under those conditions did not match the experimentally observed pressure of 1 atm. We found, however, that a very small modification of the intermolecular oxygen-oxygen potential was enough to give the correct pressure. By trial and error we decided on a correction of -0.090 35 kcal/mol to be added to this potential. In summary, the atom-atom potentials for intramolecular interactions are
VoH(r)= (1/2)(1147.65)(r - 0.9584
kcal/mol
VHH(r)= (1/2)(257.25)(r - 1.5151 A)2 kcal/mol
Swope and Andersen a 4
\
b\
* 200.K b = 303'K c = 423'K
3 C
\
000 2
1
0
2
Figure 1. Oxygen-oxygen radial distribution function for the water model used in these calculations at three temperatures and a pressure of 1 atm.
2t I ob
I
IV I
I
2
I
and the intermolecular potential is
r,
- r121) I J'l
1
6
I
I
i
8
A
Figure 2. Oxygen-hydrogen radial distribution function for the water
3
U(rll~r21rr31~r12~r22~r32)= S(lrll
I
4
Vij(lril
- r,21)
where r,, is the position of the ith atom on thejth molecule and where the subscripts on the V,, functions should be 00, OH, or HH, as is appropriate to the nature of the atoms i and j . The intermolecular V functions were the same as those of Rahman and Stillinger, with the exception of the correction of -0.090 35 kcal/mol that was added to the oxygen-oxygen function. The revised central force model is a rather good model for water and has been fit to a number of properties of the dimer and the liquid phase. It is consistent with the structure and densities of the low-pressure forms of ice studied by Morse and Rice.32 The modifications we made in the potential have no effect on the short-range forces, and thus they have no effect on the geometry of the dimer and make a trivial change in the energy of the dimer. The elimination of the long-range forces has a small effect on the liquid structure, primarily in the short-range orientational correlation functions. For the crystal, the short-ranged forces probably determine the structure even more completely than in the liquid. Thus, the potential we used can be regarded as equivalent to the revised central force model in most respects. (39) Swope, W. C.; Andersen, H. C.; Berms, P. H.; Wilson, K. R. J . Chem. Phys. 1982, 76, 631. (40) Andrea, T.A,; Swope, W. C.; Andersen, H. C. J . Chem. Phys. 1983, 79, 4516.
model at 303 K and 1 atm.
a
OO
2
6
4
r,
0
A
Figure 3. Hydrogen-hydrogen radial distribution function for the water model at 303 K and 1 atm.
Figure 1 shows the oxygen-oxygen radial distribution function for our model potential at three different temperatures for a pressure of 1 atm. Figures 2 and 3 show the OH and HH radial distribution functions at one temperature. These results were ootained by molecular dynamics runs for systems of 512 molecules that each lasted for a time of 3.75 ps after thorough equilibration.
The Journal of Physical Chemistry, Vol. 88, No. 26, 1984 6551
Inert-Gas Atoms in Water MOLAR VOLUME OF WATER
I9St
m j
Q
u and its first derivative at the inner cutoff.
I
IV. Method of the Solubility Calculation To perform solubility calculations using the method described above, it is necessary to perform simulations on a fluid containing water molecules and one solute atom. The solute atom is partially coupled to the water molecules, in the sense that forces acting between solute and water molecules are multiplied by a coupling parameter. The time average of the solute-water interaction energy must be calculated during the simulation. Such time averages, as a function of the coupling parameter, are the integrand for the integral in eq 2.1 for the solubility. The atom-water potential we use is of the form
Q
E 0
2
19.0
P
tl
a
/
Then eq 2.1 can be rewritten as
TEMPERATURE, *C
Figure 4. Molar volume of liquid water at 1 atm as a function of temperature from experiment (solid curve) and various molecular dynamics calculations (circles and crosses with error bars). Experimental data are from ref 41 and 42. The calculations represented by circles were performed for a 64-molecule system and were typically simulations of from 150 to 300 ps after lengthy equilibrations. The error bars represent f l standard deviation estimated from the volume fluctuation autocorrection function. The calculations represented by crosses were performed on 512 molecules and were typically similations of from 160 to 270 ps after lengthy equilibrations. The error bars for these calculations represent f l standard deviation estimated by breaking the long simulations into sets of 30 ps each and considering each 30-ps average volume to be statistically independent of the others. Note the consistent size dependence of the results, the molar volume being smaller for the larger system. The effect is greatest at the lowest temperature. Also note the lack of reproducibility at the lowest temperature, presumably a result of the long relaxation times of such a fluid.
The sharp first peak in the 00 function at about 2.82 8, with about 4.26 neighbors and the broader second peak at about 4.5 8,show that this model fluid has the tetrahedral structure characteristic of water. The major discrepancy between the structure of the model fluid and that of real water is that the tetrahedral structure is somewhat too pronounced in the model. Figure 4 shows the calculated and e ~ p e r i m e n t a volume l ~ ~ ~ ~as~ a function of temperature at a constant pressure of 1 atm. The temperature derivative of the model fluid's density is small but apparently nonzero near 0 OC;so the model does not have a density maximum. For the solute-water interaction, we used a Lennard-Jones potential that is a function of the solute-oxygen distance. Similar model potentials have been used in other model calculations of nonpolar solutes in water. In practice, we actually used a Lennard-Jones potential that was truncated at 6 8, and shifted upward; that is u(r) = ULj(t) - U~j(6) =O
forr>6
for r
