Thermodynamic Properties of the Williams, OPLS-AA, and MMFF94 All

2578

J. Phys. Chem. B 1998, 102, 2578-2586

Thermodynamic Properties of the Williams, OPLS-AA, and MMFF94 All-Atom Force Fields for Normal Alkanes Bin Chen, Marcus G. Martin, and J. Ilja Siepmann* Department of Chemistry, UniVersity of Minnesota, 207 Pleasant Street SE, Minneapolis, Minnesota 55455-0431 ReceiVed: NoVember 6, 1997; In Final Form: January 23, 1998

The performance of several all-atom force fields for alkanes is compared and evaluated. Configurationalbias Monte Carlo simulations in the Gibbs ensemble were carried out to calculate the vapor-liquid phase equilibria for methane, ethane, n-butane, n-pentane, and n-octane. The Williams, OPLS-AA, and MMFF94 force fields were selected as representative all-atom models for this study because they were fitted using three different strategies (Williams, crystal structures and heats of sublimation; OPLS-AA, liquid densities and heats of vaporization; MMFF94, rare gas pair potentials and quantum mechanics) and employ potentials with three different functional forms to describe nonbonded van der Waals interactions (Williams, Buckingham exp-r-6 ; OPLS-AA, Lennard-Jones 12-6; MMFF94, buffered 14-7). It is shown that seemingly small differences in the potential functions can account for very large changes in the fluid-phase behavior. The Williams and OPLS-AA force fields yield liquid densities, boiling temperatures, and critical points that are in acceptable, albeit not in quantitative agreement with experiments, whereas the fluid-phase behavior of the MMFF94 model shows very large deviations.

I. Introduction Normal alkanes have been the focus of numerous computer simulation studies for the past 20 years starting with the seminal work of Ryckaert and Bellemans.1,2 The main reasons for the great interest in alkanes are (i) alkanes are the prototypical case of molecules with an articulated structure; (ii) alkanes are of prime importance for the (petro-) chemical industries; and (iii) alkyl groups are pivotal for the function of many biological systems. The alkane models that have been used in molecular simulations can be divided into two main categories: unitedatom (UA) models and all-atom (AA) models (sometimes also referred to as explicit-hydrogen models). In the UA representation, methyl and methylene segments are treated as single pseudoatoms with their interaction sites commonly located at the position of the carbon atoms. Hydrogen atoms are treated explicitly in AA models. The primary advantage of the UA models is their computational efficiency; that is, by reducing the number of interaction sites by a factor of 3, CPU time savings of an order of magnitude can be achieved. On the other hand, the AA models give a more faithful description of the shape of real n-alkanes and are therefore considered to be more appropriate for simulations of solid or high-density (lowtemperature) liquid phases.3,4 In addition, AA models allow for the distribution of partial charges on the individual hydrogen and carbon atoms, which may be important to describe the interactions of alkanes with more polar molecules. Jorgensen and co-workers5 have recently studied the performances of three AA force fields (OPLS-AA,6 AMBER94,7 and MMFF948). They carried out Monte Carlo simulations in the isobaric-isothermal ensemble for the liquid phase of n-butane at its normal boiling point of -0.5 °C and reported that the liquid density and heat of vaporization are well reproduced by the OPLS-AA and AMBER94 force fields, whereas the MMFF94 * Corresponding author: [email protected].

force field did not yield a stable liquid phase at the selected conditions. The aim of this paper is a more thorough investigation of the ability of three AA force fields to describe the fluidphase behavior of linear alkanes. The Williams,9 OPLS-AA,6 and MMFF948 force fields were selected as representative allatom models for this study because they were fitted using three different strategies. The Williams model was parametrized using the crystal structures and heats of sublimation for n-pentane, n-hexane, and n-octane.9 For the OPLS-AA model fits to the liquid densities and heats of vaporization for ethane, propane, and n-butane were employed.6 Finally, Halgren8,10 used experimentally derived rare gas potentials and high-quality ab initio calculations to obtain the parameters for his MMFF94 force field. In addition, the three force fields employ different functional forms to describe the nonbonded van der Waals interactions (Williams, Buckingham exp-r-6; OPLS-AA, Lennard-Jones 12-6; MMFF94, buffered 14-7) and use different combining rules to determine the van der Waals interaction parameters for the interactions of unlike atoms. The remainder of this article is arranged as follows. Sections 2 and 3 give descriptions of the molecular force fields and simulation details. In the first part of section 4, the influence of the functional form used for nonbonded interactions on the fluid-phase envelope is discussed. Thereafter, thermodynamic and structural data for the three different force fields (and some minor variations thereof) are presented. Finally, section 5 emphasizes some important points. II. Force Fields A. Intramolecular Interactions. The n-alkanes are treated as semiflexible chain molecules in a way very similar to the constrained molecular dynamics model used by Ryckaert and co-workers.3 The lengths of all C-C and C-H bonds are fixed, with parameters listed in Table 1. A harmonic bond bending

S1089-5647(98)00106-0 CCC: $15.00 © 1998 American Chemical Society Published on Web 03/17/1998

All-Atom Force Fields for Alkanes

J. Phys. Chem. B, Vol. 102, No. 14, 1998 2579

TABLE 1: Bond Lengths, Equilibrium Bond Angles, and Force Constants Used for the Williams, OPLS-AA, and MMFF94 Force Fields dC-C (Å) dC-H (Å) θ0(C-C-C) (deg) kθ(C-C-C)/kB (K rad-2) θ(C-C-H) (deg) θ(H-C-H) (deg)

Williams

OPLS-AA

MMFF94

1.530 1.040 112.0 62556 112.5 106.0

1.529 1.090 112.7 58765 110.7 107.8

1.508 1.093 109.608 61637 110.549 108.836

potential is used for all C-C-C bond angles.11

