J . Phys. Chem. 1993,97, 12134-12143
12134
Theoretical Investigations of and Xe Matrices
0 3
Vibrational Relaxation and Oxygen Atom Diffusion Rates in Ar
M. Beth Ford, Angel D. Foxworthy, Gilbert J. Mains, and Lionel M. Raff Department of Chemistry, Oklahoma State University, Stillwater, Oklahoma 74078 Received: May 12, 1993; In Final Form: August 17, 1993'
The molecular dynamics of vibrationally-excited ozone isolated in Ar and Xe matrices a t 12 K are investigated using classical trajectory methods. Oxygen atom diffusion in these matrices are computed using a classical variational transition-state theory which employs a new Markov walk/damped trajectory procedure to effect convergence. The matrix model consists of a face-centered-cubic (fcc) crystal with 125 unit cells (666 atoms) in a cubic (5 X 5 X 5) arrangement. Ozone, or the oxygen atom, is held interstitially within the innermost unit cell of the crystal. The system potential is written as the sum of a lattice potential, a lattice4zone or oxygen atom interaction, and a gas-phase potential for ozone. The first two potentials are assumed to have pairwise form, while the ozone molecular potential is the one developed by Murre11 and Farantos [Mol. Phys. 1977,34, 11851. The oxygen atom/Xe pair potentials are computed a t the Hatree-Fock and Miiller-Plesset second-order perturbation level of theory for the singlet and triplet states using two different pseudopotentials for the xenon core. Vibrational relaxation is observed to be mode specific with the bending mode dominating the rate of energy transfer to the lattice. The rates over the first 0.1-0.2 ps are characterized by near-linear, first-order decay. The energy-transfer rates are found to be significantly faster in Ar than in Xe matrices. Thermal diffusion rates of oxygen atoms in fcc Xe crystals are computed for a range of possible O / X e Lennard-Jones (12,6) interaction potentials suggested by the pseudopotential calculations. In all cases, diffusion is found to be extremely slow with an activation energy greater than 3.3 kcal/mol. Quantum tunneling is shown to make a negligibly small contribution to these rates. The corresponding diffusion rates in Ar are even slower. Comparison of the results with measured diffusion coefficients indicates that almost all of the experimentally observed diffusion is occurring along lattice defects, grain boundaries, vacancies, and other lattice imperfections.
I. Introduction Energy transfer, mobility, and diffusion often play critical roles in reactions occurring under matrix isolation conditions. When the matrix is used as a means for containing high-energy density materials, the importance of these effects is obvious. They are often equally important in other systems where their influence is less obvious. The bimolecular addition of F2 to ethylene is a case in point.'" Our previous theoretical studiesb8 show that the energy-transfer rate from the vibrationally-excited 1,2difluoroethane product to the lattice phonon modes exerts a major influence on the final product ratio. Diffusion rates within the matrix cage directly affect the extent to which an atomic addition mechanism is involved in the reaction.* At present, there is very little experimental data related to diffusion, mobility, and energy transfer in matrices and even fewer theoretical treatments. Jacox9Jo has observed the reaction of atomic fluorine with both CH4 and CH3F. In particular, the reaction F
+ CH,F
-
CH,F
+HF
is presumed to occur uia F atom diffusion through the argon matrix at 14 K until it reaches a site containing CH3F. However, no measurements or estimates of such diffusion rates are reported. Feld, Kunttu, and Apkarian" have reported measurements of fluorine atom mobilities in an argon matrix subsequent to photodissociation. Alimi, Gerber, and ApkarianI2 have carried out molecular dynamics simulations of these photodissociation experiments in which the matrix is modeled as a face-centeredcubic (fcc) lattice of 365 argon atoms. No provisions are made for the dissipation of excess energy to the bulk. More recently, Lawrence and ApkarianI3 have measured the mobility of photoexcited oxygen atoms in Xe lattices. *Abstract published in Advance ACS Abstracts, November 1, 1993.
0022-3654f93 f2097-12134%04.00/0
Krueger and Weitz14 have measured oxygen atom recombination rates in a Xe matrix at 32 and 40 K. The results show that the oxygen atom lifetime is on the order of 50-60 h. Consequently, oxygen atom diffusion rates are very small. The kinetic data for recombination do not follow a simple, secondorder rate law. Initially, there is a very rapid drop of the oxygen atom concentration followed by a much slower rate of decrease. Krueger and Weitz14 have shown that an acceptable fit to the data can be obtained only if one assumes that there exists a distribution of oxygen atom diffusion coefficients in the matrix. If it is assumed that this distribution consists of precisely two distinct coefficients,their values may be uniquely extracted from the data. In the present paper, we report the results of theoretical studies of energy transfer and diffusion in rare-gas matrices at 12-80 K. Classicaltrajectory methods are employed to examine the energytransfer dynamics of vibrational-excitedozone. The dependence of the rate of internal energy transfer to the lattice phonon modes on the magnitude of the excitation, matrix composition, and the internal energy partitioning among the 0 3 vibrational modes is investigated. The results show the presence of mode specificity in the transfer rates. Oxygen atom thermal diffusion rates in face-centered-cubic (fcc) matrices are obtained using classical variational transition-state theory. Comparison of these results to the data reported by Krueger and Weitz14shows that essentially all of the diffusion observed in the experimental systemI4occurs along lattice defects grain boundaries, etc. 11. Matrix Model
In the present study, energy transfer and diffusion in facecentered-cubic rare-gas crystals are investigated. The matrix model used is the ( 5 X 5 X 5 ) fcc lattice of 125 unit cells containing 666 lattice atoms previously described by Raff.6 This model has been found to be sufficiently large to accurately represent the 0 1993 American Chemical Society
Molecular Dynamics of O3 in Ar and Xe Matrices
The Journal of Physical Chemistry, Vol. 97, No. 47, 1993 12135
density and volume expansion upon trapping of 1,Zdifluoroethanea6The heat bath effects ofthecrystal bulkarealsoadequately represented provided velocity reset functions or Langevin equations are employed on appropriate lattice atoms in the trajectory calculations.6 For computational convenience, we have utilized the velocity reset procedure developed by Riley et ~ 1 . ~To5 do so, the 666 lattice atoms are first divided into three discrete zones. This division is determined once the size of the lattice model has been chosen using density, volume expansion, or other criteria. The boundary zone (B zone) comprises the atoms located on the boundary of the crystal. There are 302 such atoms in the present case. The positions of these lattice atoms are fiied. Their presence reduces edge effects and maintains the desired lattice symmetry. The secondary zone (Q zone) contains all lattice atoms that have at least one B zone atom as a nearest neighbor. In the present model, there are 302 such atoms. The solution of Hamilton’s equations for the motion of these atoms is modified by the velocity reset functions associated with each atom in the Q zone. This procedure maintains the temperature of the lattice as energy is removed or inserted by the chemical or physical processes that areoccurring. Consequently, the model simulates an experiment conducted under constant-volume and -temperature conditions. The primary zone (P zone) comprises the remaining atoms of the crystal (62 in the present case). The motion of these atoms is affected only by the forces produced by the interaction potential. Riley et al.lf’ have shown that the velocity reset method is equivalent to a Langevin procedure in the limit of zero time between resets.
Ozone-Matrix System. The total potential, VT,of the matrixisolated ozone system is taken to be the separable sum
+ VI + V,,
where V, is the interaction potential between lattice atoms, 6 is the matrix-zone interaction, and V,, is the molecular ozone potential in thegas phase. The potential for thelatticeinteraction is assumed to have the pairwise form N
where rij is the distance between lattice atoms i and j , and N is the total number of lattice atoms in the matrix model. is taken to be a Morse potential with a cutoff radius, rij = R,, given by
K j = D[exp(-2a(rv
- ro)) - 2 exp(-a(rij - ro)]] for rij I R, for rij > R,
Vv = 0
(3)
The parameters for the Morse potential are taken from data reported by LeroyI6 for the Ar lattice and from Barker et ~1.1’ for the Xe lattice. The cutoff radius, R,, is chosen so that V,,(r,,=RJ is approximately 4% of the two-body well depth, D. There is very little experimental data relating to the interaction of rare gases with ozone in a matrix environment. Consequently, we have assumed the intermolecular term VIto also have pairwise form N
1.745 000 (A-1) 3.810 00 (A) 6.000 (A)
a r0
Rc
XeXe 0.024 214 05 (eV) 1.467 000 (A-1) 4.362 30 (A) 7.000 (A)
D a
ro Rc
Ar-O t
U
Xe-0 c 0
3
(4)
where r i k is the interatomic distance between lattice atom i and oxygen atom k. The V,k term is taken to be a Lennard-Jones
0.005 379 8 (eV) 3.213 (A)
0.007 570 2 (eV) 3.549 39 (A)
(12,6) potential given by
v,,
= e((u/rik)l2 - 2(u/rik)16 (5) The oxygen-lattice interaction parameters are derived from the data for the rare gas-rare gas and 0-0interactions using the combining rules ‘Rg-0
= (cRg-RgcM)
112
(6)
and 1
= /2[(a)Rg-Rg + (7) The potential parameter values for the ozone/matrix system are given in Table I. For the ozone intermolecular potential, V,,, we have used the global analytic function suggested by Murre11 and Farantos.18 This potential accurately reproduces harmonic force constants, Osvibrational frequencies, equilibrium bond lengths, and the dissociation energy of the molecule. Oxygen Atom-Matrix System. The functional form of the potential hypersurface used for the oxygen a atom/matrix system is similar to that employed for ozone. The total potential is the sum of a pairwise lattice potential and a pairwise oxygen-lattice interaction term. That is, (‘)Rg-4)
111. Potential Energy Surfaces
V, = V ,
TABLE I: Ozone/Lattice Pairwise Potential Parameters parameter value Ar-Ar D 0.012 228 90 (eV)
v, = v, + v,,
(8) where V, is given by eq 2. The form of the oxygen-lattice interaction, V f , is the pairwise expression given by eq 4. The O/Rg (Rg = rare gas) interaction for ozone is expected to arise principally from weak van der Waals dispersion forces since it is a closed-shell interaction. For such a system, eqs 5-7 are expected to give a reasonable description of the potential. For a bare oxygen atom/Xe pair, however, the interaction parameters are considerably less certain. The presence of vacant orbitals on oxygen opens the possibility of a much stronger coordinatecovalent interaction with a Xe atom. As a result, c may be much larger than the value predicted from eq 6. If the oxygen atom is in its ground state (3P),there are no completely vacant valence-shell orbitals present. Under these conditions, it seems likely that the O/Xe will be primarilyrepulsive since there are no unpaired electrons on the xenon atom. If, however, the oxygen atom is in an excited singlet state, there could be a bonding interaction between the vacant p orbital on oxygen and one of the electron pairs on Xe. In this case, we would expect a much stronger attractive interaction and a significantly reduced equilibrium X e O bond distance. We have carried out threesets of ab initiocalculations toobtain a better estimate of the Xe/O interaction parameters. In the first set of calculations, we employ a double-zeta basis set combined with the pseudopotential for the xenon core developed by Wadt
12136 The Journal of Physical Chemistry, Vol. 97, No. 47, 1993
Ford et ai.
TABLE II: Pseudopotential Calculations of Xe/O Pair Potential (S Dewtea Singlet State, T Denotes the Triplet State) (a4
4 (A)
0.001 25 e
2.719
0.023 79 0.OOO76
2.150 4.479
O.OO0 66 e
2.871
0.000 48 0.001 07
2.871 4.342
O(T) 0.026 25 XeO(S) XeO(T) 0.000 05 Reference 19. CReference21. dReference 21 with split s, p, and d orbitals (SPD). eNot bound. 0 Well depth in au.
2.107 4.342
atom/molecule Xe
basis set DZ-pseudocore potb
method
energy (au)
HF HF HF HF HF MP2 MP2 MP2 MP2 MP2
-15.224 33 -74.678 56 -74.803 83 -89.904 14 -90.0269 -15.234 96 -74.739 65 -74.846 92 -89.998 40 -90.082 64
HF HF HF HF HF HF HF HF HF HF MP2 MP2 MP2 MP2 MP2
-126.032 77 -74.655 81 -74.780 31 -200.689 24 -200.8 13 09 -126.032 95 -74.655 81 -74.780 31 -200.689 24 -200.814 33 -126.120 30 -74.722 43 -74.830 72 -200.868 97 -200.951 06
O(S)
O(T) XeO(S) XeO(T) Xe O(S)
O(T) XeO(S) XeO(T) Xe O(T) XeO(S) XeO(T) Xe O(S)
O(T)
XeO(S) XeO(T) Xe
O(S)
relativistic corec
relativistic core + SPW
and Hay.I9 Calculations at the HatreeFock (HF) and MBllerPlesset second-order perturbation theory (MP2) levels using GAUSSIAN 9220 have been carried out. The results are given in Table 11. The calculations show that the XeO triplet is the lowest-energy configuration of the system, as expected. This state is unbound at the H F level of theory. At the MP2 level, it is slightly bound with a well depth of 0.000 76 au (0.0208 eV) at an equilibrium Xe-0 separation of 4.479 A. As qualitative considerations would suggest, the XeO singlet state is more attractive. At the HF level, it exhibits a well depth of 0.001 25 au (0.0340 eV). The well depth increases to 0.0238 au (0.6473 eV) a t an equilibrium separation of 2.150 A at the MP2 level of theory. A second set of calculations at the HF level employs a pseudopotential for the xenon corein which relativisticcorrections have been incorporated.2I These results are also given in Table 11. With a minimal basis set, the ground XeO triplet state is again predicted to be unbound. The excited singlet is bound with a well depth of 0.0178 eV at R, = 2.871 A. If the s, p, and d orbitals are split (SPD), the HF calculations predict the triplet to be bound by 0.029 1 eV with an equilibrium separation of 4.342
A.
A third set of calculations using the relativistic core potentialzl a t the MP2 level of theory with split s, p, and d orbitals predicts that the XeO ground state is the triplet. There is a very shallow well (0.001 30 eV) at R, = 4.342 A which probably does not accommodate a bound vibrational state. The excited XeO singlet is bound by 0.714 eV at an equilibrium distance of 2.1067 A. It seems clear that in the Xe matrix environment the oxygen atom exists primarily in its ground triplet state. The 02 recombination experiments reported by Krueger and Weitz'4 are carried out over a time span of 50-60 h. Any excited singlet oxygen initially formed in the photolysis process would have decayed to the ground state on this time scale. We therefore expect the Xe/O interaction to be predominantly repulsive. We further expect the well depth characterizing the interaction to be greater than that predicted by the HF or MP2 calculations since these calculations do not contain sufficient electron correlation to properly account for the dispersion forces primarily responsible for the attractive well. In general, we would expect these dispersion forces to produce an attractive potential well somewhat deeper than that predicted by eq 6 for ozone since the oxygen
TABLE 111: Oxygen Atom/Lattice Potential Parameters for the W(12,6) Potentials Employed in the Diffusion Calculations potential no. pair Re (A) D (ev) a (A-I) &,(A) 4.362 0.024214 Xe-Xe 1.467 7.0 3.727 1 Xe-0 0.007 570 7.0 2 Xe-0 3.727 0.075 700 7.0 3 Xe-0 2.719 0.647 3 7.0 ~~~~~~
atom is not a closed-shell system. In view of this uncertainty, we have carried out several calculations using three LJ( 12,6) potentials that satisfy the foregoing expectations. The LJ( 12,6) potentials utilized in the diffusion calculations are given in Table 111. Potential 1 has an attractive well identical to that obtained from eq 6 which we used for the 03 relaxation calculations. The equilibrium distance, however, is 0.178 A greater so that this potential is slightly more repulsive at shorter Xe-O separations. We assume that this potential represents a limiting case corresponding to the smallest attractive interaction we would expect to have for the O/Xe pair. Potential 2 with a well depth 10 times that of potential 1 represents the limiting case corresponding to the largest attractive interaction we would expect to have for the XeO triplet state. This view is supported by the fact that the Xe/Xe well depth of 0.024 214 eV is only 32% that for potential 2. Since Xe has a significantly larger polarizability than oxygen, we would expect the Xe-0 triplet well depth to be less than that for Xe/Xe. The well depth of potential 2 is also 2.6 times thedeepest well obtained for the triplet state in any of the pseudopotential calculations. We have carried one set of calculations using potential 3. This potential has a well depth that matches that predicted by the MP2 results for the bound XeO singlet state using the WadtHays19 pseudopotential for the xenon core. The results obtained with this potential serve to illustrate the type of diffusion behavior to be expected if the oxygen atoms were present in the excited singlet state. IV. Methods Initial OJ/Matrix Structure. Initially, the rare-gas atoms are arranged in the (5 X 5 X 5) fcc lattice configuration. Ozone is assumed to be contained interstitially within the innermost unit cell of the crystal with its center of mass coinciding with the
Molecular Dynamics of O3in Ar and Xe Matrices
The Journal of Physical Chemistry, Vol. 97,No. 47, 1993 12137
center of the unit cell. The starting configuration for Oshas the molecule in the x-y plane with its C,, axis aligned with the z axis of the crystal. Averaging over the initial 0, configuration is achieved by executing random rotations of the molecule about the three principal axes of the crystal. That is, we take to the Oscoordinates ( x ’ J ’ , ~ ’ ) to be
(3;)
(1:)
= RX(8.J RY(8,>R,(fl,) yo
(9)
where R,(B,), Ry(By),and R,(B,) are the standard rotation matrices, ( x o , y ~ , z o are ) the starting O3coordinates, and Oa = 2aFa ( a = x , y , or z ) (10) with being a random number whose distribution is uniform on the interval [0,1]. The rare-gas matrix is now allowed to relax about the fixed 03 molecule to the nearest local minimum potential using damped trajectorymethods.6>* In this procedure, the momenta of all atoms are set to zero and the Hamiltonian equations of motion for the lattice atoms are integrated until the potential energy attains a minimum. At this point, the integration is halted, and the momentum components of all atoms are once again set to zero. This procedure is repeated until the system converges to a local minimum. In practice, the procedure is continued until convergence to at least seven significant digits is obtained. The relaxation produces a substantial expansion of the local matrix environment about ozone. The difference in potential between the ozonerigid lattice structure and the ozone-relaxed lattice structure is about -36.9 eV. Energy-Transfer Dynamics. The initial state is prepared by inserting zero-point energy into each of the threevibrational modes of ozone. Subsequently, the desired excitation energy is placed in a selected internal mode. By investigating the differences in the energy-transfer dynamics as a function of the internal mode excited, the presence or absence of mode-specific effects may be determined. A lattice temperature of 12 K is simulated by setting each initial momentum component of the lattice atoms in the P and Q zones to f[2mkbTl1/*,where m is the lattice atom mass and kb is the Boltzmann constant. In each case, the sign of the momentum component is chosen randomly. The behavior of the ozone-lattice system is calculated through the solution of Hamilton’s equations of motion for ozone and the 364 P zone and Q zone atoms moving in the potential field described by eq 1. The trajectories are integrated using a fourthorder Runge-Kutta procedure. Typically, trajectories are followed for 2.55 ps. The results from each trajectory yield the total ozone energy, E,,, as a function of time, where the ozone energy is defined to be (1 1) Eo, = To,+ voz where T,,is the Oskinetic energy and Vo,is the Murrel-Farantos18 potential. The average energy-transfer dynamics are obtained by the calculation of Eoz(t)for an ensemble of 100 trajectories at each internal Osenergy partitioning investigated. Classical Variational Transition-State Calculations. The oxygen atom is initially placed at the most stable absorption site which is often the center of the innermost unit cell of the (5 X 5 X 5) matrix. The objective of the study is to compute the diffusion rate to an adjacent adsorption site, which, in a fixed matrix, would be located at the midpoint of an edge of the center unit cell. Because of the symmetry of the fcc matrix, this point is theoretically equivalent to the initial site before relaxation of the lattice, although the finite size of the lattice model introduces a small deviation from the expected symmetry. The oxygen atom diffusion rate is calculated from the jump frequency or flux across a theoretical dividing surface separating
the two adsorption sites. If the system behaves statistically, this flux must be an upper limit to the actual diffusion rate since all diffusion events involve crossing of the dividing surface but not all crossings result in diffusion. Consequently, we seek the dividing surface that minimizes the flux. Several spherical and cubical surfaces have been examined to effect this minimization. The total energy of the O/matrix system is given by N
i= I
+
+ + v,’ (12)
+
[pdZ pyOz ~ , , , 7 / 2 m ~ V,
where Pqi ( q = x, y , z ) represents the momentum of lattice atom i in the q direction. The subscript runs over the lattice atoms, and “0” denotes the oxygen atom. It is assumed that the system behaves statistically with a canonical distribution of energies. Under these conditions, the jump frequency is proportional to the probability of the oxygen atom being on the dividing surface and to the velocity of the atom perpendicular to that surface. Thus, the flux, F( T), can be expressed by the sum of all such products averaged over the phase space of the system divided by the total available phase-space volume. That is,
where the delta function is unity when on the dividing surface and zero otherwise. The configuration space integrals cover the space corresponding to reactant conformations. The integrations over momenta in eq 13 can be done analytically for spherical or cubic dividing surfaces. Such integration yields 3N
-
i= 1
Jq 1-1
where ( u ) represents the average velocity of the oxygen atom. Since the potential being employed is separable into a lattice potential plus an oxygen-lattice interaction, eq 14 may be written in the form
The complexity of the potential precludes analytical evaluation of eq 15. We therefore utilize Metropolis sampling methods to execute the required integrations. This procedure requires that the dividing surface be replaced with a dividing “slab” of width Aw. If Aw is sufficiently small that the integrand of eq 15 is constant across the width, eq 15 becomes 3N
Jqe-’mIkT e-‘1IkT b( Aw) n d q i i= I
F(T) = [(u)/Awl
7N
(16)
whereb(Aw) isunityifaconfigurationpoint lieswithin thedividing slab and zero otherwise. In principle, eq 16 may be evaluated
12138 The Journal of Physical Chemistry, Vol. 97, No. 47, 1993
Ford et al.
using a random set of M points in the multidimensional configuration space of the system. For such a randomly selected set of points, the Monte Carlo approximation for eq 16 is M
z[e-Vm/kTe-fi/kT
2
ro Ii
I= 1
Although eq 17 yields the flux across the dividing slab, its convergence rate will be extremely slow if totally random points are selected for all atoms since virtually all points selected will correspond to highly improbable configurations. The situation may be improved by selecting the points from a Markov walk weighted by the canonical distribution function exp[-VT/kT]. In this case, the flux will given by A4
The convergence rate of eq 18 will be significantly greater than that for eq 17. However, convergence will still be very slow due to the infrequency of sampling in the regions of high potential. A more satisfactory convergence rate may be obtained by using a Markov walk weighted by the canonical distribution function for the lattice alone, exp[-V,,,/kT]. For such a selection method,
where the terms under the summations are evaluated after every Markov step on the lattice and oxygen atoms and after every damped trajectory cycle. In eq 21, Ar is the width of the dividing spherical slab and MO is the oxygen atom mass. The factor of 2 corrects for surface crossings in the wrong direction. The factor of 12 removes the degeneracy in the calculation which is present since a spherical dividing surface counts jumps to 12 equivalent diffusion sites. In practice, the Markov walk is executed by moving n randomly selected lattice atoms (7 In I20) and the oxygen atom in each step of the walk. The lattice atoms are moved according to q y = q,"ld
+ &Aq
( i = 1, 2, 3, ..., n)
(22) where qiold and qpcware the old and new x , y, and z coordinates of lattice atom i, respectively, and Aq is the Markov step size. The ti are random numbers selected from a uniform distribution on the interval [O,l]. For the oxygen atom,
M
Ax, = AQ sin 8, cos 9,
(23)
Avo = AQ sin 8, sin 9,
(24)
x [ e - v l / k T 6(Aw)]i
F(7') = OS[(u)/Aw]
A4
(19)
where the factor of 0.5 is included to correct for entries into the dividing volume from the wrong direction. We have previously used eq 19 to compute silicon and hydrogen atom diffusion rates on Si( 111) and Si(ll1)-(7X7) surfaces.21 Although eq 19 can be used to obtain diffusion rates on surfaces and in matrices, its convergence rate is still very slow. Typically, millions of Markov steps are required for surface diffusion. For matrices at cryogenic temperatures, convergence is even slower. We therefore introduce here a new method based on the assumption that the major contribution to F ( T ) in eq 16 arises from configurations in the neighborhood of the minimum-energy pathway for the diffusion process. These configurations may be conveniently located and sampled by using a combination of canonical Markov moves on the lattice atoms and totally random moves on the embedded oxygen atom followed by a series of damped trajectory cycles in which the lattice is allowed to relax toward its minimum-energy configuration in the field of a stationary oxygen atom.6-8 A Monte Carlo approximation, such as eq 17, is then expected to converge a t a much greater rate since all of the sampling is being done in statistically important regions of configuration space near the minimum-energy path. For spherical dividing surfaces, it is convenient to write eq 16 in the form
F ( r ) = [(u)/Awl X 3N
sqe-vm/kTe-v1/kT6(Aw)ndqi ro2dr, sin 8, d8, d9, i- 1
where the oxygen atom coordinates are separated and expressed in a spherical system. The Monte Carlo approximation for eq
where
The values of n and Aq are adjusted to produce a near-unit ratio between accepted and rejected moves. In most cases, n = 7 and Aq is 0.0866 A for both the lattice atoms and oxygen. The width of the dividing slab, Ar, is chosen to be equal to the maximum step size to ensure that the oxygen atom cannot traverse the dividing slab without entering its volume a t least once. Subsequent to the Markov step described above, K damped trajectory cycles are executed holding the oxygen atom stationary.6J In this procedure, the kinetic energy of each lattice atom is set to zero and the classical Hamiltonian equations of motion for the latticeatoms are integrated until the total potential energy attains a minimum. This is defined to be one damped trajectory cycle. Subsequent cycles are executed by repeating the above procedure starting with the lattice configuration achieved in the previous cycle. The actual execution of the above procedure can be greatly facilitated by making use of the fact that the system potential is independent of mass. Consequently, Hamilton's equations may beintegratedwith themassofallatomssetto 1.Oamu. Inaddition, we may employ a very large integration step size since we need not be concerned with conservation of energy. The use of these two techniques significantly reduces the computational time required to achieve convergence of the integrals in eq 20. A partial minimization of F( Z') is carried out by computing the flux through a set of spherical dividing slabs with radii of (R = 0.05jd) f o r i = 1, 2, 3,4, ..., 16, where d is the total diffusion
The Journal of Physical Chemistry, Vol. 97, No. 47, 1993 12139
Molecular Dynamics of O3 in Ar and Xe Matrices
-0.6-
0.5t A
-0.8 W
3 -1.2,
-
1
0.2 -
0
-1.4 (C)
0
-
50
100 150 Time Units
200
5
250
Figure 1. Total internal 03 energy as a function of time (1 time unit = 0.01019 ps). Each curve is the ensemble average of 100 trajectories. In
each case, 03 contains zero-point energy plus sufficientexcess excitation energy to excite the mode indicated to the Y = 3 vibrational state. (A) Asymmetric stretchexcited. (B) Symmetricstretchexcited. (C) Bending mode excited. distance measured from the initial adsorption site. Dividingslabs in the shape of cubes with the sides perpendicular to the symmetry axes of the crystal was also investigated, but these surfaces resulted in higher fluxes. The initial adsorption site is located by randomly selecting a set of 10 points within a sphere of radius 1 A centered at the geometric center of the innermost unit cell of the 666-atom matrix, which is also included in the set as an 1lth point. The oxygen atom is placed at one of these points, and the entire system is allowed to relax to the nearest potential minimum using a set of 100 damped trajectory cycles. This procedure is repeated for each of the 11 selected sites. The initial adsorption site is taken to be the point that leads to the lowest total potential. Because the ratio of oxygen atoms to the number of adsorption sites is small, the diffusion coefficient can be related to the jump frequency by where f is the fraction of vacant sites cf = 1 here), a is the dimensionality factor, which is three in the present case since diffusion within the matrix is three-dimensional, and k(7') is the jump frequency. The variational transition-state theory method assumes that k( 7') can be accurately replaced with Fmin( 7'), where Fmin(7') is the minimum flux obtained in the variational adjustment of the dividing slab.
V. Results
Ozone/Matrix Energy Transfer Rates. Energy-transfer rates have been computed for seven different initial energy states for O3in Ar matrices and for three different initial states in a Xe matrix. Thesestatescorrespond tooneofthe threeO3vibrational modes excited to either the Y = 3, 6, or 10 quantum state with zero-point energy in each of the other two modes for Ar and to v = 6 for the Xe lattice. Figure 1 shows typical results for the time variation of E,, in an argon lattice. In each case, one internal mode of Osis initially excited to an energy that corresponds to the v = 3 quantum state. Curves A, B, and C show the resulting rate of energy dissipation to the lattice for initial excitation of the asymmetric, symmetric, and O3 bending modes, respectively. In all cases, there is an initially rapid decrease of ozone energy. This rate decreases continuously until the system approaches thermal equilibrium with the lattice. Examination of the data shows that the initial energy transfer is reasonably well described by a simple firstorder rate law ln[E,,] = ln[E,,], - kt Consequently, the energy-transfer dynamics may be conveniently
10 Time Units
15
20
Figure 2. Logarithmic decay plots of the data shown in Figure 1 over the first 0.2038 ps of relaxation. The lines are linear, least-square fits to the data. 1 time unit = 0.010 19 ps. (A) Asymmetric stretch excited. (B) Symmetric stretch excited. (C) Bending mode excited.
TABLE IV Computed Rate Coefficients for 0 3 Vibrational Relaxation matrix mode excited (VI, 6, v.1 k (Ps-1) Ar asymmetricstretch (0,0,3) 1.43 Ar symmetric stretch (3,0,0) 1.48 bend Ar (0,3,0) 5.58 Ar asymmetricstretch (0,0,6) 6.93 Ar symmetricstretch (6,0,0) 2.02 Ar bend (0,6,0) 6.98 Ar
Xe Xe Xe
bend
asymmetric stretch symmetric stretch bend
(0,10,0) (0,0,6)
8.13
(6,0,0) (0,6,0)
1.76 1.45
1.53
described in terms of a rate coefficient, k,obtained by fitting eq 30 to the initial portion of the decay curves shown in Figure 1. Figure 2 shows logarithmic plots of the data given in Figure 1. In each case, the line is a linear least-squares fit to the data. Obviously, eq 30 provides an adequate description of the initial portion of the energy transfer. The negative slopes of these lines yield the respective rate coefficients for energy transfer out of the asymmetric, symmetric, and bending modes of O3when initially excited to the v = 3 state in an argon lattice. At this level of excitation, there is no statistically significant difference between the energy-transfer rates out of the asymmetric and symmetric modes. However, the rate of transfer from the bending mode is substantially larger. The ratio of rate coefficients for the asymmetric stretch and bending modes is nearly a factor of 4. The energy-transfer dynamics for matrix-isolated ozone may therefore be said to exhibit significant mode specificity. We have carried out similar calculations to those shown in Figures 1 and 2 for O3excited to the v = 6 and Y = 10 states in an argon matrix and for Y = 6 in a xenon matrix. In all cases, the resulting ensemble average decay curves are similar to those shown in Figure 1, and eq 30 is found to be an adequate description of the initial stages of the relaxation. The results are summarized in Table IV. Comparison of Y = 3 with the v = 6 results in argon shows that the rate of relaxation to the lattice phonon modes increases with increasing ozone excitation. This effect is attributed to the expected increase in potential coupling between ozone and the lattice due to the increase of the vibrational amplitudes. Mode specificity is enhanced at elevated excitation energies. When O3 isvibrationally excited to the v = 6 state, the asymmetric stretching mode is found to relax 3.4 times faster than the symmetric stretching mode. The nearly equal relaxation rates for the asymmetric and bending modes a t Y = 6 are misleading in that the two modes have different frequencies. Consequently, excitation to the Y = 6 state represents substantially more excitation energy for the stretching mode than for the bend. If the bend is excited to the Y = 10 state, the excitation energy will be nearly
Ford et al.
12140 The Journal of Physical Chemistry, Vol. 97,No. 47, 1993
equal to that for the asymmetric stretch in the u = 6 state. Under conditions of equal excitation energy, the bending mode is again found to relax faster than either of the stretching modes. We attribute the enhanced rate of energy transfer from the Osbending mode to (1) a closer match between the low-frequency bending mode (701 cm-1) and the low-frequency phonon modes of the lattice and (2) the greater amplitude of the bending motion compared to that for the stretching modes at a corresponding energy. We may also increase the coupling between lattice phonon modes and 0 9 by increasing the temperature of the lattice. The higher lattice energy increases the vibrational amplitudes of the lattice which produces a larger 03/lattice interaction. As an example of this effect, we have computed the relaxation rate for O3 in the (6,0,0)state a t an argon lattice temperature of 50 K. The result is k = 3.67 ps-I, which is a factor of 1.8 greater than the relaxation rate at a lattice temperature of 12 K. The same factors responsible for the presence of mode-specific energy-transfer rates also account for the observed differences in transfer rates in Ar and Xe matrices. The data in Table IV show that energy transfer is generally slower in Xe than in Ar. The larger lattice spacing present in the Xe matrix acts in the opposite direction from an increased vibrational amplitude. The larger matrix cage reduces the magnitude of the potential coupling between ozone and the surrounding lattice atoms. In addition, the O/Xe mass ratio is less favorable for energy transfer than the O/Ar ratio. As a result, the energy-transfer rate decreases. The apparent lack of mode specificity in the Xe lattice is again an artifact of the lower energy present in the u = 6 vibrational mode of the bend compared to the stretching modes. At comparable energies, we expect the bending mode to transfer energy more effectively to the lattice for the reasons given above. Since the oxygen/xenon interaction is due predominantly to the dispersion forces between closed-shell systems, we expect the potential given by eqs 4-7 with the parameter values listed in Table 1 to be a reasonable representation of the experimental systems. We have, however, investigated the sensitivity of the energy-transfer rates to the O/Ar potential well depth by carrying out a few calculations with e doubled to 0.010 76 eV. For this interaction potential, the relaxation rate of O3 excited to the (6,0,0) state increases from k = 2.02 to 2.99 ps-I. It appears that the relaxation rate is not an extremely sensitive function of the O/Ar well depth. Oxygen Atom Diffusion Rates. The Markov walk/damped trajectory procedure described in section IV is found to converge at a rate much faster than normally seen in the calculation of diffusion rates using classical variational transition-state theory methods.22 Figure 3 shows a typical result. After 5 X lo4 steps, each with n = 7 and K = 1, convergence to within f 1 0 % has been achieved. In contrast, we were unable to obtain convergence at 40 K using eq 19 even with lo6 Markov moves with n = 10. We are currently exploring the damped trajectory method for hydrogen atom diffusion on a Si( 11 1)-(7X7) surface.23 The potential energy of the system varies significantly as the oxygen atom moves within the unit cell. The variation in system potential for the straight-line diffusion of oxygen from one adsorption site to another in a fixed xenon lattice is shown in Figure 4 for potential 1. As expected, the potential energy along such a path attains a maximum value when the oxygen atom is located at a point midway between the two sites. The barrier for diffusionin thiscaseisabout9.4eV. For moreattractive potentials with larger values of e, the potential barrier to diffusion in a fixed lattice becomes even larger. In an argon lattice, the barrier is also considerably higher due to the smaller unit cell spacing. A diffusion barrier of 9.4 eV or higher essentially precludes the possibility of observable thermal diffusion of oxygen atoms a t temperatures characteristic of matrix experiments. However, when the lattice is permitted to relax in the potential field of the
*A.
0
10
20
30
40
50
Number of Steps11000
F i p e 3. Convergence rate for the dividing surfacelocated at the midpoint of the diffusion path (2.167 A) using potential 1. The abscissa gives the number of steps in the summationsof eq 21 where each step consists of
a Markov move of seven lattice atoms and the oxygen atom followed by one damped trajectory cycle.
Percentags of Diff&ori Eista-ce Figure4. Variation of potential energy for oxygen atom movement along a straight line connecting the two absorption sites in a fixed fcc xenon
lattice using potential 1. The abscissa gives the percent of the total diffusion distance covered.
POt.#l
-/m-
0
2
I 3
Diffusikn Distance (A) Figures. Minimum-energy profiles for oxygen atom diffusion in a xenon lattice: potential 1 (o),potential 2 (e). The plotted points are the minimum crossing potentialsobtained in the Markov/damped trajectory walk. The lines are for enhancementof visual clarity. The abscissa gives the position of the dividing surface along the diffusion path. oxygen atom and the phonon modes of the lattice are permitted tocontribute to thediffusion process, theeffective barrier generally decreases by more than an order of magnitude. We have determined minimum-energy reaction paths by recording the minimum system potential for each of the dividing slabs obtained for oxygen atom crossings observed during the Markovldamped trajectory walk. Typical examples are shown in Figure 5 for potentials 1 and 2. In each case, the barrier maximum occurs in the neighborhood of the dividing slab located at d / 2 ,as expected
Molecular Dynamics of Osin Ar and Xe Matrices
The Journal of Physical Chemistry, Vol. 97, No. 47, 1993 12141
TABLE V Oxygen Atom Diffusion Rates in a Xe Lattice Calculated from Eq 21 Using the Markov/Damped Trajectory Method Described in Section IV potential no. T (K) Ma Ebb (eV) D (cm2/s) 1 1 1 2
3
40 80 3ZC 40 40
5OOOO
41 200
50000 50 000
0.152 0.152 0.152 0.170 0.870
7.90 X 8.02 X 2.48 X 6.95 X 8.42 X 10-l1I
a Number of moves in Markov/damped trajectory walk. b Potential barrier height. Extrapolated resulted. from symmetry considerations and the results shown in Figure 4. For a LJ(12,6) potential with t = 0.007 57 eV, relaxation of the lattice reduces the diffusion barrier from 9.4 to 0.152 eV. The calculated oxygen atom diffusion rates in Xe matrices and the associated barrier heights for the three potentials described in section I11 and Table I11 are given in Table V. The results obtained with potential 1 represent an upper limit to the oxygen atom diffusion rate in a perfect fcc Xe crystal. If the rates computed at 40 and 80 K are fitted to an Arrhenius form, the result is D(7‘) = 0.819 exp[-0.143 eV/kT] cm2/s. Using this result, the extrapolated diffusion rate at 32 K is found to be 2.48 X 10-23 cmz/s. The diffusion rate at 40 K obtained with potential 2 is a factor of 11.4 slower than that predicted by potential 1, and the corresponding diffusion barrier is 11.8% higher. Potential 2 is characterized by an attractive well 10 times as deep as potential 1. As expected, the more strongly bound system exhibits the higher barrier to diffusion and the slower diffusion rate. Potential 3 is characteristic of that expected for singlet oxygen which the pseudopotential calculations indicate forms a bound state with Xe. We would expect diffusion to be extremely slow for such a bound state a t 40 K. The results of the variational transition-state calculations confirm this expectation. At 40 K, the diffusion barrier is 0.870 eV and the associated diffusion rate is 8.42 X cmz/s. In effect, a bound Xe-0 pair does not undergo diffusion at 40 K. Because of the large oxygen mass, tunneling processes would not be expected to make an important contribution to the diffusion rate. However, the classical diffusion rates given by eqs 20 and 21 are very small at low temperatures. Consequently, there exists a possibility that tunneling processes might make a significant contribution to the diffusion rate. We havecomputed an upper limit for the tunnelingcontribution by assuming that the process may be treated as one-dimensional. The tunneling jump frequency between adsorption sites at temperature T,&,(7‘), may be expressed as
where Eb is the potential barrier height, u(E) is the oxygen atom vibration frequency in the adsorption site at energy E, and Tp(7‘) is the tunneling probability at energy E. Thevibration frequency may be computed directly from the effective one-dimensional potential using
where m is the oxygen atom mass and V ( r ) is the effective onedimensional potential. A WKBJ expression is employed to compute Tp(T): (33) where w = [h/2n]-’Jbr=r. [2m(V(r) - E ) ] ’ / 2d r
(34)
with r, and rb being the two points for which V(r,) = V(rb) = E .
TABLE VI: Tunneling and Classical Diffusion Rates for Oxygen Atoms in a Xenon Matrix Using an A ~ l y t i cFit, Eq 35, to the Minimum-Energy Diffusion Path for Potential 1 temp (K) Dla (cm2/s) D t (cm2/s) 32 7.3 x 10-29 2.48 X l0-23c 40 5.1 x 10-24 7.90 X 10-19 80 2.0 x 10-14 8.02 X 10-10 a Tunneling contribution from eqs 31-35 and eq 29. Classical contribution computed from eq 21. Extrapolated result. The one-dimensional effective potential is taken to be the dynamic reaction profile shown in Figure 5 for potential 1. An acceptable fit to the computed points is given by V(r) = Ar2 exp[-B2]
for 0 5 r 5 0 . 5
=A(1-r2)exp[-E(1-r)’]
f o r 0 . 5 5 r 5 1.0 (35)
with A = 1.791 929 eV/dz and B = 4.000 d-2, where d denotes units of diffusion distance. All tunneling calculations were carried out using eq 35. The computed thermal tunneling diffusion coefficients computed from eqs 29 and 31 are given in Table VI where they are compared with the classical diffusion rates obtained from eqs 20 and 21. Clearly, the tunneling contribution to the diffusion rate is negligible at all temperatures. VI. Comparisons to Experiment and Theory RafP has previously computed the energy-transfer rates for relaxation for vibrationally-excited 1,Zdifluoroethane (DFE) into the phonon modes of Ar and Xe lattices. These rates were also found to be well described by eq 30. The average relaxation rate coefficients for DFE are 0.065 and 0.059 ps-l in argon and xenon matrices, respectively. The dynamics of ozone relaxation are similar to those for DFE in that the initial rates are adequately described by a first-order rate law. However, the average relaxation rates for ozone are about 2 orders of magnitude greater than those obtained for DFE. These results may be attributed to three factors. First, the ozone/ lattice interaction potential is larger than the corresponding one for hydrogen. In the case of DFE, the molecule-lattice interaction is dominated by the interaction of the hydrogen atoms with Ar or Xe. Hydrogen is less polarizable than the oxygen atoms of 03. Consequently, the magnitude of the interaction potential is significantlyless. This is reflected by thevaluesof the t parameters for the two systems, [ t O / t H ] A r = [ e O / e H ] X e = 3.42. In addition, the results for DFE7 show that the relaxation rates in Ar are nearly identical to those in Xe, whereas for 03,relaxation in Ar is significantly faster. This difference is due to the relative sizes of the two molecules. DFE is much larger than 03.The dimensions of DFE are comparable to the dimensions of the unit cell spacings in both Ar and Xe lattices. Consequently, DFE is equally “cramped” in both fcc matrices. This leads to a relaxation rate that is nearly the same for both lattices. In contrast, the Os dimensions are significantly less than the unit cell spacings. The larger Xe lattice therefore operates effectively to reduce the energy-transfer rates over those observed in the smaller Ar cage. The second effect involves the extent of the frequency mismatch between the relaxing modes and the phonon modes. The hydrogen atom stretching and bending frequencies in DFE are much further removed from the lattice phonon frequencies than is the O3bending frequency. The closer frequency match facilitates energy transfer. Finally, the O/Xe mass ratio is more favorable for energy transfer than is the H/Xe ratio. The optimum ratio is unity. M0lMxc is 0.122 whereas MH/hfXc is 0.007 63. Krueger and Weitz14 have extracted diffusion coefficients for oxygen atoms in xenon matrices at 32 and 40 K from the results of their recombination rate measurements. In their experiments, oxygen atoms (as N 2 0 ) were combined with xenon in a ratio of
Ford et al.
12142 The Journal of Physical Chemistry, Vol. 97, No. 47, 1993
1:729. The N 2 0 was photolyzed to produce oxygen(ID) atoms within the matrix. Relaxation of the Xe-0 exciplex follows to produce O(3P). The diffusion and subsequent reaction of oxygen atoms to form molecular oxygen were observed to occur based on the measured time variation of the O(3P) concentration obtained from the photodetection experiments. The oxygen atom decay is not well described by a second-order rate law. This indicates that a distribution of diffusion coefficients exists.24By assuming that the distribution can be adequately represented by two coefficients, Krueger and Weitz14 obtain values of 1.4 X 10-15 and 5.4 X 10-18 cm2/s at 32 K and 7.3 X 10-Is and 2.0 X 10-17 cm*/s at 40 K. Two possible interpretations of these coefficients have been suggested.14 First, the faster of the rates might be due to defects and inhomogeneities in the matrix which allow easier diffusion of some oxygen atoms while the smaller coefficient is representative of diffusion through an fcc lattice without defects. Second, the faster rate could result from recombination of oxygen atoms that are initially in close proximity, leaving the slower rate to represent atoms having to diffuse through more of the matrix before reacting. The classical variational transition-state theory calculations using potential 1 give diffusion coefficients at 32 and 40 K of 2.48 x 10-23 and 7.90 X cm2/s, respectively. At 32 K, the ratio between the two values extracted by Krueger and Weitz from their recombination measurements to those obtained from the transition-state calculations in a perfect xenon fcc lattice are 5.65 X 107 and 2.18 X 10s for the fast and slow rates, respectively. At 40 K, the corresponding results are 9.24 X lo3and 25.3. Thus, the measured ratesI4 are orders of magnitude faster than those predicted by potential 1, which we believe should yield an upper limit to the 3P oxygen atom diffusion rate. It is clear that neither of the diffusion rates obtained by Krueger and Weitz14 corresponds to diffusion through a perfect fcc xenon lattice. At 32 K, the experimental rates are between 5 and 7 orders of magnitude faster than those we compute for a perfect fcc lattice. At 40 K, the differences are between 1 and 4 orders of magnitude. The results given in Table VI show that tunneling contributions do not decrease these differences. We conclude that essentially all of the diffusion occurring in the experimental matrices is taking place along defect pathways related to the imperfections in the structure of the vapor-deposited lattices. Recent theoretical studies by RafP have shown that such defects are expected to be extensive with approximately three vacancies about each xenon lattice atom. The present results show that the experimental matrices must contain numerous imperfections, possibly even more than indicated by the Monte Carlo simulations.8 A comparison of the calculated and experimental activation energies for diffusion leads to the same conclusions. The transition-state calculations using potential 1 give an activation energy of 3.30 kcal/mol. The activation energies obtained using potentials 2 and 3 are both larger. In contrast, the temperature dependence observed by Krueger and Weitz14 yields 0.41 5 E, 5 0 . 5 3 kcal/mol. Such a low activation energy suggests that diffusion must be occurring primarily along the defect pathways present in the lattice. The transition-state theory calculations indicate that diffusion in a fcc crystal is associated with a barrier nearly an order of magnitude larger. Recent measurements by LaBrake and we it^^^ support the above conclusions. If most of the diffusion is taking place along the defects present in the lattice, we would expect to see a strong dependence of the measured results upon the procedures used to prepare the matrices. In particular, we would expect to observe a decrease in the diffusion rate as the deposition temperature for the matrix is raised since higher deposition temperatures would be expected to decrease the extent of lattice defects. LaBrake and WeitzZ5have measured hydrogen atom mobilities in xenon matrices deposited at 10 and 28 K. The results show that H
+
H recombination to form H2 occurs 5 times more rapidly in the lattice deposited at 10 K. In view of the above comparisons, we suggest that it is incorrect to conclude that the large rate coefficient extracted by Krueger and Weitzl4 from their oxygen atom recombination data represents diffusion along defects and inhomogeneities in the matrix, while the small coefficient corresponds to diffusion in a bulk fcc matrix environment. It appears that the experimental situation corresponds to one in which there is a distribution of rate coefficients present that is characteristic of the distribution of defect types in t h e l a t t i ~ e .The"s1ow" ~~ and "fast"va1ues extracted by Krueger and Weitzl4 are best viewed as representing average values over the low and high ends of this distribution, respectively.
VII. Summary The molecular dynamics of vibrationally-excited ozone isolated in Ar and Xe matrices at 12 K have been investigated using classical trajectory methods and a 666-atom model lattice with a potential energy surface that consists of the sum of pairwise lattice and lattice-ozone interactions and a three-body ozone potential. Thermal diffusion rates of oxygen atoms in a fcc xenon lattice have been computed using a classicalvariational transitionstate theory method which employs a Markov walk/damped trajectory procedure to effect convergence. Three different pairwise potentials are investigated. Average vibrational relaxation rates are determined from the ensemble results of batches of 100 trajectories at each initial energy partitioning investigated. The relaxation rates over the first 0.1-0.2 ps can be accurately characterized by a first-order rate law. Rate coefficients calculated from linear least-squares fitting of first-order decay plots are found to increase as the excess excitation energy initially present in O3increases. This effect is attributed to an increase in the ozone-lattice coupling produced by the increased vibrational amplitudes at higher levels of excitation. The variation of the rate coefficients with different initial partitioning of the excitation energy among the three vibrational modes of ozone shows the presence of pronounced mode specificity. Energy in the O3 bending mode is transferred to the phonon modes of the lattice at a rate significantly faster than for a similar level of excitation in the two stretching modes. It is suggested that this effect is primarily due to (1) the closer frequency match between the bending mode and those for the lattice phonon modes and (2) the greater vibrational amplitude of the bendingvibration compared to those for the stretching motions. At low levels of excitation, the relaxation rates of the asymmetric and symmetric stretching modes are nearly the same. However, as the excitation increases, the rate for the asymmetric mode becomes greater than that for the symmetric stretch. Ozone relaxation rates in Xe matrices are found to be less than those in Ar. This effect is primarily the result of the larger lattice dimensions of the Xe matrix which reduces the coupling between ozoneand thelatticeatoms. It isalso found that theOjrelaxation rates increase as the lattice temperature is raised. This effect is interpreted to be a result of an increased coupling to the phonon modes at higher lattice temperatures. The average relaxation rates for ozone are about 2 orders of magnitude greater than those previously reported for DFE.7 This is due to a combination of three effects: (1) the larger ozone/ lattice interaction potential, (2) a closer match between the frequencies of the03 bending mode and the lattice phonon modes, and (3) the more favorable oxygen/xenon mass ratio. Ozone relaxation in Ar is found to be significantly faster than in Xe. In contrast, the relaxation rates for DFE are nearly the same in Ar and Xe matrices.' It is argued that this difference is a consequence of the difference in the size of Os and DFE. The variational transition-state calculations show that thermal diffusion rates of oxygen atoms in fccxenon crystals are extremely
Molecular Dynamics of Osin Ar and Xe Matrices slow with an activation energy greater than 3.3 kcal/mol. The corresponding diffusion rates in Ar are found to be even slower. Comparison of the results with measured diffusion coefficientd4 indicates that nearly all of the experimentally observed diffusion is occurring along lattice defects, grain boundaries, vacancies, and other lattice imperfections. The two diffusion coefficients reported by Krueger and Weitz14 for the oxygen atom/Xe system are therefore interpreted as representing average values over the low and high ends of a diffusion coefficient distribution that is characteristicof the number and natureofthelattice imperfections present in the experimental matrices.
Acknowledgment. The authors are indebted to Professor Eric Weitz for providing preprints of his mobility measurements with Dr. D. LaBrake prior to publication and for cogent discussions concerning the general problem of atomic mobility in rare-gas matrices. We are likewise indebted to Prof. R. B. Gerber, who provided numerous cogent suggestions and comments. G.J.M. thanks the National Center for Supercomputing Applications for Grant CHE890003N and for computer time on the Cray YMP, Cornel1 Theory Center for computer time on their RS/ 6000 cluster, and Oklahoma State University for generous use of their IBM 3090. We are pleased to acknowledge financial support from theNationa1 Science Foundation under Grant CHE921 1925. Wearealsoindebted totheOklahomaStateUniversity Center for Energy Research (UCER) for financial support. References and Notes (1) Hauge, R. H.;Gransden, S.;Wang, J.; Margrave, J. L. Ber. BunsenGes. Phys. Chem. 1978,82,104;J . Am. Chem. SOC.1979,101,6950.
The Journal of Physical Chemistry, Vol. 97, No. 47, 1993 12143 (2) Frei, H.; Fredin, L.; Pimentel, G. C. J . Chem. Phys. 1981,74,397. (3) Knudsen, A. K.; Pimentel, G. C. J . Chem. Phys. 1983,78,6780. (4) Frei, H.;Pimentel, G. C. J . Chem. Phys. 1983,78, 3698. (5) Frei, H.J . Chem. Phys. 1983,79,748. (6) Raff, L. M.J . Chem. Phys. 1990, 93,3160. (7) Raff, L. M.J . Chem. Phys. 1991, 95,8901. (8) Raff, L. M. J. Chem. Phys. 1992,97,7459. (9) Jacox, M. E. Chem. Phys. 1979.42, 133. (10) Jacox, M. E. Chem. Phys. 1981.59, 199. (1 1) Feld, J.; Kunttu, H.; Apkarian. V. A. J. Chem. Phys. 1990,93,1009; 1990,92,4826;Chem. Phys. Letr. 1990, 171,425. (12) Alimi, R.; Gerber, R. B.; Apkarian, V. A. J . Chem. Phys. 1990,92, 3551. (13) Lawrence, W. G.; Apkarian, V. A. J. Chem. Phys. 1992,97,6199. (14) Krueger, H.;Weitz, E. J . Chem. Phys. 1992,96,2846. (15) Riley, M.E.;Coltrin, M. E.; Diestler, D. J. J . Chem. Phys. 1988,88, 5934. (16) LeRoy, R. J. J . Chem. Phys. 1972.57,573. (17) Barker, J. A,; Watts, R. 0;; Lee, J. K.; Schafer, T. P.; Lee, Y.T. J . Chem. Phys. 1974.61, 3081. (18) Murrell, J. N.:Farantos, S.Mol. Phvs. 1977,34, 1185. (19) Wadt, P. J.; Hay, W. R. J. Chem. Phys. 1985,82, 270,284,299. (20) Frisch, M. J.; Trucks, G. W.; Head-Gordon, M.;Gill, P. M.W.; Wong, M. W.; Foresman, J. 8.; Johnson, B. G.; Schlegel. H. B.; Robb, M. A.; Replogle, E. S.;Gomperts, R.; Andres, J. L.; Raghavachari, K.; Binklley, J.S.;Gonzalez,C.;Martin,R.L.;Fox,D. J.;Defrees,D. J.;Baker,J.;Stewart, J. J. P.; Pople, J. A. Gaussian 92,Revision A ; Gaussian, Inc.: Pittsburgh, PA, 1992. (21) LaJohn, L. A.; Christiansen, P. A.; Ross, R. B.; Atashroo, T.; Ermler, W. C. J. Chem. Phys. 1987,87,2812. (22)See, for example: (a) Rice, B. M.; Raff, L. M.; Thompson, D. L. Surf. Sci. 1988,198,360. (b) Agrawal, P. M.;Thompson, D. L.; Raff, L. M. J. Chem. Phys. 1989,91,6463. (23)Sorescu, D.; Thompson, D. L.; Raff, L. M. Unpublished work, to be submitted. (24) Spath, B. W.; Raff, L. M. J . Phys. Chem. 1992,96,2179. (25) LaBrake, D.;Weitz, E. To be published and private communication.
