J . Phys. Chem. 1986, 90, 990-992
990
Comparison of Unimolecular Rates Calculated Classically and Quantally from the Same Potential Surface L. Kenneth Cohen and Huw 0. Pritchard* Centre f o r Research in Experimental Space Science, York University, Downsview, Ontario, Canada M3J 1 P3 (Received: October 21, 1985)
Classical trajectory calculations of the specific rate constant k ( E ) for the dissociation of N 2 0 into N2 and 0 are compared with the quantal result previously derived from the same potential energy model. Above the classical threshold, the two results agree extremely well; however, in the quantal calculation, tunneling states just below the classical threshold contribute about 15% to the total rate.
A state-to-state calculation of the rate of dissociation of N 2 0 into N, and 0 has been described by Yau and Pritchard.’%*In this calculation, it was assumed that the linear vibrational motion of the N-N-O molecule could be represented as a simple harmonic motion of the two N atoms with respect to each other, coupled with a Morse oscillator motion of the 0 atom with respect to the N, pair. Further, it was assumed that the matrix elements between reactant and product states were independent of the rotational and bending quantum numbers of the reactant state. A full specification of the linear potentials, which were derived from spectroscopic data, is given in the original publication.’*3 We now report the results of a set of trajectory calculations for this reaction, using the same potential model as was used in the quantal calculations earlier. Computational Details Our original intention was to extract the nonadiabatic coupling probabilities from the original (Airy function) quantal calculations and to graft them on to the trajectory results; however, the appearance of a paper by Lefebvre4 showing good consistency between various quantal and semiclassical treatments of the curve crossing process led us to attempt a fuller treatment of the problem. We follow the development of MillerS and adopt his notation: the electronic transition probability is defined, eq 5.1 3 as
where, eq 5.39
a=--
h
(2
J::dr ( [ E- W l ( r ) ] 1 /-2 [ E - W z ( r ) ] l / 2=]
with a = V l 2 ( r o ) / [ E M ( r- 0E)0 ] ,v = ( 2 [ E M ( r o-) E 0 ] / ~ ) 1and /2, y are Eo are respectively the reduced mass for the N2-0 system and the energy of the crossing point; the subscript M refers to Morse oscillator quantities. In the limit of weak coupling, a 0 and eq 1 reduces to the elementary form (5.19) given by Miller. We compute trajectories for the coupled motion of the harmonic oscillator and the Morse oscillator at chosen values of the initial energy by standard numerical procedures.6 For each crossing of ro in the outward direction, the value of 6 was calculated by eq 1: the probability of a nonadiabatic crossing, and therefore the fraction of the flux remaining after the ith crossing at time t, is N , / N , _ , = Consequently, after n crossings, NJN, = exp(-2C:=,6J defines the rate constant kE for the decay of the coupled oscillator system at energy E . The kE value and its error were determined from the slope of an analytical fit’ of In [N,,”,] v s . t j , i = l , 2,..., n. It does not matter how the initial energy is apportioned between the two oscillators-the motion appears to become “statistical” s, in the sense that plots of E M vs. time within about 2 X look very much like Figures 1 and 6 of ref 9; in selected cases, we confirmed that kE itself was independent of the choice of initial energy distribution between the two oscillators, to within the quoted errors. For simplicity, all trajectories were started at the inner turning point of the Morse oscillator. The results of these calculations are shown in Table I for a selection of energies, beginning just above the classical threshold; notice that the calculated values of kE are not entirely monotonic with increasing -+
( 6 ) Cohen, L. K.; Pritchard, H. 0. Can. J . Chem. 1985, 63, 2374. (7) We require the slope of - xJ-126, vs. t , where i is the crossing number and n is the total number of crossings: the mean dissociation time, T = 1 /kE, and its standard error of the mean, a,, are8
”
O
n
(=I
i=i
,=I
n Z a i t , - Za, Z t ,
W , ( r ) l / l [ E- W,(r)11’2+ [ E - f+7*(r)1”21
T =
For the Landau-Zener model (eq 5.17a and 5.17b) where W , , = E , f Vl,(ro)(l- x2)1/2
and AW = 2 V I 2 ( r o ) ( -l x2)Ii2
we obtain, finally 6=
where a, = 2C;.16,, time, viz.
and
to is
a computed unbiased estimate of the starting n
to
(I) (2) (3) (4) (5)
Yau, A. W.; Pritchard, H. 0. Can J . Chem. 1979, 57, 1731. Pritchard, H . 0. J. Phys. Chem. 1985, 89, 3970. With the given potential, the crossing point occurs at ro = 1.905 A. Lefebvre, R. Chem. Phys. Lett. 1984, 107, 360. Miller, W. H. Adu. Chem. Phys. 1974, 25, 151.
0022-3654/86/2090-0990$01.50/0
=
n
. . . . n
n
n
. . . . n
nZaa,’ - ( E a , ) ’ I=I
,=I
(8) Topping, J. Errors of Observation and Their Treatment, Chapman
Hall: London, 1972; pp 97-100. (9) Sumpter, B. G.; Thompson, D. L. J. Chem. Phys. 1985, 82, 4557
0 1986 American Chemical Society
Letters
The Journal of Physical Chemistry, Vol. 90, No. 6, 1986 991
-
TABLE I: Rate Constants for Decay of N 2 0 N2 + 0, Modeled as a Pair of Linear Coupled Oscillators, as a Function of Energy E" E , kcalmol-I n, n 10-1%,, s-1 30 535 5 3.75 f 0.60 65.30 66.08 67.21 68.57 69.67 70.76 72.94 75.13 77.31 79.93 82.55 84.74 86.93 89.16 91.39 94.41 97.03 99.19 101.83 104.76 106.81
200 000 181 289 106 909 200 000 200 000 200 000 200 000 200 000 74083 200 000 138 521 15994 18 123 7 058 3 797 2 729 3228 4 185 3 292 1884
.
81 86 86 166 174 234 264 323 125 367 268 31 40 16 11 9 9 12 7 4
6.19 f 0.13 5.83 f 0.08 7.54 f 0.10 8.19 f 0.08 7.53 f 0.08 9.83 f 0.03 10.69 f 0.05 10.85 f 0.05 11.69 f 0.09 10.95 f 0.02 11.64 f 0.05 10.85 f 0.36 17.38 f 0.59 7.06 f 0.26 15.09 f 0.76 17.66 f 1.45 11.63 f 0.73 22.29 f 1.22 16.03 f 2.48 5.30 f 0.78
"Here, E is measured with respect to the zero-point energy of the s); the inteoscillator pair. n, is the number of time steps (61 = gration was terminated either at 2 X lo5&, when the integration error became too large or when (at high energies) dissociation occurred on the wrong surface. n is the number of outward crossings of ro during the trajectory.
energy, and they retain some of the oscillatory characteristics of the original quantal solution, cf. Figure 2 of ref 1 .
Energy
( k c a 1 . rno 1 -'I
Figure 1. Comparison of full quantal calculation (dots) with classical trajectory calculation (solid line). Upper curve: logarithm of p ( E ) k ( E ) 6E; lower curve p(E) k ( E ) 6E. In both cases, 6E = 0.05 kcalmol-I.
Since there are no stationary states of the Morse oscillator for E,,, 2 0, the density of states of the classical Morse oscillator is simply
Rate Constants
The simplest comparison we can make with the quantal calculations is for the value of the infinite pressure rate constant
where Q(7') = S t d E p ( E ) and p ( E ) are the partition function and density of states of the coupled oscillator system of total energy E = EM + E H + De, respectively, the subscript H denoting harmonic oscillator quantities. Given any potential V ( x ) , without rotation, E(p,x) = p 2 / 2 r V ( x ) ;hence, with (for clarity) the symbol t for energy as an integration variable, the partition function is
+
where ~ ( t is) the period of oscillation. There is no explicit expression for the partition function for a Morse oscillator Q M ( 7'): substituting VM(x) = De(1 - e-nx)2 - De into eq 3, we find that
whence
The same procedure for the harmonic oscillator, with V&) = ~ ( 2 a v ) ~ x gives ~ / 2 ,QH(T) = R T / h v as expected. Notice that, for convenience, the zero of energy for the Morse oscillator is taken to be its dissociation limit, whereas other energies are measured from the bottom of the well.
The combined density of states p ( E ) = PM(EM)PH(EH)of both oscillators can then be substituted into eq 2 to give
k, =
(1
- e-Emax/RqRT
s,"'
deM (De - tM)-l/2e-e~/RT
where the first term in the denominator appears due to truncation of the system at E,,,. The lower limit tL in the tH integration is 0 if t 5 De, or e - De if t > De, and Q M ( T ) has now been rewritten with the origin of its energy at the bottom of the well; thus ( 2 ) becomes
The energy range was divided into segments with boundaries bisecting the tabulated entries of Table I, and with k E assumed to be constant across each segment, the rate constant at 2000 K was calculatedlo to be (9.2 f 0.2) X lo4 s-I, where the uncertainty is derived from recalculation using all upper bounds or all lower bounds on kE. The full quantal result at 2000 K is 9.1 X lo4 s-], but if three tunneling states below the classical threshold are excluded, the quantally calculated rate drops to 7.9 X lo4 s-'; this relatively large change comes about because there is one state, about 1.2 kcal.mol-' below the classical threshold, with a decay rate constant of 1.8 X 10" s-]. A clearer picture of the degree of agreement between the classical and the quantal calculations can be gained by comparison (10) Notice that because k E is assumed to be independent of the bending and rotational energy of the molecule, the value of k, is the same whether or not those degrees of freedom are included in the calculation of p ( E ) , in the limit as E,,, -.
-
992
The Journal of Physical Chemistry, Vol. 90, No. 6, 1986
Letters
k(E)e-E/RT 6E. Figure 2 shows the more legitimate comparison, with the three tunneling states excluded; at high energies, where the correspondence principle might be expected to apply, the agreement is remarkable.I2 Looked at from one point of view, the spacings between the Morse oscillator levels above threshold are between 1 and 5 kcal.mol-', whereas R T at 2000 K is 4 kcal.mol-l: thus, R T is not very much greater than the level spacing; on the other hand, the first quantum state above the classical threshold is u = 24 of the Morse oscillator, whose motion is very classical indeed. Thus, on the assumption that eq 1 gives a good description of the nonadiabatic character of this reaction, the correspondence prinicple appears to hold well in this system at these energies.
60
73
Energy
EO
93
,DO
(kcal. mol -')
Figure 2. S a m e calculations as described in Figure 1, but with three tunneling states below the classical threshold excluded.
of the specific rate functions; this can most easily be done by comparing the quantally calculated* function & with its classical equivalent p ( E ) k ( E ) 6E, where k ( E ) refers to a molecule having energy E distributed among all degrees of freedom, including bending and rotation." The upper curve in Figure 1 shows the comparison of the present trajectory calculation with the original quantal calculation; the agreement is, in fact, so close that a better indication of the differences is revealed by plotting the actual contribution to the rate at each energy, Le. the function p ( E ) ( 1 1) Because the bending quanta are relatively large (1.714 kcal.mol-'), the convolution was performed by successive additions of the function listed in Table I with shifts of ME,,,and degeneracy factors' of (uknd + l), followed by a classical convolution with the rotational states; a grain width of 6E = 0.05 kcalmol-' was assumed, as in the earlier quantal calculation.
Quantal vs. Classical Calculations We have demonstrated that (at least for this reaction) classical trajectory and quantal calculations using the same potential model yield the same rate constant to within the imagined uncertainties of the two approaches. Classical calculations have significant advantages, principally that the computational methods are well-known, that they are tractable for more complicated potentials than used here, and that the uncertainty of any result is well defined. On the other hand, calculation of wave functions for complicated potentials is still the subject of much research and continues to present severe difficulties for high (Le. potentially reactive) energies; moreover, the uncertainties in the rate constant arising from approximations to the wave functions are difficult to specify.I2 However, it is clear that quantal calculations are essential whenever there are tunneling components to the reaction below the classical threshold, and also for those reactions at the small molecule limit where the assumption of ergodic behavior in the quantal calculation leads to an incorrect falloff shape.'
Acknowledgment. This work was supported by the Natural Sciences and Engineering Research Council of Canada. (12) It is not easy to assess the proper uncertainties in the quantal calculations: the wave functions were approximated by an Airy functon (cf. Appendix 1 of ref 1); our aim at that time was to obtain approximations to the wave functions which were good to 3 or 4 decimal places, in the hope that the final matrix elements might be reliable to at least 2 significant figures.
