Scaling of classical rate constants on scaled potential-energy surfaces

0 downloads 0 Views 94KB Size Report
Jun 22, 2001 - Myung Soo Kim,a) Sang Tae Park, Bong June Sung, and Jeong Hee Moon ...... 6 B. H. Yang, H. T. Gao, K. L. Han, and J. Z. H. Zhang, J. Chem.
JOURNAL OF CHEMICAL PHYSICS

VOLUME 114, NUMBER 24

22 JUNE 2001

ARTICLES Scaling of classical rate constants on scaled potential-energy surfaces Myung Soo Kim,a) Sang Tae Park, Bong June Sung, and Jeong Hee Moon National Creative Research Initiative for Control of Reaction Dynamics and School of Chemistry, Seoul National University, Seoul 151-742, Korea

共Received 18 October 2000; accepted 6 March 2001兲 The scaling relation for the classical rate constants on the scaled potential-energy surfaces has been derived using the scaling theorem in classical dynamics reported previously. This applies to the classical rate constants, both for unimolecular and for bimolecular reactions, that can be obtained by the classical trajectory method and the transition state theory. Validity of the theory has been tested for the prototype reactions, H2 CO→H2 ⫹CO and Cl⫹H2 →HCl⫹H. Exact scaling of the rate constants obtained by the classical trajectory calculations has been demonstrated. The rate-energy relations for the former reaction calculated with the statistical Rice–Ramsperger–Kassel–Marcus theory also displayed excellent scaling in the high-energy limit. The scaling relation does not hold rigorously near the reaction threshold due to the quantum mechanical zero-point energy effect. Regardless, the order of magnitude prediction of the threshold rate constant by scaling was possible even in extreme cases. The present method may allow reliable prediction of the classical rate constant by using potential energy data obtained at moderately high levels of electronic structure calculation. © 2001 American Institute of Physics. 关DOI: 10.1063/1.1374938兴

I. INTRODUCTION

Accurate calculation of the rate constant of a reaction is one of the central subjects in molecular reaction dynamics.1,2 This forms the basis for reliable analysis of experimental data, design of kinetic experiments, and prediction of the chemical fate of a system. Various theories have been developed for this purpose such as the statistical theories1–3 represented by the transition state theories and the methods based on classical or quantal dynamics calculation.4–7 When the transition state theories such as the microcanonical Rice– Ramsperger–Kassel–Marcus 共RRKM兲 theory3 are used for the calculation of a rate constant, the energy at the transition structure 共saddle point兲 relative to that of the reactants and vibrational frequencies of these structures are needed. Structure and energetics information needed for a dynamics calculation, either classical or quantal, is more extensive, which is the global potential-energy surface 共PES兲 in the configuration space spanning the reactants, transition structure, and products.4–10 Difficulty in obtaining reliable structure and energetics data is one of the main stumbling blocks in the development of theoretical methods to calculate rate constants. One often invokes chemical intuition1 to predict the magnitudes of properties which cannot be measured such as the vibrational frequencies of the transition structure or adopts an analytic PES11–13 with parameters adjusted to be compatible with the available data for the system. These introduce unwanted arbitrariness in the calculation of the rate constant. The most systematic way to obtain the information needed for the statistical or dynamical calculation of a rate a兲

Author to whom correspondence should be addressed. Electronic mail: [email protected]

0021-9606/2001/114(24)/10583/8/$18.00

constant is to use data from electronic structure calculation. Structure and vibrational frequencies at the saddle point can be calculated at moderately high levels of theory. Several schemes are available to construct a global PES from quantum chemical data.7–17 One of the successful of these, which was pioneered by Ischtwan and Collins,18–20 is to construct a PES as the weighted sum of local harmonic potentials at many configuration points obtained from quantum chemical data. Various successful cases15,16,21–30 have been reported in which global analytic PESs were constructed using the above schemes which are excellent representations of the actual quantum chemical surfaces. Even though the calculations mentioned so far can be carried out systematically, one do not expect at the moment that the rate constants calculated with quantum chemical information would be accurate enough for quantitative comparison with the experimental ones. Here the main difficulty is that one cannot obtain correct energetics data by electronic structure calculation at moderately high levels of theory. Improving the energetics accuracy a little bit by elevating the level of theory usually leads to a tremendous increase in the computation time, which may render the whole idea impractical. Recently, we reported the scaling theorem31 for classical trajectories on scaled potential-energy surfaces which states that the trajectories on each surface initiated at equivalent phase space points remain equivalent after the scaled time. When applied to the partitioning of the available energy to the products, this meant that the product mode-specific energies on scaled PESs were simply scaled. Then, it would be possible to calculate these quantities reliably on PES obtained on a moderately high level of theory if this surface resembled the real surface qualitatively, even if not quantitatively. A successful

10583

© 2001 American Institute of Physics

Downloaded 30 Jun 2001 to 147.46.44.94. Redistribution subject to AIP license or copyright, see http://ojps.aip.org/jcpo/jcpcr.jsp

10584

Kim et al.

J. Chem. Phys., Vol. 114, No. 24, 22 June 2001

