3470
J. Phys. Chem. 1985,89, 3470-3413
point contributions to the EA’S have been calculated and found to be much less than 0.1 eV except when detachment is accompanied by a change in the number of vibrational degrees of freedom. Acknowledgment. The author thanks the University of Michigan Computing Center for use of its AmdahlS860 computer.
He also thanks Professors H. B. Schlegel and W. C. Lineberger for discussions and suggestions. Registry No. H2CN-, 28892-56-0; HCNH-, 90623-30-6; CNHF, 96808-04-7; HzCN, 15845-29-1; H C N H , 15691-95-9; HZCC-, 6629536-1; H,CC, 2143-69-3; HCN-, 12334-27-9; CNH-, 96913-22-3; HCN, 74-90-8.
Monte Carlo Simulations of Alkanes in Water: Hydratlon Numbers and the Hydrophobic Eftect William L. Jorgensen,* Jiali Cao, and C. Ravimohan Department of Chemistry, Purdue University, West Lafayette, Indiana 47907 (Received: March 6, 1985)
Monte Carlo statistical mechanics simulations have been used to calculate the number of water molecules in the first hydration layer around a series of alkanes in aqueous solution. This key parameter in theories of aqueous solutions is found to have an excellent linear correlation with entropies of solution. Consequently, strong support is provided for the traditional idea that the ordering of water molecules in the first layer around a hydrocarbon is primarily responsible for the large entropy losses associated with the hydrophobic effect. However, the energetic results do not support the necessity of invoking water structure making to explain the substantial exothermic heats of hydrations for alkanes. Insights are also obtained concerning hydrophobic effects on conformational equilibria and the relative solubilities of branched and straight-chain isomers.
Introduction
The notion of “iceberg” formation around nonpolar solutes in water, introduced by Frank and Evans,‘ has become a central feature of descriptions of the hydrophobic effect.*-’ According to this idea, water molecules in the first layer around such solutes are more ordered and participate more in hydrogen bonding than bulk water. This may explain the exothermicity and substantial entropy loss that accompany the dissolution of rare gases and hydrocarbons in water.’” Furthermore, the free energy, enthalpy, and entropy of solution for the nonpolar solutes are anticipated to be proportional to the number of water molecules in the surface layer.2*4vs,7-9Indeed, roughly linear correlations have been reported for these quantities with solute cavity surface area,8,9 solute number of carbons,2 or number of hydrogens in the solute.’ Not surprisingly, the number of water molecules in the first solvation shell is also a key parameter in theoretical treatments of aqueous solutions including the significant structure theory of Eyring and co-workers7 and the Nemethy-Scheraga modeL4 However, the methods for estimating the number of water molecules in the first layer4 or cavity surface areas8v9are crude and not unique, while the correlations with numbers of hydrogenss have been criticized for having limited theoretical basis.2 Consequently, there is need for more accurate, unbiased estimation of hydration numbers for nonpolar solutes. These can be (1) Frank, H. S.; Evans, M. W. J . Chem. Phys. 1945, 13, 507. (2) Tanford, C. “The Hydrophobic Effects”, 2nd 4.; Wiley-Interscience: New York, 1980. (3) Ben-Naim, A. “Hydrophobic Interactions”; Plenum Press: New York, 1980. (4) Nemethy, G.; Scheraga, H.A. J . Chem. Phys. 1962, 36, 3401. ( 5 ) (a) Gill, S. J.; Wadso, I. Proc. Norl. Acad. Sci. U S A . 1976, 73, 2955. (b) Dec, S. F.; Gill, S. J. J. Solution Chem. 1984, 13, 27. (6) Abraham, M. H. J. Am. Chem. SOC.1982,104, 2085. (7) (a) Jhon, M. S.; Van Artsdalen, E. R.; Grosh, J.; Eyring, H. J . Chem. Phys. 1966, 44, 1465. (b) Eyring, H.; Jhon, M. S. ‘Significant Liquid Structures”; Wiley-Interscience: New York, 1969. (8) Hermann, R. B. J . Phys. Chem. 1972, 76, 2754. (9) Reynolds, J. A.; Gilbert, D. B.; Tanford, C. Proc. Norl. Acad. Sci. U.S.A. 1974, 71, 2925. (10) Leo, A.; Hansch, C.; Jow, P. Y. C. J. Med. Chem. 1976, 19, 611. ( I 1 ) Cramer, R. D., 111 J. Am. Chem. SOC.1977, 99, 5408.
0022-3654/85/2089-3470%01 SO10
used to better test the notions described above on hydrophobic effects as well as to obtain a more detailed knowledge of the structure of aqueous solutions. As described here, this has been achieved for a series of seven alkanes through Monte Carlo statistical mechanics simulations of the corresponding dilute aqueous solutions. The hydration numbers obtained in this way reflect proper thermal averaging and weighting for the various conformational isomers of the solutes. Correlations with experimental thermodynamic data and prior indices are presented. In addition, the present thermodynamic results indicate that the traditional rationalization of the exothermic heats of hydration for alkanes through enhanced hydrogen bonding is unnecessary. Insights are also obtained into the origin of the solubility differences between branched and straight-chain isomers. Computational Procedure
Monte Carlo Simulations. Statistical mechanics simulations were carried out for dilute solutions of methane, ethane, propane, butane, pentane, isobutane, and neopentane in water at 25 OC and 1 atm using the isothermal-isobaric (NFT) ensemble. Butane in water was also modeled at 60 OC and 1 atm. The systems consisted of 1 solute plus 216 water molecules in a cube with periodic boundary conditions except for methane which had been studied earlier by using 125 solvent molecules. The computational procedure was essentially the same as we have used for many other systems, including aqueous solutions of amides12and small ions,13 and an earlier study of butane in water.I4 This includes preferential sampling of the solute and near solvent molecules using I/($ + c) weighting and umbrella sampling over chopped barriers for torsional degrees of freedom. The earlier works may be consulted for details.lS In each case, 750-1000K configurations were discarded in equilibrating the systems followed typically by an additional 3000K configurations for averaging the computed properties, though for the solution of pentane 5000K configurations (12) Jorgensen, W. 1.; Swenson, C. J. J . Am. Chem. Soe. 1985,107, 1489. (13) Chandrasekhar, J.; Spellmeyer, D. C.; Jorgensen, W. L. J. Am. Chem. SOC.1984, 106, 903. (14) Jorgensen, W. L. J. Chem. Phys. 1982, 77, 5757. (15) For a review, see: Jorgensen, W. L. J . Phys. Chem. 1983,87, 5304.
0 1985 American Chemical Societv
The Journal of Physical Chemistry, Vol. 89, No. 16, 1985 3471
Monte Carlo Simulations of Alkanes in Water
~~
CH, group alkane methane
ethane propane butane isobutane pentane
neopentane
c1 20.3 11.5 10.6 10.0 9.5 10.3 8 .o
c2
YO 6o
c3
2 wuEmAu(ANEs 2 5 c
TABLE I: Computed Coordination Numbers for Alkanes in Water at 25 OC
total 20.3 23.0
6.1 5.0 1.7 4.7 0.0
21.3
29.9 4.1
30.3 34.0 32.0
were generated for the averaging. The intermolecular potential functions used in the computations have been reported previously. Specifically, the TIP4P model was adopted for water16 and the Lennard- Jones parameters for the alkanes are the OPLS set that were optimized for liquid hydroc a r b o n ~ . ’ ~These potential functions have been shown to yield excellent thermodynamic and structural results for pure water and hydrocarbons: average errors in computed densities and heats of vaporization are less than 2%.16,17The general form of the potential functions (eq 1) includes Coulomb and Lennard- Jones
terms acting between interaction sites on the two monomers a and b. The computational procedure was essentially the same as we have used for many other systems, including aqueous solutions of amidesI2 and small ions” and and earlier study of butane in water.14 The interaction sites are located on nuclei with the addition of a negatively charged site along the HOH bisector for water.I6 For the hydrocarbons, the CH, groups are represented in a united atom model with the hydrogens implicit and they are ~ n c h a r g e d . ’ ~Standard geometries based on experimental structural data were used for the water and alkane mor~omers.’~*’’ Although the intramolecular bond lengths and angles were fvred in the simulations, torsional motion was included for butane and pentane.15 Thus, when these molecules are moved in generating new configurations, not only are they translated and fully rotated but also changes in the dihedral angles are attempted. This necessitates the inclusion of appropriate torsional potential functions. These were also reported previously; they are based on MM2 calculations and are represented by three-term Fourier series plus a term for the Cl-C5 interaction in pentane.” Proximity Analyses. To calculate coordination numbers for nonspherical solutes, Mehrotra and Beveridge have defined a proximity criterion to assign first-shell solvent molecules to atoms of the solute.I8 This requires a cutoff distance for each type of solute atom. For CH, groups,we use 5.35 A for the carbon-water oxygen distance.I2 This value corresponds to the location of the first minimum in the CO radial distribution functions for methane and other alkanes in water at 25 OC. Then, if a water molecule is nearer one CH, group than any other in the solute and it is within the cutoff range, it can be counted as in the primary coordination shell for the group. Thus, the total coordination number is directly broken down on a group by group basis. This analysis also permits numerous distributions for solute-solvent energetics to be defined as we11.I8 In the present case, the proximity analyses were performed on configurations saved at 5K intervals during the Monte Carlo runs. All computations were carried out on the Harris Corp. H-80 computer in our laboratory. Results and Discussion Hydration Numbers. The computed primary coordination numbers for the alkanes in water are listed in Table I. From (16) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R.W.; Klein, M. L. J . Chem. Phys. 1983, 79, 926. (17) Jorgensen, W. L.; Madura, J. D.; Swenson, C. J. J . Am. Chem. Soc. 1984, 106, 6638. (18) Mehrotra, P. K.; Weridge, D. L. J . Am. Chem. Soc. 1980,102,4287.
-40
t *O
23
2BORO.
8.
32
35
Figure 1. Correlations between experimental thermodynamic data on the hydration of alkanes (ref 5b) and the computed hydration numbers (Table I).
separate averages for blocks of configurations the uncertainties in the total coordination numbers are estimated to be about *OS. The coordination numbers range from 20 for methane to 34 for pentane with an average increment of ca. 3 waters per carbon atom. The larger difference (4.3) between ethane and propane follows from the limited shielding of the CH2 group in propane relative to CH2 groups in higher n-alkanes due to the absence of /3 alkyl groups. Though the difference in hydration numbers for butane and isobutane is not statistically significant, the order for the pentanes is likely correct and indicates a decrease in hydration number with branching. Considering the individual group results, the shielding of the methyl groups by lengthening the n-alkane chain is apparent; however, the effect seems to be leveling off at butane and pentane which both have roughly 10 waters in the primary coordination shell for each methyl group. The influence of adjacent branching is also evident. Thus, the average number of primary water molecules per methyl group decreases from 10.0 for butane to 9.5 for isobutane to 8.0 for neopentane. The shielding of the interior groups is far greater; the C2 methylene groups for propane, butane, and pentane have 6.1,5.0, and 4.7 water molecules in the primary region, while the central methylene group for pentane has only 4.1. The C H group in isobutane is even more shielded (1.7), and the quaternary carbon of neopentane has no water molecules nearer it than other carbons. Correlations between the computed coordination numbers and experimental free energies, enthalpies, and entropies of hydration for the hydrocarbons are shown in Figure 1. The experimental data are from the recent compilation of Dec and Gill and are the standard partial molar thermodynamic quantities for transfer of the gaseous hydrocarbons into liquid water at 25 OCSSbLinear relationships are found and may be summarized in the following equations where C is the coordination number: AGO (kJ/mol) = 0.17C 22.04 (2)
+
AHo (kJ/mol) = -0.97C
+ 4.54
(3)
ASo (J/(mol deg)) = -3.83C - 59.03
(4) The equations were obtained from least-squares analyses and gave assoCiated standard deviations of 0.67 kJ/mol of AGO, 1.63 kJ/mol for pHo,and 3.41 J/(mol deg) for ASo. The linearity is particularly impressive for the entropy as shown on an expanded scale in Figure 2. This supports the notion that the ordering of water molecules in the first layer around a hydrocarbon is primarily responsible for the large entropy losses associated with the hydrophobic effect. The linearity for the enthalpy is also good and may be made even more precise by using the AHo for neopentane (-27.8 kJ/mol) reported in the review by Wilhelm et al.” In view of these correlations, linear behavior ~
(19) Wilhelm, E.; Battino, R.; Wilcock, R. J. Chem. Reu. 1977, 77, 219.
3472 The Journal of Physical Chemistry, Vol. 89, No. 16, 1985
Jorgensen et al.
t?QUEOUS RLKANES 25 C
35
RDUEOUS RLKANES 25 C
,
1
I
32
$ 29
Q8 26 -lgY -210
1
23
I
!
2o
23
%RD.
$9.
32
20
1
150
35
180
&fmE
&
I
1
270
300
Figure 4. Correlation between Hermann's solute cavity surface areas (A2) for alkanes from ref 8 and their hydration numbers. 35
-
32
-
$ 29
-
E0
8 26
23 20
TABLE II: Calculated Conformer Populations for Alkanes in the Gas Phase and Aqueous Solution alkane T,OC conformer %gas % aq butane 25 t 68.2 54 f3 31.8 46 butane 60 t 64.1 58 g 35.3 42 pentane 25 tt 46.5 61 47.1 36 5.4 1 g+g1.o 2
1
-
-.:;-" 1 ~
&+
,
II).14915920
15
BO. OF t%OROGEdi Figure 3. Correlation between the number of hydrogens for an alkane and its hydration number.
is also found for AGO, though the small variation in AGO for this series of compounds must be noted. Thus, the dependence on coordination number is slight (eq 2). Correlations were also sought between the coordination numbers and previously reported indices. For the present alkanes, the correlation with number of hydrogens is excellent as illustrated in Figure 3 and represented in eq 5. Whether the accord extends C = 1.61NH
+ 13.82
(a = 0.74 kJ/mol)
about 15% for the gauche population for butane in water as compared to the gas phase at 25 OC (Table Thus, the conformational effect cannot be attributed to a significant reduction in the number of water molecules in the first shell around the gauche conformer as compared to the trans. The solvent effect at 60 O C is diminished to a 7% increase in the gauche population, though it should be pointed out that the statistical uncertainties (&a) for all of the aqueous conformer populations in Table I1 are ca. f4%. Nevertheless, the qualitative trend is consistent with the general temperature dependence of hydrophobic effect^.^^^ In contrast, for pentane at 25 OC there is a 14.5% increase for the trans-trans conformation in water and about 9% more trans per dihedral angle. The conformational results from this comparatively long simulation appear stable and show good symmetry for the two dihedral angles. Further study of the origin of these solvent effects on conformation is desirab1e.l5 Energetics. The total energy of the solutions (eq 7 ) consists
(5)
to the other types of hydrocarbons treated by Gill et aL5 including alkynes, alkenes, and cycloalkanes remains to be determined. The correlation with Hermann's cavity surface areas* (in his units of A*) is somewhat better, as shown in Figure 4 and represented in eq 6. C = 0.104A + 3.96 (a = 0.57 kJ/mol) (6)
In closing this section, we note some additional results for aqueous butane. Coordination numbers were separately determined for the gauche and trans conformers at both 25 and 60 O C . At 25 OC, the number of primary water molecules for trans and gauche is 30.06 and 29.80, and at 60 OC, 28.24 and 28.28. Therefore, we do not find a statistically significant difference in the primary hydration numbers for the two conformers at either temperature. Of course, there are shifts in the individual group hydration numbers. The methyl groups shield each other more in the gauche form while the methylene groups are simultaneously more exposed. The primary coordination numbers for the methyl groups in the trans and gauche conformers at 25 OC are 10.35 and 9.65, while the values for the methylene groups are 4.68 and 5.25. These findings are particularly interesting since the present calculations and earlier results consistently find an increase of
of the solutesolvent energy (E& the solventsolvent energy (Es), and the intramolecular rotational energy for the solute (Eintm)in the cases of butane and pentane. The energy change on transferring the solute from the gas phase to the solution (AEO) can be expressed by eq 8-10, where E*,, is the energy for the pure AEa = ET - (E*,, + Eintrai')
(8)
liquid, Ehmig is the torsional energy for the solute in the ideal gas phase, and AE, is the solvent reorganization energy. The enthalpy of solution may then be given by eq 11, where A P is the partial AHo = AEO
+ P A P - RT
(11)
(20) (a) Pratt, L. R.; Chandler, D. J. Chem. Phys. 1977, 67, 3683. (b) Rcsenberg, R. 0.;Mikkilineni, R.; Berne, B. J. J. Am. Chem. Sac. 1982,104, 7647.
The Journal of Physical Chemistry, Vol. 89, No. Id, 1985 3473
Monte Carlo Simulations of Alkanes in Water TABLE IIk Energetic Results (kJ/mol) for Alkanes in Water at 25 'C
alkane methane ethane propane butane isobutane pentane neopentane
E,x -12.1 -20.1 -27.6 -34.7 -32.7 -41.8 -36.0
&a
-7 -59 -26 -35 -56 4 -64
L=intm
0.4 -0.5
A E O
-19 -80 -53 -69 -89 -39 -100
AHo -2 1 -82 -56 -72 -9 2 -4 1 -103
AHO(expt1)" -13.2 -19.5 -23.3 -25.9 -24.2 -28.3 -25.1
,I Reference 5b. TABLE I V Computed Alkane-Water Interaction Energies (kJ/mol) for Water Molecules in the Primary Coordination Shell at 25 OC CH, group alkane c1 c2 c3 methane -0.46 ethane -0.68 propane -0.73 -0.88 butane -0.80 -1.05 isobutane -0.81 -1.23 pentane -0.81 -1.20 -1.24 neopentane -0.85
molar volume of the solute and the last term is the PVcontribution for the solute in the gas phase. AVO is obtained as the difference in the volume of the solution and the pure solvent. The reference values for pure TIP4P water with 216 monomers at 25 OC and 1 atm were obtained previously (E*, = -9058 f 10 kJ/mol, P = 6.575 f 0.014 nm3).16 Eintr:g values for butane (2.6 kJ/mol) and pentane (5.0 kJ/mol) were obtained from Boltzmann distributions for the rotational potential functions a t 25 O C . 1 7 The computed results for AEo and its components (eq 10) are given in Table I11 along with the calculated and experimental enthalpies of solution. The E,, values are obtained with good precision in the simulations ( u = 0.2-0.3 kJ/mol). However, the uncertainties in AE, and consequently AEo and AHo are substantial. From the fluctuations in the averages over blocks of 50K configurations 0's of 12-1 6 kJ/mol are obtained, though it appears that the true uncertainties may be 2-3 times larger.21 The problem is that two large, fluctuating numbers are being subtracted. The same is true for the computations of AVO; the results overlap the experimental range but have little analytical value. The computed enthalpies of solution are all exothermic, though they do not reproduce the experimental order a t all well. In contrast, the more precise ESxvalues nicely parallel the experimental AHO's. The only change in order is for butane and neopentane; however, this difference is removed by using again the AH" for neopentane (-27.8 kJ/mol) from the review of Wilhelm et al.19 It should also be noted that with the exception of methane the computed E, values are by themselves enough to account for the experimental heats of solution for the alkanes. The E , values consist solely of the Lennard-Jones interactions between the solute and the water molecules. In addition, it should be kept in mind that the Lennard-Jones parameters that have been used yield excellent energetic results for pure water and alkane liquids. l 6 9 l 7 Thus, it appears ~ ~ e c e s ~toa invoke r y enhanced hydrogen bonding and negative solvent disruption energies (structure making) to explain the exothermic heats of hydration for This position is consistent with our earlier analyses of the hydrogen bonding in the first hydration shell for butane.14 Substantial spectroscopic evidence has also been presented that indicates, if anything, somewhat decreased hydrogen bonding upon dissolving nonpolar solutes in water.22
Further details on the solutesolvent interactions were obtained from the proximity analyses. Specifically, the interaction energies between the solute and the water molecules in the primary solvation layer for each alkyl group were computed as listed in Table IV. Several points may be noted. (1) The average alkanefirst shell water attraction is about 1 kJ/mol which is, of course, much weaker than the average water-water hydrogen bond (17 kJ/mol) for TIP4P water a t 25 O C . 1 6 (2) The data in Tables I, 111, and IV can be used to show that the interactions between the alkanes and the water molecules in the first coordination layer account for about 75% of the total solutesolvent interaction energies (Em). (3) The most important observation is the greater attractions for the water molecules nearer the centers of the alkanes than close to the methyl groups. The weaker alkane-water interactions for water molecules near methyl groups lead to the less negative E,, values for more branched alkanes than for their isomers with fewer methyl groups. This, in turn, can account for the greater solubility of the less branched alkanes (smaller AGO). That is, the experimental AAGO values for isobutane vs. butane (+1.O kJ/mol) and neopentane vs. pentane (+0.6kJ/mol) result from the A A H O contributions (+1.7, +3.2) dominating -TAAS0(-0.6, -2.7).5b These AAH" values are remarkably similar to the total interaction energy differences for the solutes with the first-shell water molecules (+1.3, +3.7) that can be computed for the two sets of isomers from the data in Tables I and IV.
Conclusion The present Monte Carlo simulations allowed calculation of the number of water molecules in the first hydration layer around a series of alkanes. This is an important quantity in theories of aqueous solutions that has not been computed previously in such a sophisticated, direct manner including proper thermal and conformational averaging. The excellent linear correlation found between the hydration numbers and entropies of solution supports the traditional view that the ordering of water molecules in the first layer around a hydrocarbon is responsible for the large entropy losses associated with the hydrophobic effect. However, the energetic results do not support the need to invoke structure-making effects to rationalize the exothermic heats of hydration for alkanes. Insights were also obtained concening the influence of the aqueous environment on the conformational equlibria for butane and pentane and concerning the relative solubility of branched and straight-chain alkanes. In particular, relatively weaker interactions between the alkanes and water molecules near methyl groups lead to the less exothermic heats of solution for more branched isomers. Acknowledgment. Gratitude is expressed to the National Institutes of Health and the National Science Foundation for support of this research. Registry No. H20, 7732-1 8-5; methane, 74-82-8; ethane, 74-84-0; propane, 74-98-6; butane, 106-97-8; isobutane, 75-28-5; pentane, 10966-0; neopentane, 463-82-1. (22) (a) For a review, see: Hertz, H. G . Angew. Chem. Inr. Ed. Engl.
(21) Jorgensen, W. L. Chem. Phys. Leff. 1980, 70, 326.
1970,9, 124. (b) Ruterjans, H. H.; Scherage, H. A. J . Chem. Phys. 1966, 45, 3296.
