J. Phys. Chem. 1995, 99, 13868-13875
13868
Critical Evaluation of Benzene Analytical Nonbonded Force Fields. Reparametrization of the MM3 Potential Jeno NagyJ Vedene H. Smith, Jr.,*,t and Donald F. Weaver"* Departments of Chemistry and Medicine, Queen's University, Kingston K7L 3N6, Ontario, Canada Received: March 28, 1995; In Final Form: July 7, 1995@
An extensive characterization of six existing nonbonded analytical 12-site potential functions for benzene has been performed by computing various properties in all three phases. The computed properties include gas-phase benzene dimer energies and structures, benzene second-virial coefficients, liquid benzene densities and vaporization enthalpies, liquid benzene C-C and ring center-center radial correlation as well as orientational functions, crystal unit cell parameters, unit cell volumes, and heat of sublimation data for six aromatic compounds. The 12-6 potential of Jorgensen and Severance emerged as the best effective twobody potential in condensed phases, while the exp-6 potentials of Williams and Starr proved to be the most consistent in all three phases. Two new parameter sets have been obtained in the framework of the MM3 potential, and one of them represents significant improvement over the original MM3 parameters and is comparable in quality to the force field of Williams and Starr. An analysis of the 12-6 and exp-6 potential forms is given, and the question of whether point charges or point dipoles should be used and the relationships between properties in different phases have also been examined.
1. Introduction Being pivotal to an understanding of more complex aromatic systems, benzene-benzene interactions have been the subject of considerable research effort at the theoretical level in the l i q ~ i d , ~ and ~ - ~solid ' states.2,'3.'4,22.30.32-39Owing to the great computational demand of ab initio methods, empirical potential functions have routinely been employed in the characterization of these interactions in molecular mechanics, molecular dyanmics, and Monte Carlo calculation^.^^-^'^^^ Although angle-dependent potentials between molecules have also been used,',6 a more elaborate description of the interaction energy is based on the site-site pair interaction formalism.40It has been well documented that description of benzene properties is unsatisfactory by 6-site modelsi4,24,25,29,31 or by models with no provision for electrostatic interaction.^.*,^^^^^,^^-^^ Site-site anisotropy can be incorporated via either bond shortening factors32,37,38,41 or directly using site-site anisotropic potent i a l ~ . ~Electrostatic ~ , ~ ~ , ~ interactions ~ can be represented in the form of single multipoles placed at the center of the benzene rings2,25,33or distributed multipoles on the interaction sites.3,'0-22.30,44 The former approach has been found inade q ~ a t eand , ~ in the latter the expansion is usually terminated at the point charges.4.5,9.'2.31,35,37.39.45 (Altematively, bond dipoles can also be used.I3) Even though the distributed point charge model can give qualitatively incorrect results in certain cases,10.44 it is widely used due to its simplicity.4.5.9.'2.31,35,37,39,45 These charges are either set to reproduce the experimental quadrupole moment of b e n ~ e n e , obtained ~ from ab initio population a n a l y ~ i s , ' or ~ . ~simply ~ adjusted to yield the best agreement between various calculated and experimental proper tie^.^'.^^ The motivation behind this study is to rigorously evaluate six commonly used analytical benzene-benzene force fields with transferable The condition of transferability is an important attribute since it renders the force fields applicable to aromatic compounds other than ben~ e n e . ~ ' , ~As~ observed % ~ ~ -by~ Shii ~ . ~and ~ Bartell,I4 only the Department of Chemistry.
* Department of Medicine. +
@Abstractpublished in Advance ACS Abstracts, August 15, 1995.
potential of Karlstrom et al. is not comprised of transferable parameters; however, it was included here due to its importance among benzene potentials and to address the need for further elucidation of its properties despite earlier extensive r e ~ e a r c h . ~ . ~ ~ All force fields are also 12-site (all atom) models incorporating point charges, except for MM3, in which charges are supplanted by bond dip01es.I~ For the latter, equivalent charges have also been used to calculate the same properties, thereby making direct comparison possible between the performance of the two charge distributions. Since the goodness of a force field can be measured by its ability to reproduce experimentally and computationally available data in all three phases with at least reasonable accuracy, the following benzene properties have been computed in the current study: gas-phase structures and energies of benzene dimers, second-virial coefficients, liquid state densities, vaporization enthalpies, liquid benzene C-C and ring center-center radial correlation as well as orientational functions, crystal unit cell parameters, unit cell volumes, and heat of sublimation data for six aromatic compounds. The recent theoretical work or Hobza et al. on the benzene dimer2Oq2I provides an improved standard for evaluating the accuracy of analytical potentials, thus improving the relevance of the current study. In addition to the assessment of the six available force fields, two new parameter sets within the MM3 potential framework have also been obtained on the basis of different criteria. MM3 has been chosen for this purpose since it possesses a very flexible two-parameter potential f u n ~ t i o n . ~ ~ , ~ ~ , ~ ~ Furthermore, as shown below, the MM3 parameter set obtained by Allinger and LiiI3 is inadequate for the characterization of the above properties. Although several of the above properties have been computed using the six previously published benzenebenzene force fields, they were recalculated in this study to reflect certain biases and approximations inherent in our procedures. We feel that this collective treatment will render the force fields more conducive to stringent comparison. In addition, similarities and differences between our results and those from other works are analyzed. Wherever possible, comparisons are made with other, previously available force
0022-3654/95/2099-13868$09.00/0 0 1995 American Chemical Society
Benzene Analytical Nonbonded Force Fields fields as well. General trends across different phases are also identified.
2. Models and Parametrization Procedures The six all-atom benzene force fields analyzed in the current work are as follows (charges refer to hydrogens, with carbons bearing the opposite charge): (1) the empirical potential of Williams and Stan fitted to 18 hydrocarbon crystals with optimized charges of 0.153,37henceforth referred to as WS77 following the notation of Shii and Bartell;I4 (2) the ab initio potential of Karlstrom et al. with charges of 0.15 and denoted here as KLWJ83;5(3) the potential of Weiner et al. with charges of 0.1545 (since this potential is incorporated in the program it will be referred to as AMBERb'); (4)MM241 with the added charges of qH = 0.15, henceforth referred to as MM2b;IZ ( 5 ) MM3 utilizing bond dipoles of 0.6 D;13338and (6) the potential of Jorgensen and Severance (termed JS90) with qH = 0.115 fitted to benzene liquid proper tie^.^' With the exception of the last two potentials, the charges give rise to the experimentally observed quadrupole moment of %-30 x C m2 of b e n ~ e n eand ~ ~are . ~approximately ~ equal to the charges derived from infrared band densities of b e n ~ e n e . ~Considering ' the dielectric constant values of 1.O and 1.5 associated with the Coulomb and Jeans formulas, re~pectively,~' and also adopted in MM3,I3 the 0.6 D in MM3 corresponds to qH = 0.0923. (This qH value is not consistent with the statement of Allinger and Lii that their point dipoles are "nearly identical" to the charges of 0.15 used by Pettersson and Li1Jef0rs.I~)To establish which of the two charge distributions was superior, the corresponding charges of 0.0923 were also used together with the MM3 van der Waals parameters (termed MM3c) to calculate all properties, except crystal data. Three of the potentials, WS77, MM2b, and MM3c, are of the exp-6 form, while the other three use the 12-6 function; the KLWJ83 potential is a more sophisticated and softer variant of the original Lennard-Jones functional form. Previous calculations with the existing force fields include crystal lattice3' and selected aromatic dimer calculations for the WS77 p ~ t e n t i a l more ; ~ recently, they showed the superiority of electrostatic potential derived charges52to Milliken charges in predicting crystal structure^.^^ WS77 was also used in a liquid Monte Carlo simulation in the NFT ensemble.30 The KLWJ83 potential is characterized by its second vinal coeff i c i e n t ~and ~ its liquid and solid properties.26 The only reason, apart from its importance, that this potential was included in our study is the somewhat unrealistic volume at which Linse performed his Monte Carlo simulations in the NVT ensemble, resulting in presumably too low internal liquid energy,26 as pointed out by Bartell et al.I4 Thorough and rigorous condensed phase analyses on these two potentials have also been performed by the Bartell g r ~ u p . ' ~They , ~ ~not only rescaled the original potentials to reconcile them with benzene crystal data at 0 K but also simplified the potential forms to make their liquid simulations computationally more fea~ib1e.I~ The remaining four force fields have also been used in dimer energy calculat i o n ~ . ' * ~ 'In~ addition, ~ ~ ' . ~ ~three crystal structures for MM3'3938 and liquid properties for the potential of Jorgensen et al. have also been ~ o m p u t e d . ~ ' Three different procedures were used to obtain new parameter sets for the MM3 potential form. In the first, parameters including point charges were fitted to 21 benzene dimer energies calculated at the second-order Mldller-Plesset perturbation theory level of Hobza et d.20.21 (see MM3bener in Table 1). As shown in Table 1, the van der Waals radius of carbon is too small, whereas that of hydrogen is too large. The optimum q H value was found to be 0.1 18e. Despite this dimer energy data
J. Phys. Chem., Vol. 99, No. 38, 1995 13869 TABLE 1: Two MM3 Benzene Parameter Sets Obtained in This Studv force field
parametersa
MM3bener MM3bmccr
0.069, 0.0380,0.021, 3.48, 3.39, 3.30,0.118
a
0.057, 0.0425,0.032, 3.85, 3.50, 3.13,0.111
Parameters are in the order 6c-c; E C - H ;
= a#,; ( A ~ ( R= ~ ~aj$; ) ) ( A ~ ~ ( R=~aiaj ~ ) ) (4) where Rij is the unit vector between the two benzene centers, and ai and aj are the c6 symmetry unit vectors of molecules i and j , respectively. Static crystal lattice summations using the PEW program55 to obtain cell parameters and heat of sublimation values were performed for the following six compounds: benzene,60naphthalene?’ ~henanthrene?~ biphenyl,@and p e ~ y l e n e . ~ ~ No zero-point corrections to the heat of formations were considered, except for benzene, for which the total value of 12.50 kcal/mol is widely regarded as the most reliable estimate.22 In the lattice energy summations, 125 unit cells with 50 8, cutoff values were used to attain proper convergence. The relationship between the heat of sublimation and the lattice energy of the crystal is given AHsub
= Egas -
1 F Xtd - 2RT
(5)
where E,,, and EXbl denote the steric energy of the isolated molecule and the crystal, respectively. Z is the number of molecules per unit cell, and the correction factor 2RT allows for pV-work and the difference between classical rotation and translation of the gas and vibration of the ~ r y s t a l . ~ ~ , ~ ’ Uncertainties in all calculations are estimated to be under 1-2%. The only exceptions are crystal heats of sublimation other than that of benzene, in which zero-point vibrations are not allowed for. 4. Results and Discussion
Dimer energies for the four configurations in Figure 1 as wellas the global minima calculated from the nine force field equations including MM3c are listed in Table 2. According to the ab initio resu1ts;O the most stable dimer, 3, is of a (skewed) parallel-displaced configuration followed by a T-shaped (vertex
J. Phys. Chem., Vol. 99, No. 38, 1995 13871
Benzene Analytical Nonbonded Force Fields TABLE 3: Calculated Benzene Second-Virial Coefficients (cm3/mol)at Five Temperatures
TABLE 4: Calculated Liquid Benzene Densities and Enthalpies of Vaporization at 298.15 K
temp/K 348.25 exptl -981.0 -925.3 WS71 -987.2 KLWJ83 AMBERb+ -1004.9 MM2b -1094.8 -729.8 MM3 -807.3 JS90 MM3bener -586.5 MM3bmccr -920.1
397.59 446.98
498.54
548.84
599.82
-724.6 -670.8 -725.8 -710.4 -793.4 -539.7 -603.8 -420.5 -670.8
-429.2 -393.8 -429.2 -403.5 -476.8 -323.7 -372.1 -239.0 -399.2
-344.5 -312.5 -349.0 -316.9 -384.9 -257.2 -301.8 -184.4 -318.8
-282.1 -252.6 -282.6 -253.1 -316.9 -207.7 -247.6 -141.9 -259.3
-552.6 -508.9 -553.2 -529.4 -606.1 -415.4 -468.6 -315.9 -512.5
into ring, 1) and a displaced T-shaped configuration (4), respectively, while another T-shaped arrangement (2) has appreciably lower stability. The KLWJ83 potential seriously mispredicts the order of configurations 2 and 3 presumably due to the absence of these arrangements in the ab initio calculations of Karlstrom et aZS5In MM2b, their order is correct, but 3 still has too repulsive an energy with respect to configurations 1 and 4, which are in turn too negative. AMBERb+ has a similar problem, although to a smaller degree. MM3 and MM3c give rise to stabilization energies that are too small, but with the exception of 1 they predict correct ordering of the configurations. As expected, the results for the two force fields agree very closely. Although none exhibit entirely correct ordering, the remainder of the potentials are more consistent in predicting the stability of the four configurations. In light of the new finding of Hobza et al. on the parallel-displaced minimum point on the dimer potential energy surface,20the SW77, AMBERb’, JS90, and K w w 8 3 force fields yield incorrect minimum energy configurations: herringbone structures in the case of the former two and slightly distorted T-shaped structures for the latter two, in agreement with earlier The significance of the minimum energy structures of these four force fields should not be overemphasized, however, given the shallow potential energy surface of the benzene dimer. The two newly parametrized force fields presented in this study offer pleasing consistency with the ab initio data, although configuration 3 is slightly less stable than desired for MM3bener, and the opposite is true for MM3bmccr. Computed and experimental6* benzene second-virial coefficients at five temperatures are given in Table 3. Clearly evident is an excellent agreement between the KLWJ83 values calculated in this study and experimental data, which is in contrast with the finding of Karlstrom et aL5 whose obtained B(T) values were 25% too After cross-checking our results using two distinct integration techniques, we are confident that our results are correct. (The integration method used by Karlstrom et al. was not mentioned in their work.5) MM2b is too attractive, and MM3 is too repulsive, in line not only with their dimer energies but with alkane properties derived from these two potentials employing similar parameters as ~ e 1 1 . Expectedly, ~ ~ , ~ ~ the MM3 and MM3c results are, again, very close to each other. The AMBERb+ results at low temperatures are promising, but the repulsive wall of the potential starts to dominate at higher temperatures, suggesting that in condensed phases AMBERb+ will substantially overestimate the AHvapof the liquid and A H s u b of crystals. The hardness of the 12-6 potential form alone cannot be responsible for this since the JS90 results are more in parallel to the experimental curve (even though the discrepancy for JS90 is substantial). WS77 and MM3bmccr provide reasonable results with slight underestimation of the experimental data, whereas MM3bener gives rise to unrealistically positive B(T) values, which indicates that the ab initio energies are too shallow.
exptl WS17 KLWJ83 AMBERb+ MM2b MM3 JS90 MM3bener MM3bmccr
AHvap(kcal/mol)
density (g/cm3)
8.09 8.40 10.67 9.53 10.07 6.83 8.06 6.35 8.30
0.874 0.849 1.062 0.921 0.966 0.808 0.872 0.832 0.867
Calculated liquid benzene densities and heats of vaporization at 298.15 are summarized in Table 4. KLWJ83, AMBERb+, and MM2b substantially overestimate the attraction between molecules, while MM3 and MM3bener underestimate it. MM3c once again yields the same results as MM3 within the estimated error limits. It is to be noted that the order between MM2b and KLWJ83 has been reversed in comparison to their virial coefficients. This phenomenon can be attributed partly to the effect the surrounding molecules in the bulk phase exert on a given benzene dimer. This “quasi many-body’’ effect is not included in the second-virial coefficients, which is an integrated property of two-body interactions alone. The problem with the KLWJ83 potential is its extremely high density, or small atomic radii, that allows the molecules to get into very close contact with each other in the liquid (as well as in the solid). The density for MM2 is also very high, but is substantially closer to the experimental value. Linse computed the properties of liquid benzene by means of an NTV Monte Carlo simulation, in which the density was fixed to the experimental value, effectively creating negative pressure in the system,29 which leads to an approximate AHvapof 8.84 kcal/mo126vs our 10.67 kcaymol. A typical difference in the characteristics of the 12-6 and exp-6 potentials can be discerned by examining the virial coefficients and AHvapdata for the AMBERb+-MM2b or JS90-WS77hlM3bmccr pairs of potentials (or any other combinations of 12-6 and exp-6 potentials). While the 12-6 potentials yield relatively more repulsive virial coefficients than their exp-6 counterparts, this difference either almost disappears or even reverses in the condensed phases. This behavior can be ascribed to the steep repulsive wall of the 12-6 functional form, which plays a diminishing role with decreasing temperature as the molecules spend more time in the potential well at lower temperatures. The same conclusion was drawn in a previous study from alkane virial coefficients and condensed phase proper tie^.^^ The best result in Table 3 is given by JS90, which was actually fitted to benzene liquid data.3’ WS77 and MM3bmccr are also in reasonable agreement with experiment, the latter giving somewhat better results. The former potential was also used in the NPT simulation of Yashonath et al., but their AHvapresult is 5% higher, and their volume/molecule is 6% lower.30 The difference may originate from the smaller size of their system and possible differences in the treatment of longrange corrections, which were not explicitly defined in their paper.30 Also, a slight discrepancy must result from their Rc-c = 1.393 A instead of 1.397 8, in the original force field.4,37 According to their calculations,their six-site ansiotropic potential yielded good agreement between experimental and computed internal energies of liquid benzene at 300 K, but it underestimated the density of the Increasing the density in line with experiment would, however, certainly shift the internal energy farther away from the experimental value. When going from the gas phase to condensed phases, the effect of many-body interactions has to be considered. For
Nagy et al.
13872 J. Phys. Chem., Vol. 99, No. 38, 1995
3.00
i
.*I,
.
t
AAAAA
MM3bmccr
2.00
1 .oo
JS90 KLWJ83 *****AMBER
6
A A A A
00 0 00
0.00 3.00
I
5.00
~
~
~
7.00
'
~
~
~
~
9.00
~
I
~
"
~
~
'
'
11.00
R/Angstrom
Figure 2. Computed center-center radial distribution functions for liquid benzene. For some force fields, abscissas are shifted by one unit for clarity.
2.00 9
* * MM2b
1.50 -
4
"
l?0
0
~
a a
ws77 0 ~ JSSO
KLWJ83 *****AMBER 0 0 000
0.50
- EXD
2.00
4.00
6.00
8.00
10.00 R/An gstro m
Figure 3. Computed C-C radial distribution functions for liquid benzene. For some force fields, abscissas are shifted by one unit for clarity. apolar systems, the leading three-body terms are the tripledipole dispersion energy and the first- and second-order exchange repulsion as inferred from calculations on the H270and CH471trimers, respectively. The overall effect of the nonadditive terms is destabilizing and is similar to atomic system^.^^,^^ Evans and Watts calculated the anisotropic threebody triple -dipole and quadrupole-induceddipole-quadrupole interactions for benzene and found that the former contributes 8-12% to the total two-body energy, while the effect of the latter is negligible.22 Considering the WS77 and MM3bmccr potentials, which give the overall best description of the gas and liquid phases, the estimated total three-body (repulsive) contribution to the liquid is in the 8-10% range of the twobody interactions and about 5-6% in the benzene crystal at 0 K (see below). Plots of the benzene center-center radial distribution functions are shown in Figure 2. Effects of too high stabilization in the liquid for MM2 and KLWJ83, and to a smaller degree
AME!ERb+, and of too low stabilization for MM3 and MM3bener are clearly visible. The JS90 and WS77 curves are almost identical, while the first peak of MM3bmccr is wider and flatter than the f i s t two. The f i s t shell integrates to 12.6 molecules, which is slightly higher than the experimentally measured 12,23 also predicted by WS77 and JS90.31 Another conspicuous feature of Figure 2 is the early ascent of the MM3 curves in relation to the others, which is attributable to the softer MM3 exchange repulsion exponent than that for the other potentials. The C-C radial distribution functions depicted in Figure 3 are rather structureless. The two peaks between 5 and 7 8, are generally well predicted by the force fields, but the plateau between 3 and 5 8, is only vaguely present in the plots, with the exception of IUWJ83, which predicts all structural features too early. The eight sets of average orientational functions in Figure 4 attest to strong qualitative similarities in the orientational structures of the liquids. Taking the JS90 force field as an
~
~
I
J. Phys. Chem., Vol. 99, No. 38, 1995 13873
Benzene Analytical Nonbonded Force Fields 1 .oo
1 .oo
0.80
0.80
3
0.60
\ \
\ \
0.60
I
0.40
I \
,
0.20 I , , , , , ,,, ,,, ,,, ,,, , 1 , , , / , , , , , 1 , , , , , , , , , , 3.60 5.d0 7.00 9.00 1 1 .oo
0.20
0.80
0.80
0.60
0.60
I
3
I
I
I
I
V
"
5.00
'
~
' ( 'I f ~ , ' , ~
7.00
I
I
I
I
,
I
I
I
I
,
,
9.00
7.60
I
,
,
I
/
,
,
,
,
I
,
1 1 .oo
,
,
,
R/Angstrom
1 .oo
jbr
5.00
AMBER, KLWJ83
1 .oo
3.00
r l l , , I , r , ( / , , , , , , , ,
R/Angstrom
JS90, WS77
0.20
1
3.00
'
'
"
I
9.00
f
'
'
'
~
'
~
1 1 .oo
'
J
0.20
r
1###
3.00
I
I
I
I
t
I
I
5.00
'
"
'
~
6
'
'
~
"
7.60
*
'
I
8
'
8
R/Angstrom
1 1 .oo
R/Angstrom
MM3bener, MM3bmccr Figure 4. Computed orientational functions in liquid benzene. Ai,
I ~ ' ' ~ ' I ' J ' I ' ' J ' '
9.00
MM2, M M 3 Aj,
example, the stacked configuration prevails at very short distances (3.5-4.5 A), where all three A's are close to unity, although the associated probability is rather low according to the center-center radial distribution function in Figure 2. (The probability for the MM3 potentials is somewhat higher owing to the early rise of their curves in Figure 2.) The stacked configuration is gradually replaced by the perpendicular structure (Ai,