ubend )

kθ (θ - θ0)2 2

(1)

where kθ and θ0 are the bending force constant and the equilibrium bond angle (the parameters are given in Table 1). The H-C-H bond angles of all units and the C-C-H bond angles containing a methyl group hydrogen are fixed, whereas the C-C-H bond angles, for which the central carbon atom belongs to a methylene segment, are determined via the same geometric construction as used in ref 3. Thus, no contributions to the intramolecular potential energy arise from H-C-H and C-C-H bond angles. The OPLS united-atom torsional potential12 governs the motion of C-C-C-C dihedral angles,

utors(C-C-C-C) ) c1[1 + cos(φ)] + c2[1 - cos(2φ)] + c3[1 + cos(3φ)] (2) with the following numerical constants: c1/kB) 355.03 K, c2/ kB) -68.19 K, and c3/kB ) 791.32 K. The torsional motion of dihedral angles involving the hydrogen atoms on a methyl group is governed by a 3-fold symmetric torsional potential,

utors(X-C-C-H) ) cX[1 - cos(3φ)]

(3)

with cC/kB ) 854 K13 and cH/kB ) 717 K.14 There are no potential energy contributions from torsional angles involving methylene group hydrogens. Nonbonded van der Waals interactions are included for all atoms that are separated by at least three C-C bonds; that is, there are no 1-4 interactions of any type. 1-5 interactions involving one hydrogen atom and 1-6 interactions for two hydrogen atoms are also not considered as nonbonded interactions, because they are accounted for by the torsional potentials. As mentioned above, our choice of intramolecular potentials is closely related to the constrained molecular dynamics model of Ryckaert and coworkers,3 who studied the rotator phase of n-triacosane using the Williams force field. This choice is also computationally very efficient for our Monte Carlo calculations (see section 3). However, it differs from the more flexible models used commonly in biomolecular simulations,15 which treat all bond angles as flexible, explicitly include torsional potentials for all dihedral angles, and often contain 1-4 nonbonded interactions. Nevertheless, it is believed that the choice of intramolecular potentials is of secondary importance for the determination of thermodynamic properties, such as vapor-liquid phase equilibria, as long as it does not result in a change of the conformational characteristics of the molecule. We have carried out some test simulations for more flexible representations, which allowed all bond angles to vary, but found the influence of this change on the fluid-phase behavior to be negligible (see section 5). B. Intermolecular Potentials. As mentioned above, the Williams, OPLS-AA, and MMFF94 force fields employ po-

Figure 1. Comparison of the nonbonded interaction potentials used in the three force fields. The left, middle, and right set of curves are for H-H, C-H, and C-C interactions, respectively. Solid, dashed, and dotted lines represent Williams VII, OPLS-AA, and MMFF94 potentials, respectively.

tentials with different functional forms for the nonbonded van der Waals interactions. The different van der Waals potentials are compared graphically in Figure 1. The Williams force field9 makes use of Buckingham exp-6 potentials:16

uvdW(rij) ) Aij rij-6 + Bij exp(-Cijrij)

(4)

where rij is the separation between two atoms of types i and j. The size parameters C are CCC ) 3.74 Å-1, CCH ) 3.67 Å-1, and CHH ) 3.60 Å-1. Eight different combinations of A and B parameters were proposed by Williams.9 Here, we used parameter sets IV and VII, because set VII had been used earlier for simulations of alkane crystals3 and set IV was recently found to reproduce the experimental liquid densities of n-alkanes over a wide range of chain lengths and temperatures.17 The parameters for set IV are ACC/kB ) -2.858 × 105 K Å6, ACH/ kB ) -6.290 × 104 K Å6, AHH/kB ) -1.374 × 104 K Å6, BCC/ kB ) 4.208 × 107 K, BCH/kB ) 4.411 × 106 K, and BHH/kB ) 1.335 × 106 K; those for set VII are ACC/kB ) -2.541 × 105 K Å6, ACH/kB ) -6.441 × 104 K Å6, AHH/kB ) -1.625 × 104 K Å6, BCC/kB ) 3.115 × 107 K, BCH/kB ) 5.536 × 106 K, and BHH/kB ) 1.323 × 106 K. The OPLS-AA force field6 utilizes Lennard-Jones 12-6 potentials18 for the van der Waals interactions:

[( ) ( ) ]

uvdW(rij) ) 4ij

σij

12

-

rij

σij

6

rij

(5)

with CC/kB ) 33.2 K, HH/kB ) 15.1 K, σCC ) 3.50 Å, and σHH ) 2.50 Å. The parameters for unlike interactions are computed using standard Lorentz-Berthelot combining rules:19,20

1 σij ) (σii + σjj) 2

ij ) xiijj

(6)

In addition, a Coulombic interaction is used in the OPLS-AA force field:

uq(rij) )

qiqj 4π0rij

(7)

where 0 is the permittivity of vacuum. The partial charge

2580 J. Phys. Chem. B, Vol. 102, No. 14, 1998

Chen et al.

carried by a hydrogen atom is qH ) +0.06e, while the charge of a carbon atom depends on the number of hydrogens bonded to it. Each methyl, methylene, or methane group is required to carry no net charge; that is, qmethyl C ) -0.18e, qmethylene C ) -0.12e, and qmethane C ) -0.24e. The MMFF94 force field is based on well-characterized interaction potentials for rare gas atoms and ab initio calculations.8 Halgren10 proposed a buffered 14-7 potential for the van der Waals interactions:

[

uvdW(rij) ) ij

][

1.07 (rij/r* ij) + 0.07

7

]

1.12 -2 7 (rij/r* ij) + 0.12

(8)

where r*ij and ij are the idealized minimum energy separation and well depth for the interaction of atoms of types i and j.21 The interactions diameters for like atoms are calculated from 1/4 r* ii ) AiRi

