Adiabatic Connection for Kinetics

Adiabatic Connection for Kinetics

(Parte 1 de 2)

Adiabatic Connection for Kinetics

Benjamin J. Lynch, Patton L. Fast, Maegan Harris, and Donald G. Truhlar*

Department of Chemistry and Supercomputer Institute, UniVersity of Minnesota, Minneapolis, Minnesota 55455-0431

ReceiVed: February 7, 2000; In Final Form: April 1, 2000

A new hybrid Hartree-Fock-density functional (HF-DF) model called the modified Perdew-Wang 1-parameter model for kinetics (MPW1K) is optimized against a database of 20 forward barrier heights, 20 reverse barrier heights, and 20 energies of reaction. The results are compared to other hybrid HF-DF methods with the 6-31+G(d,p) basis. The new method reduces the mean unsigned error in reaction barrier heights by a factor of 2.4 over MPW1PW91 and by a factor of 3 over B3LYP.

1. Introduction

Various methods, both ab initio and semiempirical, are available to calculate the barrier heights of chemical reactions. However, most of these methods have significant systematic or unsystematic errors, e.g., ab initio methods tend to overestimate barrier heights, and density functional theory (DFT) methods tend to underestimate them, whereas semiempirical NDDO methods are generally unreliable for quantitative work. The goal of this project is to devise a method that will produce increased accuracy for transition-state energies.

The approachwe will take is motivatedby the work of Becke1 on hybrid Hartree-Fock density functional theory, as justified by the adiabatic connection formula.2 This theory involves mixing various amounts of the Hartree-Fock (HF) nonlocal exchange operator3 with local DFT exchange-correlation functionals4 and gradient-corrected density functions,5-7 and the parameters have been optimized empirically against a database8 of thermochemicaldata. The originalB3PW91methodof Becke et al.,1 and the followingB3LYP methoddevelopedby Stephens et al.9 have both demonstrated remarkably high performance/ cost ratios for calculatingaccuratemolecularstructures,frequencies, and energetics.10 Experience has shown that local density functionals exhibit systematic overbinding, and gradient-corrected density functionals also exhibit systematic (albeit much smaller) overbinding, but the empirical 3-parameter hybrids eliminate most of the systematic error in binding energies. The most important parameter in these methods is the fraction of HF exchange, which is set1 equal to 20%.

The B3LYP and B3PW91 methods have also been successful for kinetics.1 However, the balance between HF exchange and DFT exchange that is needed to get barrier heights correct is apparently different from the fraction that yields the best results for stable molecules. Thus, it has empirically been found that a method called BH&HLYP gives more accurate barrier heights than B3LYP.12 The BH&HLYP method is very much like B3LYP, with the exception that the fraction of HF exchange is 50% (in fact, H&H denotes “half and half”). Local density functionalmethods tend to greatly underestimatebarrier heights; methods such as BLYP or BPW91 that have gradient-corrected density functionals still lead to obvious systematic underestimates.1 Methods such as B3PW91 and B3LYP do even better, on average, but still, it appears (on the basis of results in this paper) that they systematically underestimate barrier heights. Although BH&HLYP appears to be the best compromise for calculating barrier heights among currently available methods, there has been no systematic examination of whether one could do better. That is the first goal of the present work.

Although considerable experience has been gained with hybrid methods based on Becke’s gradient-corrected exchange functional combined with either the Lee-Yang-Parr (LYP) or

© Copyright 2000 by the American Chemical Society VOLUME 104, NUMBER 21, JUNE 1, 2000

10.1021/jp00497z C: $19.0 © 200 American Chemical Society Published on Web 05/09/2000

Perdew-Wang-1991 (PW91) gradient-corrected functionals, work in other quarters has led to improved functionals. It is important to note at this point that the choice of gradientcorrected exchange functional has a much greater effect on the thermochemical results than does the choice of gradientcorrected correlation functional. In our opinion, the most important practical advance in gradient-corrected exchange functionals is the recent work of Adamo and Barone13 who found excellent performance with a modified version of the Perdew-Wang gradient-corrected exchange functional, to be denoted MPW. When employed with 25% Hartree-Fock exchange and the Perdew-Wang gradient-correctedcorrelation functional, this yields a method called MPW1PW91. (In all cases, we adopt the name for a method by which it is known in the Gaussian9814 program.) In the present work, we shall compare calculations based on combining a variable amount of Hartree-Fock exchange with Becke’s gradient-corrected correlation functional or the MPW gradient-corrected exchange functional and with the LYP gradient-corrected correlation functionalor the PW91 gradient-correctedcorrelationfunctional.

