J. Phys. Chem. 1996, 100, 4071-4083
4071
Reaction Cross Section and Rate Constant Calculations for the D + H2(W)0,1) f HD + H Reaction on Three ab Initio Potential Energy Surfaces. A Quasiclassical Trajectory Study F. J. Aoiz,* L. Ban˜ ares, and T. Dı´ez-Rojo Departamento de Quı´mica Fı´sica, Facultad de Quı´mica, UniVersidad Complutense, 28040 Madrid, Spain
V. J. Herrero Instituto de Estructura de la Materia (CSIC), Serrano 123, 28006 Madrid, Spain
V. Sa´ ez Ra´ banos Departamento de Quı´mica General y Bioquı´mica, ETS Ingenieros de Montes, UniVersidad Polite´ cnica, 28040 Madrid, Spain ReceiVed: September 20, 1995; In Final Form: December 5, 1995X
Reaction cross sections and rate constants for the D + H2(V ) 0-1, j ) 0-7) reaction have been obtained by quasi-classical trajectory (QCT) calculations on the three ab initio potential energy surfaces (PESs) available for this system. A good agreement has been found between the QCT and quantum mechanical (QM) reaction cross sections and rate constants for the D + H2(V ) 0-1, j ) 0) reactions. Thermal rate constants for the D + n-H2(V)0) and D + n-H2(V)1) have been calculated from the excitation functions, over a wide range of temperatures. The comparison with the quantum mechanical (QM) calculations and experimental results shows that, in general, QCT thermal rate constants are smaller than their QM and experimental counterparts, and this can be traced back to a decrease in the classical reactivity in the threshold region with rotational excitation of the reagents. In addition, the analysis of the QCT results provides an explanation for the differences found in thermal rate constants calculated on the three PESs in terms of specific features of each of these potentials.
I. Introduction The rate of the hydrogen atom exchange reaction has always been of fundamental importance in the theory of chemical kinetics.1,2 Experimental investigations of the H + H2 thermal rate constant, k(T), have been carried out since the 1930s,3 coincident in time with crucial theoretical developments like the introduction of the concept of a potential energy surface4-6 (PES) for the description of nuclear motion during reactive encounters, or the transition-state theory7,8 (TST) to calculate rate constants without directly solving the dynamical problem at a microscopic level. The H + H2 system, being the simplest bimolecular reaction, was used as a testing ground for the abovementioned theoretical developments. From this time on, every major advance in the theory has used this reaction as a prototype. The first rate constants obtained from the direct solution of the nuclear equations of motion in full dimensionality were reported by Porter et al.9 The dynamical problem was solved by using classical mechanics and the relatively precise empirical PES of Porter and Karplus10 (hereafter denoted as PK). This was the seminal work of the method of quasiclassical trajectories (QCT), which is now widespread in the field of reaction dynamics. In the same work, a TST calculation on the PK PES was also performed and gave rate constants systematically lower than those from the QCT treatment, especially at the lower temperatures. QCT rate constants on the PK PES for the deuterated isotopic variants of the reaction were reported subsequently.11 Consistent experimental values of k(T) for the H + H2(V)0) f H2 + H reaction and its isotopic analogues were reported in * Corresponding author. X Abstract published in AdVance ACS Abstracts, February 1, 1996.
0022-3654/96/20100-4071$12.00/0
the 1960s and 1970s by Le Roy and co-workers12-18 and by Westenberg and de Haas19 over a wide temperature range, which in the case of the D + H2(V)0) isotopic variant extended from 167 K to about 745 K. These measurements gave smaller values of the H + H2 rate constants and its isotopomers than those predicted by the QCT calculations on the PK PES.9,11 An Arrhenius plot, i.e., a plot of log k(T) vs 1/T, of the measured rate constants showed, at low temperatures, an upward curvature clearly visible for H + H2 and less marked in the reactions involving deuterium atoms. It was conjectured that this curvature was due to tunneling. In 1976, Kuppermann et al.20 and Schatz and Kuppermann21 reported the first accurate, quantum mechanical (QM) cross sections for the H + H2 reaction in two and three dimensions, respectively. From these cross sections, obtained on the PK PES, thermal rate constants were calculated up to a temperature of 600 K. A comparison of these k(T) with the previous QCT ones by Porter et al.9 showed good agreement for temperatures higher than 400 K, being thus larger than the experimental data. For lower temperatures, the QM rate constants were appreciably higher than the QCT ones, deviating even more from the measured values but exhibiting an upward curvature qualitatively similar to the experimental one. For the PK PES, the accord of the TST rate constant calculations9 with experiment was better than that of the QCT9 or accurate QM21 treatments, but this was a fortuitous fact due to inaccuracies in the potential surface. The construction by Liu and Siegbahn22,23 of an accurate ab initio potential surface for H + H2 and the subsequent analytical fit to this surface by Truhlar and Horowitz24 (hereafter LSTH) allowed a sound appraisal of current reaction rate theories. Dynamical calculations of the rate constants of H + H2, D + © 1996 American Chemical Society
4072 J. Phys. Chem., Vol. 100, No. 10, 1996 H2, and H + D2 were performed by Mayne and Toennies25 using the QCT method on the new ab initio PES.22-24 The results for the reaction with H2(V)0) and D2(V)0) were in very good agreement with the experimental data. Only in the case of H + H2, the low-temperature experimental rate constants were significantly larger than the QCT ones. This discrepancy was attributed to tunneling. QCT calculations were also performed for the first two vibrationally excited states of the H2 molecule; in the case of the reaction with H2(V)1), the calculated k(T) were largely at variance with the experiments (see below). From a quantum mechanical point of view, the problem proved to be more difficult than anticipated and for more than a decade after the publication of the LSTH22-24 surface the only QM rate constants directly comparable with experimental data were obtained with cross sections coming from approximate calculations of varying accuracy.26-29 Good to satisfactory agreement with experiment was obtained with these methods. Transition-state theory was also applied to the H + H2 reaction on the LSTH PES and led to lower k(T) values than the ones observed in the experiments.30 The divergence between TST and experiment became notably larger at low temperatures. Several improvements in the conventional TST led to the formulation of generalized transition-state theories,31,32 among which the canonical variational transition-state theory (VTST) of Garrett and Truhlar32 has found widespread application. A major point for making TST theories comparable to experiment is the incorporation of tunneling corrections. In particular, the VTST calculations with tunneling corrections, performed for the H + H2(V)0) system on the LSTH PES, gave rate constants30,33,34 in good agreement with the existing experimental data. A major issue in the discussion of the theoretical results during the 1980s was the possible identification of tunneling effects35,36 in the calculated rate constants. A transmission coefficient κ assumed to account for tunneling was defined, in a rather phenomenological way, by dividing the rate constant calculated with a given approach by a reference rate constant with no tunneling; the reference rate constant was taken from VTST. In fact, this way of defining the transmission factor includes also other dynamical effects such as recrossing, but tunneling is expected to be the most important contribution at low temperatures. Several models, mostly based on onedimensional treatments, were also developed in an attempt to predict the tunneling factors. The good agreement of the QCT rate constants with experimental and QM values was attributed to a lucky cancelation of errors:35 the absence of zero-point energy constraints in the classical dynamics was supposed to compensate the absence of a tunneling contribution. Whereas a reasonably good agreement between experiment and theory has been found in general for the reaction with H2(V)0), a noteworthy discrepancy has long existed in the case of the reaction with H2(V)1) molecules. The rate constants for H + H2(V)1) and D + H2(V)1) measured until 1985 in various laboratories and using different techniques37-44 were substantially larger (sometimes an order of magnitude) than most theoretical rate constants, which at the time had been obtained from QCT calculations,25 from several approximate quantum dynamical methods27,28,45-50 and from variational transition-state theory.34,51 Later experiments by Wolfrum and co-workers,2,53 in which reactants and products were directly monitored by using refined laser techniques, yielded smaller rate constants in the region close to 300 K, in much better agreement with the above-mentioned rate constant theoretical predictions and much more consistent with an accurate quantal calculation of the threshold energy.52
Aoiz et al. In addition to these important new data on the V ) 1 reaction, new electron spin resonance (ESR) rate constant data for H + D2(V)0) have been reported54 in the low-temperature region and also measurements of the k(T) for the H + D2(V)0)55 and D + H2(V)0)56 up to approximately T ) 2000 K have been carried out by using a flash photolysis-shock tube technique. The high-temperature experimental rate constants were larger (by 20-40%) than the predictions of VTST,34 but approximate QM calculation using a reduced dimensionality method57 gave good agreement with the experimental rate constants for both reactions over the whole temperature range;58 the theoretical values are somewhat lower than the mean of the experimental data, but they remain within the dispersion of the experimental points. During the past decade, impressive advances in the quantum mechanical methods have produced remarkable achievements in the studies of the H + H2 reaction. Two new ab initio PESs have been released,59,60 and the pioneering work of Schatz and Kuppermann21 has been extended by different techniques (see ref 61 for a review) to higher collision energies, to internally excited states of the reactants and to other isotopic variants of the reaction, so that by the end of the 1980s the first threedimensional (3D) and fully converged reaction probabilities, cross sections, and rate constants had been calculated on ab initio PESs.62,63 In particular, the use of a theoretical method based on the S matrix version of the Kohn variational principle, allowed the calculation of accurate QM rate constants for the D + H2(V ) 0, 1; j ) 0) reaction63 from fully resolved stateto-state cross sections. Recently, further methodological progresses have allowed the performance of accurate QM calculations of the D + H2(V)0,1) f HD + H thermal rate constants64-66 and have thus enabled the direct comparison of a bulk quantity calculated from first principles with experimental measurements. Park and Light64 performed QM calculations up to 1500 K, but their results diverge gradually from the measurements, so that at 1500 K the theoretical k(T) is about one-half of the experimental one.56 Mielke et al.66 carried out an exhaustive accurate QM determination of rate constants for the ground-state reaction on the three available ab initio PESs,24,59,60 covering the range 167900 K (although extrapolations of the cumulative reaction probability allows estimates of the rate constants up to 1500 K). Their results on two of these PESs were in excellent agreement, for T < 1000 K, with the experimental results. For the reaction with H2(V)1), Auerbach and Miller65 have calculated the excitation functions for j ) 0-3, which allows the determination of thermal rate constants up to 500 K, and their results are compatible with the latest experimental measurements at temperatures close to 300 K.2,53 These works have also revealed that the extension of accurate QM calculations to systems with more than three atoms will be extremely difficult and that, in most practical cases, the recourse to approximate methods will be unavoidable for the calculations of rate constants. In this respect, the recent accurate QM data reinforce the role of the H + H2 reaction as a benchmark for different theoretical approaches, usable for larger systems. Different QM approximate methods have been recently applied to the calculation of rate constants for the D + H2 reaction and the results have been compared to the accurate QM results and to the experimental data. Usually these calculations are based on the J-shifting approximation. Takada et al.67 used the constant centrifugal potential (CCPA) method for the calculation of rate constants for D + H2(V ) 0, j ) 0) and found good agreement with the accurate values of Zhang and Miller.63 Wang and Bowman68 refined the previous reduced
The D + H2(V)0,1) f HD + H Reaction dimensionality calculations on D + H2(V)0).58 The new calculations maintain the same accord with the experimental thermal rate constants up to 2000 K and give cumulative reaction probabilities in better agreement with accurate QM results. Mielke et al.66 used a separable rotation approximation (SRA) to calculate the thermal rate constants for this system up to 1500 K. The agreement with experiment and with accuate QM was found to be good until about 900 K. In the higher temperature region, the SRA rate constants deviate from the measured values so that at 1500 K they are 30-40% lower. For the D + H2(V)1) reaction, Auerbach and Miller65 showed that the thermal rate constants obtained with the J-shifting approximation were in good agreement with the accurate QM values, but noticeable differences between the two approaches appeared in the rate constants for specific rotational states of H2. Although a detailed comparison of the fundamental theory of reaction dynamics with experiment should be performed at the level of state-resolved cross sections, most of the experimental information about the reactivity close to the threshold for reaction still comes from thermal rate constant measurements. On the other hand, these rate constants are indeed very sensitive to the shape of the PES as clearly demonstrated by Mielke et al.,66 who found appreciable differences between the low-temperature accurate QM rate constants calculated on the surface of Boothroyd et al.60 (hereafter BKMP) and those obtained on the LSTH22-24 PES and on the double many-body expansion one of Varandas et al.59 (hereafter DMBE), despite the similarity between the three mentioned PESs. Their calculations also showed that the results on the BKMP PES were in worse agreement with experiment than those obtained with the DMBE or LSTH PESs. In view of the rich information, both experimental and theoretical, gathered during the past years about the kinetics of the H + H2 reaction, a reconsideration of the QCT approach for kinetics studies seems very timely. This method has a high computational efficiency and is very flexible for an extension to higher energies and to larger systems as compared with most QM methods and has performed well even for detailed microscopic observables (state-to-state cross sections). In addition, the method allows a consideration of dynamical effects, like the influence of rotation on reactivity, in terms of easily understandable classical mechanisms. In a previous work, QCT rate constants for the D + H2(V ) 1, j ) 0) reaction were calculated69 on the LSTH PES and the results were in good agreement with experiment53 and accurate QM calculations.63 In this work, the QCT calculations have been extented to the set of D + H2(V ) 0, j ) 0-7) and D + H2(V ) 1, j ) 0-7) reactions on all the available ab initio PESs, using a new methodology that allows the calculation in a quite efficient way of the collision energy dependence of the reaction cross section, i.e., the excitation function. Using the excitation functions so obtained, we have calculated thermal rate constants k(T) for D + n-H2(V)0) from 200 and 1500 K, and D + n-H2(V)1) from 200 to 500 K. The results of reaction cross sections and thermal rate constants are compared with accurate and approximate QM calculations and with experiment. The organization of the paper is as follows. Section II describes the QCT methodology used in this work. In section III, the results of reaction cross sections and rate constants are presented and compared with QM and experimental results. Section IV is dedicated to the discussion of the results of the preceding section in terms of dynamical aspects related with the different topology of the PESs used in this study. Finally, section V contains the conclusions of this work.
J. Phys. Chem., Vol. 100, No. 10, 1996 4073 II. Method The QCT method used in the present work is similar to the one described in previous papers (see, for instance, refs 70 and 71). However, the main objective of this work is to determine the collision energy dependence of the reaction cross section, i.e., the translational excitation functions, σR(ET), for the D + H2(V ) 0, 1; j ) 0-7) f HD + H reaction using the three available ab initio potential energy surfaces: LSTH, DMBE, and BKMP.72 The usual way to do this in QCT calculations consists in running batches of trajectories, each of them at a given collision energy and at an internal quantum state of the reactive diatom. This procedure becomes cumbersome if the calculations are to be carried out on three different PESs and for 16 different initial states of the diatom. Instead, the method followed in the present work consists in running trajectories, for every initial rovibrational state V, j of the H2 reagent, whose collision energy is randomly sampled within the interval [E1, E2], such that the collisional threshold, E0 > E1. The maximum value E2 ) 1.6 eV (chosen to determine rate constants k(T) up to 1500 K) was the same for most of the calculations. Once the value of the collision energy is randomly (uniformly) sampled within ∆E ) E2 - E1, the impact parameter is obtained for every trajectory by
b ) β1/2bmax(ET)
(1)
where β is a random number in the [0,1] interval and the maximum impact parameter, bmax, at a given collision energy, ET, is given by
bmax(ET) ) D(1 - ED/ET)1/2
(2)
The values of the parameters D and ED < E1 < E0 were previously determined by fitting the values of the maximum impact parameters, found at several selected collision energies, ET, to the line-of-the-centers expression of eq 2. The resulting bmax(ET), as given by eq 2, are such that, for the selected ET, there are no reactive trajectories for impact parameters even somewhat smaller than bmax. With this kind of energydependent sampling of the impact parameter, each trajectory is weighted by wi ) bmax2/D2. The integration step size in the trajectories was chosen to be 5 × 10-17 s. This guarantees a conservation of the total energy better than 1 in 105 and 1 in 107 in the total angular momentum. The initial rovibrational energies of the H2 molecule are calculated using a Dunham expansion containing 16 terms73-75 (fifth power in V + 1/2 and third power in j(j + 1)). The classical H2 molecule rotational angular momentum is equated to j(j + 1)p2. Using this methodology, batches of 60 000 trajectories for every H2(V ) 0, 1, j ) 0-7) rovibrational state and within the collision energy range [E1, E2] have been run for each of the three PESs. For the lowest initial rotational states, j ) 0-5, the stratified sampling technique76 was used to decrease the statistical uncertainty near the threshold. To this purpose, extra batches of 6 × 104-105 trajectories were run for H2(V ) 0, j ) 0-5) in a restricted range of collision energies with E2 not higher than 0.6 eV. The whole set of trajectories was then weighted accordingly. In addition, a batch of 106 trajectories was run for the D + H2(V ) 0, j ) 0) to further check if the calculation of sets of 60 000-120 000 trajectories was enough to get reliable results for σR(ET). The overall number of trajectories calculated for this work amounts 5 × 106.
4074 J. Phys. Chem., Vol. 100, No. 10, 1996
Aoiz et al.
The σR(ET) were subsequently calculated by the method of moments expansion in Legendre polynomials (see refs 71, 76, and 77) using the reduced variable:
2ET - E2 - E1 E2 - E1
x)
(3)
The expression of σR(ET), truncated in the Mth term is given by
σR(ET) )
[
1
2R
E2 - E1 2
]
M
+ ∑cnPn(x) n)1
(4)
where R is the Monte Carlo estimate of the integral
R ) 〈σR(ET)〉∆E ) ∫E σR(ET) dET ≈ πbmax2∆ESNR/N (5) E2 1
where N is the total (reactive and nonreactive) number of trajectories, ∆E ) E2 - E1, and the sum of the weights of the reactive trajectories (wi), SNR, is given by NR
SNR ) ∑wi
(6)
i)1
The coefficients of the Legendre expansion, cn, of eq 4 are calculated as the Monte Carlo average of Legendre moments: N
2n + 1 -1 R 2n + 1 SNR ∑wiPn(xi) ) 〈Pn〉 cn ) 2 2 i)1
(7)
To calculate the errors in the reaction cross sections, thresholds, and rate constants, an estimate of the error of the coefficients is needed, and it can be easily found as
[(2n 2+ 1) 〈P 2
γn2 ≡ var[cn] ) SNR-1
]
〉 - cn2
2
n
(8)
The error in σR(ET) is given by the square root of the variance which can be evaluated in terms of γn as
var[σR(ET)] )
( ) N - SNR NSNR
[σR(ET)] + 2
() 2R
∆E
2 M
γn2Pn2(x) ∑ n)1 (9)
The Smirnov-Kolmogorov test comparing the cumulative probability distributions was used to decide when to truncate the series. Significance levels higher than 98% could be achieved using 6-10 Legendre moments, ensuring a very good convergence such that the inclusion of more terms does not produce any significant change. The translational energy threshold, E0, for every initial rovibrational state is determined by finding the roots of eq 4 by the Newton-Raphson method. Special care was paid to the analysis of the threshold in the σR(ET), which remains unaffected, within the statistical uncertainty, when the number of Legendre moments are changed in (2. The error bars correspond to (1 standard deviation, calculated according to eq 9. To check the validity of the above-described methodology, calculations of the translational excitation function for some specific initial V, j have been carried out by running different number of trajectories, as well as by comparison with the reaction cross sections obtained in calculations at definite collision energies.
Figure 1. Total reaction cross section as a function of collision energy for the D + H2(V ) 0, j ) 0) reaction calculated by QCT on the LSTH PES. The lines correspond to the calculations varying the collision energy continuously (see text). Solid and dashed lines represent calculations with 106 and 1.2 × 105 trajectories, respectively. Solid points correspond to calculations at fixed collision energy (mostly taken from ref 71), whose statistical uncertainties are always smaller than the size of the points. The error bars indicate one standard deviation of the calculation obtained with 1.2 × 105 trajectories.
Figure 1 displays the excitation function, σR(ET), for the D + H2(V ) 0, j ) 0) f HD + H reaction on the LSTH PES. The solid line corresponds to the calculation of a total number of 106 trajectories, whereas the dashed line corresponds to the same calculation but using only 120 000 trajectories. Solid points are the σR(ET) from previous work,70,71 calculated by running trajectories at fixed collision energies, which, overall, required a total number of 1.2 × 106 trajectories. As can be seen, the results obtained calculating 120 000 trajectories are coincident, within the error bars, with the results obtained by running 106 trajectories and also with the cross sections at fixed collision energies. The specific thermal rate constant from the V, j initial state of the H2 molecule is given by
k(T;V,j) ) (8kBT/πµ)1/2(kBT)-2∫0 dET ET × ∞
exp(-ET/kBT)σR(ET;V,j) (10) where µ is the reduced mass of the D, H2 system, kB the Boltzmann constant and σR(ET;V,j) represents the translational excitation function for the D + H2(V,j) reaction. In practice, the lower limit in the integral is E0 and the upper limit is E2. The thermal (averaged on initial j) rate constant from an initial vibrational state can be written as ∞
k(T;V) ) ∑pV,j(T) k(T;V,j)
(11)
j)0
where pV,j(T) are the Boltzmann’s statistical weights of the H2 rotational states (including the 3:1 nuclear spin weights), such that ∑jpV,j(T) ) 1. The sum of eq 11 is truncated in j ) 7 (the highest j state calculated in this work). At T ) 1000 K, this accounts for 99.35% of the total population and at 1500 K for 95.45%. As such, the thermal rate constant calculated from V ) 0 are also averaged on initial vibrational states. At 1500 K, the V > 0 population is just 1.96% and the j > 7 population in both V ) 0 and V ) 1 is only 2.5%. Special attention was paid to the estimate of the uncertainties of the rate constants. To this purpose, for every V, j and temperature, a large series of rate constants was calculated using
The D + H2(V)0,1) f HD + H Reaction
J. Phys. Chem., Vol. 100, No. 10, 1996 4075 TABLE 1: Rate Constants k(T) (cm3 s-1) for the D + H2(W ) 0, j ) 0) Reaction as a Function of Temperature Calculated on the LSTH PES (Numbers in Parentheses Represent Powers of 10)
a
Figure 2. Comparison of the accurate QM results of ref 63 and those obtained in the present QCT calculations for the D + H2(V ) 0, j ) 0) reaction on the LSTH PES. Top: reaction cross section as a function of the collision energy; solid points and dashed line: QM results; solid line: QCT results based on calculations with 106 trajectories. The insert displays the cross sections near the threshold. Bottom: Arrhenius plot of the corresponding rate constants; circles and dashed line are the QM results of ref 63; squares and solid line are present QCT results.
eq 10 and the σR(ET) given by eq 4, by randomly sampling the cn coefficients from a Gaussian distribution whose mean value was given by eq 7 and with a standard deviation given by the γn of eq 8. The error of the rate constant was estimated as the standard deviation of the series of calculated k(T). III. Results The comparison between present QCT results and the QM ones of ref 63 for the excitation function of the D + H2(V ) 0, j ) 0) reaction on the LSTH PES is presented in the top panel of Figure 2. Whereas the two excitation functions show a similar behavior near the threshold, for collision energies higher than 0.4 eV, the QM σR(ET) is systematically larger than the one obtained from QCT calculations, and the difference increases with collision energy. The respective rate constants, k(T; V ) 0, j ) 0), calculated from both QM63 and QCT excitation functions, are depicted as a function of the inverse of the temperature in the lower panel of Figure 2 and listed in Table 1. The high threshold for this reaction and the good accordance between the two sets of calculations for ET e 0.4 eV give rise to a remarkably good agreement between QM and QCT k(T;V ) 0, j ) 0) from 300 to 1000 K (the highest temperature for which QM calculations were reported). Only at 200 K, the QM k(T; V ) 0, j ) 0) becomes a factor of ≈2 larger than the QCT one. Although the classical threshold was found to be 0.289 ( 0.005 eV, there
T (K)
QCT
175 200 250 300 350 400 450 500 550 600 650 700 750 800 850 900 950 1000 1050 1100 1150 1200 1250 1300 1350 1400 1450 1500
4.63 ( 0.89(-20) 5.69 ( 0.89(-19) 1.98 ( 0.22(-17) 2.18 ( 0.18(-16) 1.24 ( 0.78(-15) 4.69 ( 0.23(-15) 1.34 ( 0.05(-14) 3.12 ( 0.10(-14) 6.32 ( 0.18(-14) 1.15 ( 0.03(-13) 1.91 ( 0.04(-13) 2.98 ( 0.05(-13) 4.40 ( 0.07(-13) 6.22 ( 0.09(-13) 8.46 ( 0.11(-13) 1.12 ( 0.01(-12) 1.43 ( 0.01(-12) 1.80 ( 0.02(-12) 2.22 ( 0.02(-12) 2.68 ( 0.02(-12) 3.20 ( 0.02(-12) 3.77 ( 0.02(-12) 4.38 ( 0.03(-12) 5.04 ( 0.03(-12) 5.75 ( 0.03(-12) 6.50 ( 0.03(-12) 7.29 ( 0.03(-12) 8.12 ( 0.03(-12)
QMa 1.09(-18) 2.20(-17) 2.12(-16) 1.19(-15) 4.54(-15) 1.33(-14) 3.18(-14) 6.57(-14) 1.21(-13) 2.05(-13) 3.22(-13) 4.80(-13) 6.81(-13) 9.30(-13) 1.23(-12) 1.58(-12) 1.99(-12)
Results from ref 63.
is some reaction below this collision energy in the QM calculations (see insert of Figure 2, top panel). However, the QCT σR(ET) is slightly larger than the QM one in the collision energy range between 0.30 and 0.37 eV, and this compensates, at high enough temperatures, the low-energy contribution to reaction, presumably from tunneling, in the QM calculation. At temperatures below 250 K, this latter contribution is predominant and gives rise to the observed discrepancies between classical and quantal rate constants. The comparison between QM and QCT excitation functions for the D + H2(V ) 1, j ) 0-3) reaction on the LSTH PES is displayed in Figure 3. As was already shown,69 there is a very good agreement between the σR(ET) from QCT calculations and that from exact QM results63 for V ) 1, j ) 0. In a subsequent work, Auerbach and Miller,65 using a polynomial expansion of the absorbing boundary condition Green’s function, calculated the σR(ET) for V ) 1, j ) 0-3. The results for V ) 1, j ) 0 so obtained proved to be coincident with those from the full stateto-state exact QM calculation.63 The comparison between QCT and QM63,65 σR(ET) exhibits an overall good agreement. However, whereas the classical threshold energy grows with initial j, as clearly shown in Figure 4, the QM one remains practically constant. Although these low ET discrepancies in the cross section are relatively small, they have an appreciable influence on the calculated rate constants k(T; V ) 1, j) as can be seen in Figure 5. Whereas the QCT rate constants are quite similar to the QM ones for D + H2(V ) 1, j ) 0), even at the lowest temperatures, the differences between QCT and QM k(T; V ) 1, j) become larger as the initial j increases, and for j ) 2, 3 only the QCT and QM rate constants calculated at high temperatures are in good agreement. From the above results, it is clear that the rotational excitation of the hydrogen molecule decreases the reactivity of the system and, except for the small discrepancy at threshold just commented on, the observed effect is very similar in the classical
4076 J. Phys. Chem., Vol. 100, No. 10, 1996
Aoiz et al.
Figure 3. Theoretical reaction cross sections as a function of the collision energy for the reactions D + H2(V ) 1, j ) 0-3) as calculated on the LSTH PES. Solid line and error bars (1 standard deviation): present QCT results; solid points and dashed line: QM results from ref 65, using a polynomial expansion of the absorbing boundary condition Green’s functions. The triangles of the curve for V ) 1, j ) 0, coincident within the precision of the figure with the previous QM ones, are the accurate QM results of ref 63, where the Kohn variational principle method was used.
Figure 4. Dependence of the reaction cross section on collision energy and rotational quantum number for the D + H2(V ) 1, j) system. Upper panel QM results from ref 65; lower panel QCT results from the present work in the same range of collision energies as the QM calculations.
and quantal results, in fact, when the collision energy is larger than 0.25 eV, the QM and QCT σR(ET,j) become almost
identical, irrespective of the initial j, and show a better agreement than that for V ) 0. As the collison energy increases, the negative effect of the initial rotation on reactivity seems to be overcome, and for collision energies larger than 0.6 eV the QCT σR increases with j. Table 2 contains the values of the QCT thermal rate constants k(T), averaged over initial j, for the D + n-H2(V)1) reaction from 200 to 400 K calculated on the three ab initio PESs, as well as the QM k(T) results on the LSTH PES. Since at the lowest temperatures the most populated state of the n-H2 is j ) 1, the discrepancies found for the V ) 1, j specific rate constants are reflected in the j-averaged thermal values. Thus, the differences between thermal QCT and QM k(T) for V ) 1 are caused not by strong tunneling for j ) 0 but by the increase of the classical threshold with the rotational excitation of the hydrogen molecule. At room temperature, this effect has little influence on the thermal rate constant. The most recent experimental k(T,V)1) values are also shown in Table 2.2,53,78,79 The QM and QCT theoretical results listed in the table all lie within the present experimental uncertainty. There are interesting differences in the QCT rate constants calculated on the three PESs, especially at low temperatures. The largest values for k(T) in this range of temperatures are predicted by the BKMP PES and the lowest by the DMBE one. There are also exact QM calculations available for the σR(ET) of the D + H2(V ) 0, j ) 0, 1) reaction on the DMBE PES80 that do not include the threshold region. Figure 6 shows the comparison between the QCT and QM results. As in the case of the LSTH PES, the QM σR(ET) are systematically larger than the QCT ones. The σR(ET) for j ) 1 seems to be in slightly better agreement than in the j ) 0 case. The role of initial rotation on the QCT excitation functions for the reaction with H2(V)0) is assessed in Figure 7. In this
The D + H2(V)0,1) f HD + H Reaction
J. Phys. Chem., Vol. 100, No. 10, 1996 4077
Figure 5. Arrhenius plot of the rate constants for the D + H2(V ) 1, j ) 0-3) reactions obtained by thermal averaging of the data of Figure 3. Dashed line and open circles (with error bars representing 1 standard deviation): present QCT results. Solid line and solid circles: QM results from ref 65. The open squares of the curve for V ) 1, j ) 0 are the QM results of ref 63 based on the Kohn variational principle method. The QCT-QM agreement worsens with the j value due to the discrepancies in the energy threshold region as j increases.
TABLE 2: Thermal Rate Constants k(T) (cm3 s-1) for the D + n-H2(W)1) Reaction at the Indicated Temperature and Potential Energy Surface (Numbers in Parentheses Represent Powers of 10)
a
T (K)
k(T)LSTHQCT
k(T)DMBEQCT
k(T)BKMPQCT
k(T)LSTHQMa
experimentalb
200 250 300 310 330 350 400
7.46 ( 2.02(-15) 4.06 ( 0.83(-14) 1.30 ( 0.21(-13) 1.57 ( 0.24(-13) 2.22 ( 0.32(-13) 3.04 ( 0.40(-13) 5.87 ( 0.64(-13)
5.70 ( 1.41(-15) 3.25 ( 0.61(-14) 1.08 ( 0.16(-13) 1.31 ( 0.19(-13) 1.88 ( 0.24(-13) 2.59 ( 0.31(-13) 5.10 ( 0.50(-13)
1.13 ( 0.26(-14) 5.86 ( 0.98(-14) 1.80 ( 0.24(-13) 2.16 ( 0.27(-13) 3.03 ( 0.35(-13) 4.09 ( 0.44(-13) 7.69 ( 0.69(-13)
1.23(-14) 5.56(-14) 1.57(-13) 1.87(-13) 2.59(-13) 3.48(-13) 6.54(-13)
1.0 ( 0.4(-13) 2.1 ( 0.4(-13)
Results from ref 65. b References 2, 53, 78, and 79.
figure, the σR(ET) for the D + H2(V ) 0, j ) 0, 3, 5, 7) reaction are shown from threshold up to ET ) 1.5 eV for the LSTH PES. A similar behavior (not shown) is obtained on the DMBE and BKMP surfaces. In all the three cases, the alreadymentioned negative effect of rotation on the reactivity is patent when initial j increases from j ) 0 to j ) 5 at low collision energies. For values of j > 5, the reactive yield seems to recover. However, the H2 rotation does not affect equally the σR(ET) calculated on the three different PESs. In the case of the calculations on the DMBE PES, the effect of rotation is more marked than for the other two PESs. On the other hand, at sufficiently high collision energies (0.6 eV for the BKMP PES and 0.7 eV for the LSTH PES), H2 rotation always favors the reaction. In the calculations on the DMBE PES, this only happens when ET > 0.9 eV. This dual effect of the initial j on the σR(ET) depending on the collision energy, causes the curves for different j to cross within a small interval of ET. This behavior was already noticed in the pioneering QCT calculations by Karplus et al.,9 who found that the translational threshold increased with j. Subsequently, and for different vibrational states of the H2 molecule, the effect of rotation has been studied for the hydrogen-exchange reaction on the LSTH PES in several
works.70,81 Unfortunately there are not available data from accurate quantum mechanical calculations about the evolution of the excitation function with j in the vicinity of the threshold for the reaction of D with hydrogen molecules in the ground vibrational state, but, by analogy with the V ) 1 case and from the values of the thermal rate constants (see below), one can expect a behavior similar to that described for the classical case. To clarify further the effect of the rotation of the reagents on the reactivity of the D + H2(V)0) system for the three PESs studied, Figure 8 depicts σR(ET) for initial j ) 0 and j ) 3 as calculated on the LSTH, DMBE and BKMP PESs. For H2(V ) 0, j ) 0), the reactivity grows in the sequence LSTH f BKMP f DMBE, over the whole range of collision energies. The threshold on the LSTH PES is slightly larger than in the other two PESs. In contrast, for D + H2(V ) 0, j ) 3), the reactivity on the DMBE PES is smaller than that on the BKMP PES and similar to that on the LSTH PES. Interestingly, the threshold energy for the DMBE PES seems to be in this case the largest. Thus, the initial rotation of H2(V)0) affects the system in such a way that the reactivity calculated on the DMBE PES changes from being the highest for j ) 0 to being the lowest when j ) 3 at low collision energies. Only at high enough
4078 J. Phys. Chem., Vol. 100, No. 10, 1996
Figure 6. Reaction cross section as a function of the collision energy for the D + H2(V ) 0, j ) 0) (top) and D + H2(V ) 0, j ) 1) (bottom), calculated on the DMBE PES. Solid line and error bars: present QCT results; solid circles joined by dashed line: accurate QM results of ref 80.
Figure 7. Present QCT reaction cross sections as a function of the collision energy for the D + H2(V ) 0, j ) 0, 3, 5, 7) reactions calculated on the LSTH PES. The insert displays the cross sections near the threshold.
collision energies, ET g 0.7 eV, the cross section on the DMBE PES becomes larger than that on the LSTH PES. In contrast, rotation seems to influence the reactivity on both LSTH and BKMP PESs in a very similar way, and for all the j values, the BKMP PES predicts more reactivity than the LSTH PES. Analogous results have been obtained for the reaction with H2 in V ) 1 and are shown in Figure 9 for the three ab initio PESs. The effect of the rotation of the H2(V)1) reagent on the reactivity is entirely analogous to the one shown for H2(V)0).
Aoiz et al.
Figure 8. Present QCT excitation functions, σR(ET), for the D + H2(V ) 0, j ) 0) (top) and D + H2(V ) 0, j ) 3) (bottom) calculated on the LSTH (solid line), DMBE (dotted line) and BKMP (dashed line) PESs. The inserts display the σR(ET) near the threshold to reaction.
Except, perhaps, in the vicinity of the threshold, the DMBE PES yields the highest reactivity when j ) 0. Once more, when j ) 3 the DMBE PES predicts the lowest cross sections for ET < 0.6 eV, even smaller than those for the LSTH PES. Thus, the effect of the initial rotational energy seems to be enhanced for V ) 1. This is also clearly seen in the data of k(T;V)1) of Table 2 commented on above; the k(T;V)1) predicted by the DMBE PES are always smaller than in the other two PESs and at 200 K the rate on the BKMP PES is twice that on the DMBE PES. An Arrhenius plot of the QCT thermal rate constants, k(T), averaged over initial j, for the D + n-H2 reaction calculated on the three ab initio PESs is shown in Figure 10 (right panel). The corresponding QM results66 are also shown in the same figure (left panel). Tables 3-5 contain the present QCT k(T) for the D + n-H2 reaction and the QM ones calculated by Mielke et al.66 on the three PESs, LSTH, DMBE, and BKMP, respectively, as well as the data from the three-parameter fit to the experimental results given by Michael and Fisher,56 which are represented as solid circles in the two panels of Figure 10. In addition, Table 3 includes the QM calculations of Park and Light64 on the LSTH PES and Table 4 the reduced dimensionality calculations of Wang and Bowman68 on the DMBE PES. The QCT thermal rate constants on any of the three surfaces are always smaller than the values of the mentioned threeparameter fit to the experimental data; in addition, the QCT k(T) up to 1000 K are lower than any of the corresponding QM values calculated on the same PES. The largest discrepancies between QCT and experiment are found in the low-temperature range. At 200 K, the difference is almost a factor of 3 for LSTH
The D + H2(V)0,1) f HD + H Reaction
J. Phys. Chem., Vol. 100, No. 10, 1996 4079 TABLE 4: Thermal Rate Constants k(T) (cm3 s-1) for the D + n-H2 Reaction as a Function of Temperature Calculated on the DMBE PES (Numbers in Parentheses Represent Powers of 10) T (K)
present work
200 250 300 350 400 450 500 550 600 650 700 750 800 850 900 950 1000 1050 1100 1150 1200 1250 1300 1350 1400 1450 1500
5.20 ( 6.30(-19) 1.70 ( 1.27(-17) 1.81 ( 0.92(-16) 1.01 ( 0.37(-15) 3.79 ( 1.04(-15) 1.07 ( 0.23(-14) 2.51 ( 0.42(-14) 5.10 ( 0.70(-14) 9.30 ( 1.06(-14) 1.56 ( 0.15(-13) 2.44 ( 0.20(-13) 3.63 ( 0.26(-13) 5.17 ( 0.32(-13) 7.09 ( 0.39(-13) 9.40 ( 0.46(-13) 1.21 ( 0.05(-12) 1.54 ( 0.06(-12) 1.90 ( 0.07(-12) 2.32 ( 0.08(-12) 2.78 ( 0.08(-12) 3.29 ( 0.09(-12) 3.84 ( 0.10(-12) 4.44 ( 0.11(-12) 5.08 ( 0.11(-12) 5.76 ( 0.12(-12) 6.48 ( 0.13(-12) 7.24 ( 0.13(-12)
a
Figure 9. As in Figure 8, but for the D + H2(V ) 1, j ) 0, 3) reaction.
TABLE 3: Thermal Rate Constants k(T) (cm3 s-1) for the D + n-H2 Reaction as a Function of Temperature Calculated on the LSTH PES (Numbers in Parentheses Represent Powers of 10) T (K)
present work
200 250 300 350 400 450 500 550 600 650 700 750 800 850 900 950 1000 1050 1100 1150 1200 1250 1300 1350 1400 1450 1500
5.52 ( 4.18(-19) 1.78 ( 0.88(-17) 1.87 ( 0.67(-16) 1.03 ( 0.28(-15) 3.79 ( 0.80(-15) 1.06 ( 0.18(-14) 2.47 ( 0.34(-14) 4.98 ( 0.58(-14) 9.04 ( 0.88(-14) 1.51 ( 0.13(-13) 2.36 ( 0.17(-13) 3.51 ( 0.22(-13) 4.98 ( 0.27(-13) 6.82 ( 0.33(-13) 9.06 ( 0.40(-13) 1.17 ( 0.05(-12) 1.48 ( 0.05(-12) 1.84 ( 0.06(-12) 2.24 ( 0.07(-12) 2.69 ( 0.07(-12) 3.18 ( 0.08(-12) 3.72 ( 0.08(-12) 4.30 ( 0.09(-12) 4.92 ( 0.10(-12) 5.59 ( 0.10(-12) 6.29 ( 0.11(-12) 7.03 ( 0.11(-12)
a
b
Mielke et al.a 1.76(-18) 3.19(-17) 2.76(-16)
Park and Lightb
3.2(-16)
5.03(-15) 3.17(-14)
3.3(-14)
1.14(-13) 2.92(-13)
2.8(-13)
6.09(-13) 1.10(-12)
9.7(-13)
1.80(-12) 2.72(-12)
2.2(-12)
3.89(-12) 5.32(-12)
4.1(-12)
7.02(-12) 8.99(-12)
6.6(-12)
experimentc 1.47(-18) 3.39(-17) 2.96(-16) 1.47(-15) 5.11(-15) 1.39(-14) 3.17(-14) 6.36(-14) 1.15(-13) 1.94(-13) 3.07(-13) 4.62(-13) 6.67(-13) 9.29(-13) 1.26(-12) 1.66(-12) 2.14(-12) 2.72(-12) 3.39(-12) 4.16(-12) 5.04(-12) 6.04(-12) 7.17(-12) 8.42(-12) 9.81(-12) 1.13(-11) 1.30(-11)
Wang and Bowmanb experimentc
Mielke et al.a 1.64(-18) 3.14(-17) 2.78(-16)
1.90(-18)
5.13(-15)
6.10(-15)
3.24(-14)
3.80(-14)
1.16(-13)
1.40(-13)
2.98(-13)
3.50(-13)
6.21(-13)
7.30(-13)
1.12(-12)
1.30(-12)
1.83(-12)
2.10(-12)
2.77(-12)
3.20(-12)
3.95(-12)
4.60(-12)
5.40(-12)
6.30(-12)
7.13(-12)
8.30(-12)
9.12(-12)
1.10(-11)
3.30(-16)
1.47(-18) 3.39(-17) 2.96(-16) 1.47(-15) 5.11(-15) 1.39(-14) 3.17(-14) 6.36(-14) 1.15(-13) 1.94(-13) 3.07(-13) 4.62(-13) 6.67(-13) 9.29(-13) 1.26(-12) 1.66(-12) 2.14(-12) 2.72(-12) 3.39(-12) 4.16(-12) 5.04(-12) 6.04(-12) 7.17(-12) 8.42(-12) 9.81(-12) 1.13(-11) 1.30(-11)
Reference 66. b Reference 68. c Reference 56.
TABLE 5: Thermal Rate Constants k(T) (cm3 s-1) for the D + n-H2 Reaction as a Function of Temperature Calculated on the BKMP PES (Numbers in Parentheses Represent Powers of 10) T (K)
present work
200 250 300 350 400 450 500 550 600 650 700 750 800 850 900 950 1000 1050 1100 1150 1200 1250 1300 1350 1400 1450 1500
7.80 ( 4.45(-19) 2.48 ( 0.89(-17) 2.55 ( 0.65(-16) 1.38 ( 0.27(-15) 4.98 ( 0.78(-15) 1.37 ( 0.18(-14) 3.13 ( 0.34(-14) 6.23 ( 0.58(-14) 1.11 ( 0.09(-13) 1.84 ( 0.13(-13) 2.85 ( 0.18(-13) 4.19 ( 0.24(-13) 5.90 ( 0.30(-13) 8.03 ( 0.37(-13) 1.06 ( 0.04(-12) 1.36 ( 0.05(-12) 1.71 ( 0.06(-12) 2.11 ( 0.07(-12) 2.56 ( 0.07(-12) 3.06 ( 0.08(-12) 3.61 ( 0.09(-12) 4.20 ( 0.10(-12) 4.84 ( 0.10(-12) 5.53 ( 0.11(-12) 6.25 ( 0.12(-12) 7.02 ( 0.13(-12) 7.82 ( 0.13(-12)
a
Mielke et al.a 3.21(-18) 5.24(-17) 4.21(-16) 6.93(-15) 4.11(-14) 1.42(-13) 3.54(-13) 7.23(-13) 1.29(-12) 2.07(-12) 3.11(-12) 4.40(-12) 5.98(-12) 7.84(-12) 9.97(-12)
experimentb 1.47(-18) 3.39(-17) 2.96(-16) 1.47(-15) 5.11(-15) 1.39(-14) 3.17(-14) 6.36(-14) 1.15(-13) 1.94(-13) 3.07(-13) 4.62(-13) 6.67(-13) 9.29(-13) 1.26(-12) 1.66(-12) 2.14(-12) 2.72(-12) 3.39(-12) 4.16(-12) 5.04(-12) 6.04(-12) 7.17(-12) 8.42(-12) 9.81(-12) 1.13(-11) 1.30(-11)
b
Reference 66. Reference 56.
c
Reference 66. Reference 64. Reference 56.
and DMBE surfaces, whereas it is less than a factor of 2 for the BMKP PES. These differences become much smaller for higher temperatures especially in the case of the BKMP surface.
Given the good agreement between QM results and experiment over this temperature range, an analogous effect is observed in the comparison of QCT and QM results; the largest discrepancies correspond to the lowest temperatures and the accordance
4080 J. Phys. Chem., Vol. 100, No. 10, 1996
Figure 10. Thermal rate constants k(T) for the D + n-H2 f HD + H reaction calculated on the three ab initio PESs, LSTH, DMBE, and BKMP. Left: accurate QM results from ref 66. Right: present QCT results. The solid points represent a fit to the experimental results (see ref 56).
becomes better with increasing T. The degree of agreement between quantal and classical rate constants is similar on any of the three PES. In the range between 1000 and 1500 K, the QCT k(T) agree well with the full-dimensionality QM rate constants reported by Park and Light64 (see Table 3), but both QCT and these QM results deviate gradually from the experimental data, and at 1500 K the theoretical values are about one-half of the measured ones. A better accordance with experiment in this high-temperature range is found with the calculations of Mielke et al.,66 which include accurate values of the cumulative reaction probability for total energies E < 1.2 eV and extrapolated values for higher energies. These authors point out that the previous QM calculations by Park and Light might not be fully converged at the highest temperatures. The rate constants calculated by Mielke et al. at 1500 K are 30% smaller than the experimental ones. A similar degree of agreement is obtained by the same authors with the SRA approximation66 (not shown in the table). A similar accordance between k(T) measurements and calculations over the high-temperature range from 800 to 2000 K is found with the reduced-dimensionality calculations of Wang and Bowman68 on the DMBE PES. As pointed out in ref 66, the accurate QM results obtained using the LSTH and DMBE PESs are quite similar to each other and seem to agree very well with the experimental data up to temperatures of 1000 K. In contrast, the k(T) calculated on the BKMP PES are systematically higher than those on the other two PESs and also higher than the experimental data for temperatures up to 800 K. The discrepancy between the k(T) calculated on the three PESs and especially between the results obtained on the LSTH or DMBE PESs and those on the BKMP PES is worth noting. Only at temperatures high enough (T > 800 K) the three sets of data seem to converge to similar values. The same discrepancy is found between the QCT k(T) calculated on the three PESs, as well as the same sequence in the k(T) values at low temperatures; namely, KDMBE(T) < kLSTH(T) < kBKMP(T). The results of the QCT σR(ET) for different initial j shown above, indicate that the distinct effect of rotation on reactivity obtained for each of the potential surfaces, is at the origin of the discrepancies found in the respective thermal rate constants. IV. Discussion At the lowest temperatures, the classical k(T) of the present work are smaller than those from quantum mechanics and experiment. In general, tunneling is invoked to justify this behavior. Moreover, in the cases where QM and QCT rate
Aoiz et al. constants are similar, as for j ) 0, a cancelation of two errors: absence of tunneling and absence of zero-point energy constraints in the classical calculation, has been conjectured as an explanation for the similitude.35 Nevertheless, a direct analysis of tunneling is not straightforward; most studies are based on unidimensional formulations applied to model systems,35 and it is usually not possible to infer from them the actual influence of tunneling in accurate three-dimensional calculations of dynamical observables. On the other hand, it is also not clear that quantum mechanical calculations describing the passage of the saddle-point region are totally adiabatic.32,82 Although a lucky cancelation of errors is strictly possible, the similitude in the shapes of the QM and QCT excitation functions for D + H2(V ) 0, j ) 0) and D + H2(V ) 1, j ) 0) (see Figures 2 and 3) for collision energies close to the reaction threshold (ET < 0.4 eV) suggests also the possibility of very similar quantal and classical dynamics for reactive collisions of D atoms with rotationless H2 molecules. However, as j increases, there is a neat enhancement of tunneling, that is, whereas the classical threshold increases with j for low j, the quantum dynamical threshold (i.e., the collision energy at which the cross section reaches a given, very small, value) remains practically unchanged, as clearly shown in Figure 4 for the V ) 1 case. Unfortunately, there are no data available from accurate quantum mechanical calculations about the evolution of the excitation function with j in the vicinity of the threshold for the reaction of D with hydrogen molecules in the ground vibrational state, but the analogy with the V ) 1 case and the comparison of the low-temperature k(T,j)0) and k(T) from both QM and QCT calculations strongly suggest that this rotational enhancement of tunneling is the main cause of discrepancy between the quantal and classical k(T) for the reaction of D with n-H2. There are additional indications along this direction from the analysis of the role of the quantized transition states in the accurate quantal dynamics of the H + H2 and D + H2 reactions carried out by Chatfield et al.86 The quantum dynamical thresholds for H + H2 for total angular momentum J ) 0, first decrease, for low j, and then show an oscillatory behavior with increasing j, and this interpreted in terms of the correlations between the quantal states of the reagents and those of the transition states. The authors conclude that the quantized levels of the transition state determine the state-selected chemical reactivity for the entire range of energies relevant to the thermal rate constants. In addition to the role of the rotation at the threshold commented on in the previous paragraph, which is clearly different in the quantal and classical cases, there is another remarkable effect of rotation on reactivity common to both the QCT and QM approaches; namely, a decrease of reactivity with increasing j, starting at collision energies slightly above the threshold, which can be also seen in Figure 4. The QCT σR(ET,j) values indicate that rotational excitation (j < 5) of the H2 molecules causes a decrease in the reactivity for collision energies up to 0.6-0.9 eV depending on the PES. This negative influence of rotation is commonly described as an orientation effect of the potential energy surface which has been extensively discussed in the literature in QCT calculations and several models.83-85 The QM excitation functions for V ) 1, j ) 0-3 calculated on the LSTH PES (see Figures 3 and 4) and those for V ) 0, j ) 0, 1 on the DMBE surface (see Figure 6) indicate that this orientation effect is also present in the QM results. The general idea is that the forces outside the barrier tend to steer the reactants into the cone of acceptance in the absence of rotation. Initial rotation would hinder this steering, leading thus to a decline of the reaction cross section with j.
The D + H2(V)0,1) f HD + H Reaction
J. Phys. Chem., Vol. 100, No. 10, 1996 4081
Figure 11. Potential energy versus bending D-H-H angle, R, for the LSTH, DMBE, and BKMP PESs; R ) 180° corresponds to the collinear configuration. For every angle the intenuclear distances, RAB ) RBC, are those of the saddle point at each angular configuration.
As shown in the previous section, small differences in these “orientational” forces outside the barrier, as in the cases of the three PESs here considered, can lead to appreciable differences in the dynamical observables. Even averaged magnitudes as total cross sections and thermal rate constants can differ considerably depending on the PES where the calculations have been carried out. The next step is to try to relate these observations with the topological features of the PES. The characteristics of the PES more relevant to the interpretation of the above-mentioned dynamical observables are summarized in Figures 11-13. In particular, Figure 11 shows the bending potential of the three PES, i.e., the potential energy at the saddle point as a function of the bending angle R (DHH angle), or, in other words, the location of the barrier as a function of R. The bending potentials for all the PESs are very similar. In fact, as the transition state becomes more bent, the barriers of the DMBE and BKMP PESs become almost identical. It is in the collinear configuration (R ) 180°) where the differences, albeit small, are clearer (see insert of Figure 11). The collinear barrier increases in the sequence BKMP f DMBE f LSTH (0.4137, 0.4184, and 0.4249 eV, respectively). According to this order, the expected sequence of rate constants would be kLSTH(T) < kDMBE(T) < kBKMP(T). Neither the thermal k(T) at low temperatures, nor the σR(ET, j ) 0, 3) follow this sequence. It is evident, therefore, that considerations solely based on the bending potential are insufficient. A good picture of the PES is provided by the R-γ plots,83 V(R,γ;r), which consist of contour maps of the potential as a function of the distance between the incoming D atom and the center of mass of the H2 molecule, R, and the angle, γ, between R and the intermolecular rH2 axis for a fixed value of this latter distance. The R-γ plots for each of the three PES are depicted in Figure 12. In each case, the H2 internuclear distance is fixed at its value in the collinear saddle point (0.929 77 Å for the LSTH PES; 0.928 51 Å for the DMBE PES; 0.929 76 Å for the BKMP PES). Therefore, the absolute minimum in each of the contour plots of Figure 12 corresponds to the collinear saddle point. In spite of the similitudes between the three PESs, there are interesting differences. The curvature of the contours clearly indicates that in the case of the DMBE PES the gradient of the potential in the direction of γ is larger than in the other two PESs (indicated by the thick arrows that result from the decomposition of the forces tangent to the contours). These anisotropies of the potential imply orientation forces and will tend to dynamically reorient the reactants during the
Figure 12. R-γ contour plot of the LSTH, DMBE, and BKMP potential energy surfaces at the internuclear distance, rH2, corresponding to the collinear saddle point (0.929 77 Å, LSTH; 0.928 51 Å, DMBE; 0.929 76 Å, BKMP). The thick arrows resulting from the decomposition of the forces tangent to the contours represent the gradient of the potential in the direction of γ. The contours are plotted every 0.02 eV, and the first contour (corresponding to the largest R) is 0.44 eV. The minima shown in every figure correspond to the respective collinear saddle points for each PES (0.4137 eV for the BKMP, 0.4184 eV for the DMBE, and 0.4249 eV for the LSTH).
approach motion to the barrier, therefore aligning the incoming D atom with the H-H axis. In the absence of rotation, and mainly at low collision energies, these orientation forces contribute to increase the reaction yield. The DMBE PES is the one that has larger orientation forces, and consequently the reactivity at low collision energies is significantly enhanced with respect to the other two PESs. The orientation effects are also clear in Figure 13, where the three potential energy surfaces are represented as a function of the γ angle at a fixed distance of approach, R ) 2.1 Å, and for an internuclear distance corresponding to the minimum of the H2 potential. These R, r values roughly correspond to the stage
4082 J. Phys. Chem., Vol. 100, No. 10, 1996
Figure 13. H3 potential energy as a function of γ at RD-H2 ) 2.1 Å and rH2 at the equilibrium H2 internuclear distance. Solid line: LSTH PES; short-dashed line: DMBE PES; long-dashed line: BKMP PES. Notice that for the DMBE PES the gradient of the potential with respect to γ is larger than for the other two PESs.
where the three atoms start to interact more closely. It is evident from the figure that the DMBE potential changes with the attacking angle more than the other two PESs, thus creating a somewhat larger torque, that tends to direct the incoming atom toward the minimum of the potential more efficiently than the other two PESs. This explains the order in the calculated reaction cross sections for j ) 0. However, when some rotation is present, the efficiency of this orienting force becomes negligible as compared with the rotation. At low collision energies, initial reagent rotation would hinder the entry into the cone of acceptance. This explains the negative effect of rotation on the reactivity on all the PESs. In addition, the net effect of rotation is to average out the potential to all the orientations. As a result, the “average” potential that sees the approaching atom in the case of the DMBE PES is slightly larger than that of the LSTH or BKMP PESs. Thus, the reaction cross section for j > 0 calculated on the DMBE PES is somewhat smaller than those on the LSTH and BKMP PESs. As the collision energy increases, the orientation effects become less important and the cross sections on the DMBE PES tend to the values obtained on the BKMP PES. The classical mechanics approach used throughout this work facilitates thus the investigation of the effects of rotation on reactivity in terms of the properties of the potential energy surface and provides an explanation for the relative magnitude of the thermal rate constants on each of the three PES studied. This rate constant trend is qualitatively the same as the one encountered in accurate QM calculations,66 but the classical k(T)’s are systematically shifted toward lower values (especially at low temperatures) than the QM ones due to the higher classical thresholds for reaction with rotationally excited hydrogen. Consequently, the QCT results cannot be used to discriminate between the three accurate PES currently used for this reaction. In fact, the best accord between QCT and experiment is obtained on the BKMP surface, whereas the calculations by Mielke et al.66 clearly show that the lowtemperature QM k(T) produced on the DMBE or LSTH PESs are in better agreement with the experimental results than those obtained on the BKMP PES, which, on the other hand, is supposed to be the most accurate of the three PESs and the one that fits better the ab initio points. Although the absolute values of the QCT rate constants cannot be used to distinguish the best between the very similar ab initio surfaces existing for this reaction, the insight provided by the classical study can be of help in order to design a more rigorous accuracy test. If the potential surface dependent rotational effect on reactivity stressed throughout this work is confirmed, as expected, by accurate QM calculations for the V ) 0 case and
Aoiz et al. for the different PES, further experiments at low temperatures with selected rotational states of the H2 would supply a more stringent test of the adequacy of a given ab initio surface. In fact, measurements of rate constants for the D + p-H2(V)0,1) at low temperatures, where j ) 0 is the most populated state, might be checked against the QM calculations and provide a conclusive test to elucidate which of the mentioned PES is more accurate. The QCT predictions indicate that at 300 K the k(T) for the reaction with p-H2 is 16% larger than that for n-H2 if calculations are carried out on the DMBE PES, whereas, in the case of the LSTH PES, it is about 7% smaller, and in the BKMP both rate constants are practically the same. At lower temperatures, these differences are more pronounced and thus easily contrastable with experimental measurements. Even more decisive than the measurement of rate constants would be the experimental determination of absolute cross sections for the reaction with rotationally excited hydrogen in a similar way to that reported in ref 87 but with a narrower H2 rotational distribution. V. Conclusions Exhaustive QCT calculations have been carried out in order to determine the collision energy dependence of the reaction cross section for the D + H2(V ) 0, 1; j ) 0-7) reaction from threshold up to 1.5 eV, on the three available, chemically accurate, ab initio PESs. The comparison with the available exact QM results reveals a quite good agreement. A very strong effect of initial rotation has been found for collision energies below 0.6-0.9 eV; reactivity decreases with initial j e 5. This negative effect seems also to be present in the QM calculations. Thermal state specific rate constants, k(T;V,j), as well as rate constants averaged over initial j, k(T;V), have been calculated from the excitation functions, σ(ET;V,j), over the range 2001500 K. For the reaction with rotationless H2, there is a quite good agreement between QCT and QM results for both V ) 0, 1 states of the reagents, and in both LSTH and DMBE PESs. However, the rotational state-specific rate constants for j > 0 show significant differences, at low temperatures, where the classical k(T,j) are markedly lower than the QM ones. This different rotational behavior leads to QCT thermal rate constants, k(T), that are smaller than those from quantum mechanics and experiment. The influence of the rotation can be related to the anisotropy of the potential in the way to the barrier, which originates a relatively distinct orienting character of these PESs. In the case of the DMBE the potential changes with the attacking angle more than in the other two PESs and thus is the one with more orienting tendency. This explains why the low collision energy cross sections for j ) 0 are larger than in the LSTH or BKMP PESs and why they become the smallest for j > 0. Therefore, although the three PESs are very similar and their collinear barriers do not differ by more than 0.01 eV, both QM and QCT rate constant calculations at low temperatures yield distinct results depending on the PES. This study suggests that careful measurement of lowtemperature rate constants or of cross sections for the reaction with rotationally excited hydrogen molecules, could provide invaluable data to discern between the several accurate PESs. Acknowledgment. We are most grateful to Mark Brouard for many enlightening discussions. We are indebted to S. M. Auerbach for sending us the results of ref 65 and to S. L. Mielke for letting us know their QM results prior to publication. We also acknowledge W. J. Keogh for providing us with the code of the BKMP potential and T. Dreier for his comments on the
The D + H2(V)0,1) f HD + H Reaction experiments of refs 2 and 53. T.D.-T. acknowledges the financial support through a F.P.I. fellowship of the Ministry of Education and Science of Spain. This work has been financed by the DGICYT of Spain (PB92-0219-C03). The SpanishBritish scientific exchange program Acciones Integradas is also acknowledged. References and Notes (1) The Farkas Symposium-Fifty Years of H+H2 Kinetics; Levine, R. D., Ed.; Int. J. Chem. Kinet. 1986, 18, 9. (2) Buchenau, H.; Toennies, J. P.; Arnold, J.; Wolfrum, J. Ber. BunsenGes. Phys. Chem. 1990, 94, 1231. (3) (a) Farkas, A. Z. Phys. Chem. 1930, 10, 419. (b) Farkas, A.; Farkas, L. Proc. R. Soc. London A 1935, 152, 124. (4) (a) London, F. In Probleme der Modernen Physik; Hirzel: Leipzig, 1928; p 104. (b) London, F. Z. Elektrochem. 1929, 35, 552. (5) Eyring, H.; Polanyi, M. Z. Phys. Chem. 1931, B12, 279. (6) Pelzer, H.; Wigner, E. Z. Phys. Chem. 1932, B15, 445. (7) Evans, M. G.; Polanyi, M. Trans. Faraday Soc. 1935, 31, 857. (8) Eyring, H. J. Chem. Phys. 1935, 3, 107. (9) Karplus, M.; Porter, R. N.; Sharma, R. D. J. Chem. Phys. 1965, 43, 3259. (10) Porter, R. N.; Karplus, M. J. Chem. Phys. 1964, 40, 1105. (11) Porter, R. N.; Karplus, M. Discuss. Faraday Soc. 1967, 44, 164. (12) Schulz, W. R.; Le Roy, D. J. Can. J. Chem. 1964, 42, 2480. (13) Schulz, W. R.; Le Roy, D. J. J. Chem. Phys. 1965, 42, 3869. (14) Ridley, B. A.; Schulz, W. R.; Le Roy, D. J. J. Chem. Phys. 1966, 44, 3344. (15) Le Roy, D. J.; Ridley, B. A.; Quickert, K. A. Discuss. Faraday Soc. 1967, 44, 92. (16) Quickert, K. A.; Le Roy, D. J. J. Chem. Phys. 1970, 53, 1325. (17) Quickert, K. A.; Le Roy, D. J. J. Chem. Phys. 1971, 54, 5444. (18) Mitchell, D. N.; Le Roy, D. J. J. Chem. Phys. 1973, 58, 3449. (19) Westenberg, A. A.; de Haas, N. J. J. Chem. Phys. 1967, 47, 1393. (20) Kuppermann, A.; Schatz, G. C.; Baer, M. J. Chem. Phys. 1976, 65, 4596. (21) Schatz, G. C.; Kuppermann, A. J. Chem. Phys. 1976, 65, 4668. (22) Liu, B. J. Chem. Phys. 1973, 58, 1925. (23) Siegbahn, P.; Liu, B. J. Chem. Phys. 1978, 68, 2457. (24) Truhlar, D. G.; Horowitz, C. J. J. Chem. Phys. 1978, 68, 2466; 1979, 71, 1514(E). (25) Mayne, H. R.; Toennies, J. P. J. Chem. Phys. 1981, 75, 1794. (26) Sun, J. C.; Choi, B. H.; Poe, R. T.; Tang, K. T. J. Chem. Phys. 1980, 73, 6095; 1983, 79, 5376. (27) Bowman, J. M.; Lee, K. T.; Walker, R. B. J. Chem. Phys. 1983, 79, 3742. (28) Abu Salbi, N.; Kouri, D. J.; Shima, Y.; Baer, M. J. Chem. Phys. 1985, 82, 2650. (29) Colton, M. C.; Schatz, G. C. Int. J. Chem. Kinet. 1986, 18, 961. (30) Garrett, B. C.; Truhlar, D. G. J. Chem. Phys. 1980, 72, 3460. (31) Pechukas, P. Annu. ReV. Phys. Chem. 1981, 32, 159. (32) Truhlar, D. G.; Isaacson, A. D.; Garrett, B. C. In Theory of Chemical Reaction Dynamics; Baer, M., Ed.; CRC Press: Boca Raton, FL, 1985; Vol. IV, p 65. (33) Garrett, B. C.; Truhlar, D. G. Proc. Natl. Acad. Sci. U.S.A. 1979, 76, 4755. (34) Garrett, B. C.; Truhlar, D. G.; Varandas, A. J. C.; Blais, N. C. Int. J. Chem. Kinet. 1986, 18, 1065. (35) Schatz, G. Chem. ReV. 1987, 87, 81. (36) Schatz, G. Annu. ReV. Phys. Chem. 1988, 39, 317. (37) Gordon, E. M.; Ivanov, B. I.; Perminov, A. P.; Medvedev, E. S.; Ponomarev, A. N.; Tal’Rose, V. L. Chem. Phys. 1975, 8, 147. (38) Gordon, E. B.; Ivanov, B. I.; Perminov, A. P.; Balalaev, V. G.; Ponomarev, A.; Filatov, V. V. Chem. Phys. Lett. 1978, 58, 425. (39) Kneba, M.; Wellhausen, U.; Wolfrum, J. Ber. Bunsen-Ges. Phys. Chem. 1979, 83, 940. (40) Wellhausen, U.; Wolfrum, J. Ber. Bunsen-Ges. Phys. Chem. 1982, 89, 314. (41) Glass, G. P.; Chaturvedi, B. K. J. Chem. Phys. 1982, 77, 3478. (42) Rozenshtein, V. B.; Gershenzon, Y. M.; Ivanov, A. V.; Il’in, S. D.; Kucheryavii, S. I. Chem. Phys. Lett. 1984, 105, 423. (43) Rozenshtein, V. B.; Gershenzon, Y. M.; Ivanov, A. V.; Il’in, S. D.; Kucheryavii, S. I.; Umanskii, S. Y. Chem. Phys. Lett. 1985, 121, 89. (44) Wellhausen, U.; Wolfrum, J. Ber. Bunsen-Ges. Phys. Chem. 1985, 89, 316. (45) Sun, J. C.; Choi, B. H.; Poe, R. T.; Tang, K. T. Phys. ReV. Lett. 1980, 44, 1211. (46) Pollak, E.; Wyatt, R. E. J. Chem. Phys. 1983, 78, 4464. (47) Walker, R. B.; Hayes, E. F. J. Phys. Chem. 1255, 87, 1255. (48) Abu Salbi, N.; Kouri, D. J.; Schima, Y.; Baer, M. Chem. Phys. Lett. 1984, 105, 472. (49) Walker, R. B.; Pollak, E. J. Chem. Phys. 1985, 83, 2851.
J. Phys. Chem., Vol. 100, No. 10, 1996 4083 (50) Pollak, E.; Abu Salbi, N.; Kouri, D. J. Chem. Phys. Lett. 1985, 113, 585. (51) Garrett, R. C.; Truhlar, D. G. J. Phys. Chem. 1985, 89, 2204. (52) Haug, K.; Schwenke, D. W.; Shima, Y.; Truhlar, D. G.; Zhang, J.; Kouri, D. J. J. Phys. Chem. 1986, 90, 6757. (53) Dreier, T.; Wolfrum, J. Int. J. Chem. Kinet. 1986, 18, 919. (54) Jayaweera, I. S.; Pacey, P. D. J. Phys. Chem. 1990, 94, 3614. (55) Michael, J. V. J. Chem. Phys. 1990, 92, 3394. (56) Michael, J. V.; Fisher, J. R. J. Phys. Chem. 1990, 94, 3318. (57) Bowman, J. M. AdV. Chem. Phys. 1985, 61, 115. (58) Michael, J. V.; Fisher, J. R.; Bowman, J. M.; Sun, Q. Science 1990, 24, 269. (59) Varandas, A. J. C.; Brown, F. B.; Mead, C. A.; Truhlar, D. G.; Garrett, B. C. J. Chem. Phys. 1987, 86, 6258. (60) Boothroyd, A. I.; Keogh, W. J.; Martin, P. G.; Peterson, M. R. J. Chem. Phys. 1991, 4343, 95. (61) Miller, W. H. Annu. ReV. Phys. Chem. 1990, 41, 245. (62) Mladenovic, M.; Zhao, M.; Truhlar, D. G.; Schwenke, D. W.; Sun, Y.; Kouri, D. J. J. Phys. Chem. 1988, 92, 7035. (63) Zhang, J. Z. H.; Miller, W. H. J. Chem. Phys. 1989, 91, 1528. (64) Park, T. J.; Light, J. C. J. Chem. Phys. 1991, 94, 2946. (65) Auerbach, S. M.; Miller, W. H. J. Chem. Phys. 1994, 100, 1103. (66) Mielke, S. L.; Lynch, G. C.; Truhlar, D. G.; Schwenke, D. W. J. Phys. Chem. 1994, 98, 8000. (67) Takada, S.; Ohsaki, A.; Nakamura, H. J. Chem. Phys. 1992, 96, 339. (68) Wang, D.; Bowman, J. M. J. Phys. Chem. 1994, 98, 7994. (69) Aoiz, F. J.; Buchenau, H. K.; Herrero, V. J.; Sa´ez Ra´banos, V. J. Chem. Phys. 1994, 100, 2789. (70) Aoiz, F. J.; Herrero, V. J.; Sa´ez Ra´banos, V. J. Chem. Phys. 1991, 94, 7991. (71) Aoiz, F. J.; Herrero, V. J.; Sa´ez Ra´banos, V. J. Chem. Phys. 1992, 97, 7423. (72) Throughout this work we have used the source routine of this PES provided directly by the authors (version 7.06). In the comments of this routine, the authors made clear that some few parameters (the Cbend coefficients) were changed from those originally published in ref 60 in order to avoid a spurious van der Waals well for compressed geometries of the H2. In any case, it is unlikely that these changes could affect the reactive scattering. (73) Huber, K. P.; Herzberg, G. Molecular Spectra and Molecular Structure. Part IV. Constants of Diatomic Molecules; Van Nostrand: New York, 1979. (74) Dabrowski, I.; Herzberg, G. Can. J. Phys. 1976, 54, 525. (75) Kolos, W.; Wolniewicz, L. J. Mol. Struct. 1975, 54, 303. (76) Truhlar, D. G.; Muckerman, J. T. In Atom-Molecule Collision Theory; Bernstein, R. B., Ed.; Plenum: New York, 1979. (77) Truhlar, D. G.; Blais, N. C. J. Chem. Phys. 1977, 67, 1532. (78) Dreier, T., personal communication. (79) In ref 65 the experimental rate constant at 310 K quoted from ref 53 is 1.9 × 10-13 cm3 s-1. This value, however, corresponds to the rate of total decay of H2(V)1), vibrational quenching plus reaction, whereas the actual reaction rate constant is 1.0 × 10-13 cm3 s-1. In any case, the experimental uncertainty probably includes the theoretical QM value.2,53,78 (80) Zhao, M.; Truhlar, D. G.; Schwenke, D. W.; Kouri, D. J. J. Phys. Chem. 1990, 94, 7074. (81) (a) Boonenberg, C. A.; Mayne, H. Chem. Phys. Lett. 1984, 100, 67. (b) Sathyamurthy, N.; Toennies, J. P. Chem. Phys. Lett. 1988, 143, 323. (c) Aoiz, F. J.; Herrero, V. J.; Sa´ez Ra´banos, V. Chem. Phys. Lett. 1989, 161, 270. (82) Marcus, R. A. Ber. Bunsen-Ges. Phys. Chem. 1977, 81, 190. (83) (a) Levine, R. D. J. Phys. Chem. 1990, 94, 8872. (b) Kornweitz, H.; Persky, A.; Levine, R. D. Chem. Phys. Lett. 1986, 128, 443. (c) Kornweitz, H.; Persky, A.; Schechter, I.; Levine, R. D. Chem. Phys. 1990, 169, 489. (84) Loesch, H. Chem. Phys. 1986, 104, 213; 1987, 112, 85. (85) (a) Mayne, H. R. Chem. Phys. Lett. 1986, 30, 249. (b) Harrison, J. A.; Mayne, H. R. J. Chem. Phys. 1988, 88, 7424. (c) Harrison, J. A.; Mayne, H. R. Chem. Phys. Lett. 1989, 158, 489. (86) (a) Chatfield, D. C.; Friedman, R. S.; Truhlar, D. G.; Schwenke, D. W. Faraday Discuss. Chem. Soc. 1991, 91, 289. (b) Chatfield, D. C.; Friedman, R. S.; Schwenke, D. W.; Truhlar, D. G. J. Phys. Chem. 1992, 96, 2414. (c) Chatfield, D. C.; Friedman, R. S.; Mielke, S. L.; Schwenke, D. W.; Lynch, G. C.; Allison, T. C.; Truhlar, D. G. In Dynamics of Molecules and Chemical Reactions; Wyatt, R. E., Zhang, J. Z. H., Eds.; Marcel Dekker: New York, in press. (87) Goetting, R.; Herrero, V.; Toennies, J. P.; Vodegel, M. Chem. Phys. Lett. 1987, 137, 524.
JP9527822