(9)

where R is the atomic polarizability and A is a numerical constant. The diameters for unlike atoms are obtained from

[ (

[ ( ) ])]

2 r* ii - r* jj 1 r*ij ) (r* + r* ) 1 + 0.2 1 exp -12 jj 2 ii r* ii + r* jj

(10)

and the interaction strength is calculated from

ij )

-6 181.16GiGjRiRj(r* ij)

(Ri /Ni)1/2 + (Rj /Nj)1/2

(11)

where G and N are a scale factor obtained for rare gas atoms and the Slater-Kirkwood effective number of electrons, respectively. The resulting values for hydrogen and carbon atoms are r* CC ) 3.94 Å, r* CH ) 3.60 Å, r* HH ) 2.97 Å, CC/kB ) 34.12 K, CH/kB ) 14.13 K, and HH/kB ) 10.86 K. There are no partial charges on carbons and hydrogens belonging to alkanes in the MMFF94 force field. III. Computational Details A combination of the Gibbs ensemble Monte Carlo method (GEMC)22-24 and the configurational-bias Monte Carlo (CBMC) algorithm25-29 was employed to study the vapor-liquid phase equilibria for five linear alkanes. The system sizes (number of alkane molecules and combined volume of the two simulation boxes) were chosen to yield liquid-phase simulation boxes with linear dimensions of at least 20 Å and larger vapor phases containing at least five molecules. Three hundred methane, 240 ethane, 150 n-butane, 120 n-pentane, or 90 n-octane molecules were used in these simulations. For a given alkane, the initial configuration was constructed as a layered crystal with all chains in the all-trans conformation and oriented with their long axes parallel to the z-axis of the Cartesian coordinate system. In the x-y plane, the molecules were arranged in a hexagonal structure with the backbone planes forming a herringbone pattern. Five different kinds of Monte Carlo moves were used to sample phase space: translations of the center-of-mass, rotation around the center-of-mass, conformational changes using CBMC, volume exchanges between the two boxes, and CBMC particle swaps between the two boxes. The maximum displacements used for the translational, rotational, and volume moves were adjusted to yield acceptance rates of 50%, where different maximum translational and rotational displacements were used

for the vapor and liquid boxes.30 Several thousand Monte Carlo cycles (N moves) at elevated temperature were used to “melt” the initial structures, and the system was then “cooled” to the desired temperature. During this initial phase of the equilibration procedure, only translational, rotational, and conformational moves were employed. Thereafter, the systems were allowed to equilibrate for at least 10 000 Monte Carlo cycles using all five types of moves. The frequencies of swap and volume moves were adjusted for this part and the production runs to yield approximately one accepted swap and one accepted volume move per 10 Monte Carlo cycles. The remainder of moves were equally divided among the other three move types. For each alkane, simulations were carried out for five to seven different temperatures below the critical temperature. Most simulations were started using the final configuration of a neighboring temperature as the initial configuration (avoiding the melting and cooling phases). Again, these simulations were equilibrated for at least 10 000 Monte Carlo cycles. The results reported in the next section were averaged over either 10 000 or 20 000 Monte Carlo cycles, and the statistical uncertainties were estimated by dividing the simulations into 10 blocks. A semiflexible representation of the alkanes was used for all calculations, with the exception of one simulation for methane. As mentioned earlier, all bond lengths are fixed throughout the simulations. During the conformational and swap moves only the carbon atoms are grown using the CBMC algorithm.26-29 The number of trial sites used during the CBMC growth is 4, 8, 8, and 12 for ethane, n-butane, n-pentane, and n-octane, respectively. Ten trial sites are used for the insertion of the first bead during a swap move.31,32 Once all carbon atoms are successfuly generated, the positions of the methylene hydrogens are obtained by a geometric procedure similar to that used to enforce constraints in the molecular dynamics calculations.3 Since the C-H bond length and the methylene H-C-H angle are fixed, forcing the H-C-H plane to be perpendicular to the C-C-C plane and the H-C-H bisector to lie in the C-C-C plane and point away from the methylene carbon is sufficient to assign the methylene hydrogen positions. To calculate the methyl hydrogen positions, a random C-C-C-H torsional angle is selected from the correct distribution using a Boltzmann rejection technique.33 The torsional angle then defines the positions of the methyl group hydrogens using a C3 symmetry and the fixed C-H bond length and methyl C-C-H angle. For hydrogens on methane, the positions are determined by the tetrahedral angle and C-H bond length with a randomly chosen orientation. After all hydrogen positions are known, the Lennard-Jones interactions of these are calculated and the corresponding Boltzmann weight is used in addition to the Rosenbluth weight for the acceptance of the CBMC move. Similar procedures for the anisotropic Toxvaerd model or for CBMC with two potential truncations have been described in detail elsewhere.34,35 Spherical potential truncations at rcut ) 9 Å are used for all nonbonded interactions. An atom-based truncation is used for the van der Waals interactions, whereas a group-based truncation is used for the Coulombic interactions.36 In this case, the distance between the two carbon positions is used to decide whether or not the interactions of an entire methane molecule, methylene, or methyl group should be included. An additional center-of-mass (COM)-based radial cutoff at rCOM ) rcut + dA+ dB, where dA is the maximum distances between any atom on molecule A and its COM, is used to determine whether any nonbonded interactions between molecules A and B have to be computed.37 If not noted otherwise, long-range corrections were



Figure 2. Comparison of the Lennard-Jones 12-6 (solid line), Buckingham exp-6 (dashed line), and buffered 14-7 (dotted line) potential functions.

Figure 3. Comparison of the second virial coefficients obtained for the Lennard-Jones 12-6 (solid line), Buckingham exp-6 (dashed line), and buffered 14-7 (dotted line) potentials shown in Figure 2. T* is the reduced temperature, i.e. T* ) T/(/kB).