application of the scaling theorem to the product energy partitioning was demonstrated for the photodissociation of formaldehyde in our recent report.30 Since the scaling theorem holds generally in classical dynamics, it is expected that a simple relation exists among classical rate constants on scaled PESs also. Presentation of this relation and demonstration of its validity using the following reactions as examples are the main subjects of this paper: H2CO→H2⫹CO,

共1兲

Cl⫹H2→HCl⫹H.

共2兲

It will also be demonstrated that scaling of the classical rate constant established in this work will be useful to correct the dynamics results calculated on the ab initio PESs. Compared to the previous approaches which adjust PESs to comply with the experimental data,32–35 the present method may have the advantage of being simple and systematic.

II. THEORETICAL

FIG. 1. Schematic potential-energy diagrams for reaction 共1兲, 共a兲 unscaled 共s⫽1.0兲 and 共b兲 scaled 共s兲. Energy barrier in 共b兲 is higher than in 共a兲 by a factor of s 2 . Accordingly, the total internal energy for a trajectory on 共b兲, E s , which is equivalent to that on 共a兲 is larger than the latter by the same factor, E s ⫽s 2 E. k(E) and k s (E s ) are the rate constants for the two equivalent cases.

The scaling theorem in classical dynamics reported31 previously can be summarized as follows. Let us consider two PESs related by simple scaling: V s 共 qs 兲 ⫽s 2 V 共 q兲 .

共3兲

Two phase space points, one on each surface, will be called equivalent when the following conditions are met for the positions and momenta: qs ⫽q,

共4兲

ps ⫽sp.

共5兲

Then, the total energies at these points are related by simple scaling also: E s 共 qs ,ps 兲 ⫽s 2 E 共 q,p兲 .

共6兲

It was shown previously that if the equivalence is met at time t 1 and t s1 on V and V s surfaces, respectively, the equivalence is maintained after the equivalent times

␦ t s ⫽s ⫺1 ␦ t.

N s ⫽N 0 e ⫺k s t s .

共11兲

N and N s are the same because the trajectory equivalence means that the numbers of unreactive trajectories in dynamics on scaled PESs are the same if counted at the equivalent times, namely t s ⫽t/s. Hence, comparing Eqs. 共10兲 and 共11兲, one obtains the scaling relation for the rate constant: k s ⫽s•k.

共12兲

It is to be emphasized that this relation holds when not only the PESs, Eq. 共3兲, but also the initial positions and momenta are properly scaled, Eqs. 共4兲 and 共5兲, and hence the initial internal energies also. The situation is drawn graphically in Fig. 1. In the case of bimolecular reaction, the trajectory equivalence means that the reaction cross sections on the scaled PESs are the same. Namely,

共7兲

k⫽u ␴

共13兲

k s ⫽u s ␴ .

共14兲

Namely, qs 共 t s1 ⫹ ␦ t s 兲 ⫽q共 t 1 ⫹ ␦ t 兲

共8兲

ps 共 t s1 ⫹ ␦ t s 兲 ⫽sp共 t 1 ⫹ ␦ t 兲 .

共9兲

and

Hence, the classical trajectories on scaled PESs are equivalent if initiated at the equivalent phase space points. To see the effect of the trajectory equivalence on the rate constant, let us first consider the calculation of the rate constant of a unimolecular reaction from N 0 trajectory results. Under the random lifetime assumption3,7 the rate constant 共k兲 is related to the number of unreactive trajectories 共N兲 at time t as follows: N⫽N 0 e ⫺kt , and, for trajectories on the scaled PES

共10兲

and

Here, u and u s are the relative speed of reactants. Scaling of the momentum, Eq. 共5兲, means scaling of the relative speed, namely u:u s ⫽1:s. Then, the scaling relation between the rate constants can be derived as follows, which is the same as for the unimolecular case: k s ⫽u s ␴ ⫽su ␴ ⫽s•k.

共15兲

The same scaling relation can be derived also for the transition state theory rate constant in the classical limit. As an example, let us derive the scaling relation for the microcanonical transition state theory rate constant in the classical limit which has the form36,37

Downloaded 30 Jun 2001 to 147.46.44.94. Redistribution subject to AIP license or copyright, see http://ojps.aip.org/jcpo/jcpcr.jsp

J. Chem. Phys., Vol. 114, No. 24, 22 June 2001 E⫺E 0

k⫽

兰0

Scaling of classical rate constants

‡ ‡ dE ‡1 兰 H ‡ ⫽E⫺E ‡ ⫺E 0 dq ‡2 ¯dq 3N dp ‡2 ¯d p 3N 1

兰 H⫽E dq 1 ¯dq 3N dp 1 ¯dp 3N

. 共16兲

Here E is the total energy, E 0 is the critical energy, and E ‡1 is the translational energy in the reaction coordinate 共1兲 which

k s共 E s 兲 ⫽

is assumed to be separable from the remaining degrees of freedom 共2, 3, . . . ,3N兲. The surface integral in the numerator divided by h 3N⫺1 is the density of states for the 3N⫺1 degrees of freedom orthogonal to the reaction coordinate at the transition state. Similarly, the integral in the denominator is the density of reactant states multiplied by h 3N . The rate constant on the scaled PES, V s , at the scaled total energy (E s ) can be written similarly,

