J . Phys. Chem. 1987, 91, 6083-6085
6083
Use of Statistical Perturbation Theory for Computing Solvent Effects on Molecular Conformation. Butane in Water William L. Jorgensen* and J. Kathleen Buckner Department of Chemistry, Purdue University, West Lafayette, Indiana 47907 (Received: August 26, 1987)
Statistical perturbation theory has been applied in Monte Carlo simulationsto examine the effect of hydration on the torsional energy surface for butane. The methodology is demonstrated to yield results with high precision and to have significant advantages over umbrella sampling. Consistent with traditional notions about hydrophobic effects on molecular conformation and earlier theoretical results, the population of gauche conformers for butane is found to increase by 12% upon transfer from the gas phase to aqueous solution at 25 OC.
Introduction The past 2 years have witnessed a surge in applications of statistical perturbation theory to compute free energy changes for transformations in solution. The activity was much stimulated by recognition of the potential application of the methodology to calculate relative binding affinities' and by demonstration of the high precision that could be obtained in computing relative free energies of hydration for organic solutes.2 Since that time the list of applications has expanded to include prediction of relative free energies of binding for anions to a cryptand3 and for substrates to enzymes,"-5free of hydration for protein side chains, noble gases, hydrocarbons, and small ions,&I0 the effect of hydration on the structure of an SN2 transition state," free energy profiles for separating ion pairs in solution,I2and pK,'s for weak organic acids in water.I3 The methodology has great versatility so many additional areas of application will emerge. In the present paper, the ability of such computations to provide precise results for the effect of solvation on conformational equilibria is demonstrated. The prototypical system of butane in water was chosen for this initial study. Traditional ideas about hydrophobic effects and polymer folding suggest that the more compact gauche conformer should be relatively favored over the trans form in water compared to the gas phase. Prior results from integral equation t h e ~ r i e s ' ~ , ' ~ and Monte Carlo simulation^'^'^ have indeed predicted increases in the gauche population of 12-35% upon transfer to water at 25 OC. Nevertheless, there remains concern about these findings in view of their scatter, the approximations in the integral equation (1) Tembe, B. L.; McCammon, J. A. Compur. Chem. 1984,8,281. (2) Jorgensen, W. K.; Ravimohan, C. J . Chem. Phys. 1985, 83, 3050. (3) Lybrand, T:P.; McCammon, J. A.; Wipff, G. Proc. Natl. Acad. Sci. U.S.A. 1986,83, 833. (4) Wong, C. F.; McCammon, J. A. J. A m . Chem. SOC.1986,108,3830. ( 5 ) Bash, P. A.; Singh, U. C.; Brown, F. K.; Langridge, R.; Kollman, P. A. Science 1987, 235, 574. (6)Bash, P. A.; Singh, U. C.; Langridge, R.; Kollman, P. A.Science 1987, 236, 264. (7) Straatsma, T.P.; Berendsen, H. J. C.; Postma, J. P. M. J. Chem. Phys. 1986,85, 6720. (8) Fleischman, S. H.; Brooks, C. L., 111 J . Chem. Phys., in press. (9) Lybrand, T.P.; Ghosh, I.; McCammon, J. A. J. Am. Chem. Soc. 1985, 107,7793. (10)Singh, U.C.; Brown, F. K.; Bash, P. A,; Kollman, P. A. J. Am. Chem. SOC.1987,109, 1607. (11) Jorgensen, W. L.; Buckner, J. K. J . Phys. Chem. 1986,90, 4651.
(12)Jorgensen, W. L.; Buckner, J. K.; Huston, S. E.; Rossky, P. J. J . Am. Chem. SOC.1987, 109, 1891. (13)Jorgensen, W. L.; Briggs, J. M.; Gao, J. J. Am. Chem. SOC.,in press. (14)Pratt, L. R.;Chandler, D. J . Chem. Phys. 1977,67, 3683. (15) Zichi, D.;Rossky, P. J. J . Chem. Phys. 1986,84, 1712. (16)Rosenberg, R. 0.; Mikkilineni, R.; Berne, B. J. J . Am. Chem. SOC. 1982, 104, 7647. (17)Jorgensen, W. L. J . Chem. Phys. 1982,77, 5757. (18)Jorgensen, W.L. J . Phys. Chem. 1983,87, 5304. (19)Jorgensen, W.L.;Gao, J.; Ravimohan, C. J . Phys. Chem. 1985.89, 3470.
TABLE I: Computed Free Energy Differences and Potential of Mean Force for Butane in Water" AAGaq 9, a, i-j j d i AGaq(ai) AEgas(9j) W(0,) 4.73 4.18 0 15 -0.01 f 0.02 -0.01 f 0.03 -0.55 4.21 3.66 15 30 0.09 i 0.03 -0.06 f 0.03 -0.55 2.94 2.47 30 45 0.11 i 0.03 -0.07 f 0.04 -0.47 1.61 1.23 45 60 0.01 i 0.04 -0.10 f 0.02 -0.38 0.91 0.59 60 75 0.08 i 0.02 -0.02 f 0.04 -0.32 1.13 0.85 75 90 0.01 f 0.04 -0.09 f 0.04 -0.28 2.05 1.82 90 105 0.03 f 0.04 -0.15 f 0.03 -0.23 3.01 2.87 105 120 0.11 f 0.04 -0.07 i 0.04 -0.14 3.35 3.30 120 135 0.04 f 0.03 -0.07 i 0.03 -0.05 2.80 2.81 0.01 135 150 0.10 f 0.03 0.04 f 0.02 1.63 1.67 0.04 150 165 -0.02 f 0.02 0.06 f 0.04 0.00 0.47 0.47 165 180 -0.04 i 0.04 -0.04 f 0.03 0.00 0.00 0.00 180
'Energies in kcal/mol; angles in deg. AAGaq is the incremental difference in free energy of hydration, AGaq is the cumulative free energy of hydration relative to 0 = 180°, @gas is the gas-phase torsional energy from eq 2, and Wis the potential of mean force in water (AGaq + AE*as). Error bars for AAGaq are i l a and were obtained from separate averages over blocks of lo5 configurations. theories, and the convergence of the simulations.18s20 Though umbrella sampling2' was used to promote the efficient sampling of the conformational space in the simulations, proof of convergence is difficult and the choice of surrogate (umbrella) potentials is an undesirable variable. These problems can be avoided with the perturbation approach as illustrated here in the context of providing a definitive result for the conformational shift of a model for butane in water.
Computational Procedure A series of Monte Carlo statistical mechanics simulations were carried out for a system of 1 butane molecule plus 216 water molecules in a cubic cell with periodic boundary conditions. The details for the individual calculations were identical with those described for the previous s i m ~ l a t i o nincluding '~ the use of Metropolis and preferential sampling. The TIP4P model was employed for the water-water interactions,22while the OPLS parameters were used for butane which is represented as four united atoms.23 The butane-water interactions in this model consist solely of the Lennard-Jones interactions between oxygen and the (20) Berne, B. J.; et al., to be published. (21)Valleau, J. P.; Torrie, G. In Statistical Mechanics, Part A; Berne, B. J., Ed.; Plenum: New York, 1977;p 169. (22)Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. J . Chem. Phys. 1983, 79,926. Jorgensen, W. L.; Madura, J. D. Mol. Phys. 1985,56, 1381. Tse, J. S.;Klein, M. L. J . Phys. Chem., in press.
(23) Jorgensen, W.L.;Madura, J. D.; Swenson, C. J. J . Am. Chem. SOC. 1984,106, 6638.
0022-3654/87/2091-6083$01.50/00 1987 American Chemical Society
6084
The Journal of Physical Chemistry, Vol. 91, No. 24, 1987 .q
Letters 5 .O
Effect of Solvation on the Rotational Potential Function
Functions
Rotational Potential
A
,"'
FO+". -
-. --+/
--'.I
St
1
30
0
60
90
.o
-1
1
120
150
I 0
180
60
Figure 1. Computed change in the free energy of hydration, AGaq, for butane as a function of the ClCZC3C4 dihedral angle. The two curves are for the perturbations from 0' to 180' (solid) and from 180' to 0' (dashed), both zeroed at 180'.
united CH2 and CH3 groups. Spherical cutoffs were imposed at 8 and 9 A for the water-water and butane-water interactions, respectively, based roughly on the center-of-mass separations. Statistical perturbation theory provided the change in free energy as the ClC2C3C4 dihedral angle, @, was varied in increments, A@, of 15' from @ = 0' (cis) to @ = 180' (trans). Specifically, the free energy change for mutating to aj = @i A@ is given by eq 1,24where 0 = l / k T , Ej - Ei is the energy Gj - Gi = k T In (exp[-/3(Ej (1)
120
-
-
Results and Discussion The incremental free energy changes, AAGaq, and their total, AGaq, are recorded in Table I along with the gas-phase torsional energy, AEgas, and the net potential of mean force, W, for the rotation in water. A E g " has been computed from eq 2 where the AEgas(@) =
Y2V,(1
+ cos @) + Y2V2(1 -cos
29)
+ l/ZV3(1 + cos 3@) (2)
(24) (a) Zwanzig, R. W. J . Chem. Phys. 1954,22, 1420. (b) Bennett, C H.J . Comput. Phys. 1976, 22, 245.
300
360
Figure 2. Torsional energy functions for butane. Ai?(@) is the gas-phase potential energy from eq 2, AG(@) is the change in free energy of hydration from Figure 1, and W(@) is their sum, the potential of mean force in water.
.20.
+
difference between the systems with @ j and ai, and the average Naturally, if the reference, ai, and is based on sampling for @i. perturbed, aj,systems are too disparate, convergence of the average in eq 1 will be slow. As shown below, choosing A@ = 15' was not problematic in the present case. It should also be noted that Ej - Ei involvesjust the difference in solutesolvent energies2 The simulations were run with the BOSS program on Gould 32/8750 and Microvax I1 computers in our laboratory. Operationally, the program monitors three solutes simultaneously, the reference one with and two perturbed solutes with aj = aif A@. Thus, twice A@ or 30' could be covered in one simulation. It should be noted that the perturbation was achieved by moving both CH3 groups by 7.5' rather than just one by 15' in order to minimize disturbance of the solvent. Initially, simulations were run for Oi= 15', 45', 75', 105', 135', and 165'. Then to check each computed free energy change forwards and backwards, Le. Oi aj and aj @i, simulations were run for ai= Oo, 30°, 60°, 90°, 120°, 150°, and 180'. This should generally not be necessary once the choice of A@ has been supported by running some increments in both direction. Each simulation consisted of an equilibration phase for 8 X lo5 configurations followed by averaging for 2 X lo6 configurations. Thus, a total of 3.6 X lo7 configurations were utilized, which is roughly an order of magnitude greater than in the earlier s t ~ d i e s . ' ~ -The ' ~ present calculations were run in the isothermal-isobaric ensemble at 25 ' C and 1 atm.
2110
180
Ph I Idegrees I
Phl [degrees1
I! I 1 I I
*
I 1 I 1 I 1
15-
-Gas
Phase
I
.lo.
Rqueous Phase
\, .OS
.oo 0
60
120 Ph I
180
2qO
300
:
( degrees 1
Figure 3. Population distributions for the dihedral angle about the central CC bond of butane in the gas phase and in aqueous solution at 25 'C. Units for the ordinate are mole fraction per degree.
constants were obtained by fitting to results of MM2 molecular mechanics calculation^.'^ Specifically, Vl, V2,and V, are 1.522, -0.315, and 3.207 kcal/mol. Though the MM2 results are reasonable,25the exact form of AEgaSis not critical in the present context where the emphasis is on AGa% Furthermore, the calculations have not considered intramolecular vibrations whose effect on Wshould be similar for the gas phase and solution.18s26 The simulations provided two independent evaluations of the solvent effect on the torsional energy surface from the perturbations for @i ejand @,. The average of the two results is listed as AGaq in Table I. The consistency of the results is clear from the individual AAGaq values in Table I and from their sums i series that are illustrated in Figure 1. for the i j and j Many of the AAGaq are near zero; however, the trend from the two series is unequivocal. The free energy of hydration for butane does indeed decrease in proceeding from trans through gauche (@ = 60') to cis; the gauche conformer has a lower AGaq than trans by 0.3 kcal/mol, and the cis form is better hydrated by an additional 0.2-0.3 kcal/mol. The effect on the torsional energy surface is illustrated in Figure 2; the deepening of the gauche well
-
- -
-
( 2 5 ) (a) Raghavachari, K. J . Chem. Phys. 1984,81, 1383. (b) Kanesaka, I.; Snyder, R. G.; Straws, H. L. Ibid. 1986, 84, 395. (26) Bell, C. D.; Harvey, S. C. J . Phys. Chem. 1986, 90, 6595.
J. Phys. Chem. 1987, 91, 6085-6087 and lowering of the cis barrier by hydration are apparent. The change in the torsional energy profile leads in turn to the shift in the distribution of conformer populations, S(@),shown in Figure 3. For the gas phase, S o ( @ )is obtained from a Boltzmann distribution for AEgaS(@,), . . eq- 3, while in water eq 4 "(@) = exp(-flAEgas(@))/JeXp(-flAEgas(Q.)) dQ. (3)
applies where
w(@)now replam p(@).
The two distributions
'(@) = e x p ( - B w ( @ ) ) ~ ~ e x p ( - P W O )
(4)
are interrelated by eq 5, where AGaq(Q.)is again the solvent effect
S(9) = cSo(Q.)exp(-/3AGaq(@))
(5)
illustrated in Figure 1 and c is a normalization constant. The increase in the populations for the gauche conformers and the decrease in the trans population upon hydration are apparent in Figure 3. Quantitatively, integration reveals the shift to be 12.3% from 68.2% trans in the gas phase to 55.9% in water. The uncertainity in this result is small, 1-2%, based on the difference i series. in the i j and j The computed shift is in accord with our previous result of 14% obtained via umbrella sampling with the same potential f~nctions.'~ The latter procedure entails a simulation in which 9 is allowed to vary during the sampling like any other degree of freedom. Recording the Occurrence of the different values of @ then directly yields S(@).Umbrella sampling comes in through the use of a surrogate torsional potential during the simulation in which the rotational barriers have been reduced. This permits more frequent transitions between conformers and, therefore, faster convergence for the conformer populations. The influence of the surrogate potential can then be removed analytically to yield S(@)for the true torsional potential. Though the accord for the results for butane in water seems to support the umbrella sampling method, there are several important advantages to the present perturbation procedure. Foremost, one is always concerned with umbrella sampling that a bottleneck may occur that does not allow proper sampling of the conformational space. This is very difficult to diagnose, and the only obvious, potential remedy is to run multiple, long simulations from different starting points that, however, could
- -
6085
all still encounter the same bottleneck. The choice of the optimal surrogate potential is also unclear. Any remaining torsional barriers promote the likelihood of bottlenecks, while their total removal may slow convergence by excessive sampling in regions with small S(9).Based on the present results, it seems likely that a bottleneck occurred in our simulation of pentane in water that used umbrella ~amp1ing.l~Though the averaging covered 5 X lo6 configurations, a 14.5% increase in the trans-trans conformation was computed for transfer from the gas phase to water. The problem is being reinvestigated with the present method. Another disadvantage with umbrella sampling is that all computed properties other than S(@)correspond to averages over a range of @. For ease of analysis and interpretation, it is more advantageous to have correct averages for the other properties at discrete values of @. In the perturbation procedure a surrogate potential is not used; in fact, the form of the gas-phase torsional potential is itself irrelevant since it is not evaluated during the simulations. There is also no concern about evenly covering the full conformational space since this is done by force through the sequential perturbations. Thus, the potential for the occurrence of conformational bottlenecks is eliminated. Convergence can also be readily verified by running the perturbations in both directions as shown here. The only remaining issue is the choice of A@. The present results provide some guide on this, though for new systems a few increments should be checked by perturbing in both directions. Overall, the perturbation method appears to be the most straightforward procedure for obtaining precise results on solvent-induced shifts in conformer distributions. It can be used in conjunction with either Monte Carlo or molecular dynamics simulations for which direct sampling methods and even umbrella sampling may be too slowly convergent.'8s20 Application to the present model of butane in water reveals enhancement of the gauche population in water consistent with earlier results. The microscopic origin of the shift has also been considered previo~sly.'~,''It is important to note that the effect is not just due to the fluid environment, since no conformational shift is found for pure liquid butane.18-23 Acknowledgment. Gratitude is expressed to the National Institutes of Health (GM-32136) for support of this research.
Nonunity I"Quantum Yield In the 193-nm Photodissociation of Methyl Iodide Wayne P. Hess, Ron Naaman,+ and Stephen R. Leone*$ Joint Institute for Laboratory Astrophysics, National Bureau of Standards and University of Colorado, and Department of Chemistry and Biochemistry, University of Colorado, Boulder, Colorado 80309-0440 (Received: September 8, 1987)
The I* quantum yield of CHJ at 193 nm is determined to be 70 f 4% by the diode laser gain versus absorption technique. This value is lower than the 100% I* yield determined from the analysis of molecular beam time-of-flight results. A critical evaluation of the gain versus absorption technique is presented and the possible sources of the discrepancy between the two results are discussed.
Recently, we attempted to investigate the effects of dimer formation' upon the photodissociation I* quantum yield of methyl iodide. In the course of this work we found it necessary to verify the I* quantum yield of CH31 at 193 nm (ArF excimer laser). In contrast to the 100% I* quantum yield reported by Van Veen et al. who used the molecular beam time-of-flight technique,2 we 1986-87 Visiting Fellow, Joint Institute for Laboratory Astrophysics. Permanent address: Department of Isotope Research, Weizman Institute of Science, Rehovot, 76 100, Israel. *Staff member, Quantum Physics Division, National Bureau of Standards.
0022-3654/87/209 1-6085$0 1.50/0
obtain an I* quantum yield of 70 f 4% by the diode laser gain versus absorption t e c h n i q ~ e . ~In . ~ this Letter we describe our results, evaluate the errors in the method, and consider possible sources of the discrepancy. (1) Donaldson, D. J.; Naaman, R.; Vaida, V. J . Chem. Phys. 1987, 87, 2522. (2) Van Veen, G. N. A.; Baller, T.; De Vries, A. E. Chem. Phys. 1985, 97, 179. (3) Hess, W. P.; Kohler, S. J.; Haugen, H. K.; Leone, S.R. J . Chem. Phys. 1986, 84, 2143. (4) Hess, W. P.; Leone, S.R. J . Chem. Phys. 1987.86, 3773.
0 1987 American Chemical Society