applied to the van der Waals interactions. The usual analytical tail corrections33 were used for the Lennard-Jones and Buckingham potentials. The energy and pressure tail corrections for the buffered 14-7 potential were obtained from numerical integration and are given by 3 U ULC ) 2πNFXYr* XY CXY

(12)

3 P PLC ) -(2/3)πF2XYr* XY C XY

(13)

where C UCC ) -0.024 84, C UCH ) -0.017 58, C UHH ) -0.008 377, C PCC ) 0.1696, C PCH ) 0.1204, and C PHH ) 0.055 78 for rcut ) 9 Å (see also ref 38). A long-range correction for the Coulombic interactions in the OPLS-AA force field was not included in this work. As will be demonstrated later, the small partial charges on the hydrogens and carbons have only a negligible effect on thermodynamic properties, since the architecture of the alkanes leads to an effective cancellation of the charge-charge interactions between different molecules. IV. Results and Discussion A. Fluid-Phase Equilibria. Dependence on the Functional Form of the Van der Waals Potential. To determine the influence of the functional form of the van der Waals potential on the fluid-phase behavior, we have first studied systems consisting of single-bead particles. To allow for a meaningful comparison, the parameters for the Buckingham exp-6 and buffered 14-7 potentials were scaled to yield the same position and depth of the potential energy minimum (and for the exp-6 potential also the same prefactor in the r-6 term) as the LennardJones potential. The three resulting potential functions are shown in Figure 2. At first glance, the differences between the three potential functions are very minor. The potential well of the exp-6 potential is slightly wider than the well of the Lennard-Jones 12-6 potential, where most of the difference is in the repulsive region. In contrast, the buffered 14-7 potential yields a narrower well compared to Lennard-Jones 12-6, and most of the difference is in the attractive region. To what extent can these apparently small differences affect thermodynamic properties? As should be expected, the differences become more pronounced when the second virial coefficients, B2, are compared (see Figure 3). The magnitude of B2 is approximately

Figure 4. Comparison of the vapor-liquid coexistence curves and critical points obtained for the Lennard-Jones 12-6 (solid lines and circles), scaled Buckingham exp-6 (dashed lines and diamonds), and buffered 14-7 (dotted lines and triangles) potentials. The open symbols are used for the coexistence densities obtained from the GEMC simulations, the filled symbols are the estimates of the critical points, and the lines are fits to the scaling law. T* and F* are the reduced temperature and density.