E ⫺E s0 ‡ ‡ ‡ ‡ ‡ dE s1 兰 H ‡ ⫽E s ⫺E ‡ ⫺E s0 dq s2 ¯dq s3N d p s2 ¯d p s3N s s1

兰0 s

兰 H s ⫽E s dq s1 ¯dq s3N dp s1 ¯d p s3N

10585

共17兲

.

Upon coordinate transformation 关Eqs. 共4兲 and 共5兲兴 and using the energy scaling relations 关Eqs. 共3兲 and 共6兲兴, the denominator becomes s 3N



H⫽E

dq 1 ¯dq 3N dp 1 ¯dp 3N .

共18兲

The second surface integral in the numerator is converted similarly with the multiplication factor of s 3N⫺1 . Finally, the energy scaling relation converts the first integral in the numerator into the form in Eq. 共16兲 with the multiplication factor of s 2 . Namely E⫺E 0

k s共 E s 兲 ⫽

s 2兰 0

‡ ‡ dE ‡1 •s 3N⫺1 兰 H ‡ ⫽E⫺E ‡ ⫺E 0 dq ‡2 ¯dq 3N d p ‡2 ¯d p 3N 1

s

3N

兰 H⫽E dq 1 ¯dq 3N dp 1 ¯d p 3N

This is the same relation as derived for the classical trajectory rate constant 关Eqs. 共12兲 and 共15兲兴. Derivation of the scaling relation for the canonical rate constant in the classical limit is also straightforward. In this case, the scaled temperature (T s ) must be used due to the Boltzmann factor appearing in the conversion from the microcanonical to canonical rate constant.1–3 Namely k s 共 T s 兲 ⫽sk 共 T 兲 ,

共20兲

with T s ⫽s T.

共21兲

2

2

When the potential energy is scaled by s , the vibrational frequencies are scaled by s, which results in the scaling of the zero-point energies by s. Namely, consideration of the zero-point energies disturbs the equivalence condition and hence the simple scaling relation does not hold rigorously for the quantum mechanical version of the transition state theory rate constant. The quantum effect will be especially prominent near the reaction threshold, which will be investigated using reaction 共1兲 as an example. III. COMPUTATIONAL

To corroborate the scaling relation derived in the previous section, the rate constants for the unimolecular reaction, reaction 共1兲, and the cross sections for the bimolecular reaction, reaction 共2兲, have been calculated by classical trajectory method. PESs for reaction 共1兲 were constructed by interpolation in our previous investigation of scaling of the product mode-specific energies.30 Since calculation of the rate con-

⫽sk 共 E 兲 .

共19兲

stant in the present work requires a better description of the entrance region of PES, a new PES has been constructed in this work using the coordinate system appropriate for the reactant region. For the H1H2CvO system, this includes the ⫺1 ⫺1 ⫺1 ,r CH , and r CH , inverses of three interatomic distances, r CO 1 2 two bond angles, ␪ H1CO and ␪H2CH1, and the sine function of a dihedral angle, ␾ H2H1CO,



␾H2H1CO⫽cos⫺1

共 rH2H1 ⫻rCH1 兲 • 共 rCH1 ⫻rCO兲

兩 rH2H1⫻rCH1 兩兩 rCH1 ⫻rCO兩



.

共22兲

Energy, gradient, and hessian were calculated on the Hartree–Fock 共HF兲 level using the 6-31G** basis set at 60 configuration points located on the reaction path or intrinsic reaction coordinate 共IRC兲. Limiting the electronic structure calculation to a moderate level would not cause any problem here because the purpose of the present work is just to confirm the validity of the scaling relation. Data at each point were used to obtain a local harmonic potential in the vicinity of this point using the coordinate system mentioned above. Then, the global PES was constructed as the weighted sum of these local potentials.23 Since the two hydrogen atoms are not equivalent in the transition structure, while they are in the reactant, there exist two equivalent reaction paths with equivalent transition structures. To account for this, mirror images of the above 60 local potentials were added, resulting in the interpolated PES consisting of 120 local potentials. Unlike our previous studies on reaction dynamics,23,30 no attempt was made in this work to improve the global PES by adding off-IRC local potentials for the same reason as men-

Downloaded 30 Jun 2001 to 147.46.44.94. Redistribution subject to AIP license or copyright, see http://ojps.aip.org/jcpo/jcpcr.jsp

10586

Kim et al.

J. Chem. Phys., Vol. 114, No. 24, 22 June 2001

TABLE I. Rate constants and cross sections for reactions 共1兲 and 共2兲 obtained by classical trajectory calculations on unscaled 共s⫽1兲 and scaled 共s⫽0.5 and 2兲 PESs using scaled internal energies.

Scale factor Reaction共1兲

1.0 0.5 1.5 2.0 Scale factor

Reaction共2兲

1.0 0.5 2.0

Barriera (kcal mol⫺1 )

