Direct ab initio dynamics studies of N + H2<-->NH + H reaction

Kinetics of the N+H2<-->NH+H reaction have been studied using a direct ab initio dynamics method. Potential energy surface for low electronic states have been explored at the QCISD/ cc-pVDZ level of theory. We found the ground-state reaction is N(4S)+H2-->NH(3Σ-)+H. Thermal rate constants for this reaction were calculated using the microcanonical variational transition state theory. Reaction path information was calculated at the QCISD/cc-pVDZ level of theory. Energies along the minimum energy path (MEP) were then refined at the QCISD(TQ)/ cc-pVTZ level of theory. The forward and reverse barriers of the ground-state reaction are predicted to be 29.60 and 0.53 kcal/mol, respectively. The calculated rate constants for both forward and reverse reactions are in good agreement with available experimental data. They can be expressed as k(T)52.3331014 exp(230.83 (kcal/mol)/RT) cm3 mols-1 s-1 for the forward reaction and k(T) 55.553108T1.44 exp(20.78(kcal/mol)/RT) cm3 mol-1 s-1 for the reverse reaction in the temperature range 400-2500 K.

Direct ab initio dynamics studies of N ➾ H2❫ NH ➾ H reaction Shaowen Zhang and Thanh N. Truonga) Henry Eyring Center for Theoretical Chemistry, Department of Chemistry, University of Utah, 315 S 1400 E, Room 2020, Salt Lake City, Utah 84112 ⑦ Received 1 May 2000; accepted 21 July 2000 ✦ Kinetics of the N ✶ H2 ✩ NH ✶ H reaction have been studied using a direct ab initio dynamics method. Potential energy surface for low electronic states have been explored at the QCISD/ cc-pVDZ level of theory. We found the ground-state reaction is N(4S) ✶ H2 􀀀 NH(3 ❙✷) ✶ H. Thermal rate constants for this reaction were calculated using the microcanonical variational transition state theory. Reaction path information was calculated at the QCISD/cc-pVDZ level of theory. Energies along the minimum energy path ⑦ MEP ✦ were then refined at the QCISD ⑦ TQ ✦ / cc-pVTZ level of theory. The forward and reverse barriers of the ground-state reaction are predicted to be 29.60 and 0.53 kcal/mol, respectively. The calculated rate constants for both forward and reverse reactions are in good agreement with available experimental data. They can be expressed as k(T) ✺ 2.33 ✸ 1014 exp( ✁ 30.83 ⑦ kcal/mol ✦ /RT) cm3 mol✷ 1 s✷ 1 for the forward reaction and k(T) ✺ 5.55 ✸ 108T1.44 exp( ✁ 0.78 ⑦ kcal/mol ✦ /RT) cm3 mol✷ 1 s✷ 1 for the reverse reaction in the temperature range 400-2500 K. © 2000 American Institute of Physics. ❅ S0021-9606 ⑦ 00 ✦ 30139-8 ★ I. INTRODUCTION Nitrogen-containing compounds play an important role in atmospheric chemistry and in many combustion and ex- plosion processes. The mechanisms of these processes are generally very complicated and involve hundreds of reac- tions or more.1,2 Kinetic simulations are often used to eluci- date the mechanisms of these processes.3,4 However, the ef- fectiveness of such simulations depends on the accuracy of rate constants of elementary reactions. It is difficult to mea- sure thermal rate constants of reactions involving reactive radicals due to the rapid secondary reactions of the product radicals. Thus, direct ab initio dynamics methods, which have been confirmed to be sufficiently accurate, are often effective tools for predicting the kinetics of these reactions.5-8 The N(4S) ✶ H2 􀀀 NH ✶ H reaction is a very important elementary reaction in many combustion processes of N-containing compounds, e.g., the thermal decomposition of RDX and HMX.9 However, there have been very few direct experimental studies on this reaction due to the difficulty in preparation of N atoms. Koshi et al.10 studied this reaction by the atomic resonance absorption technique in a shock tube apparatus over the temperature range of 1600-2300 K. The Arrhenius expression is k(T) ✺ (2.8 ✻ 0.2) ✸ 1014 exp( ✁ 33 ✻ 7 ⑦ kcal/mol ✦ /RT) cm3 mol✷ 1 s✷ 1. However, Davidson et al.'s experiment gave a smaller pre-exponential factor of 1.6 ✸ 1014 and a lower activation energy of 25.14 kcal/mol using the same experimental technique over the temperature range of 1950-2850 K.11 The only measured rate constant for the reverse reaction determined by Morley12 is k ✺ 1.02 ✸ 1013 cm3 mol✷ 1 s✷ 1 for the temperature range of 1790- 2200 K. This value is further recommended by Baulch et al.13 for the temperature range of 1500-2500 K. Hanson and Salimian, also based on Morley's data, suggested the rate constants to be 3.01 ✸ 1013 cm3 mol✷ 1 s✷ 1 for the tem- perature range of 1790-2200 K.14 Koshi et al. also carried out ab initio study of the N(4S) ✶ H2 􀀀 NH ✶ H reaction.10 They optimized geometries of the stationary points at the HF/4-31G level of theory and calculated the energies at the MP4/6-31G(d,p) level of theory. The calculated forward and reverse barriers for quartet state reaction are 35.3 and 4.4 kcal/mol, respectively. These values are in accordance with Koshi et al.'s experimental results, but too large compared to Davidson et al.'s results. In the present study, we employed a high-level ab initio molecular orbital method to study the reactions of N atom with H2 on several low-lying electronic states. The thermal rate constants of the N(4S) ✶ H2 􀀀 NH ✶ H reaction and its reverse reaction were calculated using the direct ab initio dynamics approach.5,7,8,15 Rate calculations are done using the microcanonical variational transition state theory ⑦♠ VT ✦ .16-21 II. METHODOLOGY A. Microcanonical variational transition state theory Microcanonical variational transition state theory ⑦♠ VT is based on the idea that by minimizing the microcanonica✦l rate constants k(E) along the MEP, one can minimize the error caused by the ‘‘recrossing'' trajectories.16-21 Within the framework of ♠ VT, the thermal rate constant at a fixed temperature T can be expressed as k✂ VT ✄ T ☎✺ ✯0❵ min ✆ NGTS ✄ E,s ☎✪ e✷ E/kBTdE hQR , ⑦ 1 ✦ where s is the statistical factor of the reaction, which is 2 for the forward and 1 for the reverse direction of the N(4S) ✶ H2 􀀀 NH ✶ H reaction, h and kB are the Planck and Boltz- a ✝Author to whom correspondence should be addressed. Electronic mail: [e-mail] JOURNAL OF CHEMICAL PHYSICS VOLUME 113, NUMBER 15 15 OCTOBER 2000 0021-9606/2000/113(15)/6149/5/$17.00 6149 © 2000 American Institute of Physics mann constants, respectively. QR is the total reactant parti- tion function, which is the product of electronic, rotational, and vibrational partition functions. The relative translational partition function was calculated classically. However, the rotational and vibrational partition functions of the reactant were calculated quantum mechanically within the rigid rotor and harmonic oscillator approximations. NGTS(E,s) is the sum of states of electronic, rotational, and vibrational mo- tions at energy E of the generalized transition state located at s, where s is the reaction coordinate. NGTS(E,s) along the MEP were also calculated quantum mechanically using the rigid rotor and harmonic oscillator approximations. The rate constant calculations were carried out employ- ing the THERATE program.15 B. Electronic structure calculations Spin contamination has been a problem in the unre- stricted Hartree-Fock wave functions of free-radical reac- tions. A recent study by Truhlar et al.22 showed that UQCISD23 and UCCSD24 methods give reasonable geom- etries even though the spin contamination can be relatively large. However, the inclusion of triple excitations in the UQCISD ⑦ T ✦ 23 or UCCDT ⑦ T ✦ 23 is sometimes necessary for a good estimation of energy. In the present study, we used the UQCISD method with Dunning's correlation-consistent double- ③ basis set25 ⑦ ccp-VDZ ✦ to perform all the geometri- cal optimizations, IRC calculations, and frequency analyses. Since the ground-state reaction is on the quartet-state poten- tial energy surface, we included corrections for triple and quadruple excitations in the QCISD method26 ❅ QCISD ⑦ TQ ✦􀀀 and used the Dunning's correlation-consistent triple- ③ basis set27 ⑦ cc-pVTZ ✦ to refine the energies of the stationary points and selected points along the IRC of this reaction. The IRC was calculated in the mass-weighted internal coordinate with the small step size of 0.02 amu1/2 bohr using the Gonzalez- Schlegel method.28 Twenty-one points were automatically selected for Hessian calculations along the MEP by a focus- ing technique.15 The information of these points was further used for calculation of variational rate constants. All elec- tronic structure calculations were performed using the GAUSSIAN 98 program.29 III. RESULTS AND DISCUSSION A. Potential energy surfaces In order to have a good understanding of the reactions between N atom and H2, we calculated reaction channels involving different electronic states of N atom. The three lowest electronic states of N atom are 4S, 2D, and 2P. We found that the ground state of N atom 4S reacts with H2 to form the ground-state NH(3 ❙ ✷) and H products along the 4 ❙ electronic state reaction path. An excited state of N atom, possibly the 2P electronic state, reacts with H2 to form the FIG. 1. Relative energies and schematic diagram for the potential energy surfaces of N atom with H2. The energies of stationary points were calcu- lated using QCISD/cc-pVDZ level of theory. FIG. 2. Geometric parameters and frequencies of the stationary points. The italic numbers are experimental data taken from Ref. 30. FIG. 3. Plot of classical potential energy along the minimum energy path vs the mass-weighted internal coordinate s. TABLE I. Relative energies ✁ kcal/mol ✂ of the stationary points on the po- tential energy surfaces of N atom with H2.a N(4S) ✶ H2 N(2P) ✶ H2 TS1 NH2(2B1) NH ✁ 3 ✄✂✶ H QCISDb 0.0 70.73 35.46 ☎ 55.56 32.92 QCISD ✁ TQ ✂ c 0.0 31.33 29.16 QCISD ✁ TQ ✂✶ ZPEd 29.60 29.07 aSee Fig. 1 for designation of TS1. bUsing cc-pVDZ basis set. cUsing cc-pVTZ basis set. dZPE were calculated at QCISD/cc-pVDZ level of theory. 6150 J. Chem. Phys., Vol. 113, No. 15, 15 October 2000 S. Zhang and T. N. Truong global minimum of the N ✶ H2 system, i.e., NH2(2B1), on the 2B1 reaction path. NH2(2B1) can further decompose to form the ground state NH(3 ❙ ✷) and H atom. We have attempted to locate the reaction between the 2D state of N atom and H2, but found the spin contamination of the wave function is rather large during the self-consistent field ⑦ SCF ✦ calculation. The results in this case may not be reliable and thus are not included here. We found the transition state for hydrogen abstraction reaction on the 4 ❙ potential energy surface which suggests the reaction path to be rectilinear. A schematic dia- gram of the reaction channels is shown in Fig. 1. Table I lists the relative energies of the stationary points of these reactions. The first excited state of the reactants is 70.73 kcal/mol higher than that of the ground state at the QCISD/cc-pVDZ level of theory. Thus, the reaction is mainly proceeding on the ground-state surface. The forward barrier of the ground state channel is 35.46 kcal/mol, which is close to the MP4/6-31G(d,p) value of 35.3 kcal/mol cal- culated by Koshi et al. However, the calculated reverse bar- rier of 2.54 kcal/mol is smaller than that of 4.4 kcal/mol from Koshi et al.'s MP4 calculation.10 The best estimates for the forward and reverse classical barrier heights calculated at the QCISD ⑦ TQ ✦ /cc-pVTZ level of theory are 31.33 and 2.17 kcal/mol, respectively. Including zero-point energy ⑦ ZPE ✦ correction reduces the forward and reverse barriers by 1.73 and 1.64 kcal/mol, respectively. Note that the ZPE-corrected forward barrier is within the experimental range of 33 ✻ 7 kcal/mol for the activation energy.10 The calculated geometries and frequencies at the station- ary points are shown in Fig. 2 along with the available ex- perimental data.30 The calculated bond lengths of equilib- rium species are generally larger than experimental values by less than 0.02 Å. The calculated bond angle of NH2 is 2.4° smaller than that of the experimental data. The differences in the harmonic vibrational frequencies between calculations and experiments are smaller than 5%. At the transition state of the 4 ❙ reaction channel, the N-H bond ⑦ 1.109 Å ✦ is shorter than the H-H bond ⑦ 1.275 Å ✦ . B. Minimum energy path of the ground-state reaction We calculated the minimum energy path for the N ✶ H2 reaction on the 4 ❙ ground-state surface at the QCISD/ cc-pVDZ level of theory. The reaction path was found to be a rectilinear one. Figures 3 and 4 display the classic potential energy profile and the variation of bond lengths along the IRC, respectively. Note that the potential energy profile was corrected by a series of QCISD ⑦ TQ ✦ /cc-pVTZ single-point calculations. The energy rises sharply in the region from 􀀀 1.0 to 0.0 amu1/2 bohr of the mass-weighted coordinate due to the breaking of the H-H bond, whereas the N-H bond is still in the process of forming. However, the rates of bond formation and dissociation are not the same. As the N(4S) atom approaches the H2 molecule, a repulsive interaction be- tween electrons in the p orbitals of the N atom with electrons in the s orbital of H2 gives rise to the first ❀ 18 kcal/mol of the barrier. At s, about 􀀀 0.8H2 begins to dissociate and the barrier continues to rise another ❀ 13 kcal/mol. At the tran- sition state, the N-H bond is already 95% formed. Figure 5 shows the plot of generalized frequencies along the reaction coordinate s. There are three vibrational modes along the reaction path; a doubly degenerate P bending vi- bration, and a ❙ stretching vibration. The bending vibrations reach the maximum at the reaction coordinate s about 􀀀 0.5 and they become the asymptotic values of zero at both reac- tant and products. The stretching shows the transformation of the H2 stretch to the NH stretch vibration. FIG. 4. Plot of bond lengths vs the mass-weighted internal coordinate s. FIG. 5. Plot of frequencies vs the mass-weighted internal coordinate s. TABLE II. Comparison of the sum of states at the transition state and the minimum of sum of state along the MEP for the forward N(4S) ✁ H2 ✂ NH(3 ✄ ) ✁ H reaction. Excess energya 5.0 20.0 50.0 100.0 125.0 N(E) at TS 2.01E3 1.06E5 2.67E6 3.62E7 8.71E7 Minimum of N(E) 1.93E3 8.77E4 1.78E6 2.05E7 4.68E7 aUnit: kcal/mol. J. Chem. Phys., Vol. 113, No. 15, 15 October 2000 N 6151 ✁ H2 reaction C. Rate constants As shown in Fig. 3, the potential energy profile of the reaction N(4S) ✶ H2 ✦ NH(3 ❙ ✷) ✶ H is rather flat in the exit channel and the reverse barrier is quite small. The tunneling effect is expected to be small in this case. However, this type of reaction often has large recrossing effects. Table II lists the sum of states at the transition state and the minimum of the sum of state along MEP of the forward reaction at a given energy E. It can be seen that the differences between them become larger with increase of excess energy. At the excess energy of 125.5 kcal/mol, they differ by a factor of almost 2. Table III shows the calculated rate constants with and without variational effects for the forward reaction. The variational effects become larger with the increase of tem- perature. At 400 K, ♠ VT rate constant is about 8.0% smaller than the nonvariational transition state theory ⑦ TST 􀀀 rate constant. At 2500 K, the difference becomes 38%. Figure 6 displays the calculated forward rate constants along with the available experimental data.10,11 It can be seen that the two series of experimental data approach each other when the temperature increases. The ❉ log k deviation of the two ex- periments is 0.9 at 1500 K and decreases to 0.44 at 2500 K. The calculated rate constants are well located between the experimental data. However, the predicted values are closer to the data of Koshi et al. than that of Davidson et al. at 1500 K and vice versa at 2500 K. The calculated rate constants can be expressed in the Arrhenius form as k(T) ✺ 2.33 ✸ 1014 exp( ✁ 30.83 ⑦ kcal/mol 􀀀 /RT) cm3 mol✷ 1 s✷ 1. The calculated rate constants along with the available experimental data for the reverse reaction are depicted in Fig. 7. No activation barrier was observed in the experiment of Morley. In our calculation, the best estimation of the barrier is 0.53 kcal/mol, which is close to the experimental result.12,13 The calculated rate constants are again in good agreement with the experimental data. The largest difference in ❉ log k between the predicted data and experiment data is only 0.47. The fitted Arrhenius expression for the reverse reaction is k(T) ✺ 5.55 ✸ 108T1.44 exp( ✁ 0.78 ⑦ kcal/mol 􀀀 /RT) cm3 mol✷ 1 s✷ 1. Since the agreement with experiment data is excellent in the high-temperature limits and there are no ex- perimental data available in the low-temperature range, the expressions obtained in the present work provide a good es- timates for the kinetics of this reaction over a wide tempera- ture range. IV. CONCLUSIONS We presented in this work a high-level direct ab initio dynamic study of the reaction between N atom and H2. Two electron state potential energy surfaces of this reac- tion, namely N(4S) ✶ H2 and N(2P) ✶ H2, were studied at the QCISD/cc-pVDZ level of theory. The reaction of N(4S) ✶ H2 ✦ NH(3 ❙ ✷) ✶ H was found to be the ground- state reaction. Energetic properties along the ground-state reaction path were further improved by carrying out a series of single point QCISD ⑦ TQ 􀀀 /cc-pVTZ calculations. Microca- nonical variational transition state theory was used to calcu- late the rate constants for both forward and reverse reac- tions. The calculated rate constants for both forward and TABLE III. Comparison of rate constants ✂ cm3 mol✄ 1 s✄ 1 ☎ using nonvariational ✂✆ TST ☎ method and variational ✂✆ VT ☎ method for the forward N(4S) ✝ H2 ✞ NH(3 ✟ ) ✝ H reaction. Temperaturea 400 600 800 1000 1500 2000 2500 ✆ TST 7.96E-27 2.20E-21 1.29E-18 6.45E-17 1.45E-14 2.55E-13 1.57E-12 ✆ VT 7.38E-27 1.94E-21 1.10E-18 5.35E-17 1.13E-14 1.91E-13 1.13E-12 aUnit: K. FIG. 6. Comparison of experimental data and calculated forward rate con- stants of the N(4S) ✝ H2 ✞ NH(3 ✟ ) ✝ H reaction vs 1000/T. FIG. 7. Comparison of experimental data and calculated reverse rate con- stants of the N(4S) ✝ H2 ✞ NH(3 ✟ ) ✝ H reaction vs 1000/T. 6152 J. Chem. Phys., Vol. 113, No. 15, 15 October 2000 S. Zhang and T. N. Truong reverse reactions are in good agreement with available ex- perimental data in the high-temperature limit. They can be expressed as k(T) ✺ 2.33 ✸ 1014 exp( ✷ 30.83 ⑦ kcal/mol ✦ / RT) cm3mol􀀀 1 s􀀀 1 for the forward and k(T) ✺ 5.55 ✸ 108T1.44 exp( ✷ 0.78 ⑦ kcal/mol ✦ /RT) cm3mol􀀀 1s􀀀 1 for the reverse reaction in the temperature range 400-2500 K. These expressions are particularly useful for estimating thermal rate constants for temperatures below 1500 K where experimen- tal data are not available. ACKNOWLEDGMENTS This work is supported by the University of Utah Center for the Simulation of Accidental Fires & Explosions, funded by the Department of Energy, Lawrence Livermore National Laboratory, under Subcontract No. B341493. An allocation of computer time from the Center for High Performance Computing is gratefully acknowledged. 1 T. B. Brill, P. E. Gongwer, and G. K. Williams, J. Phys. Chem. 98, 12242 ✁ 1994 ✂ . 2 C. F. Melius, in Chemistry and Physics of Energetics Materials, edited by S. N. Bulusu ✁ Kluwer Academic, Dordrecht, 1990 ✂ . 3 Y.-X. Zhang and S. H. Bauer, Int. J. Chem. Kinet. 31, 655 ✁ 1999 ✂ . 4 D. Chakraborty and M. C. Lin, J. Phys. Chem. A 103, 601 ✁ 1999 ✂ . 5 T. N. Truong, J. Chem. Phys. 100, 8014 ✁ 1994 ✂ . 6 K. Song and W. L. Hase, J. Chem. Phys. 110, 6198 ✁ 1999 ✂ . 7 K. K. Baldridge, M. S. Gordon, R. Steckler, and D. G. Truhlar, J. Phys. Chem. 93, 5107 ✁ 1989 ✂ . 8 D. G. Truhlar and M. S. Gordon, Science 249, 491 ✁ 1990 ✂ . 9 K. Prasad, R. A. Yetter, and M. D. Smooke, Combust. Sci. Technol. 124, 35 ✁ 1997 ✂ . 10M. Koshi, M. Yoshimura, K. Fukuda, H. Matsui, K. Saito, M. Watanabe, A. Imamura, and C. Chen, J. Chem. Phys. 93, 8703 ✁ 1990 ✂ . 11 D. F. Davidson and R. K. Hanson, Int. J. Chem. Kinet. 22, 843 ✁ 1990 ✂ . 12 C. Morley, in Proceedings of the 18th Symposium on Combustion ✁ But- terworth, London, 1981 ✂ , pp. 23. 13 D. L. Baulch, C. J. Cobos, R. A. Cox, C. Esser, P. Frank, T. Just, J. A. Kerr, M. J. Pilling, J. Troe, R. W. Walker, and J. Warnatz, J. Phys. Chem. Ref. Data 21, 411 ✁ 1992 ✂ . 14 R. K. Hanson and S. Salimian, in Combustion Chemistry, edited by W. C. Gardiner ✁ Springer, New York, 1984 ✂ . 15W. T. Duncan, R. L. Bell, and T. N. Truong, J. Comput. Chem. 19, 1039 ✁ 1998 ✂ . 16 D. G. Truhlar, B. C. Garrett, and S. J. Klippenstein, J. Phys. Chem. 100, 12771 ✁ 1996 ✂ . 17W. Hase, Acc. Chem. Res. 31, 659 ✁ 1998 ✂ . 18 B. C. Garrett and D. G. Truhlar, J. Phys. Chem. 83, 1079 ✁ 1979 ✂ . 19 B. C. Garrett, D. G. Truhlar, R. S. Grev, and A. W. Magnuson, J. Phys. Chem. 84, 1730 ✁ 1980 ✂ . 20 A. D. Issacson, M. T. Sund, S. N. Rai, and D. G. Truhlar, J. Chem. Phys. 82, 1338 ✁ 1985 ✂ . 21 D. G. Truhlar, A. D. Isaacson, and B. C. Garrett, in Theory of Chemical Reaction Dynamics, edited by M. Baer ✁ CRC, Boca Raton, FL, 1985 ✂ , Vol. 4, p. 65. 22 Y.-Y. Chuang, E. L. Coitino, and D. G. Truhlar, J. Phys. Chem. A 104, 446 ✁ 2000 ✂ . 23 J. A. Pople, M. Head-Gordon, and K. Raghavachari, J. Chem. Phys. 87, 5968 ✁ 1987 ✂ . 24 G. E. Scuseria and H. F. Schaefer III, J. Chem. Phys. 90, 3700 ✁ 1989 ✂ . 25 D. E. Woon and T. H. Dunning, Jr., J. Chem. Phys. 98, 1358 ✁ 1993 ✂ . 26 K. Raghavachari, J. A. Pople, E. S. Replogle, and M. M. Head-Gordon, J. Phys. Chem. 94, 5579 ✁ 1990 ✂ . 27 R. A. Kendall, T. H. Dunning, Jr., and R. J. Harrison, J. Chem. Phys. 96, 6796 ✁ 1992 ✂ . 28 C. Gonzalez and H. B. Schlegel, J. Phys. Chem. 94, 5523 ✁ 1990 ✂ . 29M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Robb, J. R. Cheeseman, V. G. Zakrzewski, J. Montgomery, R. E. Stratmann, J. C. Burant et al., Gaussian, Inc., Pittsburgh, PA, 1998. 30 D. R. Lide, CRC Hand Book of Chemistry and Physics, 80th ed. ✁ CRC, Boca Raton, FL, 1999 ✂ . J. Chem. Phys., Vol. 113, No. 15, 15 October 2000 N 6153 ✶ H2 reaction