Electron diffraction and Monte Carlo studies of liquids. 1

A Monte Carlo Study of Isomers and Structural Evolution in ... Monte Carlo Studies of Isomers, Structures, and Properties in .... Efforts by companies...
0 downloads 0 Views 1015KB Size
Download PDF
J. Phys. Chem. 1988, 92, 5661-5613

5667

Electron Dlffractlon and Monte Carlo Studies of Liquids. 1. Intermolecular Interactions for Benzene Xianqian Shi and Lawrence S. Bartell* Department of Chemistry, University of Michigan, Ann Arbor, Michigan 48109 (Received: January 25, 1988)

Various models of intermolecularinteractions for benzene were investigated for their suitability in simulationsof condensed-phase properties. It is proposed that, for a treatment to be considered realistic, it should satisfy the following four criteria. The potential energy function incorporated should accurately reproduce (1) the structure and (2) heat of sublimation of the crystalline phase and (3) the structure and (4) thermodynamic properties of the liquid. The present paper screens a number of model potentials for their ability to meet the first two of the four criteria. Results corroborate and amplify prior analyses of the role that charges distributed over the benzene molecule play in influencing the way the molecules pack together. A model incorporating a point-quadrupole or six point-dipolesis inadequate. Models with only six interaction sites are poor, and 12-site models without charges are mediocre. While criteria 1 and 2 impose severe constraints, they do not lead to a unique result for they can be satisfied by a variety of quite different model potentials. Electrostatic effects can be modeled well by truncated r-2interactions parametrized to reproduce electrostatic forces at atom-atom contact distances, a fact showing that the required effects are short-range discriminations. Two particularly promising models with very different properties, those of Williams and Starr and of Karlstrom et al., were rescaled and simplified to make them better suited for computer simulations.

Introduction Molecular clusters produced by homogeneous nucleation in supersonic flow can be generated in a variety of forms whose structures can be probed by electron diffraction.’ A particularly promising possibility is to nucleate and grow highly supercooled liquid clusters that are large enough to have bulk properties, yet small enough to transmit electrons with negligible multiple scattering. What is attractive about this approach is that it offers a particularly demanding test of the theory of molecular liquids. Diffraction patterns of cold, dense liquids contain more detailed information about structure and intermolecular interactions than do diffraction patterns of relatively warm liquids with freer molecular motiom2 Moreover, liquids, in principle, may reveal even more about intramolecular interactions than do crystals, for molecules in mutual contact in crystals are constrained to very special orientations with respect to each other. Experiments to date confirm the feasibility of such a program of research. Clusters of many of the substances studied have yielded diffraction patterns characteristic of liquid or liquidlike phase^,^ and cluster size can be controlled to be in the vicinity of 100 A, a favorable range for electron diffraction. The diffraction patterns, by themselves, however, cannot provide all of the information sought. For one thing, the degree of supercooling carriers a penalty as well as a reward. No direct measures of the cluster temperatures are available. For another, the required partitioning of the diffracted intensity into atomic, intramolecular, and intermolecular components cannot be made with precision without auxiliary information. In the most rigorous studies of gas-phase molecular structure this information is generated in a simultaneous refinement of the atomic/intramolecular partitioning on the one hand, and of the molecular force field (implying atomic positions and motions) on the other. An analogous cycle-by-cycle refinement of liquid diffraction patterns is out of the question at this time, of course,although an effort in this direction is desirable. What has become clear is that a full analysis of the information in electron diffraction patterns of liquid clusters cannot be made without detailed comparisons between experiment and calculations based on the statistical thermodynamics of liquids. It is equally clear that available theoretical computations for systems as complex as that of benzene, for example, are almost worthless for the purpbse. What is required is the development of realistic potential functions for implementation into whatever theoretical (1) See, for example: Bartell, L. S. Chem. Rev. 1986,86, 491. (2) Valente, E.J.; Bartell, L. S. J . Chem. Phys. 1984,80, 1451. (3) Bartell, L.S.;Harsanyi, L.; Valente, E . J. In Physics and Chemistry of Small Clusters; Jena, P., Rao, 8. K., Khanna, S. N., Eds.; Plenum: New York, 1987;pp 31-36.

0022-3654/88/2092-5667$01.50/0

framework is adopted (e&, RISM, Monte Carlo, or molecular dynamics). Moreover, the adequacy of the various frameworks has yet to be evaluated for moderately complex systems. For the present purposes it would be considered sufficient for a given interaction potential energy to account for crystal structure, heat of sublimation, liquid structure, and thermodynamic properties of liquid. The present series of papers begins with condensed benzene as a model system that is simple enough to be feasible, yet considerably more complex than previous systems subjected to extensive statistical analysis. Many model interaction energy functions for benzene have been constructed in prior work, yet each is seriously flawed, as will be discussed in this and subsequent papers. Moreover, a number of somewhat subtle points not analyzed explicitly in the previous literature should be taken into account in the construction of a realistic model potential. These will be aired in the following. In view of our incomplete knowledge of interactions between molecules, the goal of achieving a truly precise interaction law has not been attained in the present work. Nevertheless, appreciable progress has been made and the present paper describes the course of development of a function to characterize intermolecular interactions in benzene. Even if such a potential surface were known, accurately, however, there would remain problems in carrying out definitive comparisons of theoretical simulations and diffraction experiments. These are analyzed in paper 114of this series. Results of comparisons between theoretical computations and experiment are presented in paper 111.5