Internal energy (kcal mol⫺1 )

Rate constant (1011 s⫺1 )

140.0 35.0 315.0 560.0

1.78 0.89 2.67 3.56

104.6 26.15 235.35 418.4 Barrierb (kcal mol⫺1 )

Internal energyc (kcal mol⫺1 )

Cross section (Å2)

Rate constant (10⫺11 cm3 molecule⫺1 s⫺1 )

26.2 6.55 104.8

3.77 3.77 3.77

1.12 0.56 2.24

7.700 1.925 30.80

The forward barrier on the unscaled 共s⫽1兲 PES was calculated by HF/6-31G**. From LEPS PES, Ref. 44. c For the unscaled PES, this includes 20 kcal mol⫺1 of collision energy and 6.2 kcal mol⫺1 of H2 vibrational energy. a

b

tioned above. The initial conditions for trajectory runs were randomly selected at the rotational temperature of ⬃0 K by orthantlike momentum sampling38–40 with the molecular geometry fixed at the minimum energy geometry of the reactant. The classical total energy of 140 kcal mol⫺1 above the reactant energy minimum was used for trajectory calculations on the unscaled PES. For each initial condition sampled, atomic momenta were scaled according to Eq. 共5兲 to obtain the equivalent initial condition for the corresponding trajectory on the scaled PES. The fourth-order Runge– Kutta method41 was used for integration with the time step fixed at 0.01 fs. Integration of a trajectory was terminated when two fragments were 4.5 Å apart or when 10 ps elapsed. Trajectories with the energy conservation error, 兩 ⌬E 兩 /E, larger than 10⫺4 were discarded. These were 87 out of a total of 1050 trajectories. Number of reactive trajectories was counted at a constant time interval. Then, the rate constant was evaluated from the ln(N/N 0 ) vs t plot. Classical and quantal dynamics of reaction 共2兲 on the London–Eyring–Polany–Sato 共LEPS兲42–44 and LEPSinterpolated PES were reported previously24 together with the construction of PES at the QCISD/6-31G** level. The original LEPS PES has been adopted in this work for the same reason as mentioned previously. The trajectory calculations on the unscaled LEPS PES were essentially the same as those in the previous work24 which used 20 kcal mol⫺1 collision energy and 6.2 kcal mol⫺1 H2 vibrational energy. For calculations on the scaled PESs, all the atomic momenta in the initial condition for a trajectory on the unscaled PES were multiplied by s to obtain classically equivalent initial conditions. It is to be mentioned, however, that two classically equivalent initial conditions are not quantum mechanically equivalent because of the difference in the scaling factor between the energy and the vibrational frequency. Integration was carried out with the fourth-order Runge– Kutta method with the time step of 0.01 fs. The energy conservation error in each trajectory was less than 10⫺5 . 3000 trajectories were run, about half of which were reactive. The cross section was evaluated from the trajectory results using the following formula:44

2 ␴ ⫽ ␲ b max 共 N r /N 兲 .

共23兲

Here b max is the maximum impact parameter used in the trajectory calculation, N is the total number of trajectories, and N r is the number of reactive trajectories. To test the validity of the scaling relation for the transition state theory rate constants, the RRKM theory rate constants were calculated for reaction 共1兲 using the following equation without tunneling correction:3 k共 E 兲⫽␴

W ‡ 共 E⫺E 0 兲 . h␳共 E 兲

共24兲

Here, ␴ is the reaction path degeneracy, h is the Planck constant, E 0 is the forward critical energy, ␳ is the density of states of the reactant, and W ‡ is the state sum at the transition structure. The usual harmonic approximation for the oscillators was adopted. The critical energy and harmonic frequencies were taken from the HF/6-31G** results. For the calculation on the scaled PES, both E and E 0 were multiplied by s 2 while the harmonic frequencies were multiplied by s. IV. RESULT AND DISCUSSION

Before calculating rate constants, we checked the equivalence of trajectories on the scaled PESs initiated at the equivalent phase space points for reaction 共1兲. A trajectory on the original PES 共s⫽1兲 was taken as the reference. Equivalent trajectories on two scaled PESs with s⫽0.5 and 2.0 were calculated. Since the classical barrier for this reaction on the original HF/6-31G** PES is 104.6 kcal mol⫺1 , the corresponding barriers on the s⫽0.5 and 2.0 PESs are 26.15 and 418.4 kcal mol⫺1 , respectively 共Table I兲. Similarly, scaled reactant internal energies were used, for example, 35.0, 140.0, and 560.0 kcal mol⫺1 for calculations on s⫽0.5, 1, and 2 surfaces, respectively. Trajectory calculations were carried out for 20, 10, and 5 ps for s⫽0.5, 1, and 2 surfaces, respectively. This was to match the equivalent time scales. To test the equivalence of trajectories, atomic positions and momenta on these trajectories at the equivalent times were taken at every 100 integration step 共1 fs when

Downloaded 30 Jun 2001 to 147.46.44.94. Redistribution subject to AIP license or copyright, see http://ojps.aip.org/jcpo/jcpcr.jsp