20% smaller for the buffered 14-7 potential than for the exp-6 potential. The Lennard-Jones 12-6 potential falls in between, but much closer to the exp-6 potential. Gibbs ensemble simulations were carried out to determine the vapor-liquid coexistence curves for single-bead systems interacting with the three scaled potentials.39 The resulting phase diagrams are shown in Figure 4. Again, the differences are appreciable. The reduced critical temperatures are 1.401 ( 0.008, 1.297 ( 0.008, and 1.126 ( 0.006 for the exp-6, Lennard-Jones 12-6, and buffered 14-7 potential, respectively. Thus great care should be taken when we want to compare force fields that are based on different functional forms for the van der Waals interactions. The potentials (such as plotted in Figure 1) will require different well depths to yield the same thermodynamic properties. It should also be noted that this will obviously lead to problems when van der Waals potentials are fitted to ab initio energies.8,40

2582 J. Phys. Chem. B, Vol. 102, No. 14, 1998 TABLE 2: Details of the Variations of the Williams, OPLS-AA, and MMFF94 Force Fields Used in This Work

W4 W7 W7-LC OA OA-q MM MM-LC MMOPLS

vdW parameters

partial charges

vdW tail corrections

Williams IV Williams VII Williams VII OPLS-AA OPLS-AA MMFF94 MMFF94 MMFF94

no no no yes no no no no

yes yes no yes yes yes no yes

intramolecular structure Williams Williams Williams OPLS-AA OPLS-AA MMFF94 MMFF94 OPLS-AA

Figure 5. Vapor-liquid coexistence curves for variations of the Williams force field: (A) methane, ethane, and n-butane; (B) n-pentane and n-octane. The solid lines and filled circles represent the experimental coexistence densities and critical points.42,43 The calculated coexistence densities and critical points are shown as diamonds, squares, and triangles for the W4, W7, and W7-LC force fields, respectively.

Williams Force Field. Simulations were carried out for three different variations of the Williams force field9 (see also Table 2): W4, Williams parameter set IV; W7, Williams parameter set VII; W7-LC, same as W7 but without analytical tail corrections for the vdW interactions. As mentioned earlier, these two parameter sets were used in previous simulations for n-alkanes.3,17 The vapor-liquid coexistence curves (VLCC) for the Williams force field are shown in Figure 5, the corresponding numerical data are listed in Table S1, and the estimated critical data and standard boiling temperatures are given in Table 3. The critical temperature and density are determined by leastsquares fits to the scaling law (using a scaling exponent of β ) 0.32 ) and the law of rectilinear diameters24 (see ref 37 for a detailed discussion of the uncertainties in the predicted critical

Chen et al. properties that may arise from finite-size effects and from changes in the scaling exponent), and the critical pressure and boiling temperature are obtained from a Clausius-Clapeyron plot.37,41 With the exception of methane, the VLCC for W4 and W7 are very similar and yield critical points that agree to within the uncertainties of the calculations. For methane, W7 gives a critical temperature that is 5% higher than that for W4, and the saturated liquid densities for W7 are shifted to higher densities. As a result, the VLCC for W7 is in better agreement with experiment.42,43 For ethane, the critical temperatures for W4 and W7 are approximately 10% too low, which results in a shift of the VLCC and too high vapor and too low liquid densities. The saturated liquid densities of n-butane, n-pentane, and n-octane for the W4 and W7 force fields are in good agreement with the experimental results.43 However, while the estimated critical temperatures for n-butane and n-pentane are too low, that for n-octane is slightly too high. W7-LC performs much less satisfactorily for all chain lengths; that is, tail corrections should always be used with the Williams force field. Here it should be noted that the energy parameters for the Williams force field were fitted using relatively short potential truncations, but assuming that the truncated potentials would account for only 80% of the heat of sublimation.9 OPLS-AA Force Field. Simulations were carried out for two versions of the OPLS-AA force field,6 which differed only with respect to the use of partial charges (see also Table 2). Jorgensen and co-workers6 reported that the partial charges have a negligible effect on the liquid densities and heats of vaporization of pure alkanes. Thus we would like to investigate whether the same holds true for the fluid-phase behavior. The VLCC for the OPLS-AA force field are shown in Figure 6, the corresponding numerical data are listed in Table S2, and the estimated critical data and standard boiling temperatures are given in Table 3. As is evident from Figure 6 and Table 3, no significant difference was found between simulations for the OPLS-AA force field with (OA) and without (OA-q ) partial charges. For both vapor and liquid phases Coulombic interactions contribute less than 0.2% to the total nonbonded interactions; that is, the structure of the normal alkanes leads to an effective cancellation of the Coulombic interactions. A similar result has also been observed for perfluorinated alkanes, where the partial charges are much larger in magnitude.44 The OPLS-AA force field yields a very good prediction of the critical and boiling temperatures of methane. However, the saturated liquid densities and the critical density are overestimated by 10% and 5%, respectively. For all other n-alkanes, the saturated liquid densities agree well with the experimental data for the lower temperatures, which were used in the fitting of the OPLS-AA force field.6 In contrast, the critical temperatures (and correspondingly also the liquid densities at higher temperatures) are underestimated by as much as 6% for ethane and 2% for n-octane. It might be of interest to note that our calculations predict 270 K as the standard boiling temperature for n-butane, which is 3 K lower than the experimental result. Thus the stable phase of OPLS-AA n-butane at the experimental boiling point is the vapor phase and not the liquid phase (see ref 5). Two additional tests were carried out for the OPLS-AA force field. First, we have tested whether the limited flexibility used in our simulations might be the cause for the disagreement in the VLCC of methane. To this extent, simulations were carried out for a methane model with flexible H-C-H bond angles governed by a harmonic bending potential (see eq 1) with kθ/



TABLE 3: Thermodynamic Properties of n-alkanes for the Williams, OPLS-AA, and MMFF94 All-Atom Force Fields (Details of the Force Fields As Described in Table 2) and the TraPPE United-Atom Force Field. Tb, Tc, Gc, and Pc Are the Normal Boiling Point, the Critical Temperature, the Critical Density, and the Critical Pressure, respectively. The Experimental Data Are Taken from Ref 42. The Subscripts Give the Statistical Accuracy of the Last Decimal(s) W4

W7

W7-LC

Tb/K Tc/K Fc/(g cm)-3 Pc/MPa

1001 1722 0.1656 4.08

1022 1802 0.1625 4.04

951 1571 0.1726 3.96


1673 2744 0.2059 4.713

1692 2774 0.19510 3.812


2673 4105 0.23211 4.115


3012 4634 0.23311 3.57


OA-q

MM

MM-LC

TraPPE

expt

Methane 1141 1912 0.1706 4.69

1131 1912 0.1696 4.49

791 1232 0.1528 2.811

721 1142 0.1469 2.04

1121 1911 0.1603 4.511

112 191 0.162 4.6

1521 2413 0.2138 4.55

n-Ethane 1782 2852 0.2027 4.713

1732 2862 0.2057 4.211

1221 1821 0.2116 3.37

1151 1692 0.19410 2.910

1771 3042 0.2063 5.14

185 305 0.205 4.9

2662 4114 0.23510 3.99

2224 3455 0.23413 3.417

n-Butane 2702 4095 0.22411 3.89

2681 4074 0.22610 3.94

1852 2683 0.22313 2.410

1673 2383 0.24312 2.517

2612 4234 0.2316 4.14

273 425 0.228 3.8

3015 4614 0.23812 3.415

2614 3775 0.24216 3.523

n-Pentane 3014 4585 0.22514 3.313

3025 4564 0.22612 3.317

2103 2943 0.22713 2.110

1955 2623 0.23515 2.020

2962 4702 0.2384 3.71

309 470 0.230 3.4

n-Octane 3945 5588 0.23316 2.917

3992 5546 0.21714 2.54

3862 5702 0.2392 2.61

399 569 0.232 2.5

4023 57611 0.24521 2.911

OA

kB ) 33 234 K rad-2. In this case, hydrogens take part in CBMC conformational changes and swap moves. The VLCC for this model is also shown in Figure 6, and it is obvious that this added H-C-H flexibility has little effect on the phase behavior of methane. Secondly, we have investigated the influence of the potential truncation. Jorgensen and co-workers5,6,36 used a spherical potential truncations at 11 Å, i.e., 2 Å larger than the value used here. A set of calculations was carried out for n-pentane using rcut ) 12 Å. The change in potential truncation was found to cause no significant change in the VLCC (see Figure 6B) and thermodynamic properties (see Table 4). Thus, the analytical tail corrections do adequately account for the interactions in the range from 9 to 12 Å (see also discussion of the liquid structures). MMFF94 Force Field. Simulations were carried out for the MMFF94 force field with and without tail corrections (see eqs 12 and 13). The calculated VLCC are shown in Figure 7, the corresponding numerical data are listed in Table S3, and the estimated critical data and standard boiling temperatures are given in Table 3. As can be expected, the tail corrections for the buffered 14-7 potential have a smaller, albeit still significant, effect than those for the longer ranged Lennard-Jones 12-6 and Buckingham exp-6 potentials (see Table 3). The disagreement between VLCC obtained for the MMFF94 force field and the experimental results is striking. The predicted critical temperature for ethane is lower than the experimental result for methane, and the calculated VLCC for n-butane and n-pentane fall below the experimental VLCC for ethane. For example, T ) -0.5 °C, the experimental boiling point of n-butane which was used by Kaminski and Jorgensen in their comparison of all-atom force fields, is substantially above the critical temperature of MM-LC butane. What causes the dramatic failure of the MMFF94 force field in describing the fluid-phase behavior of alkanes? Kaminski

