9624
J. Phys. Chem. 1996, 100, 9624-9630
Comparative Study of Force Fields for Benzene Grant D. Smith* Department of Chemical Engineering, UniVersity of Missouri-Columbia, Columbia, Missouri 65211
Richard L. Jaffe NASA Ames Research Center, Moffett Field, California 94035 ReceiVed: NoVember 28, 1995; In Final Form: March 27, 1996X
We have examined the ability of several atomic force fields of different forms to reproduce the structure and binding energy of benzene dimer as determined from quantum chemistry calculations, the experimental gasphase second virial coefficient of benzene, and thermodynamic and structural properties of liquid benzene. The force fields investigated were a united atom Lennard-Jones potential, a united atom Lennard-Jones potential with a point quadrupole, an anisotropic united atom model, and explicit (all) atom force fields with partial atomic charges and without partial atomic charges. These force fields which do not include electrostatic interactions predicted the planar sandwich structure (D6h) to be the most stable dimer structure, in contrast to quantum chemistry calculations that indicate that the planar sandwich structure is a saddle point between stable parallel displaced (C2h) geometries. Only the explicit atom + partial charge model predicted benzene dimer structures and energies that are in qualitative agreement with those from quantum chemistry calculations. A force field of this form for benzene is presented that accurately reproduces the gas-phase second virial coefficient of benzene and the thermodynamic properties of liquid benzene and does a reasonable job in reproducing the structure of liquid benzene and properties of benzene dimer as determined from quantum chemistry calculations.
Introduction
Quantum Chemistry Study of Benzene Dimer
Most atomic force fields currently being used in molecular dynamics simulations of organic molecules and polymers consist of two parts: (1) specific bonded or valence terms to describe the molecular geometry, vibrations, and conformational energies; (2) general nonbonded terms to describe van der Waals exchange repulsion/dispersion forces and electrostatic effects. In practice, great care is taken to ensure that the valence terms accurately reflect the molecular structure, but although it is well known that intermolecular interactions play a critical role in determining the properties of matter in condensed phases, much less attention is often paid to the parameters for these interactions. Many organic molecules and polymers of scientific and technological interest contain aromatic groups. In our opinion, the ability of existing force fields to accurately reproduce aromatic interactions has not been established. In this work we investigate nonbonded parameters for benzene. Numerous force fields have been proposed to describe the interactions between simple aromatic molecules in the liquid and crystalline states (e.g., refs 1-9). We have examined several of these with the hope of identifying a simple and accurate force field that can be used in molecular dynamics simulations of organic molecules and polymers containing aromatic groups. For this purpose we have investigated the ability of various force fields to represent properties of benzene. We adopt the following three criteria for judging the quality of force fields: (1) agreement between calculated and measured gas-phase second virial coefficients for benzene, (2) agreement between thermodynamic and structural properties of liquid benzene from molecular dynamics simulations and experiment, and (3) how well the force field reproduces our recent quantum chemistry results for the structure and binding energy of benzene dimer.10
The anisotropic nature of benzene means that those force fields which lead to calculated liquid benzene thermodynamic properties and gas-phase second virial coefficients in good agreement with experimental data can yield qualitatively different predictions of the stable structure of benzene dimer. Thus, some additional information, such as the preferred dimer geometry, is needed to validate an empirical force field. We have therefore conducted an extensive quantum chemistry study of benzene dimer10 in order to obtain accurate data for benzene dimer geometries and energies. The salient results of this study are summarized here. The interaction energy of the benzene dimer arises from the sum of the steric repulsion, dispersion, and electrostatic interactions. Dispersion is always attractive, but the electrostatic interactions can be either attractive or repulsive depending on the dimer geometry. Although a fairly accurate estimate of electrostatic energy can be obtained from ab initio quantum chemistry calculations at the Hartree-Fock level using a moderate size basis set, calculation of dispersion energies requires inclusion of electron correlation effects using a large polarized basis set. Such calculations are difficult to carry out for a system with as many atoms and electrons as the benzene dimer. In addition, calculations of dimer interactions must be corrected for the effect of basis set superposition error (BSSE). The best previously published calculations for benzene dimer are by Tsuzuki et al.11 They carried out second-order MøllerPlesset perturbation theory (MP2) calculations with large basis sets (up to 6-311G(3d,3p)) for the planar sandwich structure and found it to be bound by nearly 3 kcal/mol after correcting for BSSE. The planar sandwich and other proposed low-energy structures for benzene dimer are illustrated in Figure 1. At the SCF level, the dimer is completely unbound, reflecting repulsive electrostatic interactions for the planar sandwich structure. Even at the 6-31G** MP2 level of calculation, the dimer is only
X
Abstract published in AdVance ACS Abstracts, May 1, 1996.
S0022-3654(95)03519-2 CCC: $12.00
© 1996 American Chemical Society
Force Fields for Benzene
J. Phys. Chem., Vol. 100, No. 23, 1996 9625 On the other hand, we have found that the MP2 method overestimates the differential binding energy by as much as 1.00-1.25 kcal/mol10 compared to MP4 level calculations including triple excitations [MP4(SDTQ)]. The MP4(SDTQ) calculations should provide a much more accurate estimate of electron correlation effects for a given basis set. For a series of basis sets we found the ratio of MP4 to MP2 correlation energy contributions ∆Ecorr to the dimer binding to be fairly constant,10 with
∆Ecorr MP4(SDTQ)/ ∆Ecorr MP2 ) 0.84
(1)
Here, ∆Ecorr for a given geometry and basis set is given by
∆Ecorr MPi ) [EMPi - ESCF] (dimer) - 2[EMPi - ESCF] (monomer) (2)
Figure 1. Schematic representation of suggested low-energy structures for benzene dimer.
weakly bound (0.2 kcal/mol). It is only for 6-311G(2d,p) basis sets, or larger, that significant attractive dimer interaction occurs. Also notable is the work of Schlag and co-workers12 who carried out MP2 calculations using several basis sets and obtained a value of 2.2 for the binding of the T-shaped dimer. They found the parallel displaced structure to be lowest in energy with a binding of 2.3 kcal/mol. We have carried out MP2 calculations at fixed geometries of the benzene dimer10 using a 6-311G(2d,2p) basis set (384 basis functions) using the Gaussian9213 and Gaussian9414 programs. Optimal structures were determined for planar sandwich, parallel displaced, and T-shaped orientations. The results are summarized in Table 1. The BSSE correction, which was between 1.5 and 2.0 kcal/mol for most cases, has already been included in the tabulated data. At the most favorable intermolecular separation (3.8 Å), the interaction energy of the planar sandwich structure is -2.13 kcal/mol at the MP2 level. The planar sandwich structure is not a minimum with respect to horizontal displacement of one of the molecules. The actual minimum, the parallel displaced structure, is found at a horizontal displacement of 1.54 Å (slightly more than one C-C bond length), where the vertical separation of the rings is reduced to 3.4 Å and the interaction energy is increased to -3.33 kcal/mol. Tilting one of the ring planes at these geometries does not lower the energy further. For the T-shaped dimer, the structure with a C-H bond directed at the center of the other ring (point-face) is more stable than one with a C-C bond directed at the center of the other ring (edge-face). The pointface structure is bound by 2.84 kcal/mol. However, the T-shaped structures are not minima with respect to the pivoting of one molecule along a path leading to the lower energy parallel displaced structure. The result that the parallel displaced structure is the most stable form of benzene dimer is in agreement with the experimental results of Bernstein et al.15 A comparison of quantum chemistry and experimental data for benzene dimer can be found in ref 10. To assess the accuracy of these calculations, we have computed the dimer energy as a function of basis set and electron correlation method.10 Adding f functions and diffuse s and p functions to the carbon basis results in larger dimer binding energies (up to 3.95 kcal/mol for the parallel displaced structure at the MP2 level with a 6-311+G(2df,p) basis set).
where MPi indicates the level of electron correction treatment and EMPi and ESCF are energies with and without electron correlation, respectively. As a result, large basis set MP4(SDTQ) dimer-binding energies extrapolated from MP2 energies are 3.24, 2.64, and 2.06 kcal/mol for the parallel displaced, T-shaped, and sandwich geometries, respectively. These values are close to the 6-311G(2d,2p) MP2 results, indicating that basis set incompleteness for 6-311G(2d,2p) MP2 calculations approximately cancels the error due to the MP2 overshoot, making the results for this level of calculation close to the extrapolated large basis set MP4(SDTQ) energies. Therefore, our 6-311G(2d,2p) MP2 calculations, summarized in Table 1, appear to be appropriate for studying the benzene dimer potential energy surface. Force Fields Atomic force fields for benzene can be classified based upon the number of force centers included and the treatment of electrostatic interactions. The force fields are usually either a united atom with 6 force centers located at the carbon atom positions, or along the C-H bonds (anisotropic united atom), or an explicit (all) atom with 12 force centers located at the carbon and hydrogen atom positions. The force fields either include electrostatic interactions, through partial atomic charges or a point quadrupole, or neglect electrostatic effects. It was not our intention to exhaustively investigate every published force field for benzene. Rather, we have investigated the behavior of what we considered to be representative force fields of each type in order to determine what general characteristics of the force fields are important for meeting the criteria listed above. Specifically, we have considered two explicit atom aromatic force fields of Williams and Starr2 (WS1 and WS2), who derived parameters for C-C and H-H exponential-6 plus Coulomb pair potentials to best reproduce the crystal structures of a set of mostly aromatic hydrocarbons. Force field WS1 was derived with partial atomic charges and WS2 without charges. We also consider the force field of Evans and Watts1 (EW), who derived their force field to reproduce the interaction of two benzene quadrupoles and represented it by an anisotropic united atom Lennard-Jones potential. The EW force field does not contain any electrostatic terms (the displacement of the force center away from the carbon atoms is supposed to compensate for the neglect of electrostatic effects). Finally, two force fields due to Claessens et al.3 were considered. They adjusted parameters for a united atom Lennard-Jones potential with a point quadrupole located in the center of the benzene ring (LJ + Q) and without the point quadrupole (LJ) to bring the results of a benzene MD simulation into agreement with experimental pressure and energy of vaporization. For convenience, the force
9626 J. Phys. Chem., Vol. 100, No. 23, 1996
Smith and Jaffe
TABLE 1: Benzene Dimer Stationary Point Energies and Geometries From Quantum Chemistry and Various Force Fields planar sandwich
T (point-face)
T (edge-face)
parallel displaced
method
∆E
Rx
Rza
∆E
Rx
Rz
∆E
Rx
Rz
∆E
Rx
Rz
quantum chemistry WS1 WS2 EW LJ LJ + Q this work
-2.13
0.0
3.8
-2.84
0.0
4.9
-2.51
0.0
5.0
-3.33
1.54
3.4
-1.36 -3.71 -3.78 -3.33
0.0 0.0 0.0 0.0 b 0.0
3.8 3.6 3.7 4.0
-2.34 -1.65 -2.10 -1.83 -2.58 -2.37
0.0 0.0 0.0 0.0 0.0 0.0
5.2 5.7 5.1 5.1 5.0 5.0
-2.43 -1.87
0.0 0.0
5.1 5.0
-2.57
3.3
-1.92
0.0
5.1
-2.46
0.0
4.9
3.4 b 0.4 b 2.8 2.1
-2.23
3.7
-3.79 -2.66 -2.49
3.7 3.7 3.5
a
Benzene A lies in the xy plane. Rh is the horizontal displacement of the center of mass of benzene B with respect to that of benzene A. Rz is the vertical displacement of the center of mass of benzene B with respect to that of benzene A. b No stable dimer was found.
TABLE 2: Force Field Parameters for Benzene force field parameter
WS1
WS2
reference Accb Bcc Ccc Ahh Bhh Chh Ach Bch Cch qc Ac C Q
16 87 775 3.6 577.0 2790.9 3.74 32.51 15 651.3 3.67 136.95 -0.15
16 71 702 3.6 511.47 2170.2 3.74 24.38 8508.6 3.67 111.62 0.00
EWa 14
LJ 15
LJ + Q
this work
15 78 998 3.6 519.3 2384.6 3.74 24.624 3888.0 3.415 124.416 -0.11
2 068 000 1125.0 0.00
3 087 600 1165.0 0.00
2 908 500 1052.5 -1.81
a The EW force field is an anisotropic force field. Force centers are at positions given by an effective C-C bond length of 1.756 Å. b Exp-6 potential: Vij (kcal/mol) ) Aij exp(-BijRij) - Cij/Rij6 + 332.08qiqj/Rij. Units: A (kcal/mol), B (Å-1), C (kcal Å6/mol), q (electrons). c Lennard-Jones potential: Vij (kcal/mol) ) Aij/Rij12 - Cij/Rij6 + F(Q), where F(Q) is given by ref 3. Units: A (kcal Å12/mol), C (kcal Å6/mol), Q (electrons Å2).
of temperature in good agreement with experimental results. This agreement indicates that the force fields accurately represent the average magnitude of the interaction between two benzene molecules. However, because of the anisotropy of the benzene molecule, the interaction energy between two benzene molecules is strongly dependent upon the relative orientation of the molecules. Since the gas-phase second virial coefficient is an average over all configurations, good agreement with experimental results does not necessarily indicate that a force field accurately represents energies of particular dimer structures, as will be shown below. Molecular Dynamics Simulations
Figure 2. Gas-phase second virial coefficient of benzene as a function of temperature, predicted using various force fields. Experimental data are from ref 16.
field parameters for the various force fields considered are summarized in Table 2. Gas-Phase Second Virial Coefficient The calculated gas-phase second virial coefficients for benzene using the various force fields are compared with those from experiment16 in Figure 2. Second virial coefficients were calculated in the manner as described elsewhere.17 The gasphase second virial coefficient is a measure of the interaction energy between two benzene molecules averaged over all orientations and separations. All the force fields investigated yielded the second virial coefficient of benzene as a function
Molecular dynamics simulations of liquid benzene were performed on an ensemble of 150 molecules at 298 K. A constant temperature and volume algorithm described elsewhere17 was used, with the volume set to correspond to the experimental density.18 This results in a periodic cubic cell with a length of 28.132 Å on a side. Bond lengths were constrained using the standard SHAKE algorithm.19 For force fields containing explicit hydrogen atoms, the hydrogen positions were constrained. An integration time step of 3 fs was employed. A spherical cutoff of 12 Å was utilized, and the pressure and energy were corrected for truncation effected due to the dispersion/repulsion energy. For systems where partial atomic charges were used, the electrostatic forces were scaled smoothly to zero at 12 Å using a spline17 that begins at 10.5 Å. No truncation correction for the electrostatic truncation was included, but since the lowest order electric moment in benzene is the quadrupole moment, truncation effects are not large. The initial ensemble consisted of the benzene molecules on a low-
Force Fields for Benzene
J. Phys. Chem., Vol. 100, No. 23, 1996 9627
TABLE 3: Liquid Benzene Properties at 298 K from Molecular Dynamics Simulations Using Various Force Fields
a
force field
P (atm)
Uvap (kcal/mol)a
WS1 WS2 EW LJ LJ + Q this work experimental
1500 750 2300 -70 10 -70 1
8.65 6.45 8.05 7.46 7.46 7.58 7.50b
Energy of vaporization. b From ref 18.
density lattice (about 10% of the experimental density). The ensemble was allowed to relax at 600 K at the reduced density for 500 ps. The density was then increased to the experimental value and the temperature decreased to 298 K over 500 ps. The resulting ensemble was equilibrated for an additional 500 ps. The ensembles were equilibrated for an additional 150 ps whenever the force field was modified. Production runs were also 150 ps in length. The calculated pressures and vaporization energies for liquid benzene using the various force fields are given in Table 3. For the LJ, EW, and LJ + Q force fields, values were taken from ref 3. Our code yields results for the EW and LJ force fields that are in good agreement with the results of ref 3. Baranyai and Evans8 dispute the accuracy of the pressure calculations for the liquid benzene simulations of Claessens et al.3 using the EW force field, citing a possible neglect of truncation corrections. However, with a longer truncation distance, and including truncation corrections, our pressure and energy of vaporization are in agreement with results of ref 3. Despite the fact that all force fields yielded good agreement of the gas-phase second virial coefficient of benzene with experiment, only the LJ and LJ + Q force fields, which were parametrized to reproduce the P-V-T behavior of liquid benzene, accurately reproduce the P-V-T properties of liquid benzene. Benzene Dimer Energies and Geometries The energies and geometries for the benzene dimer stationary points from the various force fields are compared with those from quantum chemistry calculations in Table 1. The energy as a function of vertical center of mass displacement Rz for paths that bring the monomers from large separation to the planar sandwich, T-shaped point-face, and T-shaped edge-face stationary points are shown in Figures 3-5, respectively. For the sake of clarity we do not show the EW paths, which qualitatively show the same behavior as the LJ paths. Examination of Table 1 and Figures 3-5 reveals that none of these force fields does a good job in representing the quantum chemistry dimer energies and geometries. Planar Sandwich. For the planar sandwich, the WS2, EW, and LJ force fields, which do not include representation of electrostatic interactions, yield structures that are much too tightly bound compared to the quantum chemistry results. The stability of the planar sandwich structure for these force fields arises from the fact that this geometry maximizes attractive dispersion interactions. For the WS1 and LJ + Q force fields, attractive dispersion interactions are offset by repulsive electrostatic terms. For the LJ + Q force field, where the quadrupole moment was set equal to the experimental value,3 the electrostatic repulsion is such that the planar sandwich is not a stationary point. This behavior indicates that representing electrostatic interactions between benzene molecules in terms of quadrupole interactions is not reasonable when the intermolecular separation is less than the size of the molecules. It is well known that multipole expansions do not give an accurate
Figure 3. Energy of the planar sandwich benzene dimer structure as a function of the vertical center of mass separation Rz determined from quantum chemistry and various force fields.
Figure 4. Energy of the T-shaped point-face benzene dimer structure as a function of the vertical center of mass separation Rz determined from quantum chemistry and various force fields.
Figure 5. Energy of the T-shaped edge-face benzene dimer structure as a function of the vertical center of mass separation Rz determined from quantum chemistry and various force fields.
representation of electrostatic interactions at these distances. For larger separations, representing electrostatic interactions in terms of quadrupole interactions is probably valid. The WS1 force field, which associates partial atomic charges with each atom, although not accurately predicting the planar sandwich energy, does predict a saddle point that is not too strongly bound. All force fields show the staggered planar sandwich structure
9628 J. Phys. Chem., Vol. 100, No. 23, 1996 (rotation of one ring about the C6 axis by 30°) to be a few hundredths of a kcal/mol lower in energy than the eclipsed structure, in agreement with quantum chemistry results. These comparisons of planar sandwich energies with quantum chemistry calculations indicate that electrostatic interactions must be considered in an accurate benzene force field, consistent with our quantum chemistry calculations that show electrostatic interactions to be important for the benzene dimer.10 A point quadrupole of the magnitude of the experimental quadrupole moment of benzene does not suffice, however, in representation of these interactions for small separations. T-Shaped Structures. The WS2, EW, and LJ force fields yield T-shaped structures that are not strongly enough bound compared to those from quantum chemistry calculations, in contrast to the planar sandwich, where the predicted structures for these force fields were too tightly bound. The WS1 and LJ + Q force fields at least qualitatively reproduce the T-shaped structure energies and paths shown in Figures 4 and 5. The low values for the binding energy of the T-shaped structures by force fields not containing electrostatic interactions are consistent with the conclusion from our quantum chemistry study that electrostatic interactions are much more favorable for these structures than for the parallel structures.10 Parallel Displaced. The WS2 and LJ force fields do not yield a stable parallel displaced structure. The EW force field yields a stable parallel displaced structure, but the geometry of this structure is very different from that of the quantum chemistry structure (see Table 1). In the EW structure, the molecules are only slightly displaced from the planar sandwich structure and the energy is essentially the same as that of the planar sandwich. The WS1 and LJ + Q force fields predict the existence of a low-energy parallel displaced structure, although in both cases the energy difference between the T-shaped saddle points and the parallel displaced minimum is much less than that seen in our quantum chemistry calculations, i.e., the path for interconversion of parallel displaced geometries is much flatter for the WS1 and LJ + Q force fields than that seen in the quantum chemistry calculations. Since only the WS1 and LJ + Q force fields qualitatively describe the parallel displaced structure, it appears that even the existence of a stable parallel displaced structure necessitates the inclusion electrostatic interactions in the force field. Alternative Force Field for Benzene Since none of the force fields we investigated for benzene are acceptable according to our criteria, we have made an attempt to devise an improved force field for benzene. From the discussion above, it is apparent that of all the forms of the force field considered (all atom, all atom + partial charge, united atom, anisotropic united atom, and united atom + point quadrupole), the all atom + partial charge model shows the most promise. We started with the WS1 C-C exp-6 parameters and charges and the C-H and H-H exp-6 parameters derived by Sorensen et al. for hydrocarbons20 and used previously by us. In contrast to most of the other force fields available for benzene, this force field does not treat C-H parameters as mean values of the C-C and H-H parameters. We then systematically reduced the partial charge from the WS1 value and uniformly scaled all the pair energies in order to improve agreement of predicted benzene liquid and dimer properties with those from experiment and quantum chemistry. The final parameters are given in Table 2. Our force field yields the gas-phase second virial coefficient of benzene and liquid benzene P-V-T properties that are in good agreement with experimental results, as shown in Figure
Smith and Jaffe
Figure 6. C-C intermolecular pair distribution function for benzene: present calculation (symbols); experimental results (ref 21) (smooth line).
2 and Table 3. From Table 1 and Figures 3-5, it is apparent that our model does a better job representing the dimer structures and energies determined from quantum chemistry calculations than the other force fields investigated. The energy of the planar sandwich structure and the corresponding path (Figure 3) are accurately reproduced. The force field does a fair job reproducing the T-shaped structures and paths (Figures 3 and 4). For the parallel displaced structure, the force field predicts a geometry that is in reasonable agreement with quantum chemistry results. The binding energy for the parallel displaced structure, however, is less than that predicted from quantum chemistry and is only slightly greater than that for the T-shaped dimer structures. As with the WS1 and LJ + Q force fields, the path for interconversion of parallel displaced structure through the T-shaped saddle points is much flatter than predicted by quantum chemistry. The energy surface predicted by the force field is very flat with respect to displacements and tilts from the parallel displaced geometry, which can lead to a small decrease in energy from the parallel displaced structure. The minimum energy configuration for the force field is -2.51 kcal/ mol with Rz ) 1.8 Å, Rz ) 4.3 Å, and a tilt angle of 40°. Our quantum chemistry calculations indicate that the energy surface is flat in this region of configuration space, but not as flat as predicted by the force field. A second, shallow minimum can be seen in the quantum chemistry path with a geometry quite similar to the minimum predicted by the force field.10 The C-C intermolecular pair distribution function for liquid benzene from MD simulations for our force field is compared with experimental results21 in Figure 6. The structure of liquid benzene predicted by our force field is in better agreement with experimental results than that predicted by the other force fields we investigated (see ref 3 for the LJ, LJ + Q, and EW force fields). Our model somewhat underestimates the number of close C-C pairs (the leading peak) probably because the population of the parallel displaced structure, which has close C-C approaches, is underrepresented. It is interesting to note that the quadrupole moment predicted by our partial atomic charges is in good agreement with experimental results. The out-of-plane component of the second moment of the charge distribution is zero for our model, in contrast to quantum chemistry predictions, which indicate that Szz, where S is the second moment tensor, is about 50% larger than Sxx or Syy. However, our charge distribution gives the wrong sign for the Sxx and Syy components, resulting in a quadrupole moment
Q ) 2/3(Szz - (Sxx + Syy)/2)
(3)
Force Fields for Benzene
J. Phys. Chem., Vol. 100, No. 23, 1996 9629 Conclusions
Figure 7. C-H nonbonded potential from several benzene force fields.
nearly equal to that from experiment. This implies that electrostatic interactions between benzene molecules will be accurately represented by our model at intermediate and greater distances. At short distances, our force field predicts that electrostatic interactions are more favorable for the T-shaped structures (around -0.5 kcal/mol) than the parallel structures (around 1.0 kcal/mol), in qualitative agreement with quantum chemistry calculations. In addition to the representative force fields discussed above, we have considered the force fields of Karlsto¨m et al.4 (KLWJ) and Jorgensen and Severance (JS),9 which have characteristics in common with our force field. Like our force field, the KLWJ force field was parametrized to reproduce benzene dimer quantum chemistry data. The level of theory employed in the dimer study was lower than that used in our study and dispersion effects appear to be significantly underestimated, resulting in quite different dimer binding energies than were found in our study.10 The gas-phase second virial coefficient for benzene as calculated by Karlsto¨m et al.4 using the KLWJ force field was not in nearly as good agreement with experiment as the force fields investigated in this work. The structure of liquid benzene predicted from Monte Carlo simulations22 using this force field also shows poorer agreement with experiment than that yielded by our force field, with the short-range structure (r < 4 Å) being slightly better by the KLWJ model but with the intermediate range order (4 Å < r < 6 Å) being better represented by our model. The JS force field is an explicit atom force field with partial atomic charges and was parametrized to reproduce thermodynamic properties of liquid benzene. The reported energies for dimer structures using the JS force field9 are in poorer agreement with our quantum chemistry values than those yielded by our force field. Like the other force fields investigated in this work, the JS force field does a poorer job than our force field in reproducing the C-C pair distribution function in liquid benzene. The ability of our force field to represent the structure of benzene liquid better than the other force fields considered may be due in part to the manner in which C-H interactions are treated; our C-H nonbonded parameters are not mean values of the C-C and H-H parameters. The dispersion/repulsion energy for the C-H interactions for the WS1, WS2, JS, and our force field is shown in Figure 7. Our C-H nonbonded potential is softer than the other potentials considered. Finally, in addition to static properties, we have examined the self diffusion of benzene predicted by our force field. We obtain a self diffusion coefficient of 1.9 × 10-5 cm2/s at 298 K, which compares well with experimental values of (2.1-2.3) × 10-5 cm2/s.23
We have demonstrated that atomic force fields for benzene of the united atom form, the anisotropic united atom, united atom + point quadrupole form, and the all atom without partial atomic charges form fail to qualitatively reproduce the structure and energies of benzene dimers predicted by high-level quantum chemistry calculations. These comparisons revealed the importance of accurately representing electrostatic interactions in modeling benzene dimer. We have presented an all atom benzene force field with partial atomic charges that accurately reproduces the gas-phase second virial coefficient of benzene and P-V-T properties of liquid benzene and does a respectable job in reproducing the C-C intermolecular pair distribution function in liquid benzene. Our force field does a better job than the other force fields investigated in reproducing energetics and geometries of benzene dimer structures predicted by quantum chemistry but certainly appears to leave room for improvement. Along these lines, we attempted to improve agreement of the shape of the benzene dimer potential energy surface with quantum chemistry predictions by adjusting parameters in the force field using a standard nonlinear least-squares method and by including additional terms with different dependencies on the interatomic separation (e.g., r-4, r-8, r-10, etc.). In several cases we were able to improve the quantitative representative of the low-energy benzene dimer energies at the expense of the liquid benzene P-V-T and gasphase second virial coefficient data. The stronger binding required to better reproduce the parallel displaced energy from quantum chemistry calculations resulted in overall interactions between benzene molecules that were too strong. A possible explanation for this behavior is that the quantum chemistry calculations yield binding energies that are somewhat too large (for example, the binding energies are not corrected for zeropoint vibrational effects) or that the force field energy surface is too flat with respect to distortions about the stable conformations. For all modifications investigated, the force field remained very flat along the path between the parallel displaced and T-shaped saddle points. It should also be noted that differences between the quantum chemistry and force field benzene dimer energies are within estimated uncertainties in the quantum chemistry energies resulting from limitations in the basis sets and electron correlation treatment used in the calculations.10 Therefore, further modifications to the force field that improve agreement with the quantum chemistry results but result in poorer agreement with liquid benzene properties do not appear to be justified. Indeed, future quantum chemistry calculations on benzene dimer at a higher level of theory could lead to quantitatively different dimer energies and geometries (although we do not expect such calculations to lead to qualitatively different conclusions regarding the stable structures of benzene dimer10). References and Notes (1) Evans, D. J.; Watts, R. O. Mol. Phys. 1975, 29, 777. (2) Williams, D. E.; Starr, T. L. Comput. Chem. 1977, 1, 173. (3) Claessens, M.; Ferrario, M.; Ryckaert, J. P. Mol. Phys. 1983, 50, 217. (4) Karlstro¨m, G.; Linse, P.; Wallqvist, A.; Jo¨nsson, B. J. Am. Chem. Soc. 1983, 105, 3777. (5) Anderson, J.; Ullo, J. J.; Yip, S. J. Chem. Phys. 1987, 86, 4078. (6) Yashonath, S.; Price, S. L.; McDonald, I. R. Mol. Phys. 1988, 64, 361. (7) Hoheisel, C. J. Chem. Phys. 1988, 89, 7457. (8) Baranyai, A.; Evans, D. J. Mol. Phys. 1990, 70, 53. (9) Jorgensen, W. L.; Severance, D. L. J. Am. Chem. Soc. 1990, 112, 4768. (10) Jaffe, R. L.; Smith, G. D. J. Chem. Phys., in press.
9630 J. Phys. Chem., Vol. 100, No. 23, 1996 (11) Tsuzuki, S.; Uchimaru, T.; Tanabe, K. J. Mol. Struct.: THEOCHEM. 1994, 307, 107. (12) Hobza, P. Selzle, H. L.; Schlag, E. W. J. Chem. Phys. 1990, 93 , 5893; J. Am. Chem. Soc. 1994, 116, 3500. (13) Frisch, M. J.; Trucks, G. W.; Head-Gordon, M.; Gill, P. M. W.; Wong, M. W.; Foresman, J. B.; Johnson, B. G.; Schlegel, H. B.; Robb, M. A.; Replogle, E. S.; Gomperts, R.; Andres, J. L.; Raghavachari, K.; Binkley, J. S.; Gonzalez, C.; Martin, R. L.; Fox, D. J.; Defrees, D. J.; Baker, J.; Steward, J. J. P.; Pople, J. A. GAUSSIAN 92, Revision E.2.; Gaussian, Inc: Pittsburgh, PA, 1992. (14) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Gill, P. M. W.; Johnson, B. G.; Robb, M. A.; Cheeseman, J. R.; Keith, T.; Petersson, G. A.; Montgomery, J. A.; Raghavachari, K.; Al-Laham, M. A.; Zakrzewski, V. G.; Ortiz, J. V.; Foresman, J. B.; Cioslowski, J.; Stefanov, B. B.; Nanayakkara, A.; Challacombe, M.; Peng, C. Y.; Ayala, P. Y.; Chen, W.; Wong, M. W.; Andres, J. L.; Replogle, E. S.; Gomperts, R.; Martin, R. L.; Fox, D. J.; Binkley, J. S.; Defrees, D. J.; Baker, J.; Stewart, J. P.; HeadGordon, M.; Gonzalez, C.; Pople, J. A. Gaussian 94, Revision B.1; Gaussian, Inc.: Pittsburgh, PA, 1995.
Smith and Jaffe (15) Law, K. S.; Schauer, M.; Bernstein, E. R. J. Chem. Phys. 1984, 81, 4871 (16) Dymond, J. H.; Smith, E. B. The Virial Coefficients of Pure Gases and Mixtures; Clarendon: Oxford, 1980. (17) Smith, G. D.; Jaffe, R. L.; Yoon, D. Y. Macromolecules 1993, 26, 298. (18) Selected Values of Physical and Thermodynamic Properties of Hydrocarbons and Related Compounds; American Petroleum Institute Research Project 44; Carnegie Press: Pittsburgh, PA, 1953. (19) Ryckaert, J. P.; Ciccotti, G.; Berendsen, H. J. C. J. Comput. Phys. 1977, 23, 327. (20) Sorensen, R. A.; Liau, W. B.; Kesner, L.; Boyd, R. H. Macromolecules 1988, 21, 200. (21) Narten, A. H. J. Chem. Phys. 1977, 67, 2102. (22) Linse, P. J. Am. Chem. Soc. 1984, 106, 5425. (23) Falcone, D. R.; Douglass, D. C.; McCall, D. W. J. Phys. Chem. 1967, 71, 2744.
JP9535194