J. Chem. Phys., Vol. 114, No. 24, 22 June 2001

Scaling of classical rate constants

10587

TABLE II. Molecular parameters used for the calculation of the RRKM rate constant for reaction 共1兲.a Vibrational frequencies

Classical barrier

E0

104.6

98.8

3199 2025 1377

3123 1668 1335

3275 1430 797

2126 1057

Millerc

93.6

87.8

2843 1746 1247

2766 1501 1164

2769 1137 697

1654 941

Troed

88.5

83.5

3009 1764 1288

2944 1563 1191

3125 1523 839

1830 936

HF/6-31G**

b

Reactant

TS

Energies in kcal mol⫺1 . Vibrational frequencies in cm⫺1 . Forward critical energy taking into account the zero-point energies at the reactant and TS. c Reference 45. d Reference 46. a

b

FIG. 2. The semilog plots (¯) of the fractions (N/N 0 ) of the unreactive trajectories vs time for reaction 共1兲. 共a兲 is for the trajectories on the interpolated PES constructed at the HF/6-31G** level. The same surface was scaled to obtain the other results. The s factors used were 共b兲 0.5, 共c兲 2, and 共d兲 1.5. Lines 共—兲 are the least squares fits to Eqs. 共10兲 and 共11兲. The rate constants evaluated with the slopes of the lines are indicated.

s⫽1兲 and 兩 q1 ⫺qs 兩 / 兩 q1 兩 and 兩 p1 ⫺ps /s 兩 / 兩 p1 兩 were calculated. Here, q1 and p1 are the atomic positions and momenta for the trajectory on the original PES and qs and ps are those on the scaled PES. These were exactly zero throughout the trajectory regardless of the computational precision adopted. This is the computational proof of the scaling theorem reported previously.31 Trajectory equivalence was not maintained for a prolonged time, however, when the s factor other than 2 n was used. We think that this was caused by chaos originating from the numerical roundoff error in computation. Breakup of the trajectory equivalence due to the numerical chaos measured in terms of 兩 q1 ⫺qs 兩 / 兩 q1 兩 and 兩 p1 ⫺ps /s 兩 / 兩 p1 兩 got serious (⬎10⫺5 ) after several picosecond even though the energy conservation was maintained. This will be discussed in a later publication. The semilog plots of the number of unreactive trajectories vs time for reaction 共1兲 on the unscaled 共s⫽1兲 and scaled 共s⫽0.5 and 2兲 PESs are nearly straight lines as shown in Fig. 2. The rate constant on the unscaled surface at the classical energy of 140 kcal mol⫺1 determined from the slope of the plot was 1.78⫻1011 s⫺1 . The corresponding values on the s⫽0.5 and 2 surfaces were 0.89⫻1011 and 3.56⫻1011 s⫺1 , respectively. k s⫽1 :k s⫽0.5 :k s⫽2 of 1:0.5:2 is exactly as predicted from the scaling relation, Eq. 共12兲. The results are summarized in Table I. We also calculated the rate constant on the s⫽1.5 surface. The rate constant of 2.67⫻1011 as determined in Fig. 2 satisfies the scaling relation

k s⫽1 :k s⫽1.5⫽1:1.5 even though only up to the third significant digit in this case. Namely, use of many trajectories virtually eliminated the effect of the numerical chaos in the rate constant determination. The cross sections for reaction 共2兲 on s⫽1, 0.5, and 2 surfaces obtained from the trajectory results are listed in Table I together with the reaction barriers on these PESs and internal energies used in the calculation. Since equivalent trajectories on the three surfaces were initiated at the equivalent phase space points, these cross sections are exactly the same as was mentioned in the theoretical section. The rate constants listed in the same table display the s scaling 关Eq. 共15兲兴 due to the scaling of the relative speed. Persky44 carried out the classical trajectory study of this reaction using the same LEPS PES as adopted in this work and reported the cross sections in the center-of-mass collision energy range of 0–12 kcal mol⫺1 . Since 20 kcal mol⫺1 collision energy adopted in this work lies outside of this range, we cannot compare the present result with Persky’s quantitatively. Judging from the general trend of the collision energy dependence of the cross sections in Persky’s work, the cross section of 3.77 Å2 obtained in this work seems to be in good agreement with Persky’s results. The RRKM rate constants for reaction 共1兲 were also calculated on the HF/6-31G** and its scaled PESs. The molecular parameters at the HF/6-31G** level needed for the calculation are listed in Table II. With scaling, the classical forward barrier, which is the energy difference between the energy minima of the saddle point and reactant, is multiplied by s 2 . Such a simple scaling relation does not hold for the quantum-mechanical barrier which is the energy difference between the zero-point levels because vibrational frequencies are scaled by s. In the usual use of Eq. 共24兲, E 0 is the quantum-mechanical forward barrier and E is the internal energy above the zero-point level of the reactant. Namely, the quantum-mechanical zero-point energy effect is inherent in Eq. 共24兲. Figure 3共a兲 shows the rate-energy data calculated with this equation on the original 共s⫽1兲 and scaled 共s⫽1.2 and 2兲 surfaces. As usual, the abscissa in this figure is the

