+
+
Langmuir 1996, 12, 2495-2500
2495
Molecular Dynamics Simulation of Benzene on Graphite: 1. Phase Behavior of an Adsorbed Monolayer Mark Alan Matties and Reinhard Hentschke* Max-Planck-Institut fu¨ r Polymerforschung, Ackermannweg 10, 55128 Mainz, Germany Received September 14, 1995. In Final Form: February 13, 1996X We present the results of a classical molecular dynamics simulation for a film of benzene adsorbed onto the basal plane of graphite at monolayer coverage over a range of temperatures. The phase behavior of the monolayer is examined in terms of the temperature dependence of several static and dynamic properties. The results indicate a liquid-solid transition around 140 K which compares well with the published experimental phase behavior.
1. Introduction Organic molecules adsorbed onto solid surfaces exhibit unique properties, especially as the thickness of the adsorbed film becomes small. The configuration of benzene adsorbed onto the graphite (0001) plane has been of particular interest in computational and experimental studies. In theoretical and computer simulation studies, Battezzati and co-workers1 used Steele’s (classical) method2 of potential expansion to find equilibrium configurations of single molecules of various small hydrocarbons, including benzene, interacting with a graphite surface. They found that the center of the benzene ring sits over a carbon atom in the graphite substrate in the minimum energy configuration. Vidal-Madjar and Bekassy-Molnar used an anisotropic model to study benzene/graphite interactions and confirmed the structure predicted previously.1 Bondi and co-workers4 used pair potentials to calculate the equilibrium structure for a benzene monolayer on graphite. Here, it was found that benzene adsorbed on the basal plane of graphite forms an incommensurate monolayer, i.e., with no correlation between the location of the adsorbing molecule and the underlying surface lattice, but the benzene molecules were found to lie flat on the surface. Sordo and co-workers5 employed the semiempirical method of Fraga6 for several hydrocarbon single molecules adsorbing onto graphite. For benzene, they found an equilibrium structure in which the benzene carbon atoms directly overlapped the graphite carbon atoms with a tilt of about 10° of the benzene plane relative to the surface, in contrast to previous studies. They reported another (local) minimum in which the benzene plane is tilted 70° from parallel with respect to the surface. These two minima were reported to be similar in energy, but the relation of the benzene carbon atoms to the underlying lattice in the second minimum was not clarified. In a study employing molecular dynamics simulations also utilizing the method of Steele,2 Hentschke and Schu¨rmann7 reported that, for a monolayer of benzene X
Abstract published in Advance ACS Abstracts, April 15, 1996.
(1) Battezzati, L.; Pisani, C.; Ricca, F. J. Chem. Soc., Faraday Trans. 2 1975, 71, 1629. (2) Steele, W. A. Surf. Sci. 1973, 36, 317. (3) Vidal-Madjar, C.; Bekassy-Molnar, E. J. Phys. Chem. 1984, 88, 232. (4) Bondi, C.; Baglioni, P.; Taddei, G. Chem. Phys. 1985, 96, 277. (5) Sordo, T. L.; Sordo, J. A.; Flo´rez, R. J. Comput. Chem. 1990, 11 (3), 291. (6) Fraga, S. J. Comput. Chem. 1982, 3, 329. (7) Hentschke, R.; Schu¨rmann, B. L. Surf. Sci. 1992, 262, 180.
on graphite, the molecules at the surface are oriented parallel to the surface, but a significant number were also found in a nearly perpendicular orientation. In experimental studies, Stockmayer and Monkenbusch8 used neutron scattering to characterize a monolayer of benzene on microcrystalline graphite and found a commensurate (x7 × x7) structure. Gameson and Rayment9 confirmed the formation of a commensurate layer using X-ray diffraction. They observed a (x7 × x7) structure for benzene on graphite at somewhat less (0.7) than monolayer coverage and reported a melting temperature of 130 K. Bardi and co-workers10 used LEED (low-energy electron diffraction) to study a benzene monolayer on graphite at low temperatures. They reported finding only one ordered structure and that only a (x7 × x)R19° model allowed correct fitting of the observed intermolecular H-H distances, if molecules were assumed to lie flat on the surface. The actual adsorption site, however, could not be obtained from the data. Tabony and co-workers11 conducted NMR studies of benzene monoand multilayers on graphite at coverages which ranged from 0.4 to 2.0 and confirmed that, at low temperatures and low coverages, the benzene ring lies flat on the surface, but as either temperature or coverage increased, the orientation moved toward a perpendicular configuration. They further presented a phase diagram depiciting the phase behavior as a function of temperature and coverage of up to 3 layers. At monolayer coverage, the boundary between the liquid and solid phases is not sharp, but rather a broad solid-liquid coexistence region was reported to exist between about 130 and 160 K, with 2-D12 solid below 130 K and 2-D liquid above 160 K. Also, at monolayer coverage near the bulk melting point, the ring plane is not aligned parallel to the surface. Finally, the upper limit for monolayer capacity for a flat orientation with no ring rotation is at slightly greater than monolayer coverage (1.08). In this paper, we present results on the molecular dynamics (MD) simulations of a benzene monolayer adsorbed onto graphite at a series of temperatures and discuss the results in terms of the aforementioned phase diagram.11 In the following paper, we extend the study to include multilayers of benzene on graphite. (8) Stockmayer, R.; Monkenbusch, M. Vib. Surf. 1982, 483. (9) Gameson, I.; Rayment, T. Chem. Phys. Lett. 1986, 123 (3), 150. (10) Bardi, U.; Magnanelli, S.; Rovida, G. Langmuir 1987, 3, 159. (11) Tabony, J.; White, J. W.; Delachaume, J. C.; Coulon, M. Surf. Sci. 1980, 95, L282. (12) Meehan, P.; et al. J. Chem. Soc., Faraday Trans. 1 1980, 76, 2011.
+
+
2496
Langmuir, Vol. 12, No. 10, 1996
Matties and Hentschke
2. Theoretical Model The simulations were performed using the classical molecular dynamics technique with a modified version of the AMBER 3.013 collection of programs, as was used previously in this group.7,14-17 The force field used is given by eq 1, where the first five terms
∑ f (r - r
V)
2 eq)
r
bond
∑
∑
+
fθ(θ - θeq)2 +
valence
fn(1 + cos(nφ - γ)) +
∑ x
i j
( )}
dihedral
2
ij
req [Å]
C-C C-H
1960 1420
1.40 1.08
valence
fθ[kJ/(mol rad2)]
θeq
C-C-H C-C-C
146 356
2π/3 2π/3
dihedrals
fn[kJ/mol]
γ
n
x-C-C-x x-C-C-H
22 8
π π
2 2
where r here lies in the xy plane (i.e., parallel to the surface) at a certain height from the surface and
H(‚‚‚) )
{
1
0
)〉
∆r ∆r , rij, r + 2 2
(6)
(13) Weiner, S. J.; et al. J. Am. Chem. Soc. 1984, 106, 765. (14) Winkler, R. G.; Hentschke, R. J. Chem. Phys. 1993, 99, 5405. (15) Winkler, R. G.; Hentschke, R. J. Chem. Phys. 1994, 100, 3930. (16) Hentschke, R.; Winkler, R. G. J. Chem. Phys. 1993, 99, 5528. (17) Kotelyanskii, M. J.; Hentschke, R. Phys. Rev. E 1994, 49 (1), 910.
∆r ∆r e rij e r + 2 2 otherwise if r -
(7)
The diffusion coefficient for linear motion parallel to the surface is
Dxy )
d lim 〈[rxy(t + τ) - rxy(t)]2〉 τf∞ 4 dτ 1
(8)
Perpendicular to the surface, we calculate the apparent diffusivity from
Dzapp )
Here, σ and are the Lennard-Jones parameters of the adsorbate-substrate interactions (the index labeling the different type of adsorbate-substrate pairs is omitted), asa2 is the area of the graphite unit cell, and d is the distance separating the layers in the substrate. Note that V0 only depends on zi, the height above the surface. The corrugation in the surface potential, i.e., the variation in the xy plane at any give height, is then contained in the second term in eq 2 as 12
fr [kJ/(mol Å2)]
i
describe the usual intra- and intermolecular interactions for atoms in the system, in this case, the atoms comprising the benzene molecules. The last term then describes the interaction of any atom (i) of the adsorbate located at some point b ri ) (xi,yi,zi) with the surface in terms of a static potential after the method of Steele.2 Thus, the dynamics of the surface are not included and calculation of the adsorbate/surface interactions is reduced from a calculation going as NaNs to one as Na, where Na is the number of adsorbate atoms and Ns is the number of substrate atoms. The surface potential is given by
bi) ) V0 + Vsurf(r
bond
-
rij
qiqj
i150 K) behavior of the film was studied,7 the magnitude of the first peak height is about 20% larger and the second peak 20% smaller. This difference is attributed to the change in nonbonded parameters, especially the partial changes, which allows more molecules to lie flat on the surface. The change in the configuration of the adsorbed molecules with temperature is not constant. The insert to Figure 1a shows two distinct regimes with the break between each regime occurring at 140 K for all peaks. Another indication that the change in behavior with temperature is not smooth lies in the shift observed in the third peak. From 240 to 180 K, the position of the third peak changes little. There is a small change at 160 and 140 K, with the locations of the peak centers at these two temperatures being similar. Then, at the lowest temperatures, 120 to 60 K, the peak positions shift again. These sudden shifts indicate a sudden change in the structure of the film for temperatures which coincide with the experimentally observed solid-liquid coexistence. The orientation of the molecules in the layers at each of the temperatures is depicted in Figure 2 in terms of the tilt angle distribution, f(θ), where θ is the angle that the
Matties and Hentschke
benzene ring normal makes with the surface normal. In the first peak of Figure 1 (first layer), the orientation of molecules at the surface becomes more flat as temperature decreases. Again, the amount of change in f(θ) at θ ≈ 0 with temperature is not constant. The change in the fraction is relatively small and constant at the higher temperatures, 240 to 180 K, but at 160 K, the change is greater and increases as the temperature is lowered further, thus showing a large change in orientation of the flat layer at the surface between 180 and 160 K. The second peak of Figure 1 (also in the first layer) is the benzene molecules adsorbed with the ring perpendicular to the surface. If the flat, (x7 × x7)R19° structure is the perfect two-dimensional crystal, then the molecules which are adsorbed in a perpendicular orientation represent defects in the crystal. From 240 to 120 K, the orientation is very similar with a rise beginning at about 25°, a maximum at about 57°, and then a decline from 60 to 90°. However, a sudden change in orientation is observed between 120 and 100 K. At 100 K and lower temperatures, the rise begins at about 40°, then reaches a maximum, and remains relatively constant from 60 to 90°. Thus, the average angle that the benzene ring normal makes with the surface normal is about 57° until the temperature drops to 100 K, where the orientation shifts more toward perpendicular, indicating some sudden change in the configuration of the first layer which makes it less able to accommodate the molecules which are more strongly tilted. In the third peak (second layer) of Figure 1, the orientation is also basically flat with the fraction sitting flat, increasing as the temperature decreases. However, it is more difficult to discern any sudden changes in orientation due to the small number of molecules in the second layer. These results confirm the statement in ref 11 that molecules at the surface tilt away from it as temperature increases. The two-dimensional pair correlation function, g2(r), for molecules adsorbed flat on the surface (the first peak of Figure 1) is shown in Figure 3. The curve for the perfect two-dimensional (x7 × x7)R19° crystal appears at the top of each plot in Figure 3. The curve for the perfect structure was obtained from simulation of the monolayer at 0.5 K. The first peak is that of the nearest neighbors (N), the second and third represent the next-nearest and third-nearest neighbors (NN), and the fourth and fifth, the fourth- and fifth-nearest neighbors (NNN). At the higher temperatures, five broad peaks appear which are shifted relative to the peaks observed for the perfect twodimensional crystal with the NN and NNN each blurred into single peaks showing some loss of short-range order. Over the range of 240 to 180 K, the peak heights increase slightly with no shifting of position, including the first peak. At 160 K, shoulders begin to appear in the second (from 10 to 16 Å) and third (from 16 to 21 Å) peaks. At 140 K, the shoulders begin to resolve into separate peaks. As temperature continues to decrease, the heights of the peaks increase, indicating that more molecules are found at those distances and the broad peaks continue to sharpen into the separate peaks observed for the perfect 2-D crystal. Again, as seen in previous figures, the rate of increase in the peak intensities is not constant with temperature. The rate of increase is greater at the lower temperatures than at the higher. At 60 K, the peak centers are shifted slightly versus the perfect crystal, indicating an adsorbate unit cell which is slightly expanded in comparison to the (x7 × x7)R19° structure. A plot of the constant volume heat capacity, ∆Etotal/∆T, is shown in Figure 4. The heat capacity is relatively constant from 240 to 160 K, but then increases suddenly
+
MD Simulation of Benzene on Graphite
Figure 3. Two-dimensional pair correlation function, g2(r), vs r in angstroms. The curve at the top of each plot is that for the perfect two dimensional crystal: (a) upper temperature range and (b) lower temperature range. Temperatures at right are in kelvin.
Figure 4. Constant volume heat capacity vs temperature in kelvin expressed as [E(Ti) - E(Ti-1)]/[Ti - Ti-1], where i indicates the ith temperature.
from 140 to 120 K with a sudden decrease again between 120 and 100 K. B. Dynamic Behavior. The temperature dependence of the center of mass diffusion coefficient parallel to the surface (Dxy) and the apparent diffusivity perpendicular (Dz) to the surface (with τ ) 50 ps) for all molecules in the system are shown in Figure 5a,b. The diffusion coefficient parallel to the surface is relatively constant from 240 to 180 K. Beyond 180 K, it decreases with no great jumps as the temperature is lowered through the transition region. However, for the apparent diffusivity perpendicular to the surface the transition is reflected in a more is relatively constant pronounced behavior. Again, Dapp z from 240 to 180 K but then exhibits a pronounced peak at 160 K, after which it decreases nearly linearly. Note that if the molecules were strictly confined to a layer parallel to the surface, then 〈∆rτ2〉1/2 would be bounded by the layer thickness and the diffusion coefficient Dz would be zero. Here, the diffusion in the z-direction is very slow. Also, it would be difficult to model the exchange with the gas phase because gas phase/surface collisions are very rare events occurring on the time scale of the simulation
+
Langmuir, Vol. 12, No. 10, 1996 2499
Figure 5. Diffusion constant Dxy and apparent diffusivity Dapp (with τ ) 50 ps) temperature (K) for motion resolved (a) z parallel (Dxy) and (b) perpendicular (Dapp z ) to the surface.
Figure 6. Linear velocity orientational autocorrelation function, Z(t), for the benzene molecule center of mass vs correlation time, τ (ps).
for the experimentally relevant gas phase pressures,7 and, thus, we do not observe the development of a plateau for 〈∆rτ2〉1/2. Nevertheless, because the short time behavior of 〈∆rτ2〉1/2 shows a marked dependence on temperature, we choose to define the above apparent diffusivity. It is worth noting that the maximum is not due to any artifact in the calculation of Dapp z , e.g., arising from an isolated molecule desorbing from the surface. The orientational autocorrelation function for the benzene molecule center of mass, shown in Figure 6, provides most striking indications of the transition from liquid to solid phase in the monolayer. At the highest temperatures, 240 to 180 K, the OVAF decays from one to zero monotonically without crossing zero; i.e., the OVAF is positive for all τ. Such behavior in the OVAF is typical of weakly interacting liquids for which the direction of motion is quickly randomized through thermal motion and for which no significant solvent cage formation occurs. At 160 K, the decay is similar to that at the higher temperatures, but a very slight zero crossing occurs at τ ≈ 1.7 ps. Then, at 140 K, the zero crossing becomes more pronounced. From 120 K, the negative lobe increases as the temperature is lowered until at 60 K it is quite
+
2500
+
Langmuir, Vol. 12, No. 10, 1996
significant. Over the range of 140 to 60 K, the value of τ at the minimum of the crossing shifts to shorter times to an ultimate value of 0.4 ps at 60 K. Also at 60 K, a second negative minimum is observed at τ ≈ 1.5 ps. Negative values in the OVAF indicate an anticorrelation in the direction of motion, i.e., a reversal in the direction of motion, which is characteristic of the molecule found within some type of ordered, constraining environment. The (negative) magnitude of the minimum is characteristic of the apparent strength of the constraining cage, where the deeper the minimum, the stronger the cage and the more often on average the direction of motion of any molecule is reversed. The value of τ at the minimum is characteristic of the apparent size of the cage, where a shorter time indicates a smaller cage. While simultaneously considering the two-dimensional pair correlation functions in Figure 3, we interpret this progression of the OVAF to indicate a purely liquid structure above 180 K with no solvent cage formation, the onset of some sort of ordered structure of effective cage at 160 K, the region of 160 to 120 K being some mixed or intermediate region, and solid structure of increasing order from 100 to 60 K. Since the center first peak in g2(r) (Figure 3) does not shift in position with a change in temperature, the nearest neighbor distance remains constant. As it is, the nearest neighbors which would form any sort of cage, the shift to shorter times of the negative minimum in the OVAF cannot be due simply to the constraining cage formed being reduced in size. It must then be due to the cage forming molecules being more strongly adsorbed and less willing to yield when the constrained molecule collides with a cage molecule. Thus, the molecules are more strongly adsorbed, and the cage becomes less elastic as the temperature is lowered. This interpretation agrees well with the phase diagram11 at monolayer coverage where the system is liquid above 160 K, exhibits mixed behavior between 160 and 130 K, and then forms a solid below 130 K. Finally, the orientational autocorrelation function for the angular velocity of the benzene ring spinning is presented in Figure 7. From 240 to 140 K, the OSAF decays quickly to zero with the rate of decay decreasing somewhat as the temperature decreases, indicating random rotational diffusion with the time required to randomize increasing with decreasing temperature. Then, the difference from 140 to 120 K is rather large with the curves at 120 and 100 K being comparable. Finally, there is a significant decrease in the rate of decay of the OSAF at 80 and 60 K with a significant correlation (about 0.4) remaining at long times (≈30 ps). This behavior is not an artifact of the OSAF calculation but is the contribution
Matties and Hentschke
Figure 7. Angular velocity (spinning) orientational autocorrelation function, C(τ), for the benzene molecule center of mass vs correlation time, τ (ps).
due to the very little rotational motion of the rings adsorbed on the surface with perpendicular orientation at the lowest temperatures, which accounts for the apparent high correlation at long times. Thus, the ring spinning seems to be somewhat discontinuous between 140 and 120 K and quite discontinuous between 100 and 80 K. In the experimental phase diagram,11 one line, labeled “IV”, is supposed to indicate a two-dimensional solid without rotational motion. However, it is not clear whether this line is a boundary separating regime. For example, no rotational motion is possible below line IV. This line IV lies at slightly greater than monolayer coverage, at a coverage of 1.08, so that at monolayer coverage and below about 130 K, no rotation should be possible. The OSAF’s calculated agree fairly well, where, at 140 K and above, the system exists as a liquid or mixed liquid/solid with rapid randomization of the direction of benzene ring spinning and at 120 K and below, the rotational motion becomes frozen in. 5. Conclusions We have presented the results of MD simulations for a monolayer of benzene adsorbed onto graphite at a series of temperatures. This study confirms earlier reports that the majority of the rings in contact with the surface assume an orientation parallel to the surface with the orientation becoming increasingly flat as the temperature decreases. The data presented reveal pronounced changes in properties between 180 to 160 K and 140 to 120 K in agreement with the published phase diagram. The most compelling example is provided by the orientational autocorrelation function for the linear velocity of the benzene ring center of mass, for which a clear progression from liquid to solid structure with an intermediate, mixed region is observed. LA950766Q