and Jorgensen5 point to the different ratios of the C-C, C-H, and H-H well depths between OPLS-AA and MMFF94. However, the corresponding ratios differ by an even larger extent for the Williams VII parameter set (see Figure 1), which performs rather well. Thus there seems to be no magic ratio for the interaction strengths (see also next section). OPLS-AA and MMFF94 differ also to a large degree in C-C bond length and C-C-C bond angle (see Table 1). To evaluate the influence of the intramolecular structure on the phase behavior, a set of calculations (MMOPLS) was carried out for n-pentane using the MMFF94 parameters for the nonbonded interactions and the OPLS-AA parameters for the intramolecular structure. The longer bonds and wider angles of OPLS-AA lead to a larger molar volume (lower density) for MMOPLS, but the boiling and critical temperatures remain unchanged (see Figure 7B and Table 4). Since neither the ratio between interaction strengths nor the intramolecular structure seems to be responsible for the disagreement between MMFF94 and experiment, could the problem be with the rare gas pair potentials that are the starting point in the development of the MMFF94 force field? Figure 8 shows a comparison of the experimental phase diagram for argon and the calculated phase diagram using the buffered 14-7 potential with the parameters listed in ref 10. In striking contrast to the VLCC for the linear alkanes, the VLCC for argon is shifted to higher temperatures and the buffered 14-7 potential yields an overestimate of the critical temperature by approximately 8%. However, this is exactly what we would like to see for a true pair potential because adding the unfavorable three-body Axilrod-Teller interactions45,46 should shift the calculated coexistence curve into good agreement with the experimental results. Thus judging from this test for argon, we feel that the rare gas pair potentials are not responsible for the failure of the MMFF94 force field. Here it should also be noted that the original buffered 14-7 parameters for hydrogen, which were

2584 J. Phys. Chem. B, Vol. 102, No. 14, 1998

Chen et al.

Figure 6. Vapor-liquid coexistence curves for variations of the OPLSAA force field: (A) methane and ethane; (B) n-butane, n-pentane, and n-octane. Experimental data are shown as in Figure 6. The calculated coexistence densities and critical points are shown as diamonds and squares for the OA and OA-q force fields, respectively. Triangles are used to show the results for methane using the fully flexible OA model and for n-pentane using rcut ) 12 Å.

Figure 7. Vapor-liquid coexistence curves for variations of the MMFF94 force field: (A) methane and ethane; (B) n-butane and n-pentane. Experimental data are shown as in Figure 6. The calculated coexistence densities and critical points are shown as diamonds, squares, and triangles for MM, MM-LC, and MMOPLS, respectively.

TABLE 4: Thermodynamic Properties of n-Pentane for the OPLS-AA Force Field Using Two Different Potential Truncations and for the MMFF94 Force Field Using Two Different Intramolecular Structures. The Subscripts Give the Statistical Accuracy of the Last Decimal(s) Tb/K Tc/K Fc/(g cm)-3 Pc/MPa

OA (rcut ) 9 Å)

OA (rcut ) 12 Å)

MM

MMOPLS

3014 4585 0.22514 3.313

2964 4546 0.23014 3.414

2103 2943 0.22713 2.110

2141 2904 0.22414 2.36

derived by Halgren10 using an estimated atomic polarizability and Slater-Kirkwood effective number of electrons, have a well depth that is approximately 25% smaller in magnitude than the final MMFF94 parameters, which were obtained using ab initio data.8 Chain-Length Dependence of the Thermodynamic Properties. Comparisons between the experimentally observed critical temperatures, Tc, critical densities, Fc, and standard boiling points, Tb, and the corresponding calculated quantities for the three force fields are shown as functions of chain length in Figures 9-11. The MMFF94 force field (with tail corrections) underestimates Tc and Tb by approximately 30%. The Williams VII force field underestimates Tc of methane and of ethane by 5% and 10%, respectively. For the medium-length alkanes the agreement is better, but there is obviously a trend in Figures 9 and 11 that suggests that the Williams force field will yield too

Figure 8. Comparison of the vapor-liquid coexistence curve and critical points for argon obtained for the buffered 14-7 potential (dotted lines and triangles) and from experiment (solid lines and filled circle).47 The dotted lines are fits to the scaling law and the law of rectilinear diameters.

high Tc and Tb for longer alkanes (nC > 8 ). Thus the interactions of methylene segments are too strong in the Williams force field. The OPLS-AA force field gives very good agreement for Tc and Tb for methane, but underestimates Tc for


Figure 9. Ratios of simulated versus experimental critical temperatures. The following symbols are used for the different force fields: W7 (open diamonds), W7-LC (filled diamonds), OA (open squares), OA-q (filled squares), MM (open triangles), and MM-LC (filled triangles). For comparison, the results for the TraPPE united-atom force field are shown as crosses.

Figure 10. Ratios of simulated versus experimental critical densities. Symbols as in Figure 8.

Figure 11. Ratios of simulated versus experimental normal boiling temperatures. Symbols as in Figure 8.

ethane. The agreement between experiment and simulation seems to improve with increasing length of the alkanes, which suggests that the OPLS-AA force field accurately reproduces the strength of the methylene interactions. The results for Fc are harder to interpret, since there is more scatter and the relative error on calculated Fc data is much larger than for the