A longstandingdifficulty in comparing or validating theoretical methods for the prediction of transition state barrier heights has been that so few experimental barrier heights are accurately known. The literature is filled with tables listing “experimental” values of the barrier heights in which the value actually listed is an Arrhenius activation energy. Experience with the most reliable available dynamical methods for predicting Arrhenius activation energies from given potential energy surfaces has shown that the two quantitiesoften differ by 2 kcal/molor more. In fact, the Arrhenius energy of activation is intrinsically temperature dependent so such comparisons even depend on the temperature range over which experiments were carried out to yield an Arrhenius activation energy. In the present paper, we make an attempt to create a database of accurate transition state barrier heights. These are sure to be inaccurate to some extent, but we hope that our attempt will stimulatefurther efforts to understand just how reliable our best estimates of barrier heights for prototype reactions really are.

Despite the fact that our barrier height database is not perfect, we think it is good enough to use in a constructive way. With this in mind, we optimize a 1-parameter DFT model especially for kinetics. The resulting model will be called MPW1K (modified Perdew-Wang 1-parameter model for kinetics) or, for completeness, when used with the basis set recommended here, MPW1K/6-31+G(d,p).

Section 2 presents our kinetics database. Section 3 tests existing methods against these database. Section 4 presents our new method.

2. Database

In this section,we assemblea databaseof zero-point-exclusive

Born-Oppenheimer barrier heights, henceforth called classical barrier heights.The best estimatesof the classicalbarrier heights for the forward reactions were estimated two ways. For H + tions.15-19 For the other reactions, we used a combination of experimental rate constraints and dynamical simulations.20-41 In all cases, the classical barrier height of the reverse reaction was calculated from the forward classical barrier height and the zero-pointexclusiveenergy of reaction,¢E, which is also called the classical endoergicity. The energy of reaction may be calculated as the difference in atomization energies of the reactants and the products.42

In most cases, the best theoretical estimate of the barrier height was made by comparing the reaction rate constants at 600 K for the best theoretical calculation and the best experi- ment. (The barrier heights for the Cl + CH4,N H + C2H6, and

NH + CH4 reactions were adjusted to experiment at 500, 1000, and 10 K, respectively.)In all cases, the experimentalreaction rate constant was determined from the experimentalists’ best 2-parameter or 3-parameter Arrhenius fit to their data. In adjusting the theoretical barrier height to make a dynamics calculation agree with experiment, we assume that all of the error in the calculated reaction rate constant comes from the barrier height. The best estimate of the classical barrier is determined by where x denotes either forward (f) or reverse (r) reactions, Vxq(theory) is the theoretical classical barrier height of the

TABLE 1: Values for Calculating the Forward and Reverse Vq(best estimate) and the Values of Vq(best estimate) for the Forward and Reverse Reactionsa theory experiment best estimate reaction ref Vfq kf ref kf ¢EV fq Vr q a T ) 600 K, except as noted differently in text; energies given in kcal/mol and rate constants given in cm3molecule-1s-1 .

Vxq(best estimate) ) Vxq(theory) + ¢Vq (1)

4812 J. Phys. Chem. A, Vol. 104, No. 21, 2000 Letters

potential surface used for the dynamical calculation, and ¢Vq is the adjustmentto the theoreticalbarrierheight.The adjustment to the barrier height is calculated from

where ktheory and kexperiment are the theoretical and experimental reaction rate constants at a given temperature T, respectively,

and R is the molar gas constant. The classical barrier height for the reverse reaction is calculated by adding the reaction exoergicity to the best estimate for the forward barrier height. The exoergicity for each of the reactions, with the exception of 15, was calculated from experimental total atomization energies.42 The exoergicity of reaction 15 is taken as the average of two theoretical values.16,43 Table 1 gives the values necessary for calculatingthe forward and reverse Vxq(best estimate) and also the actual values of

Vxq(best estimate) for the forward and reverse reactions. Note that the reactions in Table 1 are listed in the direction for which theory and experiment were compared for the rate constant. However, in subsequent tables of this paper, they are all listed in the direction predicted to be exoergic by the MPW1K method.

3. Tests of Existing Methods

