J. Phys. Chem. 1993,97, 1639-1648
1639
Molecular Dynamics with the AM1 Potential: Reactions on Diamond Surfaces X. C. Zhao, C. S. Carmer, B. WeinerJ and M.Frenklach' Department of Materials Science and Engineering, The Pennsylvania State University, University Park, Pennsylvania 16802 Received: August 8, 1992
Classical molecular dynamics simulations of diamond surface reactions were performed using quantum mechanically derived forces. A semiempirical A M 1 potential from the MOPAC quantum-chemistry program package was employed in the calculations. The evolving molecular geometry, energetics, bond orders, and vibrational frequencies were monitored. The conformation and bond breaking of a surface step, examples of reactions of gaseous species with the diamond ( 111) surface, and reactions of free gaseous analogs were examined. One of the main findings is the existence of a trapped, physisorption state for an acetylene molecule colliding with a diamond surface at a temperature typical of diamond CVD. The hypothesis of chemical similarity between the rate constants for gas-surface and analogous gas-phase reactions is addressed. The results obtained indicate that although a considerable similarity indeed exists, significant quantitative differences are induced by the confinement of the surface.
1. Introduction
Methods of molecular dynamics (MD) are being steadily developed for applications to computational studies of chemical reactions occurring on solid surfaces.l-10 In this approach, Newton's equations of motions, for both gas-phase and surface atoms, are integrated numerically. The principal difference between various methods lies in the way interatomic forces are calculated. In most previous studies, the calculation of forces has been based on analyticalpotential-energyfunctions,developed by fitting either experimental data4Jl-15 or results of ab initio quantumchemical calculations. While offering computational speed and thus being suitable for building up reaction statistics for well understood systems, analytical potentials are of limited applicability in exploratory investigations where one is interested in testing unknown reactions, rearrangements, and transition states. This is especially so when the chemical reaction may involve recrossing of potential energy surfaces or when dealing with multielementsystems, for which the parameterizationof analytical potentials may become in itself overwhelming. An ideal approach in such situations would be to calculate the forces from an ab initio potential, which utilizesthe exact electronic Hamiltonian that contains a description of the potential energy due to the interaction of electrons with nuclei and other electrons and the kinetic energy of the electrons at fixed nuclear positions. Unfortunately, such calculationsusually require excessively large amounts of computational effort. Several attempts of using ab initio potentials in molecular-dynamics calculations have been limited to tight binding minimal basis sets,19 small molecular systems,20or systems that can be adequately described by the combination of pseudopotentials and plane wave bases.21 Semiempirical quantum potentialsoffer a compromisebetween high-level quantum description and computational speed. These methods use model Hamiltonians that depend on parameters which are fitted to reproduce experimentally observed quantities and thus the computational efforts involved are much less than with ab initio or density-functional methods. A large number of static semiempirical calculations have been carried out over the last 30 years, and much experience has been gained in their reliability and accura~y,~2-~3 especially for carbon-containing compounds which are of interest to the present study.
' Department of Physics, The Pennsylvania State University, DuBois, PA
15801.
0022-3654/93/2097-1639$04.00/0
The use of semiempirical potentials for quantum-mechanical molecular dynamics simulations has been recently advocated and demonstrated on a small size molecule.24 The objective here is to test the feasibility of using such an approach for gas-surface chemical reactions. The details of the method are described in the next section. Then, several examples are presented examining reactions taking place on a diamond surface. Identification of the chemical processes responsible for diamond growth under chemical vapor deposition is of growing interest at present.2sOne of the hypotheses advanced is the principal of chemical postulating that the chemical reactivity of solid carbonaceous materials (Le., diamond, graphite, soot, amorphous carbon) is localized on the carbon atom sites in a manner similar to that of the corresponding gaseous species, or in other words, that the kinetics of analogous elementary chemical reactions on a per site basis is the same for all the forms of carbon. It is not clear, however, to what degree such analogy would hold or how to quantitatively transfer the kinetic information (i.e., rate coefficients) from one system to an0ther.2~928 We pay particular attention to these questions in the examples presented below. 2. Computational Method
Quantum mechanically derived forces were used to calculate classical atomic trajectories for reactions on hydrogenated diamond (1 11) surfaces. Subroutines from the MOPAC semiempirical quantum-chemistry program package29were interfaced with a numerical integration algorithm and various utility subroutines to create a molecular dynamics code. In addition to providing a well-tested many-body potential, the MOPAC package suppliescopious information on the electronicstructureof a system being studied. When used in conjunction with molecular dynamics, the time evolution of such relevant properties as, for example, bond order, Mulliken charge population, and electronic spin state can be followed. At each time step in a simulation, Schrodinger's equation was solved within the chosen semiempirical level of approximation for the lowest energy state of the electronic structure associated with the instantaneous nuclear configuration. The potential energy responsiblefor nuclear motion thus included contributions from semiempirically determined electron-electron, electronnuclear, and nuclear-nuclear interactions and the electron kinetic energy. Forces on nuclei due to the instantaneous electronicstate were calculated from first derivativesof the potential energy with respect to nuclear position and converted to atomic accelerations (6
1993 Americaq Chemical Society
1640 The Journal of Physical Chemistry, Vol. 97, No. 8, 1993
Zhao et al.
by Newton's second law. Equations of motion were then integrated numerically, producing atomic coordinates and momenta for the next time step. This process was repeated, solving for the potential energy and interatomic forces of the latest configuration, until the desired simulation time was covered. A. Potential. The level of approximation used to solve Schrodinger's equation was the unrestricted Hartree-Fock selfconsistent field (SCF) method using the AM1 Hamiltonian with Slater-typeorbitals and minimal basis sets. UnrestrictedHartreeFock theory was employed in our study as it is more capable of describingthe formation and breaking of bonds than the closedshell restricted Hartree-Fock approximation. The AM1 performance for the energetics and geometry of Figure 1. Diagram of diamond (1 11) surface model. Atoms marked in hydrocarbon species has been extensively tested and well docsolid black are the principal 16 carbon atoms of the model surface and umented in severalsources.22J3 Typically, the heats of formation those striped are the image atoms, used for the simulation of the periodic computed with NDDO methods are within about 6 kcallmol for boundaryconditions. Theatomswithprimed anddouble primed numbers hydrocarbon molecules and within about 12 kcal/mol for are the image atoms of the corresponding nonprimed numbered atoms. hydrocarbon radicals. Much less tested are potential-energy Atoms (17), (18), (19), and (50) are hydrogen atoms saturating the top surface layer. The rectifying H atoms are omitted. reaction barriers, but from few studies that are a ~ a i l a b l eit, ~ ~ ~ carbon ~ seems that the accuracy for transition states might be comparable to that of radicals. AM1 34 and other, similar but differently parameterized, semiempirical Hamiltonian~-MINDO/3,3S,~~ MND0,30~37~38 PM3,3l ASED-M0.39 and SLAB-MIND@-have been used extensively in static quantum mechanicalcalculationselucidating the geometry of diamond surface reconstructions and energetics of diamond surface reactions. The per-atom cohesive energy of diamond calculated4' with a recent vtrsion of MOPAC-AM1 is 178 kcal/mol (7.72 eV), which is very close to ab initio results of 177-183 kcal/mol (7.69-7.94 eV), obtained with pseudapotential local-density-functionalmethods."2 For comparison, the corresponding experimental value43 is 170 kcal/mol (7.37 eV); the value obtained from an empirical potential fitted to ab initio, HartrcGFock self-consistent-field calculations on C2 and CS 4 X clusters with contracted Gaussian double (for C atoms) and triple Figure 2. Top view of diamond (1 11) model surface shown in Figure 1. (for H atoms) zeta basis sets,I7 is 154 kcal/mol (6.69 eV), and The frame marks the boundary of the molecular-dynamics cell. the results of empirical potentials are 173 kcal/mol (7.5 eV)I2 and 169 kcal/mol (7.32 eV).14 by four layers containingfour carbon atoms each. Figure 1 depicts The criterion set for obtaining a self-consistent field was that a ball-and-stick diagram of the assumed model, where the atoms the potential energy determined by twosuccessiveSCF iterations marked in solid black represent the principal 16 carbon atoms of differ by less than 10-9 eV. Preliminary microcanonical simuthe model surface. Those atoms marked in stripes will be referred lations showed that this high degree of precision was necessary to as the image atoms; they are used for the simulation of the to maintain sufficient accuracy in the calculation of trajectories. periodic boundary conditions asdiscussed below. The numbering Interatomic forces are first-order response properties of' the of atoms is consistent throughout the manuscript and is used to electronic state, and as such their precision goes as the square identify the discussed chemical transformations. root of the precision of the potential. As these forces directly The carbon atoms of the model surface were initially placed couple kinetic energy and potential energy, the error is passed on in the ideal diamond structure. During a simulation the bottom to the calculation of trajectories. layer was held fixed at ideal lattice positions to simulate the bulk. Determination of the electronic etatcat each time step was the Hydrogen atoms were used to hydrogenate the top carbon layer. most computationally intensive portion of the dynamic simulaThese atoms are shown as smaller white circles in Figure 1. tions. To minimize the time required to obtain a self-consistent Periodic boundary conditions were imposed across planes perfield at subsequent time steps, the starting point of each SCF pendicular to the (1 11) surface to define a molecular dynamics calculation,exceptfor thefirst, wasraken tobe theset ofmolecular (MD) cell. As the MOPAC program is set up to perform orbital (MO) coefficientswhich dtscribod the electronicstructure calculations on finite clusters of atoms, the usual method of at the previous time step. In contrast to simulations performed incorporating periodic boundary conditions by the minimum image with purely analytical functions for potential energy, no strict conventionMwas difficult to implement. An alternative scheme relationship was found between the number of atoms and the was devised to provide pseudoboundary conditions whereby total CPU time required for a simulation using the SCF method. additional atoms, referred to as 'image" atoms and denoted by While increasing the number of atoms did increased total CPU primed numbers in Figures 1 and 2, were included in the time, other factors were found ta influence it as well. Higher calculation of potential energy. The positions of these image temperatures cause atomic mobility to increase,resulting in larger atoms were determinedby translating, in.the appropriatediraction, chanees in atomic positions between t h e steps. As the starting coordinates of atoms within the MD cell by a distance equal to point qf the SCF calculations was based on MO's of the previous the length of the cell and were updated at each time step so that time step, more time was n d e d for convergence at higher the motion of image atoms mimicked the motion of atoms within temperatures. In addition, SCF.wnvcrgence was slowed conthe MD cell. Hydrogen atoms were used to saturate carbon siderably at points during a simulatios where bond formation or atoms at the bottom and edges of the cluster in order to maintain bond breaking processes were occurring. the proper electronic state. The position of these rectifying H atoms was determined at each time step so as to best maintain B. Model Surface. The diamond (1 11) surface was modeled
Molecular Dynamics with AM1 Potential sp3 geometry. A total of 64 atoms were required to implement this periodic boundary scheme-( 16 C + 4 H) within the MD cell 14C image atoms 30 H rectifyingatoms. This technique was used solely to ensure a smooth and periodic potential across the boundary of the MD cell and so the image atoms and the rectifying hydrogen atoms were not included in the calculation of kinetic energy and temperature of a dynamic simulation. Note that while this method of invoking pseudoperiodic boundary conditionsresults in a finite model surface, it is equivalent to the minimum image convention in conjunctionwith a potential cutoff just beyond the first nearest neighbor distance. C. Molecular Dynamics. Beeman’s third-order predictor algorithm4$was used to integrate Newton’s equations of motion. This algorithm has been found to be particularly stable and efficient in calculating trajectories of particles interacting via steep welled potentials. Moreover, this integration scheme requires only one evaluation of atomic accelerations per time step, an important consideration when the potential calculation is quite time-consuming. Beeman’s algorithm is not self-starting, so the velocity form of Verlet’s algorithm46was used for the first integration. A time step of 0.5or 1.O fs was typically used in the simulations. It was found that low frequency (-200 cm-I) vibrations/rotations and some changes in total linear momentum were introduced into the motion of the model cluster as an artifact of the above described periodic boundary conditions and fixing the bottom layer of atoms. These nonphysical motions were removed from the system at each time step by first calculating the total linear momentum and angular velocity of the system as a whole and then subtracting these quantities from the motion of each atom. Simulations presented here were carried out under two types of conditions: adiabatic and isothermal. By adiabatic we mean that atomic trajectores followed undamped Newtonian dynamics. Isothermal runs were achieved by rescaling velocities at each time step to maintain the desired temperature. The reasons for maintaining isothermal conditions were 3-fold. First, reactions at the surftce release or absorb energy resulting in large changes in the temperature of the simulation due to the small thermal mass of the model surface; second, due to the small cluster size, the temperature fluctuationsare unrealistically large, and reaction events, the subject of the present work, exhibit even larger fluctuations because of their exponential dependence on temperature; third, the rescaling of the temperature reduces the computer time required to reach thermal equilibrium. In isothermal simulations, the time required for the system to reach equilibrium was estimated by the magnitude of the velocity scale factor needed to maintain a constant temperature. D. Determination of Vibrational Frequencies. Vibrational spectra were computed by summing over all atoms the square of the Fourier transform of orthogonal components of atomic velocities as in the following equation,
+
The Journal of Physical Chemistry, Vol. 97,No.8, 1993 1641
TABLE I: Vibrational Frequencies of C2H2 in Units of cm-1 AM 1/MOPAC potential
+
where E(w) is the energy of the mode with frequency w , n the number of atoms, mi the mass of atom i, T the total time of the simulation, and uij(r) a Cartesian velocity component of atom i at time 1. In the limit r m, the second summation in eq 1 has been shown to be equivalent to the Fourier power spectrum of the velocity autocorrelation f~nction.~’The factor mi/2 in the first summation converts the power spectrum to an energy spectrum and is necessary for scaling the relative contributions to each mode made by atom of different masses. Under typical simulation conditions,thespectral resolutionattained was about 1&20cm-l.
-
3. R d t s and Diwussion A. VibrrtioarrlSpectraof Model Compounds. Before discussing reactions on the diamond surface, it is appropriate to describe the
Eo vibrational mode
ab initio
(kcal/mol) harmonic
HFu
MD ~
C-H w2 C-C w3 C-H w4 C-H w5 C-H WI
stretch (us+) stretch (us+) stretch (u,+) bend (rg) bend (xu)
5.05 3.09 4.92 1.15 1.32
3528 2185 3442 804 929
3449 2171 3317 800 923
ex&
~~~
3665 2210 3556 814 868
3491 2011 3415 624 141
Restricted Hartree-Fock c a l c ~ l a t i o n swith ~ ~ basis set (8s6p3d2f)/ [6s3pld]. Experimental data of Strey and Mills.49
results of simulations that can be directly compared to experimental findings. In this section, the vibrationalspectra of groundstate acetyleneand a hydrogenated diamond surface are presented in order to evaluate the present implementation of molecular dynamics with the AM1 potential. Simulations of acetylene were prepared by choosing initial velocities such that each vibrational mode was populated to a specified energy. Normal mode frequencies and eigenvectors were first determined by normal mode analysis of the static molecule by AM1-MOPAC. Table I lists the normal mode frequencies calculated under the harmonic approximation and the corresponding ground-state energies. Taking the minimum potential energy configuration as a starting point, initial velocities were determined such that the total kinetic energy of each mode, j , was equal to the ground-state energy, E = huj/2. To illustrate, Figure 3 (top panel) shows the vibrational energy spectrum of an MD simulation of C2H2 in which a single mode, the symmetric C-C stretch, was excited to the energy level equivalent to the quantum ground state. Five such simulations were carried out, each simulation having a different vibrational mode populated. A time step of 0.5 fs and a total simulation time of 1000 fs were employed for each simulationwhich was carried out adiabatically. Figure 3 (middle panel) is a superpositionof the spectra produced by these five simulations. The energy residing in each mode is nearly equal to the ground-energy initially placed in that mode (the lowest two frequencies,correspondingto the symmetric and antisymmetric C-H bending modes, are doubly degenerate). Because no energy transfer was possible between modes (only one mode existed in each simulation),this superpositionof spectra can be viewed as representing the initial state of C2H2in which all modes are in vibrational ground states. Another simulation, carried out adiabatically for 2000 fs with time step 0.5 fs, was performed in which all vibrational modes of CzHz were simultaneously excited to the ground state. The energy spectrum of this simulation is shown in Figure 3 (bottom panel). The peak height of the three lowest frequency modes remained virtually unchanged compared to Figure 3 (middle panel), indicating little energy transfer involving these modes. However, the two highest frequency modes, corresponding to symmetric and antisymmetric C-H stretches, showed significant intermodal transfer of energy. Approximately 3.5 kcal/mol of vibrationalenergy went from the antisymmetric to the symmetric mode over the 2OOO-fsduration of this simulation. Fourier spectra from simulations over a longer time scales indicate that energy readily moves back and forth between these two modes as allowed in classical trajectory simulations. We speculate that because the frequencies of the C-H stretching are so close together there may be a resonance process involved. In addition, the groundstate energy initially supplied to these modes is apparently large enough that the atomic motion traverses the anharmonic regions of the potential hypersurface, resulting in anharmonic coupling effects. Table I lists the frequencies observed from this M D simulation. It can be seen that the dynamically determined frequenciesfor the two C-H stretch modes are significantly lower than those calculated for the static molecule under the harmonic approximation, while for the other vibrational modes the frequencies determined by both methods are in close agreement.
Zhao et al.
1642 The Journal of Physical Chemistry, Vol. 97, No. 8, 1993
,
6 1
5 1
I
(A I
L
I
0
500
1000
1500
2000
Wavenumber (cm") Figure 4. Phonon density of states of diamond at 300 K: dotted line, experimental data fitted by Weber;'O solid line, from MD simulation of 16 carbon atom surface model shown in Figure 1, resolution 10 cm-I.
~
I
7 -
5 -
3 -
Wavenumber (cm") Figure 3. Energy spectra of C2Hz: top panel, symmetric C-C stretch mode excited to the quantum ground state; middle panel, superposition of energy spectra from five simulations, each having a different normal mode excited to the corresponding quantum ground state; bottom panel, all normal modes initially excited to the corresponding quantum ground state. The resolution of the spectra is 10 cm-I.
This result further indicates the existence of anharmonic coupling between the two high-frequency modes. Normal mode frequencies for acetylene determined by ab initio Hartrec-Fock calculations with an extended basis set4*are also listed in Table I, along with experimentally measured values.49 The dynamically determined frequencies for the C-H stretching modes agree with the experimental values to approximately 1% while the C-C stretch and C-H bending modes are overpredicted by 8 and 2 5 1 , respectively. The frequencies produced by the semiempirical AM1 Hamiltonian compare very favorably with those of the ab initio calculation, giving more accurate values for all modes except the antisymmetric C-H bending mode. A simulation of the hydrogenated diamond (1 11) surface was carried out using the periodic boundary conditions described in section 2B. Initial velocities were chosen at random and initial atomic positions corresponded to the ideal diamond structure. Velocities were rescaled at every time step to maintain the temperature a t 300 K. A total simulation time of 2500 fs was covered with a time step of 1 fs. The velocities of the carbon atoms within the MD cell were used to make a vibrational energy spectrum for the final 2000 fs of the simulation. Figure 4 compares this spectrum with an experimentally fitted photon density of
states for diamond at 300 K.5O The two largest peaks in the MD spectrum, located a t 1190 and 1280 cm-I, lie well within the broad peak from 1050 to 1330 cm-I in the experimental density of states curve which is attributed to both transverse and longitudinal optical normal modes. The small broad bump in the MD spectrum between 600 and 750cm-I corresponds to transverse acousticmodes in the density of states which extends from roughly 500 to 800 cm-I. An additional low intensity peak in the MD spectrumat 3050cm-I (notshown) isattributed toC-H stretching modes on the diamond surface. Clearly, there are modes of vibration absent from the MD simulation, in particular the longitudinal acoustic modes centered at 960 cm-1, which are present in a full diamond crystal. We attribute this lack of modes to the small dimensions of our model surface which does not samplethe full Brillouin zone. The AM 1potential does, however, seem capable of accurately describing the dynamic behavior of diamond a t 300 K, based on the position and relative intensities of the peaks that do occur. An additional simulation of the hydrogenated diamond surface a t 300 K was carried out to test what effect rescaling velocities has on the vibrational spectrum. In this test case, the top two layers of atoms were allowed to follow undamped Newtonian dynamics while the third layer was maintained as a heat bath at 300 K by the velocity rescaling. The resulting vibrational spectrum, generated fromvelocitiesof the top two layers of atoms covering 2500 fs of simulation, was essentially identical to the spectrum produced under purely isothermal conditions shown in Figure 4. This agreement between the spectra obtained for the damped and undamped cases is an indication that the velocity rescaling does not substantially influence the surface dynamics. This is reasonable because the change in velocities upon rescaling was typically less than 1%once the surface reached an equilibrated state. It is pertinent to mention that the MD calculations exhibited a negative frequency shift in the diamond spectrum with increasing temperature. Such a shift for the Raman-active mode was observed e~perimentally.~'The specific value of the shift frequency we calculate is about 20 cm-I between 500 and 1500 K,which is in close agreement with the results obtained by Wang et al.,S2 who performed molecular dynamics calculations using an empirical tight-binding Hamiltonian. The negative frequency shift phenomenon is another manifestation of the anharmonic phonon effects and the present numerical results demonstrate that AM1 potential is capable of reproducing them. B. Physisorption of CzHz. Chemical reactions of molecules on a surface of a solid can occur through two principal mechanisms: by a direct reaction or via a trapped, precursor
Molecular Dynamics with AM1 Potential
The Journal of Physical Chemistry, Vol. 97, No. 8,1993 1643
Figure 5. Initial configuration used in the run of a free C2H2 molecule colliding with a model surface obtained from the one shown in Figure 1 by replacing H(17)with a CH2 group. Only the top layer of the model surface is shown for clarity.
stateaS3In the former case, a reaction event either takes place immediately upon the collision of an incoming gaseous reactant with a surface site or, otherwise, the gaseous reactant scatters away from the surface. In the precursor-state mechanism, an adsorbate molecule is trapped in a physisorbed state, which results in the exchange of the kinetic energy between the incoming molecule and the surface. The molecule remains in the trapped physisorption state, bouncing and migrating on the surface, until the energy and position of the adsorbate molecule either allow it to become chemically bonded to the surface or results in escape from the physisorption state and thus from the surface. During the course of the present study, computer simulations were performed with an acetylenemolecule colliding with a model diamond surface. In these runs, the model surface was that described in section 2B with a CH2 radical group bonded to one of the carbon atoms of the top surface layer, as shown in Figure 5 . This model surface and an acetylene molecule were each equilibratedat a temperature of 1200K, typical of diamond CVD, for about 500 fs. The acetylene molecule was initially positioned on one side of the CH2 step roughly 3.5 A above the surface, defined by the C(13)-C( 14)-C( 15)-C(16) plane, with a velocity of 9.88 X lo4 cm/s, the average thermal velocity at 1200 K, directed toward the surface. Three cases were run with slightly differing initial configurations. In two cases, acetylene rebounded from the surface. However, in a third case carried out adiabatically, acetylene was trapped on the surface. Figure 6 presents the atomic distances of the C2H2 molecule and its center of mass from the top surface layer (Figure 6, top panel) and from the carbon atom of the CH2 step (bottom panel) as functions of time. It shows that the acetylenecarbon atoms bounce on the surface at a distance larger than 3.5 A from the surface layer and 2.5 A from theCH2radical step, which are beyond the distances for chemical bonding. Apparently, the C2H2 molecule was trapped in a van der Waals well, its orientationbeing influenced by the surface and step atoms. Unfortunately, the limited size of the present model surface prevented the continuation of the simulation for an extended period of time: at about 420 fs from the start, the C2H2molecule reached theboundaryoftheMDcel1,at which point thesimulation was halted. The latter observation indicates that in addition to the bouncing motion perpendicular to the surface, the acetylene molecule also moves in the plane parallel to the surface, Le., diffuses along the surface. The described results clearly demonstrate the behavior typical of the trapped, physisorption state. The relatively large computational demand on the computer time of the present MD method prohibit us from obtaining the statistics of this behavior; however we can estimate the average lifetime of the physisorbed acetylene bys4 t d = U-leDIRT,where u is the bouncing frequency, D the potential energy of the physisorption state, R the universal
3
f
8
4
3 2
I t
t
0
100
200
300
400
Time (fs) Figure 6. Atomic distances of CzH2 from the top surface layer top panel) and from C(17)in Figure 5 (bottom panel). The thick solid nes are the corresponding distances for the C2H2 center of mass.
gas constant, and T the absolute temperature. For the example presented above, we obtained u-' = 200 fs (see Figure 6) and D = 8 kcal/mol, typical values for hydrocarbons physisorbed on metal surface^.^^^^^ With these values of u-I and D,the average lifetime of physisorbed acetylene at 1200 K is about 5 ps. Although this is much shorter than the average time between gas-surface collisions (about 0.1-1 ps under typical diamond deposition conditions), it is much longer than the characteristic time of a single collision. This means that while trapped in such a mobile physisorbed state, the acetylene molecule has many more chances to find a low-potential-energy-barrierreaction path leading to chemisorption than by a direct-reaction mechanism. In other words, the existence of a mobile precursor state increases5S the effective reaction rate coefficient compared to the directreaction-mechanismtreatment. Also, the motion of the acetylene molecule along the deposition surface, i.e. the surface diffusion, may cause the development of crystallographic growth patterns. Thus, the possibility of a mobile precursor state for acetylene and perhaps for other hydrocarbon species should be considered in future studies of processes responsible for diamond growth. C. Hydrogen Reactions. Reactions of hydrogen atoms on a diamond surface are among the most critical processes governing the deposition of diamond.2s Two reactions, abstraction of a surface H atom to create a surface radical and combination with a surface radical, are of particular importance as they control the fraction of thesurface sites being dehydrogenated,i.e., thenumber density of the surface radical sites.26.27+56 Brenner and coworkerss' recently investigatedthe probabilities of these reactions by molecular dynamics with an empirical potential.I4 Here, we performed several trajectory calculationsusing the quantum AM 1 potential with the objective of exploring the possible differences and similarities in the dynamic behavior of H-atom attacks on
1644 The Journal of Physical Chemistry, Vol. 97, No. 8, 1993
Zhao et al.
I
Timc (fa)
Figure 7. Bond orders for the H atom reactions with the surface model of Figure 5 (top panels) and a comparable free CH3 radical (bottom panels) at three different initial positions of the H atom; the indicated RHxvaluesare projectionsoftheinitial H-Cdistancesonto the horizontal (X-Y) surface plane.
a surface radical and a comparable free CH3 radical. Figure 7 compares the results of three sets of such simulations. The top row of panels in Figure 7 depicts the variation of bond orders obtained for a hydrogen atom attack on the same model surface as was used in the previous example of acetylene (shown in Figure 5 ) , but with a free H atom instead of C2H2. The surface was equilibrated at T = 1200 K for about 1 ps before the free H was sent in. The hydrogen atom was initially placed above the top layer of the surface, at different distances and azimuthal angles from the carbon C( 17) of the CH2 surface step, with a velocity of 5.02 X 105 cm/s (average thermal velocity at 1200 K) directed toward the surface. The bottom row of panels in Figure 7 depicts the variation of bond orders along trajectories similar to those of the top row, but for the H attack on a comparable free CH3 radical. The initial position of the latter was as follows: the C atom of CH3 was placed in the exact location of C( 17) of the surface step, the two H atoms of CH3 in locations of H(21) and H(22) of the surface, and the third H atom of CH3 in direction from C(17) to C(14) at a distance of 1.11 A from C(17) corresponding to the AM1 bond length of C-H. The trajectories were computed adiabatically. Depicted in Figure 7 are bond orders between the incoming H atom and the C atom of the surface step and the C atom of the free CH3, respectively. Bond orders are defined as the sum of the squares of the density matrix elements connecting two atoms. At large separations, as demonstrated in the right panels of Figure 7, no bonds are formed in either of the two cases, the H-Csurface or H-C,,,. At short separations, as on the left panels of Figure 7, the C-H bonds are formed in both the H-Csurrae and H-C,,, cases, yet the behavior is somewhat different. After about 350 fs from the beginning of the H-Csurface run, the newly added H atom broke off and abstracted a nearby surface hydrogen from the surface, H( 19),forming a free hydrogen molecule and a surface radical at C( 15). In the intermediate cases, shown in the middlecolumn panels of Figure 7, a C-H bond is formed for the H attack on the free CH3 but not on the surface. This difference was caused due to the interaction between the incoming H atom with a surface H, H(19). The latter two setsof trajectories clearly indicate the differences between analogous gas-phase and gas-surface reactions. The differences were caused by the interaction between the reacting H atom and neighboring surface sites. To what degree these differences will affect the kinetic rate coefficients cannot be established with the present calculations. However, since in no case the trapping of the free hydrogen atom had occurred, and probably should not be expected due to the poor energy transfer as a result of the low mass of H, trajectory calculations with a properly calibrated empirical potential, as for instance those of Brenner et al.,S7can suffice for determination of the rate data for these reactions.
Figure8. TheUchair"conformationoftheC,H5stepon themodelsurface.
Only the top surface layer is shown for clarity. 14
7
12
6 5 4
3
-
3 B g
2
I
dO
-35
0
35
0
70
Angle (degree)
Figure 9. Data points represented by diamonds are the potential energy calculated by static AM1 for the boat-chair conformation. The solid line is the least-squares fit through these points. The abscissa is the dihedral angle between the C(17)-C(28)-€(27) plane and the vertical plane containing atoms C(17) and C(27) in Figure 8. Negative angles correspondto theboat conformation,positiveanglestochair. Thediagram depicting the oscillating boat-chair motion in time is the result of the isothermal M D simulation at 1200 K.
D. Boat-Chair Conformation. Elucidation of processes that govern the growth of diamond is currently of primary concern to the field of diamond CVD.*j One conceivable (prototype) surface intermediateis thought tobe thestructureshownin Figure 8: a C3H5radical step, C(27)-C(28)-C( 17), on a (1 11) surface. The formation of this species on the otherwise perfect (1 11) plane may serve a nucleus for the growth of the next layer.26 The C3Hs step, however, can take two distinct conformations, each leading to a different diamond polytype, cubic or hexagonal stacking of diam0nd.ss-5~Shown in Figure 8 is the "chair" configuration of the C3H5 step, Le., carbon atoms C(28)-C(17)-C(14)-C(ll)C( 15)-C(27) form the "chair" conformation of a saturated sixmembered ring. Further growth from this configuration leads to cubic stacking of the next diamond layer. Flipping C(28) over the C( 17)-C(27) axis results in the "boat" form of the C3H5 step, which leads to the hexagonal stacking of the next layer. Potential energy calculations under static conditions with the AMI potential resulted in a 2.4 kcallmol difference between the boat and chain conformations of the C ~ H step. S This is shown in Figure 9, where the potential energy of the step is plotted against thedihedral angle between the C( 17)-C(28)4(27)plane and the vertical plane containing atoms C( 17) and C(27). The negative values of the angle correspond to the boat configuration and positive to the chair. The potential energy barrier going from the boat to the chair is calculated to be 0.36 kcal/mol. Superimposed on the static-potential-energy diagram in Figure 9 are the results of an isothermal dynamic simulation at 1200 K using the same quantum-chemical potential. Examination of the latter indicates that the C3Hs step spends most of the time, as expected, in the lower-energy, chair conformation and only
Molecular.Dynamics with AM 1 Potential
The Journal of Physical Chemistry, Vol. 97, No. 8, 1993 1645
occasionally crossing the barrier to the boat form. What is of interest is the irregularity in the chair-boat transformations: the chair-to-boat crossings do not occur at regular time intervals but rather in small uneven bunches. Having performed static AM1 calculations, we can determine the frequency of the chair-boat transformations from the knowledge of the normal-mode vibrational frequencies and applying the formalismof a transition-state theory (TST) -60Using a simple canonical TST, without consideringzero-point energies to be consistent with the MD simulations, the characteristic frequency (the rate coefficient) at 1200 K for the chair-to-boat directionis 3.2 X 101*s-'.Thezerepoint energycorrectionwould increase the estimated rate coefficient by only 28%, too small to affect our following conclusions. h comparison with the MD results shown in Figure 9, the TST result corresponds to the frequency of the individual crossings, not of the bunches, and hence does not describe the much slower periodicityof the actual phenomenon. This discrepancy is a consequence of the interaction of the C3H5step with the phonons of the "solid", which enters into the formalcalculationsthrough the entropy term. Thus, the standard entropy change for the chair-to-boat transformation is 1.98 cal mol-' KI-, as calculated by the classical limit of the partitionfunction formalism,61
rI N-6
m = kB In
W y r
i= 1 N&
(2)
ni= Iw p r wherew values are the normal-mode frequenciesobtained in static AM1 calculations, ke is the Boltzmann constant, and N is the number of atoms; all the cell and image C atoms along with the associated rectifying H atoms were included in this calculation. Essentially thesame result is obtainedushg thequantum partitionfunction equation, AS = kB In qbht/qVchair,where qv is the vibrational partition function, because only low-frequency modes are different for the two conformations. The entropy change of 1.98 cal mol-' K-' along with the energy difference of 2.4 kcal/ mol results in the equilibrium constant of about unity at 1200 K for the chair e boat isomerization reaction: In other words, the C3H5 step is expected to spend equal amounts of time, roughly 300 fs at the conditions stated, in the two conformation states. This describes correctly the series of barrier crossings occurring between 2 and 3.5 ps in the MD simulations depicted in Figure 9. The entropy change for the chair-to-boat transformation can be obtained from dynamic simulationsusing the method suggested by Karplus and Kushik,62
k,
uchair
m=T1nF where u is the determinant of the matrix whose elements are variances and covariances of the internal coordinates, q, i.e., ~y = ((qi- (qi))(q, - (4,))) and the angle brackets indicate time averages. Fifteen cell C atoms (excludingthe bottom four C atoms that were permanently fixed) and seven H atoms were included in our analysis. thus comprising 60 internal coordinates: 21 bond lengths, 20 bond angles, and 19 dihedral angles. The entropy change computed by eq 3 was -6.82 cal mol-' K-I. The difference between the latter result and 1.98 cal mol-' K-l obtained by eq 2 leads to the 100-fold lower equilibrium constant for the chair F? boat isomerization reaction, i.e., 100-fold slower chair-to-boat transformations. This represents more accurately the overall picture obtained in the MD simulations.
lo Wavenumber (cm-1) Figure 10. Fourier transform spectra of atomicvelocities from simulations at 2370 K: top panel, from 16 principal carbon atoms of the surface model shown in Figure 1; middle panel, from the carbon atoms of the C,Hs step of the model surface in Figure 8; bottom panel, from the atoms of the free C6Hll radical. The resolutions are 20, 10, and 10 cm-1, respectively.
E. Bond Breaking. Reactions involving the formation and breaking of chemical bonds, and especially of carbon-carbon bonds, are of particular interest to elucidation of underlying mechanisms of diamond CVD. Unfortunately, at temperatures of interest to diamond CVD, i.e. at and below 1200 K, C-C and C-H bond-scission reactions occur on the time scale of a microsecond and longer and hence are inaccessible to MD simulations. One way of overcoming this difficulty is to perform simulationsat much higher temperatures, when the characteristic time of bond breaking is on the order of 100-loo0 fs. Such an approach was utilized recently by Garrison et a1.I0 in their numerical studies of the CH2 surface radical insertion into a diamond (100)-(2Xl) dimer. Here we investigate bond-breaking reactions of the CjHs-step model by subjecting it to high temperatures. Our objective is 2-fold: to make a comparison between the analogous surface and gas-phase reactionsand to comparethe MD and staticTST results. First, we examined thevibrational spectraof the C3H~-stepmodel, as they compare to those of the model surface without the step (Le., the structure shown in Figure 1) and the free cyclohexyl radical, C6Hll. Figure 10 depicts such a comparison for 2730 K. Qualitatively similar results were obtained for other temperatures. As can be seen in Figure 10, the general features of the CjHs-step model (Figure 10, middle panel) resemble more closely those of the free C6HII radical (bottom panel) rather than those of the bulk surface model (top panel). This is somewhat surprising, considering the relatively fast energy transfer in the system and the phonon interaction with the chair-boat conformations discussed in the previous subsection. Nonetheless, the results shown in Figure 10 indicate that the C3H5 step behaves overall closer to a gas-phase species than to the bulk solid. The
Zhao et al.
1646 The Journal of Physical Chemistry, Vol. 97, No. 8, 1993
TABLE II: Conditions of the Bond-Breaking MD Simulations initial equilibration period temperature (K)
bond breaking
time* (fs)
temperature increment (K/fs)
s2 s3 s4 s5 S6
2000 2500 2500 2500 2000 2SOO
1000 100 200 100 200 100
500 230 500 500 500 5
G1 G2 G3 G4 G5 G6 G7 G8 G9
2000 2000 0.04 0.30 2500 2500 0.04 3000 3000
1500 1500
500 5 5 5 500 5 4 adiabatic run adiabatic run
casea
s1
200 200
final temperature (K)
C3Hs-Step Surface Model 2500 2730 3000 3000 3000 3000 Free Cyclohexyl Radical 2500 2500 3000 3000 3000 3000 3000
integration step size (fs)
timeC (fs)
typed
1.o 1.o 1.o
3300 2350 520 250 450 480
C(27)-C( 15) C(17)-H(21) C(17)-C(14) C( 17)-C( 14) C(27)-H(29) C( I3)-H( 18 ) r
0 0
C(2)-C(3) C(S)-C(6) 1,4 H shift C(5)-C(6) C(2)-C(3) C(5)-C(6) C(5)-C(6) C(2)-C(3) 1,2 H shift
0 0 A 0 0 0 0
1.O
I .o 0.5 1.o
1 .O
1.o
1.o
1.o
1.o 1.o 0.5 0.25
2400 13365 4925 2025 100 320 760 3300 500
symbol in Figure 1 1
0 0
A X
0 A
Computations for all cases except G8 and G9 were carried out under isothermal condition by scaling the temperature of the substrate for the C3Hs-stepsurface model or the whole molecule for the free C ~ HI Iradical. Runs SI and SS had identical initial conditions, those scaled from the surface model equilibrated at 1200 K. The initial conditions for each of the runs S2, S3, S4, and S6 were the conditions attained at the end of the temperature ramp in run SI. Runs G 1 and G2 had identical initial conditions. The initial conditions for each of the runs G5 and G6 were the conditions attained at the end of the temperature ramp in run GI. Runs G8 and G9 were carried out adiabatically; the initial conditions of run G8 were those attained at the end of the temperature ramp in run G3, and those of G9-in G4. The time period during which the system was allowed to equilibrate before the temperature ramp. The reaction time-the period of time from the onset of the temperature ramp to the instant of the bond breaking. The bond broken. The numbers in parentheses refer to the numbering in Figure 8 for the C3Hs-step surface model, and those for the cyclohexyl radical refer to the regular numbering convention, starting with the radical carbon site as C(1). eSurface atom H(18) migrated from C(13) to C(28). (1
1"
n
0 C(27)-C(lS) 0 C(17)-C(14) 0 C(17)-H(21) A C(27)-H(29) x C(13)-H(18) -C-CcIsT)
10''
D
b
ei 0
0
10"
Figure 11. Arrhenius plots for bond breaking reactions: top panel, for theC3H5-stepmodel of Figure 8; bottom panel, for the free C ~ H1 Iradical. Symbols are the results of the molecular-dynamics simulations and the lines the transition-state-theory estimates without zero-point energy corrections.
molecular-dynamics results reported below are consistent with this conclusion. The results of the MD simulations for the CSHrstep model and analogous free CsHlI radical are presented in Figure 1 1. In each run, the model surface was first equilibrated at lower
temperature, then the temperature was increased in a prescribed manner, and when a desired temperature was reached, the run continued at this final temperature until a bond was broken. The temperature and duration of the equilibration period, the rate of the temperature ramp, the final temperature, and the integration time step used in the MD simulations are summarized in Table 11. Displayed in Figure 11 are the rate constants for various bond breaking events, defined as the inverse of the reaction times; this formalismfollows the Poissonian natureof first-order reaction kinetics. The reaction time is defined in our study as the time interval between theonset of the temperature rampand theinstant of the bond breaking. These time intervals are listed in Table 11. Also shown in Figure 11, for comparison, are the rateconstants of C-C bond breaking obtained in the canonicalTST calculations using the static-AM1 results. As before, the TST results do not include zero-point energies, to be consistent with the MD simulations; with this correction, the rate constants are 15 and 50% higher, for the surface and free-radical simulations, respectively, which is again too small to affect our discussion. Inspection of Figure 11 indicates a substantial scatter in the obtained rate constant values for both surface and gas-phase cases. Such large scatter is expected. Different points in the graphs represent not necessarily the same reaction events, and some variations are expected for diffcrent bonds. More importantly in our case are the 90%confidencelimits for a Poisson-distributed variable, which are a factor of 10up and down from the expected value. It is clear that such computer-intensive calculations as undertaken in the present study with a quantum potential are not feasible for establishingreaction statistics. Nonetheless,we note that all the MD points lie well above the TST results and the zero-point energies cannot reconcile this discrepancy, nor can it be caused by our method of implementing periodic boundary conditions as the discrepancy is the same for both constrained (surface) and unconstrained (gas-phase) MD simulations. Since the harmonic TST provides an upper limit for the rate constant, and possible phonon interactions should decrease not increase this estimate,63we concludethat the discrepancyobserved in this study is probably due to the limited size of the model surface and, as a result of it, unrealistic temperature control causing excessive fluctuations in the energy redistribution. Indeed, inspection of
Molecular Dynamics with AM1 Potential thedata listed in Table I1 shows that starting with identical initial conditions but using more gentle heating results in longer reaction times, Le., in lower rate constantvalues, closer to the TST estimates (compare, e.g., cases S2, S3, and S6 or G1 and G2). On the other hand, extending the equilibration period or decreasing the rate of the temperature ramp increases the uncertainty in the onset of reaction and hence in the resulting rate constant. These factors seem to limit the reliability of MD simulations for bond-breaking reactions. We finally note that the TST rate constant values for the analogous reactionsat the surface,the C( 14)-C( 17) bond cleavage in the C3H5-stepmodel, and gas phase, C(2)-C(3) bond cleavage in the cyclohexyl free radical, are close to each other, 5.4.X 1O1O and 3.5 X 1Olo s-I at 2730 K, respectively, and the two are in good agreement with the experimental data typical for the gas-phase ,!?-scission reactions of alkyl radicals. For example, the highpressure-limit rate constant of reaction
*CHZ-CH,-CH,
-.+
C H 2 4 H 2 + CH,
is644.5 X lolos-I at 2730 K, as compared with the corrected for the zero-point energy results of 5.8 X 1O1Oand 5.6 X 1O1Os-1 obtained for the surface and gas-phase reaction at 2730 K in the present study. Interestingly, at 1200 K, the TST rate constants of these two reactions, corrected for the zero-point energy, are 2.1 X lo7 and 9.2 X lo6 s-I: Le., the surface reaction is about twice as fast as its gas-phase analog, reflecting the slightly lower potential energy barrier (34 vs 37 kcal/mol) and slightly tighter transition state of the former. Thus, although the two analogous surface and gas-phase reactions consideredhere exhibit certainly a similarity, there are also differences, as in the other examples discussed above, induced by the confinement of the surface.
4. Summary and Conclusions Classical molecular dynamics simulations of diamond surface reactions were performed using quantum mechanically derived forces. A semiempirical AM1 potential from the MOPAC quantum-chemistry program package was employed in the calculations. The evolving molecular geometry, energetics, bond orders, and vibrational frequencies were monitored during the simulations. Several examples of reactions of gaseous reactants with the diamond (1 11) surface, conformationand bond breaking of a surface step, and reactions of free gaseous analogs were examined. It is demonstrated that the use of quantum-mechanical forces in classical trajectory simulations is feasible for the analysis of chemical reactions occurring on surfaces. This approach may provide results not easily obtainable by other means. One such finding in the present study is the existence of a trapped, physisorption state for an acetylene molecule colliding with a diamond surface at 1200 K, at a temperature typical of diamond CVD. The implication of this result is 2-fold: first, the estimation of the rate constants for gas-surface reactions of hydrocarbon precursors, based at present solely on the assumption of a direct, Eley-Rideal-type mechanism, should be reconsidered: second, the surface diffusion of hydrocarbon species should be taken into account as it affects the not fully understood morphology of diamond growth. The main question that we were intested in is the hypotheses of chemical similarity for the carbon-hydrogen system, postulating that the rate coefficients of gas-surface reactions can be estimated from analogous reactions taking place in the gas phase.26 The results obtained in all the cases tested in the course of the present study indicatethat although a considerable similarity indeed exists, significant quantitativedifferences are induced by the confinement of thesurface: interactions with neighboring atoms, as for example in the reaction of a hydrogen atom with a surface radical (section 3C); influence of surface phonons, as in the boat-chain conformation (section 3D); anddifferingstiffnessof the transition states,
The Journal of Physical Chemistry, Yo/.97, No. 8, 1993 1647 as in the bond-breaking reactions (section 3E). The existenceof a trapped, physisorptionstate provides an additional, important factor contributingto thedifference. Thus, althoughthesimilarity concept can serve as a source for rough, order-of-magnitude estimates, the quantitative evaluation of the rate parameters for surface and gas-surface reactions needs direct experimental, theoretical, and numerical studies.
Acknowledgment. The computations were performed using the facilities of the Pennsylvania State University Center for Academic Computing. The work was supported by Innovative Science and Technology Program of the Strategic Defense Initiative Organization (SDIO/IST) via the U S . Office of Naval Research, under Contract No. N00014-92-J- 1420, and by the Office of Exploratory Research at the Electric Power Research Institute, Project No. RP2426-5 1. References and Notes (1) Brenner, D. W.; Garrison, B. J. Adu. Chem. Phys. 1989, 76, 281. (2) Agrawal, P. M.; Raff, L. M.; Thompson, D. L. Surf. Sci. 1987,188, 402. (3) Agrawal, P. M.; Thompson, D. L.; Raff, L. M. Surf.Sci. 1988,195, 283; J . Chem. Phys. 1989, 91, 5021. (4) Stillinger, F. H.; Weber, T. A. Phys. Rev. Lett. 1989,62,2144; Weber, T. A.; Stillinger, F. H. J . Chem. Phys. 1990, 92, 6239. (5) Srivastava, D.;Garrison, B. J.; Brenner, D. W. Phys. Rev. Lett. 1989, 63, 302; Langmuir 1991, 7, 683. (6) Smith, R. Proc. R. SOC.London, A 1990, 431, 143. (7) Brenner, D. W.; Dunlap, B. I.; Mintmire, J. W.; Mowrey, R. C.; White, C. T. In Proceedings of the Second International Symposium on
Diamond Materials; Purdes, A. J., Angus, J. C., Davis, R. F., Meyenon, B. M., Spear, K. E., Yoder, M., Us.The ; Electrochemical Society: Pennington, NJ, 1991; p 39. (8) Pailthorpe, B. A. J . Appl. Phys. 1991, 70, 543. (9) Weakliem, P. C.; Carter, E. A. J. Chem. Phys. 1992, 96, 3240. (IO) Garrison, B. J.; Dawnkaski, E. J.; Srivastava, D.; Brenner, D. W. Science 1992, 255, 835. (11) Stillinger, F. H.; Weber, T. A. Phys. Reo. B 1985, 31, 5262. (12) Tersoff, J. Phys. Reo. Lett. 1986, 56, 632; Phys. Reu. B 1988, 37, 6991; Phys. Reu. Lett. 1988, 61, 2879. (13) Brenner, D. W.; Garrison, B. J. Phys. Reu. E 1986, 34, 1304. (14) Brenner, D. W. Phys. Reu. B 1990, 42,9458. (1 5) Takai, T.; Lee, C.; Halicioglu, T.; Tiller, W. A. J. Phys. Chem. 1990, 94. 4480. (16) Abell, G. C. Phys. Reu. B 1985, 31,6184. (17) Pailthorpe, B. A.; Mahon, P. Thin Solid Films 1990, 193/194, 34. Mahon, P.; Pailthorpe, B. A.; Bacskay, G.B. Philos. Mag. E 1991,63, 1419. (18) Weakliem, P. C.; Wu, C. J.; Carter, E. A. Phys. Reo. Lett. 1992,69, 200. (19) Drabold, D. A,; Wang, R.; Klemm, S.; Sankley, 0.;Dow, J. Phys. Reu. B 1991, 43, 5132. (20) Hartke, B.; Carter, E. A. Chem. Phys. Lett. 1992. 189. 358. Maluended, S. A.; Dupuis, M. In?. J . Quantum-Chem. 1992, 42, 1327. (21) Buda, F.; Chiarotti, G.L.; Car, R.; Parrinello, M. Phys. Reo. B 1991, 44,5908. Buda, F.; Chiarotti, G. L.; Car, R.; Parrinello. M. Phys. Reo. Lett. 1989, 63, 294. Stich, I.; Car, R.; Parrinello, M. Phys. Reu. Lett. 1989, 63, 2240. Hohl, D.; Jones, R. 0.;Car, R.; Parrinello, M. J . Am. Chem. Soc. 1989, 111. 826. (22) Stewart, J. J. P. In Reuiews in Computational Chemistry; Lipkowitz, K. B., Boyd, D. B., Eds.; VCH: New York, 1990; p 45. Boyd, D. B. ibid., p 321. (23) Martin, J. M.; Francois, J. P.; Gijbels, R. J. Comput. Chem. 1991, 12, 52. Coolidge, M. B.; Marlin, J. E.; Stewart, J. J. P. J . Comput. Chem. 1991, 12, 948. (24) Field, M. J.; Bash, P. A.; Karplus, M. J. Comput. Chem. 1990, 11, 700. (25) For example: Diamond and Diamond-Like Films; Dismukes, J. P.. Purdes, A. J., Spear, K. E., Meyerson, B. S., Ravi, K. V., Moustakas, T. D., Yoder, M., Eds.;TheElectrochemicalSociety:Pennington,NJ, 1989. Carbon 1990, 26 (6). J . Mater. Res. 1990 5 (11). Proceedings of the Second International Symposium on Diamond Materials; Purdes, A. J., Angus, J.
C., Davis, R. F., Meyerson, B. M., Spear, K. E., Ycder, M., Eds.; The Electrochemical Society: Pennington, NJ, 1991. Diamond and DiamondLike Films and Coatings;Clausing, R. E., Horton, L. L., Angus, J. C., Koidl, P., Eds.; Plenum: New York, 1991; NATO-AS1 Series B: Physics, Vol. 266. (26) Frenklach, M.;Spear, K. E.J. Mater. Res. 1988,3,133. Frenklach, M. J . Appl. Phys. 1989,65,5142. In Carbon in the Galaxy: Studies From Earth and Space; Tarter, J. C., Chang, S., DeFrees, D. J., Eds.; NASA Conference Publication 3061, 1990; p 259; In Diamond and Diamond-Like Films and Coatings; Clausing, R. E., Horton, L. L., Angus, J. C., Koidl, P., Eds.; Plenum: New York, 1991; NATO-AS1 Series B: Physics, Vol. 266, p 499. (27) Frenklach, M.; Wang, H. Phys. Rev. B 1991, 43, 1520. (28) Belton, D. N.; Harris, S . J. J . Chem. Phys. 1992, 96, 2371.
1648 The Journal of Physical Chemistry, Vol. 97, No.8, 1993 (29) Stewart, J. P. P. Quant. Chem. Prog. Exch. (QCPE) 1983, No. 455. (30) Huang, D.; Frenklach, M.; Maroncelli, M. J. Phys. Chem. 1988,92, 6379. Huang, D.; Frenklach, M. J. Phys. Chem. 1991, 95, 3692; J . Phys. Chem. 1992, 96, 1866. (31) Valone, S. M. In Diamond and Diamond-Like Films and Coatings;
Clawing, R. E., Horton, L. L., Angus, J. C., Koidl, P., Eds.; Plenum: New York, 1991; NATO-AS1 Series B: Physics, Vol. 266; p 525. (32) Dannenberg. J. J.; Rayez, J. C.; Rayez-Meaume, M. T.; Halvick, P. J . Mol. Srrucr. (Theochem) 1985, 123, 343. (33) K(lrtvtlyesi,T.;Seres,L.J.Mol.Srruct. (Theochem) 1991,251,123. (34) Valone, S.M.; Trkula, M.; Laia, J. R. J . Mater. Res. 1990,5, 2296. (35) Tsuda. M.;Nakajima, M.; Oikawa. S.J . Am. Chem.Soc. 1986,108, 5780; Jpn. J . Appl. Phys. 1987, 26, L527. (36) Deak. P.; Giber, J.; Oechsner, H . Surf, Sci. 1991, 250. 287. (37) Venvoerd, W. S.Surf. Sci. 1981, 108, 153. (38) Tsuda, M.;Oikawa, S.;Sekine, C.; Furukawa. S.In Proceedings of the Second International Symposium on Diamond Materials; Purdes, A. J., Angus, J. C., Davis, R. F., Meyerson, B. M., Spear, K. E.,Yoder, M., Eds.; The Electrochemical Society: Pennington, NJ, 1991; p 154. (39) Mehandru, S.P.; Anderson, A. B. J . Mater. Res. 1990,5,2286; Surf. Sci. 1991, 248, 369. (40) Zheng, X. M.; Smith, P. V. Surf. Sci. 1991.256, 1; 1992, 261, 394; 1992, 262, 219. (41) Stewart, J. Personalcommunication: Thecalculationswaeperformed with a 64-atom cell using beta version of MOPAC 7.0. (42) Bachelet, G. B.; Baraff. G. A.; Schliiter, M. Phys. Reo. B 1981, 24, 4736. Yin, M. T.; Cohen, M. L. Phys. Reo. B 1984, 29,6996. Chang, K. J.; Cohen, M. L. Phys. Reo. B 1987, 35, 8196. Bernholc, J.; Antonelli, A.;
Del Sole, T. M.; Bar-Yam, Y.; Pantelid-, S.T. Phys. Reo. Lett. 1988, 61,
2689. (43) Kittel, C. Introduction IO Solid State Physics, Wiley: New York, 1986; p 55. Lide, D. R.. Ed. Handbook of Chemistry and Physics, 71st ed.; CRC: Boca Raton, FL, 1990; p D-24. (44) For an explanation of the minimum image convention and molecular
dynamic techniques in general, see Allen, M. P.; Tildesley, D. J. Computer SimulationofLiquids; Clarendon: Oxford, 1987. Heennann, D. W. Computer Simulation Methods in Theoretical Physics; Springer-Verlag: Berlin, 1990. Gould, H.; Tobochnik, J. An Introduction to Computer Simulation Methods, Part I ; Addison-Wesley: Reading, MA, 1988. (45) Beeman, D. J. Compur. Phys. 1976,20, 130. (46) Verlet, L. Phys. Reo. 1967, 159, 98. 147) Futrelle. R. P.: McGintv. D. J. Chem. Phvs. Lett. 1971. 12. 285. (48) Simandiras, E: D.; Rice: J. E.;Lee, T. J.; Amos, R. D.; Handy, N. C. J . Chem. Phys. 1988,88, 3187.
Zhao et al. Strey, G.;Mills, I. M. J. Mol. Srruct. (Theochem.) 1976, 59, 103. Weber, W. Phys. Rev. B 1977, 10,4789. Herchen, H.; Cappelli, M. A. Phys. Reo. B 1991.43, 11740. Wang, C. Z.; Chan, C. T.; Ho, K.M. Phys. Reo. B 1990,42,11276. For example: Taylor, J. B.; Langmuir, I. Phys. Reu. 1933,44, 423. Boudart. M.; DjCga-Maradassou, G.Kinetics of Heterogeneous Catalytic Reactions; Princeton University Press: Princeton, NJ, 1984. Ce er, S. T. Annu. Reo. Phys. Chem. 1988,39,479. Interaction of Atoms and &olecules with Solid Surfaces; Bortolani. V.,March, N. H., Tosi, M. P., Eds.; Plenum: New York, 1990. Zhdanov. V. P. Elementary Physicochemical Proccsseson Solid Surface; Plenum: New York, 1991. Dynamics of Gas-Surface Interactions; Rettner, C. T., Ashfold, M. N. R., Eds.; The Royal Society of Chemistry: Cambridge, 1991. (54) Grimley, T. B. In Interaction of Atoms and Molccules with Solid Surfaces; Bortolani, V., March, N. H ,Tosi, M.P., Eds.; Plenum: New York. 1990; p 25. ( 5 5 ) Weinberg, W. H. In DynamicsofGas-SurfrrceInreracrions;Rettner, C. T., Ashfold, M. N. R., Eds.; The Royal Society of Chemistry: Cambridge, 1991; p 171. (56) Frenklach. M. Phys. R w . B 1992,45, 9455. (57) Brenner, D. W.; Robertson, D. H.; Carty, R. J.; Srivastava, D.; Garrison, B. J. Presented at the MRS Spring Meeting, San Francisco, CA, April 1992. (58) Angus, J. C.; Hayman, C. C. Science 1988,241,913. Angus, J. C.; Hoffman, R. W.; Schmidt, P. H. InScienceand TechnologyofNewDiamond; Saito, S.,Fukunaga, O., Yoshikawa, M., Eds.; Terra Scientific Publishing Co.: Tokyo, Japan, 1990; p 9. (59) Spear, K. E. J. Am. Ceram. Sci. 1989. 72, 171. Spear, K. E.; Frenklach, M. In DiamondandDiamond-Like Films; Dismukes,J. P., Purdes, A. J., Spear, K. E., Meyerson, B. S.,Ravi. K.V., Moustakas, T. D., Yoder, M., Eds.; The Electrochemical Society: Pennington, NJ, 1989; p 122. Spear, K. E.;Phelps, A. W.; White, W. B. J. Mater. Res. 1990, 5, 2277. (60)For example: Eyring, H. J . Chem. Phys. 1935,3,107. Kreevoy, M. M.; Truhlar, D. G.In Inuestigation of Rates and Mechanisms of Reactions; Bernasconi, C. F.. Ed.; Wiley: New York, 1986; Part I. p 13. Steinfeld, J. I.; Francisco, J. S.;H a s , W. L. Chemical Kinetics and Dynamics; Prentice Hall: Englewood Cliffs, NJ, 1989. (61) For example: Hill, T. L. An Introduction to Statistical Thermodynamics; Dover: New York, 1986. (62) Karplus, M.; Kushik, J. N. Macromolecules 1981, 14, 325. (63) Straus, J. B.; Voth, G.A. J . Chem. Phys. 1992, 96, 5460. (64) Tsang, W. J . Phys. Chem. Ref. Data 1988. 17, 887. (49) (50) (51) (52) (53)