Figure 12. Liquid-phase site-site radial distribution functions for n-pentane at 360 K (W4, W7, and OA) and at 220 K (MM). Solid, dashed, dotted, and dashed-dotted lines are used for the W-4, W-7, OA, and MM force fields, respectively. The radial distribution functions for C-H and C-C pairs are vertically shifted by 1 and 2 units.

temperatures. For comparison, Tc, Fc, and Tb calculated for the TraPPE force field,37 a simpler united-atom force field with Lennard-Jones parameters fitted to VLCC for alkanes, are also included in Table 3 and Figures 9-11. The TraPPE force field yields better results for Tc (but this quantity was used in the fitting procedure), whereas it performs with similar accuracy as the OPLS-AA and Williams force fields for Fc and Tb. However, because of the more accurate prediction of Tc, the TraPPE force field yields better overall agreement with the experimental VLCC.37 Thus at present, use of the computationally more demanding OPLS-AA and Williams all-atom force fields does not yield gains in accuracy for the prediction of fluidphase equilibria of alkanes. Here we should address two important observations. All three force fields give larger deviations in Tc and Tb for ethane than for methane or butane. This observation suggests that different nonbonded vdW parameters are required for the carbon atom in methane, in a methyl group, or in a methylene group for an all-atom model to yield a quantitative description of thermodynamic properties; that is, the approximation that the vdW parameters of an atom do not depend on the type of its bonded neighbors6-10 is most likely not valid. The C-C well depth of the Williams force field is much deeper than that of OPLSAA, while the H-H well is 2 times shallower (see Figure 1). Considering the great differences between these two force fields, it might be a surprise that they give relatively similar results. From this we expect that there will be several combinations of carbon and hydrogen vdW parameters that may work for allatom alkane models, particularly if we allow different carbon parameters for different groups. Thus the empirical fitting of vdW parameters for all-atom models will be a very difficult task. B. Liquid Structures. Radial distribution functions (RDFs) of liquid-phase n-pentane are compared for the three force fields in Figure 12. Since the VLCC for the MMFF94 force field is very different, the RDFs are compared for temperatures that yield similar liquid densities. At 360 K, the saturated liquid densities for W4, W7, and OA are 0.556, 0.556, and 0.552 g cm-3, and that for MM at 220 K is 0.567 g cm-3. Regarding the liquid structure, the most obvious result is that there is very little difference between the three force fields. RDFs were also

2586 J. Phys. Chem. B, Vol. 102, No. 14, 1998 calculated for other alkanes and other sets of temperatures, but the differences between the RDFs were always very small. Thus neither the different functional forms of the potential, the slightly different positions of the well depths (see Figure 1), nor the different bond lengths are sufficient to cause major differences in liquid structure (if the comparison is made for similar liquid densities, i.e. different temperatures for different models). Therefore, we feel that RDFs or the corresponding radially averaged structure factors are probably not a good means to compare force fields for alkanes. Finally, at the temperatures used in this work there is very little structure beyond the first peak, and no peak is observed around 9 Å, the range of our spherical potential truncation. V. Conclusions In this work, the performance of three all-atom force fields for n-alkanes is assessed. While the Williams force field10 is specific to n-alkanes, the OPLS-AA6 and MMFF948 force fields are designed for more general purposes. Configurational-bias Monte Carlo calculations in the Gibbs ensemble were used to determine the vapor-liquid coexistence curves of methane, ethane, n-butane, n-pentane, and n-octane. It is clearly evident that the MMFF94 force field does not yield a good description of the thermodynamic properties of the alkanes. The critical temperatures are best reproduced by the OPLS-AA force field, and in particular the agreement improves for longer alkanes, whereas the Williams force field yields slightly better results for the saturated liquid densities, in particular for methane. There is no significant difference in fluid-phase behavior for the OPLSAA model with and without partial charges. All three force fields show larger deviations from the experimental results for ethane than for methane or butane, which suggests that specific van der Waals parameters for carbon atoms belonging to methane, methyl, or methylene groups may be required. The large differences between the Williams and OPLS-AA van der Waals parameters combined with their similar fluid-phase behavior illustrate that several combinations of parameters may work for all-atom alkane models. Finally, the liquid structure is not very useful to distinguish between the three alkane force fields. Acknowledgment. We would like to acknowledge many stimulating discussions with Doug Tobias, Mike Klein, JeanPaul Ryckaert, and Ian McDonald. Financial support from the Petroleum Research Fund, administered by the American Chemical Society (Grant No. 29960-AC9), through a Camille and Henry Dreyfus New Faculty Award, and through a McKnight/Land-Grant Fellowship is gratefully acknowledged. M.G.M. would like to thank the Graduate School, University of Minnesota, for the award of a Graduate School Fellowship. Part of the computer resources were provided by the Minnesota Supercomputer Institute through the University of MinnesotaIBM Shared Research Project and NSF Grant CDA-9502979. Supporting Information Available: Tables S1, S2, and S3 list the numerical results (pressure, vapor and liquid densities) for the Williams, OPLS-AA, and MMFF94 force fields, respectively (6 pages). Ordering information is given on any current masthead page. References and Notes (1) Ryckaert, J. P.; Bellemans, A. Chem. Phys. Lett. 1975, 30, 123. (2) Ryckaert, J. P.; Bellemans, A. Faraday Discuss. Chem. Soc. 1978, 66, 95. (3) Ryckaert, J. P.; McDonald, I. R.; Klein, M. L. Mol. Phys. 1989, 67, 957.