Downloaded 30 Jun 2001 to 147.46.44.94. Redistribution subject to AIP license or copyright, see http://ojps.aip.org/jcpo/jcpcr.jsp

10588

J. Chem. Phys., Vol. 114, No. 24, 22 June 2001

FIG. 3. 共a兲 The RRKM rate-energy relations for reaction 共1兲. The solid curve 共—兲 was calculated with HF/6-31G** data in Table II. The dots (¯) and dash–dots (⫺•⫺•) were obtained using parameters scaled with the s factor of 1.2 and 2.0, respectively. 共b兲 The scaled RRKM rate-energy relations for reaction 共1兲 obtained by transforming the data in 共a兲. The abscissa is the classical internal energy referred to the reactant energy minimum divided by s 2 . The ordinate is the logarithm of the rate constant divided by s. The solid curve 共—兲, dots (¯), and dash–dots (⫺•⫺•) represent the results for s⫽1, 1.2, and 2.0, respectively.

internal energy above the zero-point level of the reactant. In the classical trajectory study of this reaction described above, the rate constant on the original surface 共s⫽1兲 was calculated at the classical energy of 140 kcal mol⫺1 which corresponds to the internal energy of 121.8 kcal mol⫺1 after the zeropoint energy correction. The RRKM rate constant at this energy read from Fig. 3共a兲 is 4.15⫻1011 s⫺1 which is larger than the classical trajectory rate constant of 1.78⫻1011 s⫺1 by a factor of 2.3. Harmonic treatment of the molecular vibrations, which reduces the denominator in Eq. 共24兲 more than the numerator, would be the main reason for this discrepancy. Under the circumstances, the rate constants calculated by the two different methods can be regarded as comparable. As is well known, rate constant is a very sensitive function of the critical energy. For example, as the latter increases by a factor 4, from s⫽1 to s⫽2 in Fig. 3共a兲, the rate constant decreases by orders of magnitude. As was mentioned in the introductory section, this makes it difficult to obtain accurate rate constants from energy and frequency data obtained from electronic structure calculation. The above rate-energy data are replotted in Fig. 3共b兲 to test the validity of the scaling relation. The abscissa in this

Kim et al.

figure is taken as the classical internal energy, namely the internal energy including the reactant zero-point energy, divided by s 2 to comply with the scaling relation 关Eq. 共6兲兴 while the ordinate is taken as the logarithm of the rate constant divided by s 关Eq. 共19兲兴. The figure displays an excellent overlap of the rate-energy data calculated on the s⫽1, 1.2, and 2 surfaces after scaling, especially in the high internal energy range. Namely, the scaling relation for the rate constant would hold rigorously in the classical limit. The overlap is not quite satisfactory in the low-internal energy range especially near the reaction threshold. As was explained earlier, this mismatch is due to the fact that the zero-point energies scale differently from the potential energy. At the reaction threshold, namely at the internal energy corresponding to the zero-point level at the saddle point, the scaled rate constants (k s /s) are 1.110⫻109 , 2.706⫻109 , and 4.695 ⫻109 s⫺1 , respectively, for calculations on the s⫽1, 1.2, and 2 surfaces. As the forward critical energies increase by factors of 1.44 and 4 共from s⫽1 to s⫽1.2 and 2兲, the scaled rate constants increase to 2.44 and 4.23, respectively, times the values on the unscaled surface. This means that if one can obtain the critical energy with decent accuracy, say within 50%, by improving the level of the electronic structure calculation, the rate constant can be estimated fairly accurately by utilizing the scaling relation in spite of the quantummechanical zero-point energy effect. Further disturbance of the rate constant scaling relation will occur, however, when the quantum-mechanical tunneling, which is another important quantum effect in reaction dynamics, is present. Finally, we would like to demonstrate the utility of the rate constant scaling relation by comparing with the RRKM rate constants for reaction 共1兲 reported by Miller45 and by Troe.46 Classical forward barriers and vibrational frequencies used by these investigators are listed in Table II. Some of these are the experimental data while others are values from various levels of ab initio calculation corrected and/or scaled according to experience. Namely, each set of data was not obtained from the quantum-chemical calculation at a single level, unlike our HF/6-31G** data. We will not be concerned with how these values were obtained or whether the rateenergy relations calculated with these values are compatible with the experimental results. These sets of data will be treated just as numbers that can be used to calculate the rate-energy relation. We would like to emphasize, however, that the energy and frequency scaling relations do not hold satisfactorily among these sets of data. For example, the scaling factor (s 2 ) for the Miller’s potential with respect to the HF/6-31G** one estimated from the classical forward barrier is 0.895 共⫽93.6/104.6兲 which results in the s value of 0.946. On the other hand, the s factor estimated from the ratios of the vibrational frequencies is ⬃0.86. In this work, the RRKM calculation without rotation was done to reproduce Miller’s rate-energy relation. Since the Troe’s results are available in the numerical form, they were taken as reported. Then, both the internal energies and rate constants were scaled to the HF/6-31G** results using the scaling factors (s 2 ) estimated from the classical forward barriers, 0.895 and 0.846, respectively, for Miller’s and Troe’s cases. The scaled rate-energy relations thus obtained are in good agreement

