Langmuir 1991, 7, 3110-3117
3110
Computer Simulations of Benzene Adsorbed on Graphite. 1. 85 K Alexei Vernovt and William A. Steele* Department of Chemistry, The Pennsylvania State University, University Park, Pennsylvania 16802 Received February 26, 1991. I n Final Form: July 2, 1991
Computer simulations of the structure and energy of benzene layers adsorbed on graphite at 85 K are presented. Models for the benzene-benzene and the benzene-solid interaction potentials are developed. These give a stable 4 7 X 4 7 commensurate film, in agreement with experiment. It is shown that the orientations of these adsorbed molecules are not entirely parallel to the surface as generally assumed. The fraction of nonparallel molecules in the monolayer is small but increases with increasing coverage. The bilayer film is also studied; this film has somewhat less stabilization energy than the bulk crystal, but does not transform into the bulk during the duration of the simulation. Also, the second layer is orientationally less well ordered than the first. Comparisons with available experimental data are made. 1. Introduction
Over the past 15years, numerous computer simulations of the properties of gases adsorbed on the basal plane of graphite have been reported.’ Results have been obtained for the structures and thermodynamic properties of both monolayer and multilayer films. In some cases, dynamical properties such as diffusion constants have also been simulated. Specific adsorbates studied include not only all the rare gases, but numerous molecular adsorbates such as oxygen,zn i t r ~ g e nmethane? ,~ e t h ~ l e n eethane,6 ,~ carbon disulfide,’ and carbon dioxide.8 One of the particular advantages of this approach is that it is capable of yielding detailed information about molecular orientation in adsorbed films and its role in determining observable properties such as phase transitions, heats of adsorption, and monolayer capacities. This information can be difficult to come by with available experimental methods, especially if one is interested in “out-of-plane” orientations-that is, in the orientation of some molecular symmetry axis relative to the surface plane of the adsorbent. In-plane orientational ordering in solid monolayers can be accessible through careful diffraction studies, but even in this case, data of adequate accuracy are sparse. Experimental studies of benzene adsorption on graphite have a relatively long history, with careful isotherm and heat of adsorption measurements at room temperature having been reported in 19619and repeated several times since Of course, the shape of benzene leads one + Permanent address: Department of Physical and Colloid Chemistry, People’s Friendship University, M. Maklaya 6,117198 Moscow,
U.S.S.R. (1) Nicholson, D.; Parsonage, N. G. Computer Simulation and the Statistical Mechanics of Adsorption; Academic Press: New York, 1982. (2) Dan, R. P.; Etters, R. D.; Kobashi, K.; Chandrasekharan, V. J . Chem. Phys. 1982,77,1035. Bhethanabotla, V.; Steele, W. A. Langmuir 1987, 3, 581; Can J. Chem. 1988, 66, 866; Phys. Reu. B 1990, 41, 9480. (3) Tildesley, D. J.; Vernov, A. V.; Steele, W. A. Carbon 1987, 25, 7. Talbot, J.;Tildesley, D. J.; Steele, W. A. Mol. Phys. 1984,51,1331. Joshi, Y. P.; Tildesley, D. J. Mol. Phys. 1985, 55, 999. Peters, C.; Klein, M. Phys. Reu. B 1985, 32, 6077. Cardini, G.; O’Shea, S. F. Surf. Sci. 1985, 154,231. Etters, R. D.; Roth, M. W.; Kuchta, B., private communication. (4) Severin, E. S.;Tildesley, D. J. Mol. Phys. 1980, 41, 1401. (5) Moller, M. A.; Klein, M. Can. J.Chem. 1988,66,774; Chem. Phys. 1989, 129, 235. Nose, S.; Klein, M. Phys. Reu. Lett. 1984, 53, 818. (6)Moller, M. A.; Klein, M. J. Chem. Phys. 1989, 90,1960. (7) Joshi, Y. P.;Tildesley, D. J.; Ayres, J. S.; Thomas, R. K. Mol. Phys. 1988, 65, 991. (8) Hammonds, K. T.; MacDonald, I. R.; Tildesley, D. J. Mol. Phys. 1990, 70, 175. (9) Isirikyan, A. A.; Kiselev, A. V. J. Phys. Chem. 1961, 65, 601. (10) Pierotti, R. A.; Smallwood, R. E. J. Colloid Interface Sci. 1966, 22, 469.
0743-746319112407-3110$02.50/0
to expect that there might be pronounced orientation effects in these adsorbed layers.”-19 However,the presence of a large quadrupole moment oriented perpendicular to the molecular plane complicates one’s intuitive picture and in fact various conflicting conclusions concerning benzene orientations in the monolayers have been published. The most striking feature of the original heat of adsorption studies was that there is an extremely small variation in heat with increasing surface density, When the substrate is nearly homogeneous, as in the graphitic adsorbents used, one expects an increasing heat of adsorption that is roughly linear in surface density at high temperature (above the 2D critical temperature) with a slope that correlates with the heat of vaporization of the bulk adsorbate. Data for spherical or nearly spherical molecules support thisz0and indicate that one is observing the lateral interaction energy of molecules in the monolayer. With this thought in mind, it was argued that the observation of a nearly zero slope for benzene (and a “normal” slope for cyclohexane) meant that the benzene molecules were lying flat on the surface, with an edge-wise interaction that was much weaker than that for other orientations. In fact, bulk solid benzene is composed primarily of nearest-neighbor pairs that are perpendicular to one another2’ an orientation that is most favorable for quadrupolar interaction. In the absence of such pairs, the average interaction should be much weaker. This interpretation of the room-temperature data was thrown into doubt by low-temperature NMR studies that seemed to indicate that the molecules were oriented perpendicular to the surface.z2 However, recent reinterpretation of these experiments has concluded that these ~
(11) Pierotti, R. A. Chem. Phys. Lett. 1968, 2, 420. (12) Pierce, C.; Ewing, B. J. Phys. Chem. 1967, 71, 3408. (13) Pierce, C. J . Phys. Chem. 1969, 73, 813. Faraday (14) Dollimore, D.; Heal, G. R.; Martin, D. R. J. Chem. SOC., Trans. 1 1972, 68, 832. (15) Katir, J.; Coulon, M.; Bonnetain, L. J. Chin. Phys. 1978, 75,789. (16) Delachame, J. C.; Coulon, M.; Bonnetain, L. Surf. Sci. 1983,133, 365. ~~
(17) Avgul, N. N.; Kiselev, A. V.; Lygina, I. A.; Poshkus, D. P. Izu.
Akad. Nauk SSSR 1959. - - - ,-1196. ---
(18) Bondi, C:;-Baglioni, P.; Taddei, G. Chem. Phys. 1985, 96,277. (19) Bondi, C.; Taddei, G. Surf. Sci. 1988, 203, 587. (20) Steele, W. A. The Interaction of Gases with Solid Surfaces; Pergamon Press: New York, 1974; Section 4.9. (21) Cox, E. G.; Cruikshank, D. W. J.; Smith, J. A. S . Proc. R. SOC. London 1958, A247, 1. (22) Boddenberg, B.; Moreno, J. A. Z. Naturforsch. 1976,3IA, 854; J. Phys. 1977, 38 (C4), 52; Ber. Bunsen-Ges. Phys. Chem. 1983,87, 83.
0 1991 American Chemical Society
Computer Simulations of Benzene Adsorbed on Graphite
Langmuir, Vol. 7, No. 12,1991 3111
data are consistent with parallel 0rientations.~3.~* Further indirect support for the hypothesis of surface-parallel orientations has been supplied by diffraction measurements25-28on the monolayer at 85 K, which show a commensurate d 7 X 4 7 ordered layer with one in seven carbon hexagons in the surface occupied. Although reasonable estimates of the benzene molecular area are consistent with this area per molecule for parallel orientations, detailed energy minimization calculations using the known benzene-benzene interaction potentials give an ordered layer of coparallel molecules with spacing somewhat larger than that for the 4 7 X d7.'8J9,27 This question is not the only point of controversy concerning the benzene-graphite system. It appears that benzene coverages past monolayer do not wet the surface at temperatures below 230 K; i.e., second or higher layer formation is preempted by the formation of bulk cryst a l l i t e ~ . ' ~Of ~ ~course, ~ . ~ ~if the first layer is made up of coparallel molecules, this provides a very poor template for further layer growth, which must eventually give films with structures approaching that of bulk crystal. Thus, the structure and energetics of multilayer films above and below the wetting temperature of 230 K becomes an interesting subject for investigation by computer simulation. Monolayer melting in the temperature interval of 130150K has been sugge~ted~~s~~*2*-one can then ask whether this melting is first order or not, whether the order of the melting transition depends on coverage, and what role orientation plays in this case. Finally, a significant temperature dependence of the heats of adsorption in the benzene monolayer has been observed at 233-243 K15J6and has been interpreted to be a consequence of molecular reorientation. The computer simulations reported in this paper will address several of these points.
of the energy plus inclusion of the quadrupolar electrostatic energy. The site-site models are often 12-sitedescriptions that take both the C and the H as spherically symmetric centers of interaction; six nonspherical sites have been suggested,33and in some cases, the position of the H sites has been optimized by taking the C-H site distance to be slightly different from the C-H bond length. As a rule, the parameters of these site-site potentials have been chosen to give crystal structure and energy in agreement with the data for solid benzene. Because of the prevalence of T-shaped and herringbone orientations for the nearestneighbor pairs in the solid, the values and separations of the minima for the potentials in the literature are generally in agreement with each other for pairs with these orientations. However, not much effort has been put into the development of models that give experimentally adequate results for edge-on pairs. To illustrate this point, we present a number of energy calculations here, both for the coparallel monolayer and for the dimer energies with various relative orientations. The models considered here are as follows: Williams32 (the numerical results for the Williams and the Allinger36 model are nearly identical), Karlstrom3' (actually, an ab initio quantum calculation), and Price.33 The Williams model is 12-site but the H sites set at 1.027 A away from the C atoms rather than at the bond length of 1.084 A; Allinger took this distance to be 1.019A. For the site-site interactions, an exponential-inverse sixth-power distance dependence was taken, with somewhat different values of the parameters in the two cases. Although both representations of the electrostatic energy were adjusted to give the correct molecular quadrupole moment, Williams assumed a distribution of point charges located on the C and H sites, and Allinger took point dipoles located at the center of each C-H bond. The ab initio quantum mechanical calculations of benzene-benzene energy are not readily describable in terms of either site-site or electrostatic terms; however, Karlstrom's results are given in sufficient detail to allow us to evaluate this potential for various dimer configurations. In the model of Price et al. six anisotropic sites are placed a t the C atoms together with six point quadrupoles oriented perpendicular to the molecular plane. In Figure 1, curves of model energy versus distance are shown for a T-shaped pair (A), two edge-on orientations (B and D), and a parallel configuration (C). All these energies vary somewhat as the molecules in the pair are rotated around their symmetry axes-see B and D, for instance. However, the numerical changes are minimal, amounting to 100 K or less in the minima. Although all four potentials (including Allinger's) are reasonably close to one another for the T pairs, significant differences are present for the coplanar pairs. In this case, the ab initio calculation of Karlstrom gives deeper minima at smaller molecule-molecule separations than the semiempirical models. When one comes to calculate the lattice energy of a triangular two-dimensional array of coplanar molecules, these differences produce marked differences in the minimum energy and in the lattice spacing corresponding to the minimum. As an example, this minimum lattice energy for the Williams model is plotted as a function of a , the nearest-neighbor distance, in Figure 2. (These energies were minimized with respect to azimuthal orientation.) The minimum energy of -4300 K occurs at a = 6.95 A, a distance that is considerably in excess of the experimental result of 6.51 8, for the d 7 X d?commensurate monolayer. Of course, the benzene-graphite energy is lower for the commensurate than for the noncommen-
2. Molecular Interactions To simulate this system, we need realistic models for both the benzene-benzene and the benzene-graphite interactions. In the benzene-benzene case, there have been rougly 20 different models proposed over the past few y e a r ~ . ~ ~ 1 ~In+ a3 6compromise between ease of use in a simulation and realism of representation of molecular shape and size, it seems that the best of these are based on a site-site representation of the nonelectrostatic part (23) Boddenberg, B.; Grosse, R. Z . Naturforsch. 1986,41A, 1361. (24) Boddenberg, B.; Grosse, R. 2.Naturforsch. 1988,43A, 497. (25) Stockmeyer, R.; Stortnik, H. Surf. Sci. 1979,81,1979. Mondenbusch, M.; Stockmeyer, R. Ber. Bunsen-Ges. Phys. Chem. 1980,84,808. (26) Meehan, M.; Rayment, T.; Thomas, R. K.; Bomchil, G.; White, J. W. J.Chem. Soc.,Faraday Trans. 1 1980,76,2011. Tabony, J.; White, J. W.; Delachaume, J. C.; Coulon, M. Surf. Sci. 1980, 95, L282. (27) Gameson, 1.; Rayment, T. Chem. Phys. Lett. 1986, 123, 150. (28) Bardi. U.: Mannanelli, S.; Rovida, G. Surf. Sci. 1986, 165, L7; Langmuir 1987,3, 155. (29) Evans, D. J.; Watts, R. 0. Mol. Phys. 1975,29,777; 1976,31,83. (30) Serrano-Adan,F.;Bafion, A,;Santamanir, J. J. Chem. Phys. 1984, 86._433. _ .-. ,
(31) Karlstrom, G.; Linse, P.; Wallquist, A,; Jbnsson, B. J.Am. Chem. SOC.1983,105, 3777. (32) Williams, D. E.; Cox, S. R. Acta Crystallogr. B 1984,40, 404. ( 3 3 ) Yashonath, S.;Price, S.L.; MacDonald, I. R. Mol. Phys. 1988,64, 361. Price, S. L.; Stone, A. J . Mol. Phys. 1984, 51, 569. (34) Gupta, S.;Sediawan, W. B.; MacLaughlin, E. Mol. Phys. 1988,65, 961. MacRury, T. B.; Steele, W. A,; Berne, B. J . J. Chem. Phys. 1976, 64. 1288. Kabadi, V. N.: Steele, W. A. Ber. Bunsen-Ges. Phys. Chem. 1985,89, 2,9. (35) Ciccotti, G.; Ferrario, M.; Ryckaert, J.-P. Mol. Phys. 1982, 47, 1253. Claessens, M.; Ferrario, M.; Ryckaert, J.-P. Mol. Phys. 1983,50, 21736) Allinger,N.L.;Yuh,Y.H.;Lii,J.-H.J.Am. Chem.Soc. 1989,111, 8551. Lii, J.-H.; Allinger, N. L. J.Am. Chem. SOC.1989, I l l , 8576. Allinger, N. L.; Lii, J.-H. J. Comput. Chem. 1987.8. 1146. Petterson, I.; Liljefors, T. J. Comput. Chem. 1987, 8, 1139.
-
Vernov and Steele
3112 Langmuir, Vol. 7, No. 12,1991
Table 1. Well Depth and Size Parameters for the Benzene-Benzene Site-Site Interaction Used in This Work
A
~~
sites C-C
H-H
D
C
Y
1
Figure 1. Model calculations of benzene-benzene interaction energy u shown as a function of center-of-massseparation r for the four different orientations indicated. Models: AllingerWilliams (- -); Price (- -); Karlstrom (- - -); this work (-).
Wi llioms This work
I
6.5
(a)
7.0
'
Figure 2. Potentialenergyof a twedimensionaltriangularlattice of coplanar benzene molecules plotted as a function of nearestneighbor distance a for the Williams model and for the model suggested in this work. The energy has been minimized with respect to azimuthal orientations in each case. The v'7 X 4 7 commensurate lattice spacing is shown by the dashed vertical line.
surate film; however, in order to stabilize the WilliamsAllinger layer a t the commensurate distance, the added gas-solid energy needed is a t least 400 K, which is roughly 1 order of magnitude larger than our calculated value (see below). Thus, we have developed a new semiempirical model for the benzene-benzene interaction having the virtues that it agrees fairly well with Karlstrom's ab initio calculation, it is simple enoughto be useful in a simulation, it does not differ greatly from other potentials for T-shaped pairs and thus should give a reasonable representation of bulk solid benzene properties, and it yields a stable lowtemperature 4 7 X 4 7 monolayer.
cs-slk.
33 12
K
us-s.
A
3.65 2.8
We adopt a model of 12 spherical sites with LennardJones 12-6 interaction function. The sites are located on the C atoms and on the C-H bond vectors at a distance 0.8 A away from the C atoms. Well depth and size parameters for these sites are listed in Table I-values for unlike C-H pairs are obtained from Lorentz-Berthelot combining rules. The electrostatic part of this potential was represented by six ideal quadrupoles located on the C atoms perpendicular to the molecular plane. Since the molecular quadrupole moment is known to be 8.5 X esu, each site quadupole was set equal to 1.4 X esu. The dimer energies given by this model are shown by the solid curves in Figure 1, where one can see that calculations for all models considered are similar for T dimers, but that this new potential is generally in good agreement with the Karlstrom calculation for coplanar pairs. It differs from the Williams-Allinger and Price models in that it allows a closer approach and a deeper well for edge-on dimers. In fact, we will show that the new potential can produce a stable 4 7 X 4 7 commensurate monolayer. The dashed curve of Figure 2 shows a minimum energy for a coplanar 2D lattice equal to -6250 K at a nearest-neighbor distance of 6.62 A. The added gas-solid energy required to stabilize the commensurate lattice, which has a spacing of 6.51 A, is only 40 K. We now develop a molecule-solid potential model that gives the requisite energy change. One of the advantages of a site-site interaction model for adsorbate molecules is that it can readily be adapted to describe the molecule-graphite interactions. In the present case, one assumes that the graphite is composed of carbon sites, each of which interacts with the sites in a benzene molecule.37 The site-site energies are taken to be Lennard-Jones 12-6 functions with well depths and size parameters estimated by applying Lorentz-Berthelot combining rules to the C-C well depth (elk = 28 K) and size (u = 3.4A) for the graphitic sites. The summations over the semiinfinite graphite solid can be transformed into a Fourier series in x , y, the benzene coordinates in the plane. Thus, the moleculeaolid potential u,(r,P)is written as20
u,(r,P) = uo(z,P) + u,(z,@)fI(x,Y)
(2.1)
where uo(z,@)is a sum over planes of the surface-averaged interation with a plane for a molecule whose axis makes an angle fl with the surface normal, ul(z,p) is a sum of Bessel functions, fl(x,y) is a periodically varying function that reflects the symmetry of the graphite lattice, and z is the distance between the benzene center of mass and the surface. In principle, these are merely the leading terms in an infinite series, but calculation shows that an adequate representation of the benzeneaolid potential is given by eq 2.1. The resulting interaction has a deep minimum in uo(z,P) that is essentially the energy of adsorption in the limit of zero coverage-see below. There is a small periodic variation in the minimum when the molecule moves across the surface. Figure 3 shows this (37) Kiselev, A. V.; Poshkus, D. P. Trans. Faraday SOC.1963,59,428. Kieelev, A. V.; Porthkus, D. P. J. Chem. Soc., Faraday Trans. 2 1976,72, 960. Kiselev, A. V.; Poshkus, D. P.; Grumadas, A. J.J. Chem. Soc.,Faraday Trans. 1 1979, 75,1281. Kiselev, A. V.; Poshkus, D. P.; Grumadas, A. J. J. Chem. Soc., Faraday Trans. 1 1975, 75, 1288.
Langmuir, Vol. 7, No. 12, 1991 3113
Computer Simulations of Benzene Adsorbed on Graphite
- 1
-140-
4.5
/i'
1
5
.-c E
i
X
" -200-220-
J
-
Figure 3. Variation in minimum adsorption energy as benzene slides across a graphite surface shown for the model developed here. The letter h denotes a benzene centered over a graphite hexagon; the letter b denotes the benzene centered over a C-C bond; and the letter c denotes a benzene over a C atom. The curve with the open squares is for the system with no surface quadrupoles, while that with the filled squares is for the model including graphite quadrupole moments.
variation in minimum energy. If one remembers that the benzene-olid energy for a molecule in an incommensurate monolayer is just equal to the minimum value of U O ( Z , @ ) , with @ = 0" for a coplanar layer, the added energy of formation of a commensurate layer is just ul(z,p) evaluated at the distance of minimum energy, multiplied by fl(x,y) for a commensurate layer. A calculation of the minimum energy for an incommensurate coplanar layer gave -5200 K, while the commensurate minimum is shown to be -5215 K in Figure 3. The resulting commensurate stabilization energy is a very small 15 K. However, this model does not include electrostatic benzene-graphite energies. We38 (and others39) have suggested that the graphite surface is quadrupolar, by analogy with the measured moments for several aromatic hydrocarbons. The Fourier expansion for the molecular quadrupole-graphite quadrupolar array has been presented elsewhere;38an evaluation of this energy at the distance of minimum total energy gives the second curve in Figure 3. The inclusion of this electrostatic interaction does not contribute to UO(Z,@),but gives a somewhat larger energy of stabilization of -40 K for the commensurate phase. Thus, we can conclude that the newly devised interaction potentials should be capable of producing the experimentally observed v'7 X v'7 phase if all molecules are parallel to each other and to the surface. It is of interest to evaluate the change in minimum energy and in the position of this minimum as an isolated benzene is reoriented. Calculated curves showing these quantities are given in Figure 4 for our potential. It is clear that strong hindrance to tumbling of the benzene plane is present in for monolayer molecules. A barrier to such motion amounting to 2600 K is obtained. In addition, a significant shift in molecule-surface distance of minimum energy occurs as the molecular plane is tilted. 3. Thermodynamic Properties at Low Coverage As a check of this proposed gas-olid potential, the Henry's law constant KH and the heats of adsorption qst(0) for an isolated benzene molecule on the graphite surface (38) Vernov, A.; Steele, W. A., to be published. (39) Cracknell, R. F. Ph.D. Thesis, Imperial College, London, 1990. Nicholson,D.;Cracknell,R. F.;Parsonage, N. G.Mol. Simul. 1990,5,307.
P
P
Figure 4. Variation in the minimum benzene-graphite interaction energy uswith orientation of the molecular plane (lefthand panel). Here, 0 is the angle between the molecular axis and the surfacenormal. The right-handpanel shows how the location of the minimum changes with 0. In this case, zminis the distance between the benzene center of mass and the surface plane when the energy is at its minimum. Table 11. Calculated Thermodynamic Properties for Adsorbed Benzene at Low Coverage
T,K
KH,m
sSt(O), kcal/mol
289 353 373
34.1 3.21 1.04
9.1 9.48 9.39
were c a l ~ u l a t e and d ~ compared ~ ~ ~ ~ ~with ~ ~the limiting zero coverage experimental results. These values were calculated from
where A is the area of a triangle (l/s of the hexagon) on the graphite surface, N A is Avogadro's number, r and w characterize the position of the benzene molecule above the triangle and its orientation, and Vis the volume of the prism with height 10 i% above the triangle. The six dimensional integrals over all positions and orientations of one molecule were evaluated by using Simpson's rule as well as Gaussian with quadratures with 16 and 32 points. The calculations gave essentially the same results. Calculations of these thermodynamic properties were performed at 298,353, and 373 K and the results are shown in Table 11. These values are in very good agreement with experimental m e a s ~ r e m e n t s , 4which ~ - ~ ~ give qst(0)= 9.42 kcal/mol in the temperature interval 343-438 K46 and 9.7-10 kcal/mol at T = 298 KS9J0Experimental values of KH from different sources are more variable, but lead to a value in the range of 3-7 bm at T = 353 K.4°v42945Table I1 shows that our calculations lies in this interval. Thus, (40) Grumadas,A.;Poshkus, D. P.Russ.J.Phys. Chem. (Engl.Transl.) 1979,53,1378. Kiselev, A. V.; Poshkus, D. P.; Sherbakova, K. D. Rues. J. Phys. Chem. (Engl. Transl.) 1986,60,797. (41) Vidal-Madjar, C.; Gonnard, M.-F.; Guiochon, G. J . Colloid Interface Sci. 1975, 52, 102. (42) Vidal-Madjar, C.; Gonnard, M.-F.; Goedert, M.; Guiochon, G. J . Phys. Chem. 1975, 79, 732. (43) Elkington, P. A.; Curthoys, G. J . Phys. Chem. 1969, 73, 2321. (44) Ross, S.; Saelens, J. K.; Olivier, J. P. J . Phys. Chem. 1962,66,696.
(45)Kalaschnikova, E.V.; Kiselev, A. V.; Petrova, R. S.;Scherbakova, K. D. Chromatographia 1971,4,495. (46)Battezatti, L.; Pisani, C.;Ricca, F. J.Chem. SOC.,Faraday Trans. 2 1975, 71, 1629. (47) Vernov, A. V.; Gorelov, D. S.;Sapata, H. V. Zh. Fir. Khim. 1989,
63, 696.
Vernov and Steele
3114 Langmuir, Vol. 7, No. 12, 1991
50 1
I
1
0
5
10
1
Figure 5. Time dependence of the average temperatures obtained for benzene molecules adsorbed on graphite at 298 K and low coveragefrom the Evans isokineticalgorithm. The lower fluctuating plot is (muZ2)/k, where vz is the velocity normal to the surface, m is the molecular mass, and the brackets denote an average over all benzenes at a given time t. The upper plot is the sum ( mvX2+ mu,Z)/k and the constant line denotes of all velocity components.
6
4
time lpicosec)
2
(4)
Figure 6. Calculations of the density of adsorbed benzene molecules at a given distance z between the center of mass and the solid surface. The dotted line is given by Henry’s law p( z ) / p , = (exp - u,(z,B)/kT)and the solid curve is obtained from a simulation at very low surface density. Here, T = 298 K and pm denotes the density of the dl X dl commensurate phase.
momentum P have the form
we have some confidence that our semiempirical potential will produce realistic simulations of benzene adsorbed on graphite. 1
4. Simulation Details
In previous studies of physisorbed layers, it was found convenient to use the constant-temperature algorithm proposed by Evans and Morriss.48 In this algorithm the equations of motion are solved subject to the constraint that the total kinetic energy of translation and rotation of the molecules remains constant. The same algorithm was used here for the benzene simulation, but during the simulation it became clear that an incorrect distribution of kinetic energy among the degrees of freedom was being obtained at low surface density. These effects are illustrated in Figure 5 for the benzene-graphite system at 298 K and at a two-dimensional density p* = 0.098 (p* = p / p , where pmis the density corresponding to the commensurate 4 7 X 4 7 phase). As can be seen, the average kinetic temperature is properly constrained to remain constant during the total time of simulation. However, more detailed evaluation shows that the total is only the mean value between a fluctuating high temperature corresponding to movement of molecules in the X-Y plane (solid line) and a considerably lower temperature corresponding to the Zcomponent of the kinetic energy (broken line). In the limiting case of zero coverage (Le., one molecule on the surface), the standard Evan’s method of simulation also gives an incorrect distribution of molecules above the surface, compared to a Henry’s law evaluation. Figure 6 showsthe two distribution functions obtained in this limit. Here, the solid line represents the distribution of centers of mass, Boltzmann-averaged over all orientations and locations in the X-Y plane. The broken line in the figure shows the estimation of the same distribution obtained from the simulation. The observed differences are primarily due to the existence of different temperatures for the Z and for the X, Y degrees of translational freedom. To overcome these problems we have modified the constant-temperature algorithm to constrain the kinetic energy of nonequivalent degrees of freedom separately. In this case, the equations of motion for the translational ~
(48)Evans, D. J.; Morriss, G. P. Comput. Phys. Rep. 1983,%A, 433.
i
(4.4)
Here, f i is the force on molecule i. In addition, separate constraints were applied to the nonequivalent rotational degrees of freedom of the benzene molecule. For example, the constrained equation for the 2 component of angular velocity was written as
here Iix, l i y , and l i , are the components of the moment of inertia tensor, wix,wiy,and wit are the Componentsof angular velocity, and Ni, is the z component of torque. All variables are written in the principal axis frame for the molecule. Similar expressions were used for the X, Y rotational motion. Equation of motions for the orientational degrees of freedom were expressed in terms of quaternionsP9 which are widely used now in computer simulation of molecular systems. Periodic boundary conditions and the minimum image convention were used in the X-Y plane. A time step of length 2.16 X ps was used in all calculations. A fifth-order predictor-corrector algorithm was employed for the integration of the equation of the translation -
-
(49) Allen, M.;Tildesley,D. J. Computer Simulation ofliquidu;Oxford University Press: New York, 1987.
Langmuir, Vol. 7, No. 12,1991 3115
Computer Simulations of Benzene Adsorbed on Graphite
Figure 7. A top view of the minimum energy configuration of a benzene monolayer adsorbed on graphite. This lattice was obtained by running the molecular dynamics simulation at a very low temperature (10 K) and shows the sample in its periodically replicated computer cell.
..
Figure 9. Same as Figure 8, but for a complete monolayer and a duration of 39 ps. p*:OO98
,
p*:0196
p*=050
-_
-
p*.100
~
5
IO I
p"
= 1.14
I
.
.
(HI
p* = 1.29 - ,
5
IO
p* = 1.50
- -
Figure 8. Top and side views of the computer-generatedcenterof-masstrajectoriesare shown for l / 2 layer of benzene on graphite at 85 K. The computer cell is indicated by the rectangle. The time duration of these trajectories is 78 picoseconds. motion and a fourth-order predictor-corrector method was used for the integration of the equation of rotational motion. Equilibration runs of 3000 time steps (or more) were performed; data were then taken from runs of an additional 6000 time steps (or more). Four series of simulations were performed: (1) over a temperature interval of 85-298 K for a coverage p / p , = 0.5;(2) in the same temperature interval at monolayer coveragepm = 1;(3)for a fixed temperature of 85 K, variable coverage ranging up to 2 monolayers; (4) for a fixed temperature of 298 K, variable coverage ranging up to 2 monolayers. Reorientational time-correlation functions and twodimensional diffusion constants were evaluated together with thermodynamic and structural properties. Much of this work willbe discussed elsewhere. Here, we concentrate on the low-temperature simulations.
5
IO
5
IO
5 I
IO
5
IO
(1,
Figure 10. Center-of-massdensities in arbitrary units as a function of the benzene-graphite surface distance for various values of the surface coverage p* in units of the monolayer d 7 X
4 7 value.
5. Results at 85 K A top view of one monolayer of benzene at ita minimum energy configuration in the simulation cell is shown in Figure 7. The azimuthal or in-plane orientation angle is defined so that each molecule in Figure 7 has 9 = 40".The 2D lattice shown in this figure corresponds to a 4 7 X 4 7 monolayer phase. However, our simulations at 85 K were run over a range of coverages from -0.1 to 2.0 commen-
3116 Langmuir, Vol. 7, No. 12, 1991 I
-
Vernov and Steele
I
surface, but more importantly, because the orientations of these second-layermolecules are distributed over a range of values. Note that the quadrupolar energy is repulsive for pairs of molecules that are parallel, giving rise to a weak minimum in the total interaction relative to that for T orientations, as shown in panels C and A of Figure 1. Experimentally, it has been argued that second-layer formation at 85 K is unstable relative to formation of bulk crystal. Although we have not confirmed this, it should be noted that the average potential energiesof the benzenegraphite system that were evaluated during the simulation indicate that the two-layer system has less stabilization energy than the bulk crystal. This is shown in Figure 11, where the energy is plotted as a function of surface coverage. Of course, the thermodynamic stability of this system depends upon free energy rather than the stabilization energy, but one expects the entropy of the multilayer film to approach that of the bulk. In this case, the simulations would not be inconsistent with the experimental observations. However, entropies or free energies are notoriously difficult to calculate from molecular dynamics simulation data. The average potential energy Uis maieup of two terms: the average molecule-solid interaction UmS and the average molecule-molecule energy Umm. All three energies are plotted versus coverage (or surface density) in Figure 11. At very low coverages, the film is gaslike and a steep rise in U m m with increasing density is observed. Subsequently, a 2D crystal forms and increasing coverage merely increases the fraction of the surface covered by this crystal. In this region, the lateral interaction energy is nearly constant and equal to --1100 K per molecule. At coverages greater than monolayer, -U" gradually increases as the number of neighboring benzenes per adsorbed benzene increases. It is interesting to note that the average molecule-solid energy on this uniform surface decreases gradually with increasing coverage. For coverages greater than monolayer, this clearly reflects the smaller interaction for molecules adsorbed at larger distances from the surface (see the density curves in Figure 10). However,the change at densitiesless than monolayer is primarily due to changes in the fraction of surface-parallel molecules. As the layer builds up, there is a gradual decrease in this fraction, which accounts for the energy change. Evidently, there is considerable cancellation between the changes in Omsand Ummso that the total energy U is nearly constant for all
I
1
h 3.0
0.oy
I 0.5
I
I 1.5
1.0
I
2.0
P
Figure 11. Simulation values of the average potential energy of an adsorbedbenzene molecule at 85 K (circles) together with the (triangles)and calculated energy d-ue to benzene-benzene benzene-graphite U,, (squares) interaction. The dashed line indicates the experimental value of the crystal energy of solid benzene at 85 K.
om,
surate layers. At coverages greater than -0.2, a condensed 2D solid phase formed. Pictures of the center-of-mass trajectories generated in the simulations for l / 2 and 1 monolayer are shown in Figures 8 and 9. The top views show 4 7 X d 7 commensurate layers in both cases, with some minor drift of the entire crystal along paths of minimum energy,which is around the periphery of a carbon hexagon in our model. The side views of the trajectories indicate that a small fraction of the molecules in the condensed phase are not parallel to the surface. This can be seem more clearly in the curves of density versus molecular center-of-mass distance that are given in Figure 10. As one might expect from the data in Figure 3, the small shoulders seen on the peaks in density are due to nonparallel molecules whose centers move out from the surface as they reorient away from surface-parallel. (Of course, the main peaks are due to surface-parallel molecules.) For coverages greater than monolayer, second layer peaks at z e 7.5 8, are visible in the plots of Figure 10. These peaks are rather broad in part because the underlying molecules are not all at the same distance from the p * = 0.50
pf.100
r
-
2.2 -
3.0
1.8 -
2 5-
1.0 -
'("
0.6 -
0.2
I
I
20-
1.4 -
P(+)
1
1.0 1.5
-
0.5 I
I
0
20
+
I
I
40
60
00-1
0
I
I
I-
20
40
60
Figure 12. Plots of the distribution of azimuthal (or in-plane) orientation angle @ for the l / 2 and the monolayer benzene films. The perfect 2D crystal lattice found by energy minimization and pictured in Figure 7 has all azimuthal angles equal to 40'. The 6-fold in-plane symmetry of benzene limits the range of 4 to 60°, as shown.
Computer Simulations of Benzene Adsorbed on Graphite
Langmuir, Vol. 7, No. 12, 1991 3117
coverages greater than the 2D condensation point of p l p , rr, 0.2, for 85 K. Plots of the distributions of azimuthal angle for two surface coverages are shown in Figure 12 and indicate that the thermal motion at 85 K is insufficient to destroy the in-plane orientational order of this crystal.
average benzene-benzene interaction. As a consequence, the calculated net energy of adsorption hardly varies with coverage over a range extending from 0.2 to 2.0 layers. These results are in agreement with the available experimental data. Although the simulations of bilayer benzene films show no instability over the time span of the computer run, there are some indications that the bulk crystal may be more stable than the films. Of course, the translational diffusive motions required to transform a thick film into a monolayer plus bulk are extremely slow a t 85 K. The energies indicate that such a transformation could well occur in a simulation, in agreement with experiment, if run times could be made sufficiently long. Extensive simulations of the melting behavior, the translational and orientational dynamical properties, and the structures of the ' 1 2 layer and the monolayer films at temperatures ranging from 90 to 300 K will be reported elsewhere, as will a comparison of simulations of the coverage-dependent film properties with experimental data a t 300 K.
6. Conclusions By making minor alterations in a 12-site model of the benzene-benzene potential, an interaction law has been developed that should give a good calculation of the properties of bulk solid benzene while producing the d 7 X d7 commensurate 2D solid that has been observed at low temperature on graphite. In doing so, inclusion of the benzene-surface quadrupolar energy was found to help stabilize this phase, but does not appear to have a large effect on the total energy of the adsorbed phase. The use of this potential in computer simulations of the coverage dependence of the structure and energy of the adsorbed phase indicates that there are significantnumbers (- 10%) of molecules that are more perpendicular to the surface than parallel, at monolayer coverage. In this system, the increasing presence of nonsurface parallel benzene molecules causes a decrease in the magnitude of the benzene solid energy, which is largely canceled by changes in the
Acknowledgment. This work was supported by Grant DMR-8718771from the Division of Materials Research of the N.S.F. Registry No. Benzene, 71-43-2; graphite, 1782-42-5.
