JOURNAL O F T H E AMERICAN CHEMICAL SOCIETY 0 Copyright 1983 by the American Chemical Sociefy
March 23, 1983
Volume 105, Number 6
Solvation and Conformation of Methanol in Water’ William L. Jorgensen* and Jeffry D. Madura Contribution from the Department of Chemistry, Purdue University, West Lafayette, Indiana 47907. Received June 11, 1982
Abstract: Monte Carlo statistical mechanics simulations have been carried out for a dilute solution of methanol in water at 25 OC and 1 atm. The water-water interactions were described with the TIP4P potential, while the methanol-water interactions were represented by a modified TIPS potential that incorporates the methyl hydrogens explicitly. The latter feature allowed the internal rotation for the methyl group to be included in the simulation. Statistics for the conformational process were enhanced by umbrella sampling over chopped rotational barriers. The computed heat of solution (-17 f 5 kcal/mol) is in better accord with experiment (-10.8 kcal/mol) than earlier theoretical results. Furthermore, the preferred conformation of methanol in water is demonstrated to be staggered in contrast to the computed findings of Bolis, Corongiu, and Clementi. Detailed structural results have also been obtained and are discussed. Water molecules form a cage around the methyl group, while two to three water molecules are hydrogen bonded to the hydroxyl end of methanol. Further analyses of the hydrogen bonding reveal that the hydrogen bonds for the cage molecules have normal strengths. The exothermic heat of solution is due to favorable solute-solvent interactions and not enhanced solvent-solvent hydrogen bonding.
I. Introduction An important application for computer simulations of fluids is the study of the hydration and dynamics of biomolecules2 As a complement it is desirable to simulate aqueous solutions of monofunctional organic molecules such as alcohols, amines, carboxylic acids, and amides. Since experimental thermodynamic data, particularly heats and volumes of solution, are available for the hydration of these species, valuable insights could be obtained on the reliability of the theoretical methods and the intermolecular potential functions that are utilized. Such studies are also intrinsically important for obtaining a better understanding of the solvation organic molecules including solvent effects on conformational equilibria and for developing intermolecular potential functions. Nevertheless, there has been little theoretical work in this area. Although aqueous solutions of nonpolar solutes, especially methane,3 have been modeled, the only monofunctional organic systems that have been treated are f ~ r m a l d e h y d e , ~ m e t h a n ~ l ,and ~ , ~ethanol.’ The difficulties in obtaining viable intermolecular potential functions for these systems are apparent in the results for the alcohols. In two cases the computed energy of solution is much too In the other study the energy
of solution is not reported; however, evidence is presented which indicates methanol may strongly prefer to be eclipsed in water rather than staggered as it is in the gas phase.6 This would be a remarkable finding, particularly since experiment and traditional theory concur that solvent effects on conformational equilibria are usually slight unless there is a significant change in polarity between conformers.* However, hydrogen-bonded solvents can be anomalous and hydrophobic effects can cause conformational changes. In the present paper the results of our own initial simultations of the hydration of an organic molecule are presented. Methanol has again been chosen as the solute so comparison can be made with the previous ~ o r k . It~ is, ~clearly important to readdress the question of the preferred conformation of methanol in water. Another key issue is the ability of the potential functions that we use to yield reasonable thermodynamic results for the solution process. Along these lines we have devoted much effort to the creation of simple, transferable intermolecular potential functions (TIPS) suitable for fluid simulations. The TIPS yield good thermodynamic and structural results for a variety of pure liquids including ~ a t e r , ~ . ’ ~ alcohols,’3-1sethers,I6 and alkyl
(1) Quantum and Statistical Studies of Liquids. 25. (2) See, for example: (a) Rossky, P. J.; Karplus, M. J . Am. Chem. SOC. 1979,101, 1913. (b) Hagler, A. T.; Osguthorpe, D. J.; Robson, B. Science 1980,208, 599. (c) Clementi, E.; Corongiu, G . Biopolymers 1979, 18, 2431. (3) (a) Owicki, J. C.; Scheraga, H. A. J. Am. Chem. SOC.19177.99.7413. (b) Swaminathan, S.;Harrison, S. W.; Beveridge, D. L. Ibid. 1978, 100, 5705. (c) Okazaki, S.; Nakanishi, K.; Touhara, H.; Watanabe, N. J . Chem. Phys. 1981, 74, 5863. (d) Bolis, G.; Clementi, E. Chem. Phys. Lett. 1981.82, 147. (e) Rapaport, D. C.; Scheraga, H. A. J. Phys. Chem. 1982, 86, 873. (4) Mehrotra, P. K.; Beveridge, D. L. J . Am. Chem. SOC.1980, 102,4287. ( 5 ) Nakanishi, K.; Okazaki, S.;Ikari, K.; Touhara, H. Chem. Phys. Lett. 1981, 84, 428. (6) Bolis, G . ;Corongiu, G.;Clementi, E. Chem. Phys. Lett. 1982,86, 299. (7) Alagona, G.; Tani, A. Chem. Phys. Left. 1982,87, 337.
(8) For a review, see: Abraham, R. J.; Bretschneider, E. In “Internal Rotation in Molecules”; Orville-Thomas, W. J., Ed.; Wiley: London, 1974; Chapter 13. (9),Jorgensen, W. L. J . Chem. Phys. 1982, 77, 4156. The computed densities reported in this paper are too high by about 7%. This is the only significant correction needed. (!O) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. M., submitted for publication. (1 1) Jorgensen, W. L.; Binning, R. C.; Bigot, B. J . Am. Chem. SOC.1981, 103, 4393. (12) Jorgensen, W. L. J . Am. Chem. SOC.1981, 103, 4721. (13) Jorgensen, W. L. J . Am. Chem. SOC.1981, 103, 341. (14) Jorgensen, W. L. J. Am. Chem. SOC.1981, 103, 345. (15) Jorgensen, W. L.; Ibrahim, M. J. Am. Chem. SOC.1982, 104, 373.
0002-7863/83/1505-1407$01.50/0
0 1983 American Chemical Society
1408 J. Am. Chem.Soc.. Vol. 105,No. 6, 1983
Jorgensen and Madura
chlorides"~" under a broad range of temperatures and pressures. In particular, two similar potential functions, TIPS2 and TIP4P, have been developed for water.'JO The average errors in the energy and density for water at 25 "C and 1 atm are 1 and 7% for TII?Q9 and 3 and 0% for TIP4P;'O however, the TIPS2 potential yields a slightly better position for the first peak in the 00 radial distribution function. Several studies of dilute solutions have also been completed using the TIPS including simulations of sodium and methoxide ions in methanolI8 and atomic ions in water and tetrahydrof~ran.'~ The computed thermodynamic and structural results were again found to be in good agreement with available experimental data, though the statistical uncertainties for some computed properties are s u b ~ t a n t i a l . ~The ~ ~ present '~ study is a natural extension of this work and our previous investigations of pure methanol,13J5 water?JO and internal rotation in Besides the conformational results, detailed analyses of the thermodynamics and structure for methanol in water are presented in the following. 11. Statistical Mechanics Calculations (a) Monte Carlo Simulations. The Monte Carlo calculations were carried out for a sysem of one methanol molecule and 125 water molecules in a cube with periodic boundary conditions. The correspondence to an infinitely dilute solution is not exact in view of the additional images of the solute. However, the solute molecules only interact with the solute in the central cell, and edge effects do not appear to be significant based on the results of the hydrogen bonding analyses described below. The same boundary conditions have been found suitable for dilute solution simulations in our earlier studies18-20aand ~ t h e r s . ~ - ~ * ~ J The present simulation was performed in the isothermal-isobaric (NPT) ensemble at 25 O C and 1 atm. This ensemble is particularly appropriate for such studies because it yields directly heats and volumes of solution which may be compared with experimental data. The general details of the computational procedure have been presented p r e v i o ~ s l y ; ~however, ~ ~ ' ~ * ~some ~ specific points should be mentioned. First, in order to improve the statistics for the solute and its near neighbors, the Metropolis algorithm has been modified to include the preferential sampling procedure of Owicki.22 As in earlier ~ o r k , ' ~theJ ~probability of moving a solvent molecule is made proportional to l / ( r z c), where r is the distance from the solute to the solvent molecule and c is an adjustable constant. For the present case c was chosen to be 90 A2 which causes the nearest neighbors to be moved two to three times more often than the most distant solvent molecules. In addition, the sampling of the solute was enhanced by a factor of 3 over random sampling by attempting to move it every 40 configurations. Another special sampling procedure, umbrella sampling, was employed for the internal rotation of the s o l ~ t e . All ' ~ ~atoms ~~ in the solute are represented in the intermolecular potential functions as described in the next section. Consequently, the internal rotation of the methyl group can be included in the simulation using a threefold rotational potential, V(4)= Vo(1 cos 34)/2, where Vois taken as the experimental barrier height (1.07 k c a l / m ~ l ) . ~Umbrella ~ sampling permits the use of a surrogate rotational potential, V'(+),that is designed to allow more facile barrier crossing. Averages for a property 0 of the true system can then be obtained from eq 1 where ( )w indicates a configu-
+
+
(16) (a) Jorgensen, W. L.; Ibrahim, M. J. Am. Chem. SOC.1981, 103, 3976. (b) Chandrasekhar, J.; Jorgensen, W. L. J . Chem. Phys. 1982,77,5073. (17) Jorgensen, W. L.; Bigot, B. J . Phys. Chem. 1982, 86, 2867. (18) Jorgensen, W. L.; Bigot, B.; Chandrasekhar, J. J. Am. Chem. SOC. 1982, 104, 4584. (19) Chandrasekhar, J.; Jorgensen, W. L. J . Chem. Phys. 1982, 77, 5080.
Chandrasekhar, J.; Spellmeyer, D.; Jorgensen, W. L., to be submitted for publica tion. (20) (a) Bigot, B.; Jorgensen, W. L. J . Chem. Phys. 1981, 75, 1944. (b) Rebertus, D. W.; Berne, B. J.; Chandler, D.Ibid. 1979, 70, 3395. (21) Mezei, M.; Beveridge, D. L. J . Chem. Phys. 1981, 74, 6902. (22) Owicki, J. C. ACS Symp. Ser. 1978, No. 86, 159. (23) Lees, R. M.; Baker, J. G. J . Chem. Phys. 1968, 48, 5299.
Table I. TIPS Parameters for Water and Methanola ~
site 0 in H 2 0 ~b in H,O H in H,O 0 in CH,OH H g in CH,OH C in CH, OH HC in CH,OH
(7
0.000
-1.040 (0.520) -0.620
0.430
1 0 - 3 2~
c 2
600 0 0 560 0 1811 7
610
0 0 600
0
532 33 a Units: q , electrons;A2,kcal A"/mol; C2, kcal A6/mol. e 2 in eq 2 is 332.18 kcal A/mol. Values in parentheses determined by neutrality of the monomers. M is a point 0.15 A from oxygen toward the hydrogens on the bisector of the HOH angle. (-0.110) 0.100
rational average obtained by sampling over the surrogate system and w = exp(P(V(4) - V'(4)). We have found the procedure to be an effective way to enhance the sampling of the full configuration space and, therefore, convergence for systems with conformational degrees of freedom. Previous applications include Monte Carlo simulations of n-butane in a Lennard-Jones solvent and pure liquid 1,2-di~hloropropane.~~~~~ As in these cases, an appropriate choice for V'(4)can be made by chopping the barriers in V(4). So, for methanol the rotational barriers between staggered conformations have been chopped at 0.5 kcal/mol. The Monte Carlo simulation was then run with sampling over V'(4), and all averages were corrected according to eq 1. This combination of preferential and umbrella sampling provides a general means for the efficient simulation of dilute solutions involving a wide range of organic and biochemical solutes. Further enhancement may be obtainable by the inclusion of force-bias sampling; however, it is not clear that this provides any greater gain in convergence than additional preferential sampling for the same amount of computer time.24 The initial configuration for methanol in water was obtained by modifying a configuration from the N P T simulations of pure water.1° The latter calculations also provide a basis for comparison of the dilute solution results with the pure solvent.1° Equilibration involved over lOOOK configurations and an additional 2000K were used to obtain the final averages for methanol in water. Runs of this length have been found to be appropriate for similar simulations with preferential ~ a m p l i n g . ~ New ~ J ~ configurations ,~~ were generated by translating the selected monomer in all three Cartesian directions, by rotating it about one randomly selected axis, and for methanol by performing the internal rotation. The volume moves involved scaling all the intermolecular distances and were attempted on every 1000th configuration. The ranges for the motions were chosen to yield an acceptance rate of 4040% for new configurations. In addition, spherical cutoffs at 7.5 A were invoked in evaluating the intermolecular potential functions which included interactions with a monomer's 50-60 nearest neighbors. The computations were run on a Harris Corp. H-80 computer in our laboratory. (b) Intermolecular Potential Functions. The intermolecular interactions were described in the TIPS format. The TIP4P parameters were used for water as summarized in Table I.I0 Each water monomer is represented by four sites located at the three nuclei and a t a point M on the bisector of the HOH angle 0.15 A from oxygen toward the hydrogens. The bond length and angle are fixed at experimental values: r ( 0 H ) = 0.9572 A, LHOH = 104.52°.25 Interaction energies, e,,,,,, are determined by the intermolecular interactions between the sites including LennardJones and Coulomb terms (eq 2). For water, the q, A , and C
parameters were chosen to give good thermodynamic and struc(24) Mehrotra, P. K.; Mezei, M.; Beveridge, D. L. J. Chem. Phys., in press. (25) Benedict, W. S.;Gailar, N.; Plyler, E. K. J . Chem. Phys. 1956, 24, 1139.
Solvation and Conformation of Methanol in Water
J . Am. Chem. SOC.,Vol. 105, No. 6, 1983
Table 11. Results for Dimers and Complexesa
Table 111. Thermodynamic Results for the Hydration of Solutes at 25 "C and 1 atma
-Mr(O0) e (TIPS) (6-31G*)b complex 2.79 39 5.64 5.55 CH,OH.,.OH, 2.76 30 6.21 5.13 HOH...OHCH, 2.79 26 5.88 5.66 CH,OH...OHCH, cyclic (CH,OH), 2.84 46 4.00 2.15 46 6.24 5.64 HOH...OH, cyclic (H,O), 2.19 42 4.11 Units: r, A ; 0 , deg; AE,kcal/mol. Results are for linear complexes except as noted. Data from ref 27.
tural results for the liquid and reasonable descriptions of gas-phase dimersto It is noted that the only Lennard-Jones term is between oxygens for water-water interactions (Table I). Although TIPS parameters are available for alcohols, modification was required to allow explicit representation of the methyl hydrogens for methanol. Normally alkyl hydrogens are implicit in the TIPS so the internal rotation of a methyl group would not be considered. This approximation was found to be acceptable in a comparison of simulations of liquid methanol with the TIPS and with a similar potential called MHL (modified Hagler-Lifson) which includes the methyl hydrogens e ~ p l i c i t l y . ' ~The M H L function involves an f 9repulsion instead of r-12 so it could not be directly adopted in the present study. However, Lifson, Hagler, and Dauber (LHD) also reported Lennard-Jones 12-6 parameters and charges for alkyl carbons and hydrogens from their work on carboxylic acid and amide crystals.26 With these values and the TIPS for alcohols as a basis, a six-site model for methanol was developed with an interaction site on each atom. The experimental geometry for the monomer was again assumed; the OH, CH, and CO bond lengths are 0.945 1, 1.0936, and 1.4246 A, and the COH and H C O angles are 108.533 and 110.297°.23 The only parameters that were adjusted from the LHD or TIPS values were the charges for oxygen and hydroxyl hydrogen and the A for oxygen. These parameters were selected to yield reasonable interaction energies and geometries for methanol-water complexes and methanol dimers. The final values are recorded in Table I, and results of geometry optimizations for the complexes are shown in Table 11. The 0-0 distances in the hydrogenbonded complexes are all near 2.80 8,which is desirable for liquid simulation^,^-'^ though the best gas-phase data are 0.1-0.2 8, l ~ n g e r . ~As~ usual, , ~ ~ 0 is the angle between the HOX bisector of the hydrogen bond acceptor and the hydrogen bond vector for linear complexes, while it is the 0 4 - H angle for cyclic dimers. The optimized values are typical for TIPS and gas-phase res u l t ~ . * ' , The ~ ~ optimal interaction energies are also given in Table I1 and are compared with the ab initio results of Tse, Newton, and Allen using the 6-31G* basis set.27 These are the most sophisticated ab initio results that have been obtained so far for this series. The TIPS results compare favorably with the ab initio findings, though the complexes are somewhat more bound with the TIPS. This is again desirable for effective pair potentials in order to obtain good energies for hydrogen-bonded liquids.1° It should also be noted that the relative energies for the complexes with methanol are in the same order as the ab initio data. In particular, methanol is a somewhat better hydrogen bond acceptor than donor with water. 111. Results and Discussion (a) Thermodynamics and h r g y Distributions. The total energy of the dilute solution ( E T )is composed of solvent-solvent (Ess) and solute-solvent (Esx) terms and the intramolecular rotational energy of the solute (E,,,,,). The energy of solution of the solute from the ideal gas phase is then given by u s 0 1
= ESS + ESX+ Eintra - E&
- Ejfjtra
1409
(3)
solute property
CH,OH
ERX -_Eintra EtFitra
-18.0 f 0.3 0.33 0.03 0.32 -1269 f 2 -1211 * 3 2 i 5 -16 i 5 3864 i 1 1 3745 i 1 3 119 t 24 (6 3 . 5 ) d
EQQ E,&
Uss Msol V V* AVsoi MSOl
*
-1li5 (- 10.8)d
Na' -195 -1199 -1271 12 -122 3795 3145
H,OC
*1 i i
2 3
f
5
- 20 -1261 -1271
10
i 5 i8
-10
3175
13 50 i 21 (-11) -123 % 5 (-106)
3145 30.0 f 0 . 1 (30.0) - 1 0 . 8 f 0.0 (-10.5)
f
Energies and enthalpy in kcal/mol; volumes in A'. Experimental values in parentheses. Results from ref 19. Results from ref 10. Experimental data from ref 29. where E& is the solvent-solvent energy for the pure solvent and E:fftrais the intramolecular rotational energy for the solute in the ideal gas. Equation 3 can be rearranged to
AEsol = MSS+ AEintra + ESX
(4)
where the energy of solution is expressed ,as the sum of the solvent reorganization energy (AEss = Ess E,,), the change in intramolecular energy (AEintra= Eintra- E&a) and the solute-solvent energy. Similarly, the enthalpy of solution is the enthalpy difference between the solution and the pure liquid plus the solute in the ideal gas:
AH,,, = ET
+ PV - (E& + P P ) - (E$tra+ RT)
(5)
This can also be expressed as
AH,,I = AEs01+ PAV,,l
- RT
(6)
where the volume of solution, AV,,, is the difference between the volume of the solution (V) and of the pure liquid (P). The computed results for these quantities are shown in the second column of Table 111 for methanol in water. For comparison, the corresponding results for dissolving water in water and sodium ion in water from our previous studies with the TIP4P potential are also g , i ~ e n . l ~A, ' number ~ of points need to be noted. Some details are E#,,, for methanol is,obtained from a Boltzmann distribution over V(4)and Es and Ess include cutoff corrections of 0.1 kcal/mol per molecule for the Lennard-Jones interactions neglected beyond 7.5 8,. Also, for the solutions of methanol and sodium ion the computed heats and volumes of solution involve taking differences between large numbers, so the statistical uncertainty is substantial, particularly for AVsol.The error bars for the computed quantities are fu and were obtained via separate averages over each 50K configurations. These values are dependent on the size of the blocks; analyses of a long (4000K) NPT simulation of pure water indicate the computed uncertainties in Table I11 are probably within a factor of 2 of the true values.30 In view of the nature of the simulations and the error bars, the computed heats and volumes of solution are in reasonable agreement with the experimental values. For methanol the computed AHsol (-17 f 5 kcal/mol) is in better accord with the experimental data (-10.8 k c a l / m 0 1 ) ~than ~ ~ in the previous simulations of aqueous alcohols. For ethanol Alagona and Tani obtained a AHsol44 f 7 kcal/mol below the experimental value (-13 k ~ a l / m o l )while , ~ Nakanishi et al. only state their result for methanol is unsatisfactorily too n e g a t i ~ eand , ~ Bolis et al. provide no information on the point: The results are clearly very sensitive to the choice of potential functions. Many ~ t ~ d i e ~ have ~ ~ , ~ , ~ * ~ -
(26) Lifson, S.;Hagler, A. T.; Dauber, P. J. J . Am. Chem. SOC.1979,101, 5111.
(27) Tse, Y.-C.;Newton, M. D.; Allen, L. C. Chem. Phys. Lett. 1980, 75, 350.
(28) Dyke, T. R.; Muenter, J. S. J . Chem. Phys. 1974, 60, 2929.
~
~~
_ _ _ _ _ ~
(29) (a) Alexander, D. M.; Hill, D. J. T.Aust. J . Chem. 1969, 22, 347. (b) Jolicoeur, C.; Lacroix, G. Can. J . Chem. 1976, 54, 624. (30) Jorgensen, W. L., Chem. Phys. Left. 1982, 92, 405.
Jorgensen and Madura
1410 J. Am. Chem. SOC.,Vol. 105, No. 6, 1983 .18
ENERGY PAIR DISTRIBUTIONS
BONDING ENERGY DISTRIBUTIONS
.oo -30
-25
-20
-15
-10
-5
0
o i -8
BONDING ENERGY
Figure 1. Distribution of total solutmolvent bonding energies (kcal/mol) for methanol in water. Units for the ordinate are mole fraction per
kcal/mol. employed the MCY-CI potential for water which is known to underestimate the density of liquid water by 24%.31 Consequently, when fixed volume simulations consistent with experimental densities are performed with the MCY-CI potential, the system is under very high pressure. The ramifications of this effect warrant investigation. The computed AV,, for methanol (1 19 f 24 A3) overestimates the experimental value (63.5 A3),29bthough the discrepancy amounts to only 0.4 A3 per molecule. The difference might diminish in longer simulations; however, it would be impractical to test this much more since the error bars decrease only as l/dZ where M is the number of configurations and the simulations have already been run for 3-4000K configurations. On the other hand, the differences in the computed and experimental heats and volumes of solution may be real. Possible sources for error include the limited system size, boundary conditions, and deficiencies in the intermolecular potential functions including the neglect of specific three-body effects. Perusal of the energy components in Table 111 reveals that in all three cases the negative enthalpies of solution result from favorable solute-solvent interactions outweighing unfavorable solvent reorganization energies. A simple analysis can be given for the results for pure water: roughly two hydrogen bonds between solvent molecules must be broken (A& = +10 kcal/mol) in order to have sites for forming about four hydrogen bonds with the water solute (Esx = -20 kcal/mol). For sodium ion the solvent disruption is much greater (72 kcal/mol) due mainly to the reduced hydrogen bonding for the six waters in the first shell; however, it is dramatically offset by the strong ion-solvent attractions (-195 kcal/mol). The components for methanol in water are similar to those for water in water. Considering this and results of hydrogen bonding analyses presented below, it appears that the principal effect for methanol is the loss of one to two solventsolvent hydrogen bonds and the gain of two to three methanol-water hydrogen bonds. Overall these results for polar systems concur with the findings for formaldehyde and an alanine dipeptide in water which also have the solute-solvent attraction dominating the unfavorable solvent reorganization energyeZas4In contrast, simulations of the hydration of methane and other apolar solutes unananimously attribute the exothermic heats of solution predominantly to negative solvent reorganization energies or "structure making".3*3z It may be useful to refer to these two cases as solute dominant and solvent dominant hydration, respectively. It would be interesting to investigate how long an alkyl chain would be needed for the hydration of an alcohol to become solvent dominant. Various energy distributions were also obtained for the solute and solvent molecules. The distribution of total bonding energies for the methanol solute is shown in Figure 1. The solute experiences a continuum of energetic environments covering a 15kcal/mol range. Note that the average of the distribution is Esx. Greater insight into the environment of the methanol molecule can be obtained from the distribution of interaction energies between it and solvent molecules (Figure 2). As usual with (31) Owicki, J. C.; Scheraga, H. A. J. Am. Chem. SOC.1977, 99, 7403. (32) Pangali, C.; Rao, M.; Berne, B. J. J . Chem. Phys. 1979, 71, 2982.
r
I
-2
-'i
-6
0
2
Y
INTERRCTION ENERGY
Figure 2. Distribution of solutesolvent interaction energies (kcal/mol) for methanol in water. Units for the ordinate are number of solvent
molecules per kcal/mol. 15
BONDING ENERGY DISTRIBUTIONS
.oo -36
-31
-26
-21
-16
-11
-6
BONDING ENERGY
Distributions of total solvent-solvent bonding energies (kcal/mol) for methanol in water (dashed line) and pure liquid water (solid line). Units are the same as for Figure 1.
Figure 3.
ENERGY PRIR DISTRIBUTIONS
3
4i
O
-8
-6
-4
-2
0
2
INTERACTION ENERGY
Figure 4. Distributions of solventsolventinteraction energies (kcal/mol) for methanol in water (dashed line) and pure liquid water (solid line). Units are the same as for Figure 2.
hydrogen bonded systems, the distribution is bimodal. The band at low energy corresponds to the hydrogen bonded neighbors and the spike near 0 kcal/mol is due to the weak interactions with the many distant solvent molecules in the bulk. The minimum near -2.25 kcal/mol suggests an energetic criterion for hydrogen bonding; integration to that point indicates the methanol molecule participates in an average of 2.3 hydrogen bonds with solvent molecules. The corresponding distributions for the solvent molecules are shown in Figures 3 and 4. In both cases only water-water interactions are included. For comparison, the results for pure liquid waterlo are given by the solid curves in the figures. It is apparent that the introduction of one methanol molecule has little effect on the average energetics for the 125 water molecules. Integration of the energy pair distribution to -2.25 kcal/mol indicates each water molecule participates in an average of 3.6 water-water hydrogen bonds in both cases. It may also be noted that the greater smoothness in the solvent distributions than those for the solute is due to the statistical factor of having roughly 40 solvent moves for each solute move. (b) Structure. The eight possible solute-solvent radial distribution functions (rdfs) between different atoms are shown in Figures 5 and 6. In each instance the solid and dashed curves involve the oxygen and hydrogen in the water molecules, respectively. The solute-solvent hydrogen bonding is reflected in the first peaks of the rdfs in Figure 5. The fist peak in the 0-0 rdf occurs at 2.8 A (cf. Table 11) and its integral indicates there
Solvation and Conformation of Methanol in Water
J . Am. Chem. Soc., Vol. 105, No. 6, 1983 141 1 Table IV. Hydrogen Bonding Analyses for Aqueous MethanolQ
RRDIRL DISTRIBUTION FUNCTIONS
13-XW 0-HW
I
0
H-XW RRDIRL DISTRIBUTION FUNCTIONS H-HW
H-OW
0
I
’-,
1
2
/
I
3
I
I
5
tl
6
I
R
Figure 5. Solute-solvent radial distribution functions for methanol in water. The oxygen and hydrogen in water are designated OW and HW, while 0 and H refer to the oxygen and hydroxyl hydrogen in methanol. Distances are in A for all radial distribution functions. C-XH RRDIRL DISTRIBUTION FUNCTIONS 2
I
0
-
HC-XW
2i I 0 1 1
I
I
RRDIRL DISTRIBUTION FUNCTIONS
I
HC-OW
I
I
2
3
5
L)
6
I
R
Figure 6. Solute-solvent radial distribution functions. C and HC refer to the carbon and alkyl hydrogens of methanol. For other designations see Figure 5.
are 2.9 water molecules within 3.35 8,of the methanol’s oxygen. The second band in this distribution out to about 5 8, can be assigned to water molecules in the first layer nearer the methyl group and to the second shell around the hydroxyl group. Turning to the 0-H, (methanol @water H) rdf, a small f m t peak is found from 1.6 to 2.6 8,. Its integral is 1.9 which can be assigned to hydrogen bond donating waters. The second peak centered near 3.3 8, should contain the other hydrogens of the donors and hydrogens of waters acting as hydrogen bond acceptors. The latter species are highlighted in the first peak of the H-O, rdf in Figure 5. This band ranges from 1.6 to 2.7 8, and contains 1.0 water molecule.33 Thus, from this analysis it appears that on the average the methanol molecule could be in three hydrogen bonds, one as donor and two as acceptor. Figure 6 illustrates the rdfs for the methyl end of the solute. The peaks are broader, and less distinct structure can be recognized than near the hydroxyl group. The most interesting point is the similarity of the C-0, rdf to the same rdf for methane in water.3b In both cases the first peak ranges from about 3-5.5 8, and integrates to ca. 20 water molecules. By analogy to the results for methane3 this implies a cage exists around the methyl group of methanol. However, on the basis of the lower height of the first peak in the (2-0, rdf for methanol, the cage is not as distinctly defined as in the simulations of methane.3b Clathrate formation was also found in the previous theoretical studies of methanol in water.ss6 Furthermore, the caging is apparent in stereoplots of configurations from the present simulation. The examples in (33) The fact that the first peak in the 0-H, rdf integrates to twice the value of the first peak in the H-O, rdf is due to the normalization of the rdfs which takes account of the density of water hydrogens being twice that of oxygens.
-5.56 1.37 -4.19 0-3.5 2.92 3.5-4.5 3.41 -4.18 -5.58 1.40 -5.52 1.35 -4.16 4.5-5.5 3.52 -5.49 1.34 3.55 -4.15 5.5-6.5 -5.52 1.36 6.5-7.5 3.60 -4.17 -5.53 1.35 3.51 -4.18 1.5-8.5 -5.52 1.35 -4.17 oureH-0 3.57 E’Sin kcal/mol. (E(HB)) is the average hydrogen bond energy which can be decomposed into coulomb (E(cou1omb))and Lennard-Jones (E(LJ))contributions. Only water-water hydrogen bonds which are defined by an interaction energy of -2.25 kcali mol or less are included in the analyses. r ( O 0 ) is the methanoloxygen to water-oxygen distance. Figure 7 only include the solvent molecules within 5.5 A of the methyl carbon for clarity. Two different configurations are displayed; the upper one gives a better view of the region near the methyl group and the lower one focuses on the hydrogen bonding with the hydroxyl group. It is clear from the plots that the cage is flexible in view of its irregularity and the variety of hydrogen bond geometries. In addition, the hydrogen bonding to the solute blends nicely into the cage structure. One final observation is that hydrogens from water in the cage do not typically point in toward the methyl group. Instead 0-H bonds prefer to lie parallel to the surface of the cage. This effect is also reflected in the C-0, and C-H, rdfs (Figure 6) in which the major band ranges from 3 to 5.5 8, for both distributions. The contribution from 2 to 3 A in the C-H, rdf integrates to two hydrogens which can be attributed to the hydrogen bond donating waters near the hydroxyl group. The solvent-solvent rdfs are compared with the pure water results in Figure 8. As found before for the energy distributions (Figures 3-4), it is evident that the environment for the majority of the 125 solvent molecules is very similar to that in pure water. Since analyses of the structure of liquid water have been presented p r e v i o ~ s l y , ~the J ~discussion will not be repeated here. However, one point worth noting is the basic similarity of the prediction for the 0-0 rdf and the X-ray findings, though as with most two-body potential functions the first peak is too high.’O~~~ (c) Hydrogen Bonding Analyses. In order to obtain further details on the effect of the solute on the solvent’s structure, analyses of the water-water hydrogen bonding were carried out. An energetic criterion was used to define a hydrogen bond based on the location of the minimum in the energy pair distribution. Namely, any pair of molecules with an interaction energy of -2.25 kcal/mol or less is considered to be hydrogen bonded. The hydrogen bonding for the water molecules in shells around the methanol was then analyzed. The results are summarized in Table IV where the shells are determined by the separation of the oxygens of methanol and water. The first shell from 0 to 3.5 8, contains 3.4 water molecules which participate in an average of 2.9 water-water hydrogen bonds. These molecules are also responsible for the 2.3 hydrogen bonds with methanol discussed near the end of section 1II.a above. Consequently, they participate in a total of 3.6 hydrogen bonds each, which is exactly the same as for monomers in pure water. The average hydrogen bond energies and their Coulomb and Lennard-Jones components are also given in Table IV. The hydrogen bond analyses are made from configurations saved at 5K intervals during the simulation. By analysis of different sets of configurations the error bars on the hydrogen bond numbers and energies are estimated as h0.02 and i0.03 kcal/mol. Therefore, the average hydrogen bond energies for all regions in the dilute solution are the same as for pure liquid water. The average number of hydrogen bonds is also very close to the bulk value throughout. The only exceptions are the 0-3.5-A region discussed above and possibly the 3.5-5.5A shell which has a (34) Narten, A. H.; Levy, H. A. J . Chem. Phys. 1971, 55, 2263
1412 J. Am. Chem. Soc., Vol. 105, No. 6,1983
Jorgensen and Madura
r' I
'1-1
n
I
/
Figure 7. Stereoplots of two configurations from the simulation of methanol in water. Only water molecules within 5.5 A of the carbon in methanol are displayed. DIHEDRAL RNGLE DISTRIBUTIONS
00 RRDIRL DISTRIBUTION FUNCTIONS 3
3
- PURE WATER
-- --
G 2
SOLUTION X-RFIT
1 I
0
I
I
I
I
I
OH RRDIRL DISTRIBUTION FUNCTIONS
W I
'10
60
80
100
120
(DEG.1
Figure 9. Dihedral angle distributionsfor methanol in water (solid line) and in the ideal gas (dashed line). Staggered conformations have $ , = 60, 180, and 300O; eclipsed conformations have C#J = 0, 120, and 240O. Units for the ordinate are mole fraction per degree X lo-'. I
I
I
I
I
HH RRDIRL DISTRIBUTION FUNCTIONS
G 1 0 2
20
PHI
G l 0
0
3
5
L)
6
I
R
Figure 8. The solventsolvent radial distribution functions for methanol in water (dashed lines) and pure liquid water (solid lines). The X-ray
results of Narten and Levy (ref 34) for pure liquid water at 25 OC and 1 atm are also shown for the 00 distribution (shortest dashes). slightly lower average number of hydrogen bonds than the bulk. The latter area includes the cage around the methyl group which in simulations of aqueous methane shows enhanced solventsolvent b ~ n d i n g .This ~ region was specifically probed by analyzing the water molecules within 4.5 A of the methyl carbon excluding the molecules within 3.5 A of the oxygen in methanol. Indeed, these water molecules have a slightly lower average number of hydrogen bonds (3.39) than bulk water (3.57), but their average hydrogen bond strength (-4.15 kcal/mol) is the same as the bulk value within the statistical limits. Although stronger hydrogen bonding is not apparent for the cage molecules, what is remarkable is the ability of the cage molecules to participate in a near normal
number of hydrogen bonds even though they border an excluded region. Similar observations were made by Rossky and Karplus for water molecules in the nonpolar regions surrounding an alanine dipeptide.2a (a) Internal Rotation. By including the internal rotation in the simulation and with the umbrella sampling, good statistics could be obtained for the conformational problem. The computed intramolecular energy of 0.33 f 0.03 kcal/mol (Table 111) is the same as the ideal gas value (0.32 kcal/mol) obtained from a Boltzmann distribution for V ( 4 ) . Thus, methanol prefers to be staggered (4 = 60°) rather than eclipsed (4 = Oo) in aqueous solution just as it does in the ideal gas. This is confirmed by the computed distributions for the dihedral angle shown in Figure 9. The differences are not pronounced, so the fundamental conclusion is that there is essentially no solvent effect on the conformational preference of methanol. This is not surprising in view of our previous identical finding for pure methanolI3and analogous results for methyl rotations in alanine dipeptide.2a Furthermore, since the methyl group is not in close contact with the solvent because of the caging (Figure 7), the intramolecular potential should be relatively unperturbed. As mentioned in the Introduction, solvent effects on conformational equilibria are not normally expected to be large unless there are significant differences in polarity between conformers.8 For example, experiments and simulations
J. A m . Chem. SOC.1983, 105, 1413-1419 with the TIPS find a dramatic increase in the gauche-trans ratio for 1,2-dichloroethane in going from the gas phase to the neat liquid." In this instance the more polar gauche conformer (p = 2.6 D) is preferentially stabilized in polar media in comparison to the trans conformer ( p = 0). However, some folding of nalkanes in water has been predicted35 and found in Monte Carlo calculation^.^^*^^ For n-butane, the theory of Pratt and Chandler and a simulation both yield a ca. 20% increase in the gauche population upon transfer from the ideal gas phase to aqueous sol~tion.~~J~ A few comments can be made on the discrepancy between the present findings and those of Bolis et aL6 All the evidence from their simulations indicates methanol prefers to be eclipsed in water; however, they are aware of the limitations of their computations and say the conformational preference is "an open problem". There are several specific aspects of their study that provide possible sources of error. Again the MCY-CI potential has been employed for the water interactions, and physically unusual boundary conditions were used: the 198 solvent molecules plus the solute were constrained to a sphere with a constant volume commensurate with the density of liquid water. In addition, the water-methanol potential was obtained by fitting to modified ab initio results with a minimal basis set.6 From Figure 4 in their paper it is apparent that staggered methanol prefers to be a hydrogen bond donor rather than acceptor with water by about 3 kcal/mol. With minimal basis set calculations this difference is normally only ca. 1 kcal/mol and with larger basis sets the preference reverses (Table II).27 The authors also indicate the eclipsed form is a better hydrogen bond acceptor than the staggered by about 1 kcal/mol. This is surprisingly large for as remote an effect; the present methanol-water potential predicts the conformation of methanol to make less than a 0.1 kcal/mol difference on any hydrogen bond energy. Also, our own experience has been that it is difficult to obtain good intermolecular potential functions for fluid simulations by adding dispersion and other corrections to minimal basis set results.38 As discussed elsewhere,
1413
even high quality ab initio calculations are not necessarily particularly useful in this regard.39 In closing this section some final technical details may be noted. There were 862 barrier crossings between different staggered conformers during the final 2000K configurations. Many of the crossings reverted quickly; the number of transitions such that the new well was explored extensively was about 60. The methyl group made several complete circuits, though the tagged hydrogen spent most of its time in the wells for two of the three identical conformers. In this case the umbrella sampling increased the number of barrier crossings by a factor of ca. 2-3. It is apparent that for solutes wih higher rotational barriers umbrella sampling would be very valuable for obtaining adequate sampling of the conformational space.20
IV. Conclusion The present work demonstrates the utility of Monte Carlo simulations in the NPT ensemble with preferential and umbrella sampling for modeling dilute solutions containing organic solutes. This appears to be as efficient an approach as is currently available that includes the solvent molecules explicitly. Many detailed structural insights were obtained for methanol in water. The methyl group is surrounded by a flexible cage of water molecules that blends in to accomodate two to three methanol-water hydrogen bonds. The energies and numbers of hydrogen bonds for the cage molecules are remarkably near normal for bulk water. It was also found that the preferred conformation of methanol in water is staggered as in the gas phase. The reasonable thermodynamic results obtained with the TIPS support the validity of the observations. Acknowledgment. Gratitude is expressed to the National Science Foundation (CHE80-20466) for support of this study. The authors are also grateful to Professor B. J. Berne for a preprint of ref 36. Dr. J. Chandrasekhar and Professor M. C. R. Symons provided helpful discussions.
Registry No. Methanol, 67-56-1. (35) Pratt, L. R.; Chandler, D. J . Chem. Phys. 1977, 67, 3683. (36) Rosenberg, R. 0.; Mikkilineni, R.; Berne, B. J. J. Am. Chem. SOC. 1982, 104, 7647. (37) Jorgensen, W. L. J . Chem. Phys. 1982, 77, 5757.
(38) Jorgensen, W. L. J . Am. Chem. SOC.1980, 102, 543. (39) Jorgensen, W. L. J . Chem. Phys. 1981, 75, 2026. .
Absorption Spectra and Photochemical Rearrangements of Toluene, Cycloheptatriene, and Norbornadiene Cations to Methylenecyclohexadiene Cation in Solid Argon Benuel J. Kelsall* and Lester Andrews Contribution from the Department of Chemistry, University of Virginia, Charlottesuille, Virginia 22901. Received August 20, 1982 Abstract: Toluene, cycloheptatriene, and norbornadiene cations have been produced and isolated by matrix photoionization
methods and their absorption spectra recorded. Absorptions have been identified in agreement with photodissociation and photoelectron spectra. In addition, the photochemistry of these cations has been investigated, showing that they all rearrange into a common methylenecyclohexadiene cation, the McLafferty rearrangement product from gas-phase mass spectrometry studies.
Introduction
The toluene cation and its structural isomers are among the most thoroughly studied cations in mass Of (1) Bursey, J. T.; Bursey, M. M.; Kingston, D. F. I. Chem. Rev. 1973, 73, 191. (2) Dunbar, R. C. J. Am. Chem. Soc. 1973, 95, 472. (3) Baldwin, M. A,; McLafferty, F. W.; Jerina, D. M. J . Am. Chem. Soc. 1975, 97, 6169.
0002-7863/83/1505-1413$01.50/0
particular interest are the rearrangements among the C7HBf. structural isomers that precede decomposition. McLafferty et. a1 have proposed four rearrangement pathways that account for the observed isotopic ~crambling.~ A common theme from these (4) Jackson, J.-A. A.; Lias, S. G.; Ausloos, P. J . Am. Chem. SOC.1977, 99, 7515. (5) McLafferty, F. W.; Bockhoff, F. M. J . Am. Chem. SOC.1979, 101, 1783.
0 1983 American Chemical Society