Downloaded 30 Jun 2001 to 147.46.44.94. Redistribution subject to AIP license or copyright, see http://ojps.aip.org/jcpo/jcpcr.jsp

J. Chem. Phys., Vol. 114, No. 24, 22 June 2001

FIG. 4. The scaled RRKM rate-energy relations for reaction 共1兲 obtained using Miller’s (¯) and Troe’s (⫺•⫺•) parameters are compared with the HF/6-31G** result 共—兲. The formers were scaled using the scaling factors (s 2 ) of 0.895 and 0.846, respectively.

with the HF/6-31G** one as shown in Fig. 4. For example, the scaled rate constants on Miller’s, HF/6-31G**, and Troe’s surfaces at the classical internal energy of 125 kcal mol⫺1 are 5.99⫻1010, 3.64⫻1010, and 3.86⫻1010 s⫺1 , respectively. The agreement is all the more remarkable considering that vibrational frequencies do not follow the scaling relations dictated by the differences in the classical forward barrier. This indicates that RRKM rate constants with reasonable accuracy can be obtained by using the scaling factor determined from the classical forward barriers. This is of practical importance because it is not likely, even though one would hope, that the energy barrier and vibrational frequencies obtained by electronic structure calculation scale with the actual values simultaneously, the former by the scaling factor of s 2 and the latter by s. Also to be noted is that the above agreement is achieved using the rate-energy relation obtained at HF/6-31G** as the reference which is only a moderate level of the electronic structure calculation. V. SUMMARY AND CONCLUSION

Reliable prediction of the rate constant for an energized system of molecule共s兲 is critical to determine its chemical fate. However, accurate theoretical calculation of a rate constant is a difficult task, one of the main reason being its dramatic change with the changes in the energy content of the system and the reaction critical energy. In this work, we have derived the scaling relation for the classical rate constants of reactions occurring on scaled potential-energy surfaces. Validity of this relation has been demonstrated for the rate constants of prototype unimolecular and bimolecular reactions obtained by the classical trajectory calculation and the RRKM theory. The results from the RRKM calculation have also shown that a reasonable estimation of the rate constant is possible even at the reaction threshold where the scaling relation does not hold rigorously due to the quantummechanical zero-point energy effect. The present results im-

Scaling of classical rate constants

10589

ply that a decent estimation of a rate constant is possible by calculating it on an ab initio PES and scaling it. Namely, even when PES with chemical accuracy is not available from the ab initio calculation, a rate constant may be estimated with reasonable accuracy by using the scaling relation. Theoretically, rigorous scaling of classical rate constants holds only when PESs are scaled exactly. Namely, exact scaling of an ab initio PES to the real one is the prerequisite for the exact scaling of classical rate constants. Exact scaling over the global surface may be too much to hope for, especially with the currently available electronic structure methods. However, the fact that the rate-energy relation at the HF/631G** level scaled well with the presumably reasonable data reported by Miller and by Troe indicates that the rate constant scaling holds well between data on less than exactly scaled PESs. Then, a reliable prediction of the rate-energy relation would be possible by calculating on a raw ab initio PES and scaling the result when the reaction critical energy is known experimentally. The present scheme may not be adequate, however, to treat systems with multiple reaction pathways. This is because the ab initio critical energies for all the barriers are not expected to scale uniformly with real ones, especially at the moderately high levels of the electronic structure theory. This will certainly limit the applicability of the present method. One needs to check, however, whether the global scalability of the ab initio PES will improve as the level of the theory improves. ACKNOWLEDGMENTS

This work was supported financially by CRI, the Ministry of Science and Technology, Republic of Korea. S.T.P., J.H.M., and B.J.S. thank the Ministry of Education for the Brain Korea 21 fellowship. S. W. Benson, Thermochemical Kinetics 共Wiley-Interscience, New York, 1976兲. 2 M. J. Pilling and I. W. H. Smith, Modern Gas Kinetics 共Blackwell Scientific, London, 1987兲. 3 P. J. Robinson and K. A. Holbrook, Unimolecular Reactions 共WileyInterscience, New York, 1972兲. 4 D. H. Zhang and J. Z. H. Zhang, J. Chem. Phys. 99, 6624 共1993兲. 5 D. H. Zhang and J. C. Light, J. Chem. Phys. 104, 4544 共1996兲. 6 B. H. Yang, H. T. Gao, K. L. Han, and J. Z. H. Zhang, J. Chem. Phys. 113, 1434 共2000兲. 7 L. M. Raff and D. L. Thompson, Theory of Chemical Reaction Dynamics, edited by M. Baer 共CRC, Boca Raton, 1985兲, Vol. III. 8 G. C. Schatz, Rev. Mod. Phys. 61, 669 共1989兲. 9 D. G. Truhlar, R. Steckler, and M. S. Gordon, Chem. Rev. 87, 217 共1987兲. 10 D. G. Truhlar, Potential Energy Surfaces and Dynamics Calculations 共Plenum, New York, 1981兲. 11 O. L. Polyansky, P. Jensen, and J. Tennyson, J. Chem. Phys. 105, 6490 共1996兲. 12 E. Martı´nez-Nu´˜nez and S. A. Va´zquez, Chem. Phys. Lett. 316, 471 共2000兲. 13 N. C. Handy and S. Carter, Chem. Phys. Lett. 79, 118 共1981兲. 14 T. Hollebeek, T. S. Ho, and H. Rabitz, Annu. Rev. Phys. Chem. 50, 537 共1999兲. 15 T. S. Ho, T. Hollebeek, H. Rabitz, L. B. Harding, and G. C. Schatz, J. Chem. Phys. 105, 10472 共1996兲. 16 G. C. Schatz, A. Papaioannou, L. A. Pederson, L. B. Harding, T. Hollebeek, T. S. Ho, and H. Rabitz, J. Chem. Phys. 107, 2340 共1997兲. 17 R. D. Levine and R. B. Bernstein, Molecular Reactions Dynamics and Chemical Reactivity 共Oxford University Press, New York, 1987兲. 18 J. Ischtwan and M. A. Collins, J. Chem. Phys. 100, 8080 共1994兲. 1