Table 2 lists the forward and reverse barrier heights and classical exoergicity for the 20 test reactions as calculated by theBH&HLYP,B3LYP,andMPW1PW91methods,in allcases using the 6-31+G(d,p)4 basis set. We found that this basis set is significantly more accurate than its 6-31G(d,p) subset for the present problem. For Table 2 and throughout the whole paper, all of the saddle points were verified to have one and only one imaginaryfrequency.Beforecomparingto experiment,we added spin-orbit coupling using standard experimental values45 for reactant and product radicals. The first three rows of Table 3 give the mean signed error (MSE), mean unsignederror (MUE), and root-mean-square error (RMSE) averaged over all 60 data values (40 classical barrier heights and 20 classical energies of reaction).All of the mean errorswere computedfrom unrounded results containing spin-orbit effects.

The MPW1K rows of Table 3 will be explained in a later section, in which we will also discuss the other results in Table 3. The only point we need to make in the present section concernsthe resultsfor the 20 energiesof reaction.The MPW1PW91 method has much lower mean unsigned errors and root-meansquare errors for energies of reaction than does either of the methods based on Becke’s exchange functional and the LYP correlation function. We conclude from this that MPW1PW91 is a more balanced functional, and for this reason, we choose it as a starting point for further work.

4. Parametrization of New Method

The one-parameter hybrid Fock-Kohn-Sham operator can be written as follows:46 where FH is the Hartree operator (i.e., the nonexchange part of the Hartree-Fock operator),FHFE is the Hartree-Fock exchange operator, X is the fraction of Hartree-Fock exchange, FSE is Slater’s local density functional for exchange, FGCE is the gradient correction for the exchange functional, and FC is the total correlation functional including both local and gradientcorrected parts. Setting X ) 0.25 yields the MPW1PW91

TABLE 2: Barrier Heights and Energies of Reaction (kcal/mol) BH&HLYP B3LYP MPW1PW91 reaction Vforq Vrevq ¢EV forq Vrevq ¢EV forq Vrevq ¢E

TABLE 3: Mean Errorsa (kcal/mol) fraction HF MSE MUE RMSE a All theoretical values include spin-orbit effects.

Letters J. Phys. Chem. A, Vol. 104, No. 21, 2000 4813

method. Instead, we will choose the value that minimizes the root-mean-square error of the 60 data in the database.

Note that we explored various “2K” and “3K” models (i.e., models with 2 or 3 parameters, instead of 1), but varying the coefficients of other terms in (3) is either nonphysical or does not significantly improve the error.

The optimization process was carried out iteratively. We startedwith MPW1PW91geometriesfor the reactants,products, and transition states and found the optimum value of X. Then we re-optimizedthe geometrieswith this value of X and so forth until the method converged to X ) 0.428. All parameter optimizations were carried out with inclusion of spin-orbit effects.

The keywordsrequiredto carryout MPW1Kcalculationswith

Gaussian98 are: #mpwpw91 IOp(5/45)10000428) IOp(5/46)05720572) IOp(5/47)10001000) opt

5. Results and Discussion

Table 4 gives the results with the optimized value of X, and

Table 3 compares mean values of the errors to the previous methods. (Tables 2 and 4 give the results without spin-orbit terms, but the spin-orbit effects were included before calculating the errors in Table 3.) The barriers tend to be somewhat underestimatedby the new method,but the mean unsignederrors and root-mean-square errors in barrier heights (40 data) and overall (60 data) are significantly better than previous methods. The improvement in barrier heights compared to MPW1PW91 is a factor of 2-3, at the cost of making energies of reaction about 30% worse. The mean unsigned error in barrier heights is a factor of 3 smaller than B3LYP, which has, to date, been by far the widely used density-functional-based method for kinetics.1

Table 3 shows a very obvious correlation of the mean signed error in barrier heights with the fraction X of Hartree-Fock exchange. The value of X that minimizes the mean unsigned error and root-mean-squarederrors does not yield a mean signed error of zero, but it would nevertheless be the preferred value for applications to new problems where the classical barrier height is not known.

Table 3 does not directly test the transition state geometries.

However, as a consequence of Hammond’s postulate,47 the quality of transition state geometries should be correlated with the quality of the reaction energies and barrier heights.

We note that the 6-31+G(d,p) basis set is quite affordable even for very large molecules. Thus, the new method should be the option of choice for a wide variety of applications in kinetics.

6. Concluding Remarks

This paper providesan attemptto create a systematicdatabase of classical barrier heights that can be used for testing electronic structure methods for use in predicting chemical reactivity. It also attempts to quantify the systematic errors in density functional theory for predicting barrier heights. Finally, we propose a new hybrid Hartree-Fock-density functional theory involving a single parameter that predicts significantly more accurate barrier heights, on average, than other widely used methods.

(Parte 1 de 2)