Chen et al. (4) Moller, M. A.; Tildesley, D. J.; Kim, K. S.; Quirke, N. J. Chem. Phys. 1991, 94, 8390. (5) Kaminski, G.; Jorgensen, W. L. J. Phys. Chem. 1996, 100, 18010. (6) Jorgensen, W. L.; Maxwell, D. S.; Tirado-Rives, J. J. Am. Chem. Soc. 1996, 118, 11225. (7) Cornell, W. D.; Cieplak, P.; Bayly, C.; Gould, I. R.; Merz, K. M.; Ferguson, D. M.; Spellmeyer, D. C.; Fox, T.; Caldwell, J. W.; Kollman, P. A. J. Am. Chem. Soc. 1995, 117, 5179. (8) Halgren, T. A. J. Comput. Chem. 1996, 17, 490, 520, 553, 616. (9) Williams, D. E. J. Chem. Phys. 1967, 47, 4680. (10) Halgren, T. A. J. Am. Chem. Soc. 1992, 114, 7827. (11) Van der Ploeg, P.; Berendsen, A. J. Chem. Phys. 1978, 76, 3271. (12) Jorgensen, W. L.; Madura, J. D.; Swenson, C. J. J. Am. Chem. Soc. 1984, 106, 813. (13) Scott, R. A.; Scheraga, H. A. J. Chem. Phys. 1966, 44, 3954. (14) McQuarrie, D. A. Statistical Mechanics; Harper and Row: New York, 1976; p 140. (15) Kollman, P. A. Acc. Chem. Res. 1996, 29, 461. (16) Buckingham, R. A. Proc. R. Soc. London Ser. A 1938, 168, 264, 378. (17) Tobias, D. J.; Tu, K.; Klein, M. L. J. Chim. Phys. Phys.-Chim. Biol. 1997, 94, 1482. (18) Lennard-Jones, J. E. Proc. R. Soc. London Ser. A 1924, 106, 463. (19) Lorentz, H. A. Ann. Phys. 1881, 12, 127. (20) Berthelot, D. C. R. Hebd. Se´ anc. Acad. Sci., Paris 1898, 126, 1703. (21) For this specific buffered potential the exact minimum energy separation and well depth are 0.996r*ij and -1.0006ij, respectively. (22) Panagiotopoulos, A. Z. Mol. Phys. 1987, 61, 813. (23) Panagiotopoulos, A. Z.; Quirke, N.; Stapleton, M.; Tildesley, D. J. Mol. Phys. 1988, 63, 527. (24) Smit, B.; de Smedt, P.; Frenkel, D. Mol. Phys. 1989, 68, 931. (25) Siepmann, J. I. Mol. Phys. 1990, 70, 1145. (26) Siepmann, J. I.; Frenkel, D. Mol. Phys. 1992, 75, 59. (27) Frenkel, D.; Mooij, G. C. A. M.; Smit, B. J. Phys.: Condens. Matter 1992, 4, 3053. (28) de Pablo, J. J.; Laso, M.; Suter, U. W. J. Chem. Phys. 1992, 96, 2395. (29) Siepmann, J. I. In Computer Simulation of Biomolecular Systems: Theoretical and Experimental Applications; van Gunsteren, W. F., Weiner, P. K., Wilkinson, A. J., Eds.; ESCOM Science: Leiden, 1993; Vol. 2, pp 249-264. (30) An upper limit of 2rcut ) 18 Å and π was used for the maximum translational and rotational displacements. These upper limits are usually reached for low-temperature/low-density vapor phases, and the acceptance rates can thus be higher than 50% in these cases. Using the OPLS model for n-pentane at 400 K, the maximum translational and rotational displacements in the liquid phase were 0.8 Å and 0.5 rad, respectively, the maximum translational displacement in the vapor phase was 15 Å, and the the acceptance rate for rotations in the vapor phase (using the upper limit for the maximum displacement) was 75%. The maximum volume displacement for this system was 700 Å3. (31) Esselink, K.; Loyens, L. D. J. C.; Smit, B. Phys. ReV. E 1995, 51, 1560. (32) Mackie, A. D.; Tavitian, B.; Boutin, A.; Fuchs, A. H. Mol. Sim. 1997, 19, 1. (33) Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids; Oxford University Press: Oxford, 1987. (34) Smit, B.; Karaborni, S.; Siepmann, J. I. J. Chem. Phys. 1995, 102, 2126. (35) Vlugt, T. J. H.; Martin, M. G.; Smit, B.; Siepmann, J. I. Mol. Phys., in press. (36) In contrast, Kaminski and Jorgensen5 employed a potential truncation at 11 Å for n-butane, which is molecule-based and is determined solely on the separation between C2 atoms. In addition, the potentials were smoothed over the last 0.5 Å inside the truncation. (37) Martin, M. G.; Siepmann, J. I. J. Phys. Chem. 1998, 102, 2569. (38) Kaminski and Jorgensen5 did not include tail corrections for the MMFF94 force field in their study. (39) All simulations were carried out for 500 particles with rcut ) 3.75σ using the appropriate tail corrections. (40) Yin, D.; MacKerrell, A. D. J. Phys. Chem. 1996, 100, 2588. (41) Martin M. G.; Siepmann, J. I. J. Am. Chem. Soc. 1997, 119, 8921. (42) Teja, A. S.; Lee, R. J.; Rosenthal, D.; Anselme, M. Fluid Phase Equilib. 1990, 56, 153. (43) Smith, B. D.; Srivastava, R. Thermodynamic Data for Pure Compounds: Part A Hydrocarbons and Ketones; Elsevier: Amsterdam, 1986. (44) Chen, B.; Siepmann, J. I. Unpublished results. (45) Axilrod, B. M.; Teller, E. J. Chem. Phys. 1943, 11, 299. (46) Stenschke, H. J. Chem. Phys. 1994, 100, 4704. (47) Younglove, B. A. J. Phys. Chem. Ref. Data 1982, 11, 12.

Thermodynamic Properties of the Williams, OPLS-AA, and MMFF94 All

Recommend Documents