Downloaded 30 Jun 2001 to 147.46.44.94. Redistribution subject to AIP license or copyright, see http://ojps.aip.org/jcpo/jcpcr.jsp

10590 19

Kim et al.

J. Chem. Phys., Vol. 114, No. 24, 22 June 2001

K. C. Thompson and M. A. Collins, J. Chem. Soc., Faraday Trans. 93, 871 共1997兲. 20 K. C. Thompson, M. J. T. Jordan, and M. A. Collins, J. Chem. Phys. 108, 8302 共1998兲. 21 T. Yamamoto and S. Kato, J. Chem. Phys. 107, 6114 共1997兲. 22 T. G. Lee, S. C. Park, and M. S. Kim, J. Chem. Phys. 104, 4517 共1996兲. 23 Y. M. Rhee, T. G. Lee, S. C. Park, and M. S. Kim, J. Chem. Phys. 106, 1003 共1997兲. 24 S. Y. Lin, S. C. Park, and M. S. Kim, J. Chem. Phys. 111, 3787 共1999兲. 25 T. Ishida and G. C. Schatz, Chem. Phys. Lett. 314, 369 共1999兲. 26 K. A. Nguyen, I. Rossi, and D. G. Truhlar, J. Chem. Phys. 103, 5522 共1995兲. 27 K. C. Thompson and T. J. Martinez, J. Chem. Phys. 110, 1376 共1999兲. 28 T. Ishida and G. C. Schatz, J. Chem. Phys. 107, 3558 共1997兲. 29 T. Ishida and G. C. Schatz, Chem. Phys. Lett. 298, 285 共1998兲. 30 B. J. Sung and M. S. Kim, J. Chem. Phys. 113, 3098 共2000兲. 31 J. H. Moon, S. T. Park, and M. S. Kim, J. Chem. Phys. 110, 972 共1999兲. 32 M. Meuwly and J. M. Hutson, J. Chem. Phys. 110, 8338 共1999兲. 33 J. M. Bowman and B. Gazdy, J. Chem. Phys. 94, 816 共1991兲.

34

K. Higgins, F-M. Tao, and W. Klemperer, J. Chem. Phys. 109, 3048 共1998兲. 35 B. Gazdy and J. M. Bowman, J. Chem. Phys. 95, 6309 共1991兲. 36 J. I. Steinfeld, J. S. Francisco, and W. L. Hase, Chemical Kinetics and Dynamics 共Prentice Hall, New Jersey, 1989兲. 37 I. W. M. Smith, Kinetics and Dynamics of Elementary Gas Reactions 共Butterworths, London, 1980兲. 38 T. G. Lee, M. S. Kim, and S. C. Park, J. Chem. Phys. 104, 5472 共1996兲. 39 D. L. Bunker and W. L. Hase, J. Chem. Phys. 59, 4621 共1973兲. 40 W. L. Hase and D. G. Buckowski, Chem. Phys. Lett. 74, 284 共1980兲. 41 W. H. Press, S. A. Tenkolsky, W. T. Vettering, and B. P. Flannery, Numerical Recipes in FORTRAN, 2nd ed. 共Cambridge University Press, Cambridge, 1992兲. 42 T. Ishida and G. C. Schatz, Chem. Phys. Lett. 298, 285 共1998兲. 43 T. Takata, T. Taketsugu, K. Hirao, and M. S. Gordon, J. Chem. Phys. 109, 4281 共1998兲. 44 A. Persky, J. Chem. Phys. 66, 2932 共1977兲. 45 W. H. Miller, J. Am. Chem. Soc. 101, 6810 共1979兲. 46 J. Troe, J. Phys. Chem. 88, 4375 共1984兲.

Downloaded 30 Jun 2001 to 147.46.44.94. Redistribution subject to AIP license or copyright, see http://ojps.aip.org/jcpo/jcpcr.jsp