Prior Investigations No attempt will be made to review exhaustively the large body of prior research on interactions between benzene molecules in condensed phases. Many aspects of the problem together with work prior to 1975 have been reviewed elsewhere.6 The potential functions of concern in the present investigation are those which were introduced to account for (1) crystal structure and, in some cases, lattice modes, utilizing packing analyses via potential energy minimization, and/or (2) thermodynamic properties of the liquid. Models fall into two categories, namely those with six interaction sites per molecule (6-site) and those with twelve (12-site). In some cases electrostatic interactions associated with the molecule’s known quadrupole moment have been included as well. As first shown by Williams and co-~orkers,’-~and reinforced by Dzy(4)Bartell, L.S. Acta Chem. Scand., in press. ( 5 ) Bartell, L. S.; Sharkey, L. R.; Shi,X . J. Am. Chem. Soc., in press. (6)MacRury, T.B.; Stbele, W. A.; Berne, B. J. J. Chem. Phys. 1976, 64, 1288. 0 1988 American Chemical Society

5668

The Journal of Physical Chemistry, Vol. 92, No. 20, 19'88

Shi and Bartell

abchenko,I0-" such interactions are crucial. 6-Site Models. Of the 6-site models in the literature that have been used in Monte Carlo (MC) or molecular dynamics (MD) simulations of liquid benzene, we will consider those due to Evans and Watts,I2 Jorgensen et al.,I3 and Claessens et al.I4 The latter authors introduced two intermolecular potential functions, one with six Lennard-Jones interaction sites situated on the carbon atoms in benzene (designated the LJ6 potential) and another (the LJ6+Q potential) which in addition placed a point quadrupole at the molecular center. Jorgensen's model,13 with which we carry out Monte Carlo calculations in paper I11 of this series5 for purposes of comparison, is very close to the LJ6 model of Claessens et al.I4 Evans and Watts12 moved the interaction sites away from carbon atoms toward the hydrogens to get a more realistic compromise for the shape of benzene. Steinha~ser'~ adopted the model of Evans and Watts in computations of dynamic properties of liquid benzene. Although the six-site potentials give a good account of the equilibrium thermodynamic properties of liquid benzene, they are deficient in several respects to be discussed presently. They were introduced as a serviceable approximation because thermodynamic properties are not sensitive to details on the molecular scale and because an increase in the complexity of the model potential enormously increases the computation time. Published results for benzene allowed no useful comparison to be made with the newly available electron diffraction patterns of liquid microdrops216 in which the molecular scatterers have 12 scattering sites. This distinction between interaction sites and scattering sites does not in itself, however, exclude 6-site models from comparisons with electron scattering experiments because it is straightforward, given coordinates of the interaction sites in a series of instantaneous configurations encountered during an M C or M D run, to graft scattering sites onto the molecular frames to permit the computation of diffraction patternse5 Nevertheless, as will become evident in this and later papers in this series, a satisfactory modeling of the molecular packing in solid and liquid benzene requires a better representation of molecules than can be achieved with six interaction sites. Therefore, we shall concentrate mainly on twelve-site models. 12-Site Models. A review of a number of 12-site models for benzene that were proposed before 1983 has been written by Dzyabchenko" who assessed their ability to represent properly the crystal structure at low pressure. Of all of the potentials tested, only that of Williams and Starr9 (WS77) reproduced experiment. Interaction sites of WS77 were situated on the carbon nuclei and at positions slightly inside the protons toward the carbons (to correct for the asphericity of hydrogen's electron distribution), and partial charges were placed at each interaction site to reproduce the quadrupole moment. For each of the other potential functions, crystal-packing analyses (CPA) yielded a global minimum a t zero pressure corresponding to the space group P 2 , / c observed in the laboratory at 25 kbar pressure, or to yet another space group, so far not seen for benzene. If constraints were then imposed to force upon the CPA the Pb, orthorhombic symmetry found for benzene at 1 atm, several minima were discovered, none of which accurately reflected the observed structure in molecular orientations and lattice constants. In contrast, the potential of Williams and Starr yielded the appropriate structures both at 1 atm and 25 kbar? Moreover, the parameters of the Williams and Starr potential satisfied the commonly invoked combining rules (C-H from C-C and Ha-H) and were successfully transferable

to other compounds than benzene, as well. Such attractive features warranted special attention in simulations of liquid properties. Recently, an altogether different 12-site model was constructed by Karlstrom et a1.I' and subjected by Lime1*to MC simulations of both liquid and crystalline benzene. It was derived from a b initio quantum calculations. Hartree-Fock SCF-MO calculations established the overlap repulsion contribution and a perturbation procedure gave the dispersion energy. Numerical results were then represented analytically by a 12-interaction-site model with each interaction of the form Ar-' Br4 Cr4 D f 9 Lime's representation will hereafter be referred to as the (12-9-6-4-1) potential. Even though the potential energy function is modeled by a 12-site representation, the form is strikingly different from that of Williams and stat^.^ The "combining rules" are far from being satisfied and the site-site interactions are not accurately transferable to other molecules. Furthermore, the potential wells of deepest interaction (if the electrostatic terms are neglected) are for the H-H, not the C-C pairs. While none of these heterodox aspects make the interaction potential violate any fundamental principles, they are counterintuitive and reduce the attractiveness to those who seek simple, transferable models for application in crystallography. They are of special interest in the present study, however, because of certain promising properties. Among others, the (12-9-6-4-1) potential yields excellent molecular orientations and lattice constant ratios for crystalline benzene. It also predicts that the stable configuration of the benzene dimer is T-shaped, the shape proposed for the dimer somewhat speculatively by Janda et al.19 on the basis of its observed dipole moment. The WS77 benzene dimer is also (mildly) polar but not T-shaped.zo

(7) Williams, D. E. Acta Crysrallogr. 1974, A30, 71. (8) Hall, D.; Williams, D. E. Acta Crystallogr. 1975, A31, 56. (9) Williams, D. E.; Starr, T. L. Compur. Chem. 1977, I , 173. (10) Dzyabchenko, A. V. Zh. Strukt. Khim. 1984, 25(3), 85. (1 1 ) Dzyabchenko, A. V. Z h . Strukt. Khim. 1984, 25(4), 57. (12) Evans, D. J.; Watts, R. 0. Mol. Phys. 1976, 32, 93. (13) Jorgensen, W. L.; Modura, J. D.; Swenson, C. J. J . Am. Chem. SOC. 1984, 106, 6638. (14) Claessens, M.; Ferrario, M.; Ryckeart, J.-P. Mol. Phys. 1983,50,217. (15) Steinhauser, 0. Chem. Phys. 1982, 73, 155. (16) Heenan, R. K.; Valente, E. J.; Bartell, L. S . J . Chem. Phys. 1983, 78, 243.

(17) Karlstrom, G.; Linse, P.;Wallqvist, A.; JBnsson, B. J . Am. Chem. Soc. 1983, 105, 3777. (18) Linse, P.J . Am. Chem. Soc. 1984, 106, 5425. (19) Janda, K. C.; Hemminger, J. C.; Winn, J. S.; Novick, S . E.; Harris, S. J.; Klemperer, W. J . Chem. Phys. 1975, 63, 1419. (20) Williams, D. E. Acta Crystallogr. 1980, A36, 715. (21) Evans, D. J.; Watts, R. 0. Mol. Phys. 1975, 29, 777; 1976, 31, 83. (22) Kuharski, R. A.; Rossky, P. J. J . Chem. Phys. 1985, 82, 5164. (23) Kuharski, R. A,; Rossky, P. J. J . Chem. Phys. 1985,82, 5289. (24) Bartell, L. S.; Roskos, R. R. J . Chem. Phys. 1966, 44, 457. (25) Barker, J. A.; Fisher, R. A.; Watts, R. 0. Mol. Phys. 1971, 21, 657. (26) Barker, J. A.; Klein, M. L. Phys. Rev. 1973, B7, 4707.

+

+

+

+

Preliminary Considerations in Modeling Form of Potential Functions. Although some work has been done on many-body interactions in condensed benzene2' and on quantum correctionsz2and effects of intramolecular vibration^^^*^^ in simulations of condensed phases of polyatomic molecules, the inclusion of such effects would increase the cost and complexity of computations to a degree unwarranted in the present investigation. It is believed that many-body, quantum, and vibrational effects at ordinary temperatures are fairly small and that, up to a point, they can be absorbed in the model two-body potentia1.25*26 Enough poorly understood aspects remain to be explored that it is premature to undertake a final fine-tuning that would include many-body effects and quantum effects beyond those discussed below. Interactions between benzene molecules, then, will be taken as pairwise additive and appropriate for computations involving condensed phases rather than for calculation of the second virial coefficient of the vapor. Considerations of economy also encourage a severe simplification in the form of individual sitesite potentials, especially when as many as 12 interaction sites per molecule are invoked. Therefore, (exp6) van der Waals interactions of the sort adopted by Williams in his successful treatment of organic crystals7+' will be simplified for M C computations to n - 6 interactions where n signifies an r-n term. Furthermore, square roots to get Coulomb terms and f 9terms (as in Karlstrom's potential)I7 are replaced by terms with even powers of l / r as discussed in a later section. The reason this is important in M C runs even with today's computers is that the principal consumption of time is in the computation of potential energy terms. In a run of lo6 cycles with 128 12-site molecules there are over loloatom-atom contributions. On a Cray X-MP/48 computer it takes over 5 times as long,

Intermolecular Interactions for Benzene

1'.. .

.............................

1 1 0 -

0

100

200

T, K Figure 1. Volume per molecule as a function of temperature. Solid line, behavior of real crystal; dashed line, behavior in classical simulation corresponding to the same potential energy function. Values approximately those for benzene.

starting with 9,to compute a term of form (exp6) as it does with form (12-6-2). Confusion of Parameters by Effects of Zero-Point Motions. One of the criteria to be met by a realistic treatment is that the model potential correctly reproduce the crystal structure and heat of sublimation. All prior M C and M D studies of benzene have used as reference points the estimated volume per molecule, Vo, and heat of sublimation, -E,", of the material at 0 K. We regard these reference points as useful but troublesome, for the following reasons. The utility is partly due to the fact that Vo and the packing energy E," can be calculated routinely in crystal-packing analysis (CPA). Such analyses via, for example, the Williams PCK a l g ~ r i t h m ~are ~ *efficient ~* and effective ways to test a potential. But neither the CPA nor the limit at 0 K of a classical M C computation yields a directly observable property. The quantity E," = -52.3 kJ/mol for benzene29 was established by extrapolating the observed heat of sublimation to absolute zero and then subtracting from it the zero-point energy of the lattice computed from the phonon spectrum.30 This, presumably, would yield the potential energy of packing of molecules in a lattice corresponding to molecular volume Von, the observed (quantum) volume at 0 K. The trouble is that Von is inflated beyond Voby zero-point motions, and that, if a lattice could be freed of zeropoint motions, as it relaxed to its static state the packing energy would relax to EcW,dropping by an amount very nearly equal to

E," - ECm = (Vo" - Vo)(aE,"/av),,

(1) Unfortunately, the correction ( Von - Vo)is not known reliably , ~ be estimated reaeven though the derivative ( 1 3 E , " / d y ) ~can sonably well from the model potential. A helpful way to look at the problem and its partial resolution is in Figure 1. According to the third law of thermodynamics, the slope of the experimental curve V(T)must fall to zero at 0 K. In a classical M C or M D simulation, on the other hand, the slope does not, and the fictional intercept Vodiffers from the larger Vozp by just the inflationary effect of quantum zero-point motions. One potentially useful way to estimate the desired correction (Vozp - Vo) is through Griineissen's r e l a t i ~ n ~ ~ , ~ ~ Y = .V/tsc,, (2) where CY and 0 are the coefficients of thermal expansion and compressibility, respectively, and y is Griineissen's parameter which is very nearly independent of temperature according to experiment and to theory, with certain assumption^.^'-^^ What ~~~

~

The Journal of Physical Chemistry, Vol. 92, No. 20, 1988 5669 is needed to establish Vois a "classical" coefficient of thermal expansion ac,(7') to guide the "classical extrapolation" of V(T ) to 0 K. Because the derivation of eq 2 does not involve quantum mechanics, a way to find the needed slope Vat,( 7') is to use the classical law of Dulong and Petit (or its counterpart for molecular crystals) for C,. This works well and correctly reproduces V(T ) calculated for crystalline argon in M C runs.33 The problem for benzene is that no satisfactory low-temperature values of a,/3, or y had been established for benzene in the initial stages of the present research. Therefore, an extrapolation of the observed lattice constants of benzene to 0 K was made subjectively, using a gently curved slope of classical aspect. Because of considerable scatter in the data at the lower temperatures,3636 the result was little better than a guess. Because the results for Vo (and the lattice constants and molecular orientations) were close to those already published by Evans and Watts, who seem to have used a linear extrapolation of the more reliable lattice constants, and since these results have been used by all subsequent investigators, we adopted the Evans and Watts values as references.12 Previous authors, however, have not to our knowledge explicitly discussed the above sources of trouble or allowed for the differences between E," and EcWor the error in Vo. A rough estimate of (E," - ECW)/E,"was 0.02, a value that is small but not insignificant. This estimate, however, was based on the Evans and Watts Vo and a guess of Von and was considered to be too crude to apply to the M C calculations. A reassessment of the problem during the final stages of the research led to a promising treatment of eq 2 in its quantum and classical aspects. A plausible function @(T ) / @ ( O ) inferred from g(T)/[V(T) - VO]of Ar,25 using T/T(mp) to scale results to benzene, was adopted together with C,for the lattice modes from from the X-ray and ref 30 and the approximate values of CUV neutron measurements.3636 This established a reasonable value for y@(O). By fixing yj3(0), then, a function a(T ) was constructed to generate a curve V(T ) that nicely cleaved the scattered experimental values and extrapolated to Vo* = 114.6 A3/molecule. Redrawing a(T) to give it a smooth classical extrapolation to 0 K (not allowing it to follow the steep quantum drop to zero) seemed to establish Vo = 112 A3/molecule to within a fraction of a percent (cf. 110.8 A3/molecule according to Evans and Watts). These values reduced (E," - EcW)/E," to about 0.01. Since this value is small and speculative, results in the following will adopt E," instead of EcW as a reference. Construction of Potential Functions 6-Site Models. Despite the known deficiency of existing 6-site potentials in reproducing the crystal structure, it seemed worthwhile to reexamine the problem. Previous 6-site models were based on Lennard-Jones 12-6 interactions, and site positions of most were not optimally selected. Some evidence in the literature suggests that an m - 6 potential with m less than 12 is superior to 12-6 in interatomic interaction^.^'.^^ Therefore, 8-6 and 10-6 functions were initially considered as well as 12-6, computational considerations indicating that m be an even integer. A series of M C runs at different temperatures with 8-6 functions gave thermodynamic properties comparable to those of 12-6 runs. Finding the optimum site position for PCA results could be reduced to a one-dimensional problem of optimizing the ratio u/rs, where r, is the site radius, and then rescaling to obtain the correct cell volume. Presumably the most favorable site radius which governs the shape of the molecule (flat for small u/rs,spherical for large u/rs)would correspond to the best representation of the

~

(27) Williams, D.E. Quantum Chemistry Program Exchange; Indiana University; Bloomington, IN 47405, 1979; Program 373. (28) Williams, D. E. Quantum Chemistry Program Exchange; Indiana University; Bloomington, IN 47405, 1983; Program 481. (29) Oliver, G.D.;Eaton, M.; Huffman, H. M. J . Am. Chem. SOC.1948, 70, 1502. (30) Nakamura, M.;Mujazawa, T. J . Chem. Phys. 1969,51, 3146. (31) Griineissen, E. Handb. Phys. 1926, 10, 43. (32) Mott, N. F.;Jones, H. The Theory of the Properties of Metals and Alloys; Dover: New York, 1958; pp 19-21.

(33) Bartell, L.S.;Sharkey, L. R., unpublished research. (34) Kohzin, B. M.Zh. Fir. Khim. 1954, 28, 566. Kohzin, B. M.;Kitaigorcdskii, A. I. Zh. Fiz. Khim. 1955, 29, 2074. (35) Cox, E. G.; Cruickshank, D. W. J.; Smith, J. A. C. Proc. R. Soc. London A 1958, A247, 1 . (36) Bacon, G.E.; Curry, N. A.; Wilson, S. A. Proc. R. SOC.London A 1964, A279, 98. (37) Pauling, L . The Nature of the Chemical Bond, 2nd ed.; Cornel1 University Press: Ithaca, NY, 1945; p 339. (38) Warshel, A.; Lifson, S.J. Chem. Phys. 1970, 53, 582.

5670 The Journal of Physical Chemistry, Vol. 92, No. 20, 1988

Shi and Bartell TABLE I: Parameters in Potential Functions’ EAJ,,P

type, ij

AI2

Ai0

AI

A2

A6

AO

(1 2-6)

cc CH HH

3 439 400 451 930 29 100

-2804.2 -552.5 -9 1.6

2 899 300 369 090 21 010

-2291 -433.9 -66.2

2 293 900 384 200 24 740

-2383.8 -469.7 -77.9

35.10 -35.10 35.10

-1.40 1.40 -1.40

6023210 -2937.5 612474 -69047.2 23490.1 -33 862

-1971.5 187.6 -984.4

32.00 -32.00 32.00

-1.267 1.267 -1.267

(1 2-6-1)

cc

-R

CH HH

----

Figure 2. Configurations of benzene dimer adopted in sampling of tential energy. See text for details.

30.13 -30.13 30.13

( 12-6-2) PO-

cell constant ratios b / a and c/a, and molecular orientations. It turned out, however, that for an LJ6 model with no electrostatic components, molecular orientations and the c / a ratio were insensitive to u/rs. Moreover, the best b / a ratio corresponded to ~ of Jorgensen et a value for u/r, close to that of C l a e ~ s e n s ’and a1.I3 This places the interaction site very close to the carbon atoms, a somewhat counterintuitive result in view of the expected flattening effect of the ring due to hydrogens protruding well beyond the carbons. Optimization of 6-site models with no electrostatic component yielded such small benefits in PCA and M C runs, however, that optimizations for m = 8 and 10 were not pursued and attention was thereafter focused upon models incorporating some form of quadrupole or upon 12-site models. Furthermore, although little is known about details of benzenebenzene interactions, rare gas paradigm^^^.^^ indicated that however mediocre m = 12 is in modeling, a smaller value of m is worse (a smaller value gives a much better representation at very short distances but a much worse representation at distances that contribute heavily to binding). Inclusion of a point quadrupole of the experimentally observed value failed to improve the molecular orientations in the crystal appreciably. Moreover, it destabilized face-to-face dimer configurations considerably more than did the Williams-Starr9 or Karlstrom et a1.I’ models (with point charges) that gave excellent molecular orientations. Adding six point-dipoles onto the six interaction sites in an attempt to give an improved, delocalized representation of the quadrupole moment was a step in the right direction in crystal and dimer structures, but was not good enough.41 Placing six negative charges at carbon sites and six positive charges on the hydrogens while retaining six LennardJones interaction sites yielded quite good crystal and dimer orie n t a t i o n ~ . It ~ ~was not seriously considered for Monte Carlo simulations, however, because it required a minimum of 12 sites and, when optimized, 18 sites, and therefore offered no real computational advantage over full 12-site models. 12-Site Models. Although a variety of 12-site models for benzene were explored during the course of the research, most attention was directed to two. Because of the special virtues of the (exp-6-1) potential of Williams and Starr (WS77)9 and the (12-9-6-4-1) potential of Karlstrom et al.” as described above, an attempt was made to modify each to a form more appropriate for the present purposes while retaining, if possible, the essential characters of each. Although the WS77 potential yielded an excellent value of E? in CPA (benzene was a calibration point for WS77), it was parametrized to yield the experimental crystal structure at 138 K.36 Therefore it had to be rescaled to give the smaller classically extrapolated volume Vo. The (12-9-6-4-1) potential” gave a value for lE2l that was substantially too large and a volume that was appreciably too small. In addition, it was expedient to simplify the form of the potentials to shorten computation, as discussed above. As a guide in rescaling, a set of seven dimer configurations illustrated in Figure 2 was c o n s t r ~ c t e d . ~In~ the “standard” (39) Siska, P. E.; Parson, J. M.; Schafer, T. P.; Lee, Y. T. J . Chem. Phys. 1971, 55, 5762; 1972, 56, 1511. (40) Maitland, G. C.; Rigby, M.; Smith, E. B.; Wakeham, W. A. Intermolecular Forces; Clarendon: Oxford, U.K., 198 1; Chapter 9. (41) For details, see: Shi, X.Ph.D. Thesis, University of Michigan, 1988.

cc

CH HH ( 12-1 0-6-2)

cc

CH HH

“Units: r in A, energy in kJ/mol. The (12-6) parameters were derived from the (12-6-2) values with the r‘l term deleted, rescaled to give the same crystal binding energy for the same (no longer optimal) structure. Parameters are rescaled as discussed following eq 3 in text to yield V, = 112.0 A3/molecule for all but the (12-6) potental. For benzene, interaction site positions (also rescaled by factor fl are rcc = 1.401 A, rCH = 1.031 A, the latter value incorporating the Williams foreshortening to compensate for asphericity of hydrogen atoms.

I-

,

2

3

4

5

6

7

DISTANCE BETWEEN CENTERS OF MASS (A)

Figure 3. Potential energy curves for benzene dimer according to the (exp-6-1) model of ref 9. Configurations of Figure 2.

configurations whose results will be presented subsequently, rings were in eclipsed orientations relative to each other in the faceto-face configuration 1, and hydrogens pointed toward hydrogens in configuration 7. In other tests, one ring was rotated 30’ to give a staggered configuration 1 and enmeshed configuration 7. Results were for the most part so similar to those of the standard configurations that little reference will be made to them. For all configurations, dimer intermolecular potential energies were calculated as a function of R , the separation between benzene centers of mass. Simplifying the Buckingham plus Coulomb form of the WS77 potential meant adjusting constants in the less flexible form (12-6-1) to yield the reference values for E? and V, while preserving the excellent molecular orientations and lattice constant ratios of WS77. When this proved to be feasible, attention was turned to simplifying even further (in order to avoid square roots in M C computations) to a (12-6-2) form with the r-* term adjusted to yield the correct Coulomb forces at molecular contact distances and with the term truncated at 5 A. This stratagem, which has been incorporated into certain MD simulations reported rests in the present work on the assumption that the most important role of electrostatic interactions is that expressing the short-range preferences of hydrogens and carbons to associate with each other while avoiding contacts with like atoms. It is supposed that the long-range Coulomb effects are unimportant (42) See, for example: Brooks, B. R.; Broccholeri, R. E.; Olafson, B. D.; States, D. J.; Swaminathan, S.; Karplus, M.J. Comput. Chem. 1983.4, 187. The representation of Coulomb potentials by r-2 functions has been used in biophysical simulations for well over a decade.

The Journal of Physical Chemistry, Vol. 92, No. 20, 1988 5671

Intermolecular Interactions for Benzene

4 h

w 1

0

z

%

o

& >-

u

g

-4

W Z J

5 c

z

W

-8

E -12

- , 2

3

5

4

6

7

-1

3

2

8

4

5

6

I

a

7

DISTANCE BETWEEN CENTERS OF MASS (A)

DISTANCE BETWEEN CENTERS OF MASS (A)

Figure 4. Potential energy curves for benzene dimer according to the (12-6-2) model of Table I. Configurations of Figure 2.

Figure 6. Potential energy curves for benzene dimer according to the (12-10-6-2) model of Table I. Configurations of Figure 2.

I +--.--+-\ 9 g

z

-4

1

4

-8 i

-12

:

I

,

‘ii

2

--I0[ 20 7

DISTANCE BETWEEN CENTERS OF MASS (A)

Figure 5. Potential energy curves for benzene dimer according to the (12-9-6-4-1) model of ref 18. Configurations of Figure 2.

because of extensive cancellations and dielectric polarization. The ineffectiveness of point quadrupoles, which do correctly give the long-range electrostatic behavior, further supports the insubstantial role of long-range forces. Results of CPA corroborated these assumptions. Finally, the (12-6-2) potential was further simplified to a 12-6 potential, rescaling the coefficients to yield the same E: for the same (no longer optimal) structure, to see whether discarding the “electrostatic” terms had as deleterious an effect on the liquid structure as on the solid. Final parameters for the model potentials are listed in Table I. Comparisons between profiles of dimer interaction energies for the WS77 and (12-6-2) models are given in Figures 3 and 4. Simplifying and rescaling the (12-9-6-4-1) potential to a (12-1 0-6-2) form involved more compromises but parameters were found by trial and error that made the features of the simplified model mimic those of Karlstrom’s model where there were departures from WS77. Moreover, at the same time, these parameters yielded the adopted reference values of E: and V, while retaining excellent molecular orientations and lattice constant ratios. Profiles of Karlstrom and (12-10-6-2) dimer interaction energies are illustrated in Figures 5 and 6. It should be noted that although the truncated quasi-electrostatic terms in the (12-6-2) and (12-10-6-2) potentials of form

where r,,, = 5.055 A, are continuous functions of r, their derivatives are not. This mildly nonphysical behavior was adopted rather arbitrarily in the interests of economy. An artifact of this discontinuity can be seen in Figures 4 and 6 where the curves for configuration 1 exhibit shallow double minima. A larger value of r,,, than 5 A would be appropriate. In molecular dynamics simulations, the discontinuity of derivatives caused by the excessive truncation would be particularly disadvantageous. To help convey the important role of electrostatic interactions in the packing of benzene molecules, the errors in the three

J

p...

Figure 7. Shift from the observed equilibrium Eulerian angles of reference molecule in benzene (I) crystal as a function of partial changes on C and H sites for (12-6-1) model of Table I. TABLE I1 Crystal-Packing Analyses for Various Model Potentials” model - E , O ~ V, a b c e 6 6-site 40.7 120.8 7.18 9.58 7.03 47.1 0.1 29.8 12-6‘

+

12-site 12-6d*‘ 12-6- 1‘ 12-6-2‘ 12-10-6-2‘ 12-9-6-4-11

(exp-6)* expth

53.8 52.3 52.3 52.3 62.5 52.8 52.3

112.4 112.0 112.0 112.0 99.8 120.3 112.0

6.34 7.37 7.32 7.34 7.41 7.40 7.33

9.17 9.13 9.21 9.10 8.47 9.36 9.23

7.73 6.65 6.64 6.71 6.55 6.91 6.62

60.4 43.0 43.5 43.3 44.2 44.9 40.9

0.1 90.0 20.7 3.6 20.8 3.1 20.0 3.5 22.9 2.3 21.4 3.4 22.8 1.1

“Distances in A; angles in deg; Eulerian convention of reference 43. bCrystal binding energy in kJ/mol, corrected for zero-point vibrations. CParametersof Jorgensen et al., ref 13. dParameters of the 12-6-2 potential, rescaled after f Zterm deleted. See Table I for details. eThis research. fKarlstrom et al., ref 27. gWilliams and Starr, ref 9, parametrized to fit crystal at 135 K, not at 0 K. *Data of ref 34-36 extrapolated “classically”to 0 K and rescaled as outlined in text. Individual lattice constants reliable to about 1%. Eulerian angles of the molecular orientation in the crystal are plotted in Figure 7 as a function of the fractional charge assigned to sites in a 12-6-1 Results of crystal-packing analyses for the principal model potentials constructed are compared in Table I1 with those for the WS77 and Karlstrom potentials. The latter are the only two models in the literature that successfully yielded the correct space group and molecular orientations in the crystal at zero external pressure. Results for the pressure dependence of the WS77 and the present (12-10-6-2) models are listed in Table 111. It can be seen that the (12-10-6-2) model also gives the correct orthorhombic crystal structure at atmospheric pressure and a smaller molar volume for the high-pressure phase 11. This potential, (43) Note that the convention for Eulerian angles adopted in the present work is not the British convention but, rather, the xyz convention of ref 10.

5612 The Journal of Physical Chemistry, Vol. 92, No. 20, 1988 TABLE III: Crystal-Packing Analyses for the Two Crystalline Phases of Benzene at Two Pressures According to Model Potentials@ ( 12-1 0-6-2)

uhase I

model* phase I1

(exp-6-1) model' uhase I uhase I1

P=O

C

-52.28 7.30 9.05 6.68

Bd P

66.44

E,O a

b

-49.03 5.30 5.25 8.38 (1100) 65.88

-52.82 7.40 9.36 6.91 72.05

-52.91 5.68 5.45 7.96 (1100) 69.05

P = 25 kbar EC0 a

b

>

-50.20 7.23 9.00 6.45

-48.93 5.32 5.30 8.19 (1100) 65.34

-49.23 7.29 9.14 6.59

-50.06 5.75 5.42 7.56 (1 100) 66.66

P 63.19 66.11 a Energy in kJ/mol, cell constants in A, molar volume in cm3/mol. bThis work, scaled to V, of hypothetical vibrationless crystal. See text. 'Model of Williams and Starr, ref 9, scaled so that vibrationless crystal has molar volume close to that of real crystal at 138 K. dFollowing Williams, we constrain @ to be 110O. however, exaggerates the stability of phase I in comparison with the results of WS77 and also gives a smaller decrease in volume upon phase change. Therefore, it falls just short of making phase I1 stable at 25 kbar, even though it implies an effect of pressure on [E,"(II) - E,"(I)] that is somewhat larger than given by WS77. Because the (12-10-6-2) potential shares the well-known deficiency of 12-6 potentials in making the repulsive interaction much too steep at small r, it cannot be expected that it will model high-pressure properties as faithfully as a well-chosen (exp6-1) potential which has an extra parameter to control the behavior at small r. In the light of Monte Carlo calculations to be described in paper 1115in this series, it is evident that the value of Vothat was selected by Evans and Wattsl2 and adopted a t the outset of the present work, is too low by about 1%. Therefore, there is some merit in rescaling the present potentials to achieve a somewhat improved result. Available data are not yet sufficiently definitive to establish Vozpwith enough certainty to justify using E," of eq 1 rather than E," as a reference packing energy. Rescaling the energy, of course, can be done by multiplying all coefficients in the model potentials by a constant. Corresponding temperatures in a given prior simulation would be scaled in proportion. Rescaling the effective size of molecules without changing the binding energy is accomplished simply, as follows: Given a site-site interaction energy of form (3) if it is wished to rescale volumes by a factor off3 while maintaining all geometric shape factors (including molecular orientations and lattice constant ratios in CPA), then the coefficients Ak in eq 3 should be multiplied by fk.Of course, the distances of the interaction sites from the molecular center have to be rescaled by the same factor f,as well. Tentatively, as discussed above, we propose that a better value of the classically extrapolated molecular volume Vo is 112 A3/molecule, corresponding tof3 = 1.01 1, and this shift has been incorporated into the rescaled potentials (12-6-2) and (12-10-6-2) in Table I.

-

Discussion The principal aim of this paper has been to develop potential energy functions for intermolecular interactions that are, on the one hand, more realistic than those used in prior simulations of condensed benzene yet, on the other hand, of simple enough form to be economical in long simulations. Every one of the potentials proposed in the literature was known to be deficient in one or more aspects of crystal structure or dimer structure although two, those of Williams and Starr and of Karlstrom et al., exhibited certain very desirable properties. Because it was not yet known whether

Shi and Bartell the requirements for satisfactory simulations of liquid behavior were as stringent as those for crystals, and because it was apparent that none of the 6-site models had been fully optimized, some attention was directed toward 6-site models. As will be shown in paper 111: 6-site models, while able to fit liquid thermodynamic properties, are no better at accounting for liquid structure than solid structure and are mediocre at representing solid-state thermodynamic properties if no electrostatic interactions are included. Twelve-site models, being more flexible than 6-site models, are, of course, superior representations of benzene. Not even they yield adequate dimer and crystal structures unless electrostatic charges are distributed over the molecular framework to reproduce the observed molecular quadrupole moment." A point quadrupole added to the molecular center does not suffice. Whatever the true physical interpretation of the interactions modeled by the distributed charges, it can be concluded that it is not the long-range Coulombic aspect that is essential for success. When the longrange electrostatic r-l potentials were replaced by r-2 potentials adjusted to yield the electrostatic forces at atom-atom contact distances, then truncated at 5 A, the structural results were just as good. They are also considerably cheaper to use than r-' potentials in long simulations. Simplification from (exp-6-1) potentials to the computationally less demanding (12-6-2) representation, or from the (12-p6-4-1) form to (1 2-10-6-2), resulted in little degradation in crystal structure or energy if parameters were judiciously adjusted. It is clear, then, that, while a model potential which cannot reproduce the crystal and dimer structures cannot be considered to be realistic, the ability to reproduce the structure and packing energy is no guarantee that the model potential is correct. Figures 3-6 convey at once the influence of electrostatic charges in destabilizing the face-to-face configuration 1 that would otherwise be favored by dispersion forces. The figures also illustrate the characteristic differences between the WS77 and Karlstrom models. They show that the crystallographer's (WS77) model makes benzene a little flatter and broader and a little less reluctant to tolerate face-to-face approaches than the Karlstrom model. A comparison of Figures 3 and 4, and of 5 and 6, confirms that the simplified representations preserve these differences. One very important property of the (exp-6-1) potential of Williams and Starr is that, even though it is mediocre in modeling the benzene dimer and (as will be seen5) liquid benzene, as well, it is transferable to other organic compounds than benzene. It has proven to be excellent in predicting crystal structures of a variety of substances. It is reasonable to enquire, then, whether the present (12-10-6-2) potential, which is superior to the WS77 potential in modeling certain properties of b e n ~ e n e ,is~ also transferable to a useful degree. The packing analyses summarized for several compounds in Table IV indicate that it is far inferior to WS77 in this respect. The inferior transferability of the (12-10-6-2) potential is probably due in part to chance elements in the nonunique analytical representation of an array of numerical quantum results. The original (12-9-6-4-1) representation was derived solely to fit, approximately, the benzene data while the WS77 representation was based on crystallographic data including many hydrocarbons. It is not unlikely that an alternative representation of the Karlstrom data could be devised that would exhibit better transferability. After the foregoing discussion was written new MD simulations of liquid and solid benzene came to our attention. In this on a 32-molecule system, individual molecules were given all 36 atomic degrees of freedom in order to examine vibrational frequencies in condensed phases.46 Intermolecular interactions were based on a model potential constructed by Powell et aL4' to re(44) Battaglia, M. R.; Buckingham, A. D.; Williams, J. H. Chem. Phys. Left. 1981, 78, 421. (45) Anderson, J.; Ullo, J. J.; Yip, S. J . Chem. Phys. 1987, 86, 4078. (46) ,Because the molecular vibrations were modeled as classical and harmonic, they would be unable to account for the isotope effects in liquid benzene discussed in ref 24.

The Journal of Physical Chemistry, Vol. 92, No. 20, 1988 5673

Intermolecular Interactions for Benzene

TABLE I V Comparison between Experimental Properties for Selected Hydrocarbon Crystals and Properties Calculated bv Model Potentials’ ~~

n-pentane property

expt

Ab

-E,O

41.5d 4.16 9.08 14.86

105.2 4.12 8.81 14.80

a

b C

a

B 7 lattice

orthorhombic

n-hexane

Be

4.18 8.96 14.83

expt

Ab

52.6d 4.179 4.70 8.57 96.6 87.2 105.0

136.9 3.81 4.68 8.53 102.6 75.4 120.2 triclinic

naphthalene

Be

4.19 4.53 8.67 96.5 88.0 103.3

anthracene

expt

Ab

Be

expt

Ab

Be

72.4# 8.24* 6.00 8.66

64.6 8.20 5.81 8.35

8.15 5.92 8.50

102.1d 8.56’ 6.04 11.16

80.0 8.64 5.93 10.96

8.54 5.89 11.15

122.9

122.1

122.5

124.7

124.7

123.4

monoclinic

monoclinic

‘Energies in kJ/mol, distances in A, angles in deg. b(12-10-6-2) potential, this work, charges assigned as in ref 9. C(Exp-6-l)potential of Williams and Starr, from ref 9. dU.S. National Bureau of Standards: Circular C461,US.Department of Commerce. cBradley, R. S.; Cleasby, T. G.J . Chem. Soc. 1952, 1690. ’Mattison, H.;Norman, N . Acta Chem. Scand. 1967, 21, 127. #Norman, N.Acta Chem. Scand. 1961, 15, 1755. *Cruickshank,D. W.J. Acta Crystallogr. 1957, 10, 504. ‘Cruickshank, D.W.J. Acta Crystallogr. 1956, 9, 915; 1957, 10, 470.

produce phonon dispersion curves of crystalline benzene derived by inelastic neutron diffraction. The Powell model, in t u p , was constructed by starting with Williams (exp-6) potential terms devoid of electrostatic interactions and refining coefficients of the attractive and repulsive terms until calculated dispersion curves agreed optimally with experiment. We had initially disregarded this potential as artificial because it deviated in the extreme from the conventional combining rules for atomic potentials and seemed unphysical in relying entirely on C-H interactions for,stabilization. The C-C and H-H interactions were repulsive over the entire range of internuclear distances corresponding to atomic contacts. Because this counterintuitive apportionment of a t o m a t o m contributions to the cohesive interactions between benzene molecules is precisely that which occurs in the (12-9-64-1) Linse model and our own (12-10-6-2) as a result of the “electrostatic interactions”, an examination of the Powell (exp6) model seemed appropriate. It was found that the individual C 2 and increasing 1A21 to 58.68 kJ A2/mol (whence A. = -A2/144 A2),the CPA runs were degraded by only about 1% in cell constants and by 1 or 2 O in angles. An optimization of all A , parameters would recover a better fit. On the other hand, modifying the (12-10-6-2) potential to the form (12-10-6-1) gave a mediocre representation of the crystal if the A, were frozen for n > 1. After this paper was written, two more model potentials for benzene appeared [Pettersson, I.; Liljefors, T. J . Comput. Chem. 1987,8, 1139; Allinger, N . L.; Lii, J.-H. J. Comput. Chem., 1987,8, 11461. Neither of these was subjected to a rigorous CPA. Acknowledgment. This research was supported by the National Science Foundation under Grant No. CHE-8419389. We thank Prof. D. E. Williams for helpful remarks. We gratefully acknowledge a generous allocation of computing time from the University of Michigan Computing Center. Registry No. C6H,, 71-43-2.