13077
J. Phys. Chem. 1994, 98, 13077-13082
Free Energies of Hydration and Pure Liquid Properties of Hydrocarbons from the OPLS All-Atom Model George Kaminski, Erin M. Duffy, Tooru Matsui,? and William L. Jorgensen* Department of Chemistry, Yale University, New Haven, Connecticut 0651 I Received: July 7, 1994; In Final Form: August 26, 1994@
The OPLS force field for organic and biomolecular systems has been expanded to include an all-atom representation for saturated hydrocarbons. The model was formulated and tested via Monte Carlo simulations for both the pure liquids and dilute aqueous solutions of methane, ethane, propane, and butane and as such leaves no ambiguity in terms of its performance in either hydrophilic or hydrophobic environments. To evaluate accuracy, key comparisons with available experimental thermodynamic and structural data including heats of vaporization, molecular volumes, fluctuational properties, and radial distribution functions for the liquids as well as free energies of hydration were made and found to be in reasonable accord. Revisitation of the united-atom model showed that the all-atom force field performs with overall superiority (e.g., average deviations from experimentally measured free energies of hydration were reduced from 0.36 to 0.13 kcall mol). Comparisons of the results with those from other force fields are also provided.
Introduction Molecular dynamics and Monte Carlo statistical mechanics simulations have become powerful tools for studying organic and biomolecular systems.' They depend critically on potential functions to describe the intramolecular and intermolecular energetics.* In 1984, the OPLS force field for hydrocarbons was r e p ~ r t e dwith , ~ several attractive features: (i) the simplicity of the Coulomb plus Lennard-Jones 6-12 format, (ii) the computational efficiency of an united-atom representation for the CH,, groups, and (iii) the excellent (average errors of 2%) reproduction of experimental energies and densities of the 15 liquids tested . At that time, a continual evolution of force fields paralleling growth and availability of computing resources was forecast. Since then, several research groups have initiated programs targeting increased complexity and improved performance of the potential functions, culminating in reports of new molecular mechanics force field^.^ Though nonbonded contributions to the potential energy are emphasized, different functional forms are used. Furthermore, the bases for the design of the force fields vary; some were devised with the intention of expanding to a large repertoire of atom types which reproduce experimental while others aimed for a limited set of common organic functional groups, transferable to biological systems of interest!a-b In terms of improvements over the OPLS model, these efforts advanced an all-atom representation of the molecules; the benefits of such models have been discussed p r e v i o u ~ l y .However, ~ ~ ~ ~ only one group attempted to address the problem of intermolecular interactions in condensed phases; even then only relative free energies of hydration for methane, ethane, and propane were although they note that their nonbonded parameters were determined by empirically fitting pure liquid properties. Since knowledge of how a force field behaves in both hydrophilic and hydrophobic environments is of obvious importance to the study of biomolecular systems, further testing is needed on condensed-phase systems. Herein the development of an all-atom force field that is compatible with the OPLS potential functions is reported for t Permanent addrss: Matsushita Electric Industrial Co., Energy Research Laboratories, Moriguchi, Osaka 570, Japan. Abstract published in Advance ACS Abstracts, November 1, 1994. @
0022-365419412098- 13077$04.50/0
simple hydrocarbons and tested on methane, ethane, propane, andbuthe through Monte Carlo simulations of the pure liquids and dilute aqueous solutions. The force field is simple in form and yields accurate results for thermodynamic and structural properties of the pure liquids as well as relative and absolute free energies of hydration. The focus here is on intermolecular interactions. A subsequent paper dealing with a greater variety of systems including cycloalkanes will address the intramolecular energetics including conformational equilibria.6 Computational Details Intermolecular Potential Functions. While an united-atom representation for CH, groups is computationallyappealing, this approximation is known to become increasingly problematic in polar molecule^.“^+^^^^ Consequently, an all-atom version of the OPLS force field for organic and biomolecular systems* has been under development, beginning with aromatic molec u l e ~ and ~ ~ ,alkanes. ~ The intermolecular interactions are described in the usual Coulomb plus Lennard-Jones format (eq 1). n
h
As such, the interaction energy between molecules a and b is described by a summation over all pairwise interactions between sites i and j on the two molecules. Moreover, the A's and C ' s in eq 1 are functions of the Lennard-Jones u's and E ' S , such that Aii = 4 ~ i c 7 and , ~ ~ Cii = 4ciui6. Furthermore, standard combining rules are used, where Au = (Ai&)"* and Cu = (CiiCj)l/z.The intermolecular interactions are evaluated within a spherical cutoff radius and, for the aqueous solutions, are quadratically feathered via a switching function to zero over the last 0.5 A.8 In order to derive parameters for alkanes, numerous Monte Carlo simulations were performed for the pure liquids. In each case, the object of this procedure was a unique set of parameters which accurately reproduced primarily the experimental heat of vaporization,molecular volume, or, equivalently,density, and C-C radial distribution functions (rdf). For the simplest hydrocarbon, methane, which has one carbon and four equiva-
0 1994 American Chemical Society
13078 J. Phys. Chem., Vol. 98,No. 49, I994 TABLE 1: OPLS Parameters for Alkanes model atom q(e) a(A) all-atom C(&) -0.24 3.500 -0.18 3.500 C(H3) -0.12 3.500 C(H2) +0.06 2.500 WC) united-atom CH4 0.00 3.730 0.00 3.775 CH3 (ethane) 0.00 3.905 CH3 0.00 3.905 CH2
E
Kaminski et al.
(kcdmol) 0.066 0.066 0.066 0.030 0.294 0.207 0.175
ref
0.118
3 3 3 3
TABLE 2: Nonbonded Parameters for Alkanes from Other Force Fields model atom q (e) o(A) E (kcdmol) ref -0.27 3.207 4a 0.0903 -0.18 3.207 4a 0.0903 +0.09 0.0045 2.616 4a 9 O.Oo0 3.473 0.0951 O.Oo0 3.473 0.0951 9 O.Oo0 2.846 0.0152 9 -0.27 0.080 3.6705 10b -0.18 3.8754 0.055 10b +0.09 0.022 2.3520 10b lent hydrogens, there is a maximum of five adjustable parameters: qc = - 4 q ~plus u and E for both C and H. However, except for minor adjustments in the well depths, the knnardJones Parameters were required to be consistent with other values?bv8 In addition, earlier parametrizationssteered the initial guess for q H in the neighborhood of 0.05-0.10.5b.8 The resultant parameters were tried in pure liquid simulationsof the remaining hydrocarbons to test the transferability; of course, the charges on the CH3 and CH2 carbons were modified to maintain neutrality. The optimized values, as well as the OPLS unitedatom hydrocarbon parameters3 that were used in verifying the potential functions, are given in Table 1. To compare this all-atom model with some other recently reported alternatives, sample pure liquid and dilute aqueous solution simulations were carried out using the and DREIDING9 force fields. As both employ the same Coulomb plus Lennard-Jones 6- 12 functional form to evaluate the intermolecular interactions, the associated nonbonded parameters were easily incorporated into the OPLS functional form. Though bond stretching and angle bending may also be sampled, the emphasis here was placed on the intermolecular part of the potential functions. Finally, a set of unpublished CHARMM all-atom parameters that was used in work recently reported by Dunbrack and Karpluslo was also tested." The parameters from these force fields are given in Table 2. It may be noted that the parameters in these cases are readily transferable to larger alkanes, while for the approach of Sun et al. an ab initio calculation is needed to obtain the partial charges for each Intramolecular Potential Functions. Torsional motion was included in the simulations for ethane, propane, and butane. Specifically, rigid-body rotation about each C-C bond in the alkanes was performed. This requires torsional potentials, V(4), for the dihedral angles 4, which are represented by Fourier series (eq 2). V(4) = l/zvl(l
+ cos 4) + V2V*(1 - cos 24) + Y2V3(1+ cos 34) (2)
The rotation of methyl groups required a single 3-fold term with a barrier height V3 = 3.00 kcal/mol,12 while the rotation about the central C-C bond in butane necessitated all three terms. In the latter case, the previously reported Fourier coefficients (Vl = 1.522, V;!= -0.315, V3 = 3.027), which came from MM2 calculations, were used.3
Fluid Simulations. Monte Carlo statistical mechanics simulations were performed for liquid methane, ethane, propane, and butane as well as for the single hydrocarbon molecules in water. Fixed standard bond lengths and angles were adopted for the monomers;13 specifically, rCH = 1.09 A, rcc = 1.53 A (1.54 8, for ethane), LCCH = LHCH = 109.5' (for C-CH3), LHCH = 107.0' (where C is bonded to two other carbons), and LCCC = 112.0'. In all cases, the simulations were performed in the isothermal-isobaric ("T)ensemble at 1 atm and at either the boiling points for the pure liquids14 or 25 'C for the solutions. Additionally, standard procedures including periodic boundary conditions as well as Metropolis sampling,15 which was augmented by preferential sampling16 for the solutions, were emp10yed.l~ Further, the ranges for translations and rotations of the monomers were chosen to yield new configurations with a ca. 40% acceptance rate. Finally, the simulations were run in 2 x lo5 configuration batches to monitor convergence. The uncertainties in the computed properties ( f l u ) were obtained via the batch means procedure.17 All calculations were run on Silicon Graphics 4D/35 and W4000 workstations in our laboratory. For the pure liquids, each system consisted of 128 hydrocarbon molecules in a cubic cell, ca. 20-27 8, on a side. In all cases, the intermolecular interactions were spherically truncated at 11 8, based on the distances between central atoms, and a standard cutoff correction to the total energy beyond that radius was made.3 Additionally, changes in the dihedral angles were adjusted to obtain reasonable acceptance rates. Lengthy equilibrations were performed during the parametrization. Subsequently, simulations involving 1 x IO6 configurations of equilibration followed by averaging over 8 x lo6 configurations were performed to ensure adequate convergence of the thermodynamic properties. As a point of comparison, simulations employing the OPLS united-atom representation were then executed over the same number of configurations for each of the hydrocarbons studied. Additionally, all-atom simulations were carried out for methane and ethane using only the LennardJones interactions with the atomic charges set to zero in order to assess the role of the Coulombic component of the intermolecular potential functions. Simulations of butane were also performed with the CHARMM92, DREUIING, and CHARMM94 nonbonded parameters covering 1 x lo6 configurations of equilibration and 2 x lo6 configurations of averaging. The key structural and thermodynamic properties, namely, the heat of vaporization (AHvap),molecular volume (V), heat capacity (Cp), coefficient of thermal expansion (a), isothermal compressibility (K), and C-C radial distribution functions, were computed in the usual manner.3 It is worth noting that the heats of vaporization were computed from eq 3 where Einmand Ei are the average torsional and intermolecular energies. Ehwa(g)was determined from Monte Carlo simulations of an isolated monomer and was essentially identical to the total torsional energy computed in the pure liquid ~imulation.~ All pure liquid simulations were performed with the MCLiquid and BOSS programs.18 To test the potential functions in aqueous media, calculations were performed to determine the free energies of hydration of the hydrocarbons in TIP4P water.19 In the absence of intramolecular degrees of freedom, the absolute free energy of hydration is readily obtained through gradual disappearance of the solute; in fact, the value for united-atom methane was computed in this manner. For all-atom methane and ethane, mutations in which the Coulombic and Lennard-Jones components were decoupled were performed in order to avoid the production of unshielded charges in aqueous solution. For the larger solutes,
Pure Liquid Properties of Hydrocarbons
J. Phys. Chem., Vol. 98, No. 49, 1994 13079
TABLE 3: Computed and Experimental Thermodynamic Properties of Liquid Methand system all-atom 2M 4M 6M 8M L-Jall-atom 2M 4M 6M 8M
V
e
AHvaD
CD
a
K
57.6 f 0.1 57.3 f 0.1 57.2 f 0.1 57.2 f 0.1
0.462 f 0.001 0.465 f 0.001 0.466 f 0.001 0.466 f 0.001
2.175 f 0.006 2.184 f 0.006 2.187 f 0.004 2.187 f 0.004
12.4 f 0.6 14.4 f 0.9 14.2 f 0.7 13.8 f 0.5
328 f 28 414 f 38 399 f 30 386 f 23
160 f 13 192 f 17 184 f 12 180 f 10
57.2 f 0.2 57.3 f 0.1 57.3 f 0.1 57.3 f 0.1
0.466 f 0.002 0.465 f 0,001 0.465 f 0.001 0.465 f 0.001
2.189 f 0.007 2.184 f 0.005 2.186 f 0.004 2.186 f 0.004
12.9 f 0.8 13.2 f 0.6 13.5 f 0.6 13.8 f 0.6
337 f 35 349 f 27 369 f 27 384 f 25
156 f 14 162 f 11 172 f 12 180f 11
64.0 f 0.1 64.3 f 0.1 64.0 f 0.1 64.1 f 0.1 62.8
0.416 f 0.001 0.414 f 0.001 0.416 f 0.001 0.416 f 0.001 0.424
united-atom
2M 4M 6M 8M experimentb
1.924 f 0.004 13.5 f 0.9 356 f 49 237 f 30 1.923 f 0.002 13.3 f 0.5 348 f 29 238 f 17 1.925 f 0.002 13.3 f 0.4 352 f 24 241 f 15 1.924 f 0.002 13.5 f 0.4 362 f 20 248 f 13 1.96 13.2 349 157 Units are as follows: V, A3;e, g/cm3; AHvap,kcaumol; Cp,cal/(mol K); a,1 x atm. All properties measured at the boiling K; K , 1 x point, -161.49 O C . The reported statistical uncertainties are f l u as obtained from the fluctuations in separate averages over batches of 2 x lo5 configurations. See refs 3 and 25. TABLE 4: Computed and Experimental Thermodynamic Properties of Liauid Ethane system V e m v a p CP all-atom 93.4 f 0.5 0.535 f 0.003 2M 3.41 f 0.02 19.7 f 2.5 92.2 f 0.4 4M 0.542 f 0.002 22.1 f 2.1 3.45 f 0.02 92.5 f 0.3 6M 0.540 f 0.002 2.15 f 1.5 3.44 f 0.01 92.5 f 0.2 0.540 f 0.001 8M 3.44 f 0.01 20.6 f 1.5 L-J all-atom 2M 92.6 f 0.4 0.539 f 0.002 3.43 f 0.01 15.4 f 1.0 4M 92.4 f 0.3 0.540 f 0.002 3.44 f 0.01 17.7 f 1.1 92.5 f 0.2 6M 0.540 f 0.001 3.46 f 0.01 17.7 f 0.9 92.6 f 0.2 8M 0.539 f 0.001 3.43 f 0.01 17.5 f 0.7 united-atom 2M 91.0 f 0.7 0.549 f 0.004 3.54 f 0.03 27.4 f 4.6 4M 91.2 f 0.3 0.548 f 0.002 3.53 f 0.01 10.9 f 1.7 91.2 f 0.2 6M 0.548 f 0.001 3.53 f 0.01 19.0 f 1.0 91.2 f 0.2 8M 0.548 f 0.001 3.53 f 0.01 18.0 f 0.8 91.5 experimentb 0.546 3.62 17.6 Details as in Table 3. All properties measured at the boiling point, -88.63 "C. See ref 3. relative free energies of hydration were computed through the following thermodynamic cycle and connected to methane to obtain the absolute free energies of hydration. Here, A and B refer to two homologues (e.g., ethane and methane); the difference in free energy (AG) between them is computed through statistical perturbation theory (SPT).1,2*20As ethane, propane, and butane were mutated into the smaller alkanes, those atoms that disappear were converted to dummy atoms (q = u = E = 0) with shortened bond lengths. However, the torsional potentials are retained, and since the torsional energy is then the same for the reference and perturbed solutes, the torsional energy did not contribute to the free energy changes.21 By only including the intermolecular contributions to the free energy changes, the gas-phase AG is zero in each case, and the relative free energies of hydration (AhGyd) are given by -AG(,)(AB).
a
K
323 f 62 389 f 58 368 f 40 344 f 30
242 f 45 283 f 39 265 f 28 247 f 20
200 f 30 264 f 31 263 f 24 260 f 19
150 f 22 190 f 22 189 f 17 188 f 14
513 f 114 354 f 42 308 f 26 284 f 19 231
356 f 70 262 f 26 235 f 17 219 f 12 N/A
solute-solute and solute-solvent cutoffs of 9.5 A, based on roughly the center-of-mass separations. Attempted changes to the dihedral angles were made in ranges of f20". Each mutation involved four simulation steps with double-wide sampling22such that eight incremental free energy changes were calculated. As a final check of the potential functions, SET was used to compute the free energy change for rotation about the central bond for all-atom butane in water. Specifically, the C-C-C-C angle was varied in increments of 15" from the cis to the trans conformation; again, double-wide sampling was employed such that 12 incremental AG's were obtained. Earlier computations with the united-atom model employed this procedure and can be consulted for further details.23 All SPT computations involved an equilibration phase of 1 x lo6 configurations followed by 2 x lo6 configurations of averaging. These solution calculations were performed with the BOSS program.24 Results and Discussion
Additionally, the CHARMM92, CHARMM94, and DREIDING models for ethane were mutated to OPLS ethane to compare their performance in water. In all cases, one monomer plus 265 water molecules (ca. 20 x 20 x 20 A) were used, with
Pure Liquids. Tables 3-6 contain the computed thermodynamic properties of liquid methane, ethane, propane, and butane from the OPLS all-atom model. The results are reported in increments of 2 million configurations over the full averaging period of 8 million configurations in order to illustrate the convergence. Using the all-atom model, the heats of vaporization and molecular volumes converge rapidly and demonstrate good agreement with experiment; the average percent error is 3.5% overall or 1.2% excluding the results for methane. The
13080 J. Phys. Chem., Vol. 98, No. 49, 1994
Kaminski et al.
TABLE 5: Computed and Experimental Thermodynamic Properties of Liquid Propand system V e m v a P CP all-atom 2M 125.4 f 0.3 0.584 f 0.001 4.55 f 0.03 28.1 f 3.0 0.582 f 0.002 4.53 f 0.02 28.3 f 2.2 4M 125.9 f 0.5 6M 125.4 f 0.4 0.584 f 0.002 4.55 f 0.02 26.8 f 1.5 0.585 f 0.001 4.55 f 0.01 25.2 f 1.1 8M 125.2 f 0.3 muted-atom 2M 133.1 f 0.2 0.550 f 0.001 4.15 f 0.01 21.1 f 0.8 0.551 f 0.001 4.16 f 0.01 21.5 f 0.5 4M 132.9 f 0.3 0.550 f 0.001 4.15 f 0.01 21.8 f 0.5 6M 133.1 f 0.2 8M 133.1 f 0.2 0.550 f 0.001 4.15 f 0.01 21.8 f 0.4 experimentb 126.0 0.581 4.49 23 a Details as in Table 3. All properties measured at the boiling point, -42.1 "C. See ref 3. TABLE 6: Computed and Experimental Thermodynamic Properties of Liquid Butaid system V e m v a CP all-atom 2M 162.8 f 1.1 0.592 f 0.004 5.38 f 0.04 40.3 f 4.4 4M 161.9 f 0.7 0.596 f 0.003 5.40 f 0.02 37.1 f 2.3 6M 161.5 f 0.5 0.598 f 0.002 5.41 f 0.02 36.4 f 1.7 8M 161.3 f 0.4 0.598 f 0.001 5.43 f 0.02 37.4 f 1.6 united-atom 2M 162.7 f 0.6 0.593 f 0.002 5.38 f 0.03 30.2 f 1.5 4M 162.1 f 0.5 0.596 f 0.002 5.38 f 0.02 30.1 f 1.1 6M 162.5 f 0.3 0.594 f0.001 5.36 f 0.01 30.3 f 0.9 8M 162.3 f 0.3 0.595 f 0.001 5.37 f 0.01 30.2 f 0.8 experimentb 160.3 0.602 5.35 31.8 (I
K
234 f46 236 f 33 232 f 27 221 f 21
292 f 61 271 f 38 244 f 27 215 f 19
188 f 19 200 f 14 207 f 12 206 f 11 115
223 f 22 239 f 16 264 f 14 245 f 13 NIA
a
K
318 f 77 258 f 38 247 f 29 264 f 27
352 f 82 296 f 44 273 f 31 282 f 28
171 f 28 163 f 20 173 f 17 172 f 15 176
204 f 28 203 f 24 222 f 21 216 f 17 NIA
Details as in Table 3. All properties measured at the boiling point, -0.5 "C. See ref 3.
TABLE 7: Comparison of Thermodynamic Properties of All-Atom Liquid Butane from OPLS, CHARMM92, DREIDING,and CHARMM94 Force Field9 proOPLS CHARMM92 DREIDING CHARMM94 perty V 162.8 f 1.1 430.5 f 13.5 190.4 f 1.5 155.4 f 0.8 e 0.593 f 0.004 0.224 f 0.007 0.507 f 0.004 0.621 f 0.003 1.33 f 0 . 0 3 4.76f0.04 5.59f 0.05 AHva 5.38 f 0 . 0 4 25.8 f 1.7 38.4 f 4 . 0 29.0f2.8 C, 40.3 f 4 . 4 a 318f77 682f223 309f71 262f58 K 352f82 NIA 525 f 124 232 f 45 Details as in Tables 3 and 6. Results represent averages over 2 x lo6 configurations.
'
a
fluctuational properties, Cp,a,and K , converge more slowly, as duly noted b e f ~ r e . ~ Moreover, J ~ - ~ ~ it is clear that attainment of high precision for these properties, particularly a and K , would require further averaging beyond the 8 million configurations performed here. Nevertheless, the computed heat capacities of the alkanes are in reasonable agreement (average errors of 12%) with the experimental values. Within the OPLS framework, the computed thermodynamic properties of the pure liquids were found to be essentially independent of the model used; Le., the results were mostly the same, given the statistical uncertainties, for the united-atom, all-atom, and all-atom Lennard-Jones representations (Tables 3-6). The lack of impact of the Coulombic terms is notable though not unexpected in view of the small charges on hydrogen and the zero or negligible dipole moments for these alkanes. However, the CHARMM92 and DREIDING all-atom force fields gave substantially different results for the thermodynamic properties of liquid butane, as can be seen in Table 7. Notably, both parameter sets underestimated the density by 16% (DREDING) and 63% (CHARMM92). Moreover, the CHARMM92 system expanded over the entire 3 x lo6 configurations that were sampled. The problem with the DREIDING parameters appears to be a too large u and a too small E for hydrogen (Table 2), while the extreme expansion with CHARMM92 is most likely due to the small magnitude of the hydrogen e. This makes the system resemble a gas of weakly interacting molecules. To
verify this hypothesis, a simulation of OPLS butane was performed with a 10-fold decrease in the CH (0.003 kcdmol), which led to a similar volume explosion (V = 441.6 f 20.7 A) after the equilibration phase of 1 x lo6 configurations and 2 x lo6 configurations of averaging. The CHARMM94 parameters, which appear to have been tested to some degree on condensedphase systems, reproduced the experimental density and heat of vaporization with average errors of 3.2% and 4.5%, respectively. Finally, unpublished results with the all-atom AMBER force field for the pure liquids have yielded average errors for the densities and heats of vaporization of 3.8% or 4.4% excluding the results for methane.27 The computed structural properties, on the other hand, show some model dependence. Figures 1-3 indicate that, although it is difficult to discern any differences between the all-atom and L-J all-atom models, both lead to more liquid structure than their united-atom counterparts. Further, the all-atom models better reproduce the experimental C-C radial distribution functions for methane and ethane (Figure 1),28as they describe the alkanes with greater detail. Only the all-atom models were able to reproduce the experimentally observed minimum at 4.8 A and second maximum at 5.2 8, in the ethane C-C rdf.28b Additional structural details may be noted. Figures 2 and 3 illustrate the effects of chain length and molecular position of the carbon atoms on the C-C radial distribution functions. For instance, it is easily seen in Figure 2 that a longer carbon chain suppresses the peaks in the terminal C-C rdf; this is attributable to the greater variety of intermolecular contacts available to larger molecules and to shielding of the terminal carbons by the rest of the molecule. Similarly, the CHZ-CHZ peaks in Figure 3 are situated at greater interatomic distances than the CH3-CH3 maxima, presumably because of shielding provided by the external methyl groups. Finally, it should be mentioned that the radial distribution functions converge rapidly; the results after 2 x lo6 configurations of averaging are indistinguishable from those after 8 x lo6 configurations. Free Energies of Hydration. The results from the various mutations in dilute aqueous solution (TIP4P water) are given in Table 8. The average deviation (0.13 kcal/mol) from the
Pure Liquid Properties of Hydrocarbons 4.0
m All- A I ~ _ . All. VDW
-
- Eapnmnt
3.0
2.0
_ _ Unitcd Atom
- - UnitcdA
CH,
J. Phys. Chem., Vol. 98, No. 49, 1994 13081
All-Aim All. VDW EXpnmni
2 1.0
h
Q 2.0
1.0
0.0 C,H,
I
- CH,-CH,. Unitcd Atan '1.1 Cq-CH,. United Atan CH,-CH,. All-Aim CH,-CH,. All-Atom
3.0
-
Q 2.0
0.0 - CH,-CH,. Uniied Atom
C,H,,
... CH,-CH,. United Atom
- - CH,-CH,.
. .
All-Atom CH,-CH,. All-Atom
Q
1.0
1.0
0.0 3.0
4.0
5.0
.
6.0
3.0
4.0
..,A *
5.0
..
6.0
7.0
4.0
- Methane United-Atom Model
Ethane Propane, CH,-CH, Butane. CH,-CH,
3.0 -
1
4.0
5.0
6.0
7.0
3.0
4.0
r, A
Figure 1. C-C radial distribution functions for methane, ethane, propane, and butane from the pure liquid simulations.
g 2.0
0.0 3.0
..,A
I -
.
I \
5.0 6.0 r, A
7.0
8.0
Figure 3. Comparison of C-C rdf for the methyl and methylene groups in liquid propane and butane.
TABLE 8: Free Energies of Hydration Calculated with the OPLS Force Field; Intermediates in Parentheses Are for the All-Atom Simulations Only process
-
AG(all-atom) AG(united-atom) AG(exptl)b
C& (C&",) 0 -2.20 f 0.22 CZHS cH4 0.25 f 0.06 CZHS~(C&~~)-O -1.67f0.18 C3HS c2H6 -0.22 f 0.08 Ca10 C3H8 -0.29 f 0.11 4
-2.46 f 0.36 0.13 f 0.20 -2.26f0.33 -0.98 f 0.10 -0.27 f 0.17
-2.00 0.17 -1.83 -0.13 -0.12
Units are kcal/mol. Experimental values are from ref 29. C4HIo 00
2.47 2.08
- Methane 30
-
All-Atom Model
__
Ethane Propane, CH,-CH,
- - Butane,CH,-CH,
-0.24 C3H8 2.18 1.96
- - -0.22
C2H6
1.95 1.83
f0.25
CH4
2.20 2.00
-2.20
0 OPLS exptl
.
- -.-._ --_ _ _ --- -
experimentally measured AG's is smaller for the all-atom model, although in most cases the united-atom model also affords reasonable results. Connection of the individual A G s , as shown below, yields computed absolute free energies of hydration of the hydrocarbons that have the correct magnitudes and ordering as experimentally observed.29 Notably, ethane is the most watersoluble.
Moreover, it should be noted that both the direct (by full annihilation) and indirect (by conversion to methane) methods for computing the absolute free energy of hydration of ethane yield values in accord with experiment (1.67 f 0.18 and 1.95 .f 0.23 kcdmol). The comparison studies between the OPLS all-atom model and other all-atom force fields in aqueous solution revealed some interesting differences. The models with the small E on hydrogen, CHARMM92 and the contrived OPLS/lO (10-fold reduction in EH), were markedly more hydrophilic than the OPLS all-atom model. Though a change in E does not alter the position of the minimum or intercept for a Lennard-Jones potential, the slope (dEldr) is proportional to E . Thus, a 10-fold reduction in E makes the repulsive wall much softer, which presumably leads to the reduced hydrophobicity. Mutation of the OPLS all-atom ethane to the CHARMM92 ethane in water yielded a free energy difference of -2.67 f 0.08 kcdmol, and the mutation to OPLS/ 10 ethane gave - 1.22 f0.02 kcdmol. The DREIDING model, on the other hand, was excessively hydrophobic; mutation of the OPLS all-atom ethane to the DREIDING ethane favored the former by 0.76 f 0.03 kcdmol. However, the CHARMM94
Kaminski et al.
13082 J. Phys. Chem., Vol. 98, No. 49, 1994
MacKerell for the CHARMM94 parameters that were used in ref 10. E.M.D. is grateful for fellowship support from the Organic Division of the American Chemical Society sponsored by the Rohm and Haas Co. References and Notes (1) (a) Beveridge, D. L.; DiCapua, F. M. Annu. Rev. Biophys. Biophys. Chem. 1989,18,431. (b) Jorgensen, W. L. Acc. Chem. Res. 1989,22, 184. (c) Kollman, P. A.; Men, Jr., K. M. Acc. Chem.Res. 1990, 23, 246. (d) Jorgensen, W. L. Chemtracts: Org. Chem. 1991, 4, 91. (e) Straatsma, T. P.; McCammon, J. A. Annu. Rev. Phys. Chem. 1992,43,407. (f) Kollman, P. A. Chem. Rev. 1993, 93, 2395. (2) See, for example: (a) Brooks, B. R.; Bruccoleri, R. E.; Olafson, B. D.; States, D. J.; Swamhathan, S.;Karplus, M. J. Compur. Chem. 1983, 4, 187. (b) Weiner, S. J.; Kollman, P. A.; Nguyen, D. T.; Case, D. A. J. Comput. Chem. 1986, 7, 230. (c) Levitt, M. J. Mol. Biol. 1983, 168, 595. (3) Jorgensen, W. L.; Madura, J. D.; Swenson, C. J. J. Am. Chem.
1
J
- 0 8 0.0
30.0
60.0
90.0
120.0
150.0
1800
C-C-C-C (dcg)
Figure 4. Computed hydration effects on the barrier to rotation about the central C-C bond in butane at 25 "C. The simulations were performed in increments of 15 deg. The results for the united-atom model are from ref 23.
model like the all-atom OPLS model performed well for the free energy of hydration of ethane; mutation of the OPLS allatom ethane to the CHARMM94 ethane yielded a free energy difference 0.18 f 0.01 kcdmol. With the AMBER all-atom model and TIP3P water, methane to ethane and ethane to propane perturbations yielded reasonable free energy changes of 0.16 f 0.04 and 0.20 f 0.04 kcdmol, though the sign is incorrect in the first Another ethane to propane mutation with a united-atom model and TIP3P water also yielded the wrong sign and a 1.0 kcdmol e ~ ~ othough r , ~ convergence ~ of the results was probably not achieved in this earlier work based on current procedures. The all-atom OPLS model was also found to furnish the expected effects of hydration on the butane torsional profile, as depicted in Figure 4. Results from the earlier united-atom are also shown. Though there are some differences between the two models, the trend is unambiguous. The gauche conformer is better hydrated than trans butane by 0.3 kcdmol, and the cis form is still lower by nearly the same amount. This result reflects hydrophobic stabilization of the more compact gauche rotamer and is consistent with prior results from integral equation theories and simulations using united-atom model^.^^*^^ Conclusion Studies of both the pure liquids and dilute aqueous solutions for methane, ethane, propane, and butane were undertaken to develop and test the OPLS all-atom potential functions for alkanes. In general, the computed thermodynamicand structural properties of the pure liquids were shown to be in very good agreement with available experimental data. Moreover, the allatom model definitely improves the representation of the structure in the pure liquids and of the free energetics for hydrocarbons in aqueous media over the simpler united-atom model. An important observation that emerged from comparisons with other force field parameters is that the Lennard-Jones E ' S do significantly affect the description of the alkanes in condensed phases, emphasizing once more the need for thorough testing in these media. Acknowledgment. Gratitude is expressed to the National Institutes of Health for support of this work and to Dr. Alex
SOC. 1984, 106, 6638. (4) (a) Smith, J. C.; Karplus, M. J. Am. Chem. SOC. 1992, 114, 801. (b) Sun, Y.;Spellmeyer, D.; Pearlman, D. A.; Kollman, P. J. Am. Chem. SOC. 1992, 114, 6798. (c) Halgren, T. A. J. Am. Chem. SOC. 1992, 114, 7827. (d) Amcdeo, P.; Barone, V. J. Am. Chem. SOC.1992,114,9085. (e) Rap@, A. K.; Casewit, C. J.; Colwell, K. S.;Goddard, W. A. Skiff, W. M. J. Am. Chem. SOC.1992, 114, 10024. (5) (a) Jorgensen, W. L.; Briggs, J. M. J. Am. Chem. SOC. 1989, 111, 4190. (b) Jorgensen, W. L.; Severance, D. L. J. Am. Chem. SOC. 1990,112, 4768. (c) Pranata, J.; Wierschke, S . G.; Jorgensen, W. L. J. Am. Chem. SOC. 1991, 113, 2810. (6) Maxwell, D. S.;Tirado-Rives, J.; Jorgensen, W. L. To be published. (7) Brooks,C. L. Karplus, M.; Pettitt, B. M. Proteins: A %oretical
m,
Perspective of Dynamics, Structure, and Thermodynamics; Wiley: New York, 1988. (8) Jorgensen, W. L.; Tirado-Rives, J. J. Am. Chem. SOC. 1988, 110, 1657. (9) Mayo, S. L.; Olafson, B. D.; Goddard, W. A. ID J. Phys. Chem. 1990, 94, 8897. (10) (a) Dunbrack, R.; Karplus, M. Struct. Biol. 1994, 1, 334. (b) MacKerell, A. Private communication. (11) For clarity, the two CHARMM parameter sets will be denoted by the year in which they appeared in the literature, CHARMM92 and cHARMM94. (12) Lowe, J. P. Barriers to Internal Rotation; Streitwieser, A., Jr., Taft, R. W., Eds.;Progress in Physical Organic Chemistry; Interscience: New York, 1968; Vol. 6. (13) (a) Pople, J. A.; Gordon, M. J. Am. Chem. SOC. 1967, 89, 4253. (b) Hehre, W. J.; Ditchfield, R.; Stewart, R. F.; Pople, J. A. J. Chem. Phys. 1970, 52, 2769. (14) CRC Handbook of Chemistry and Physics, 73rd ed.; Lide, D. R., Ed.; CRC: Boca Raton, E,1992c. (15) Metropolis, N.; Rosenbluth, A. W.; Teller, A. H.; Teller, E. J. Chem Phys. 1953, 21, 1087. (16) Owicki, J. C. ACS Symp. Ser. 1978, No. 86, 159. (17) Allen, M. P.; Tildesley, D. J. Computer Simulations of Liquids; Clarendon: Oxford, England, 1987. (18) Jorgensen, W. L.; Matsui, T. MCLiquid, User's Manual, Yale University, New Haven, CT, 1991. Jorgensen, W. L. BOSS, Version 3.5, Yale University, New Haven, CT, 1994. (19) Jorgensen,W. L.; Chandrasekhar,J.; Madura, J. D.; Impey, R. W.; Klein, M. L. J. Chem. Phys. 1983, 79, 926. (20) Zwanzig, R. W. J. Chem. Phys. 1954, 22, 1420. (21) Jorgensen, W. L.; Nguyen, T. B. J. Comput. Chem. 1993,14, 195. (22) Jorgensen,W. L.; Ravimohan, C. J. J. Chem. Phys. 1985,83,3050. (23) Jorgensen. W. L.; Buckner, J. K. J. Phys. Chem. 1987, 91, 6083. (24) Jorgensen, W. L. BOSS Version 3.2, Yale University, New Haven, CT, 1992. (25) Sychev, V. V. Thermodynamic Properties of Methane, Hemisphere: Washington, DC,1987. (26) Jorgensen, W. L. Chem. Phys. Len. 1982, 92, 405. (27) Kollman, P. A. Personal communication. (28) (a) Habenschuss, A.; Johnson, E.; Narten, A. H. J. Chem. Phys. 1981, 74, 5234. (b) Sandler, S. I.; Lombardo, M. G.; Wong, D. S.-H.; Habenschuss, A.; N e n , A. H. J. Chem. Phys. 1982, 77, 2144. (29) Ben-Naim, A.; Marcus, Y.J. Chem. Phys. 1984, 81, 2016. (30) Fleischman, S. H.; Brooks, C. L., IlI. J. Chem. Phys, 1987, 87, 3029. (31) (a) Ran, L. R.; Chandler, D. J. Chem. Phys. 1977, 67, 3683. (b) Mikkilineni, R.; Beme, B. J. J. Am. Chem. SOC. 1982, Rosenberg, R. 0.; 104,7647. (c) Jorgensen, W. L. J . Chem. Phys. 1982, 77,5757. (d) Zichi, D.; Rossky, P. J. J. Chem. Phys. 1986, 84, 1712.
