Revisited with Multireference Methods - American Chemical Society

522

J. Phys. Chem. A 2008, 112, 522-532

Kinetics of CH + N2 Revisited with Multireference Methods Lawrence B. Harding* and Stephen J. Klippenstein Chemical Sciences and Engineering DiVision, Argonne National Laboratory, 9700 South Cass AVenue, Argonne, Illinois 60439

James A. Miller Combustion Research Facility, Sandia National Laboratories, LiVermore, California 94551-0969 ReceiVed: September 18, 2007; In Final Form: October 26, 2007

The potential energy surface for the CH + N2 reaction was reexamined with multireference ab initio electronic structure methods employing basis sets up to aug-cc-pvqz. Comparisons with related CCSD(T) calculations were also made. The multireference ab initio calculations indicate significant shortcomings in single reference based methods for two key rate-limiting transition states. Transition state theory calculations incorporating the revised best estimates for the transition state properties provide order of magnitude changes in the predicted rate coefficient in the temperature range of importance to the mechanism for prompt NO formation. At higher temperatures, two distinct pathways make a significant contribution to the kinetics. A key part of the transition state analysis involves a variable reaction coordinate transition state theory treatment for the formation of H + NCN from HNCN. The present predictions for the rate coefficients resolve the discrepancy between prior theory and very recent experimental measurements.

1. Introduction Prompt NO has fascinated combustion chemists for more than three decades. The phenomenon was first discovered by Fenimore1 and reported at the Thirteenth International Combustion Symposium held in Salt Lake City in August of 1970. Fenimore’s paper was the first one ever presented at a session of the Combustion Symposium organized around combustiongenerated air pollution. His results were controversial from the beginning. Fenimore was studying NO formation in the post-flame gases of a series of flames fueled by H2, CO, and various hydrocarbons at atmospheric and higher pressures. In his flames, the primary reaction zone, immediately adjacent to the burner surface, was largely unresolvable. Fenimore observed in all the flames that the NO increased with distance above the burner through the thermal, or Zeldovich, mechanism

Fenimore suggested two possibilities for the prompt- NO reaction

CH + N2 f HCN + N

(R1a)

C2 + N2 f CN + CN

(R2)

and

In 1970, both reactions (R1a and R2) were considered unlikely. R1a is spin-forbidden, and R2 is four-centered. There ensued a search for alternative explanations. The search was pursued by many people, and several other explanations were proposed.2-4 However, in 1989, Miller and Bowman4 examined the possible sources of prompt NO from a number of perspectives and concluded that indeed R1a must be responsible for the phenomenon, with perhaps a minor contribution at very high temperatures from the reaction of carbon atoms with N2

O + N2 f NO + N

C + N2 f CN + N

N + O2 (OH) f NO + O (H)

For a variety of reasons, direct experimental determinations of k1 have been limited to very high temperatures5,6 (T > 2300 K), although there have been some indirect inferences at lower temperatures.3,4,7 At the same time, theoretical investigations of R1 have been hampered by the problem of handling the spin change in R1a, including finding the intersection that must occur between the doublet and the quartet potential energy surfaces and calculating the probability of a collision complex hopping from one surface to another. Many theorists have contributed to our understanding of these issues.8-17 However, two investigations were pivotal. Cui et al.16 showed conclusively that the spin-orbit coupling responsible for the spin flip in R1a was not strong enough to account for experimental observations. This

However, the NO profiles extrapolated to zero at the burner surface in the H2 and CO flames, whereas there was a nonzero intercept in all the hydrocarbon flames. This intercept NO, or prompt NO, increased with an increasing equivalence ratio up to about Φ ) 1.4, prompting Fenimore to attribute the phenomenon to a reaction in the primary reaction zone between a hydrocarbon free radical and molecular nitrogen, the products of which presumably react further to form NO. * Corresponding author. E-mail: [email protected]; fax: (630) 252-9292.

10.1021/jp077526r CCC: $40.75 © 2008 American Chemical Society Published on Web 01/03/2008

(R3)

Kinetics of CH + N2 with Multireference Methods

J. Phys. Chem. A, Vol. 112, No. 3, 2008 523

prompted Moskaleva et al.17,18 to look for alternative product channels. They showed that indeed another set of products

CH + N2 f HNCN f NCN + H

(R1b)

could account for the experimental results available at the time, even though R1b is significantly endothermic. Subsequently, NCN has been detected in flames,19 and photodissociation experiments on the HNCN20 and DNCN radical21 support the mechanism of Moskaleva and Lin. Additionally, a CASPT2 study has been reported22 confirming the existence of the R1b pathway. The controversy concerning R1 thus appeared to have been completely resolved15 until very recently, when two independent studies called the results into question. A modeling study23 of natural gas flames, using the rate calculated by Moskaleva and Lin, led to an underestimation, by more than a factor of 6, of the amount of prompt NO. Furthermore, Hanson and co-workers have extended the temperature range of their k1(T) measurements to below 2000 K.24 These measurements indicate that the Moskaleva and Lin prediction of the rate coefficient is an order of magnitude too low at 2000 K and probably even more so in the practically important temperature range, 1000 K < T < 2000 K. The present investigation was initiated to see if this discrepancy was due simply to an inadequacy of the theoretical methodology employed by Moskaleva and Lin or whether some new unidentified mechanism might be at work. Indeed, recently Berman et al.25 have reported a new study of the potential surface suggesting that HNNC may play a role in the kinetics of this reaction. In this paper, we extend the previous theoretical studies in several significant ways. Geometries and frequencies for all of the stationary points were obtained using a multireference method (CASPT2/aug-cc-pvtz). These calculations demonstrate that previous single reference approaches are inadequate for one of the key transition states. A second key extension of previous work is the use of variable reaction coordinate transition state theory to calculate the rate for the final dissociation of HNCN to H + NCN, a reaction for which there is no saddle point. Finally, we also briefly consider the possible contribution from the quartet state of CH. 2. Methods 2.1. Electronic Structure Methods. In this work, we made use of three ab initio electronic structure methods: (i) multireference second-order perturbation theory, CASPT2, (ii) internally contracted multireference configuration interaction with, CAS+1+2+QC, and without, CAS+1+2, the Davidson correction, and (iii) coupled-cluster theory restricted to single and double excitations with perturbative triple excitations, CCSD(T). All of the ab initio calculations employ the Dunning26-28 correlation consistent basis sets. All of the multireference calculations and most of the CCSD(T) calculations were performed using the MOLPRO program package.29-34 For comparison to previous work, we also report some CCSD(T) calculations performed with the Gaussian 98 program.35 The CASPT2 calculations can be performed either with or without internal contraction. For the largest active spaces, with no symmetry, it was only possible to perform the contracted CASPT2 calculations. For smaller active spaces and points of higher symmetry, the CASPT2 calculations were performed both ways. Note that the Gaussian 98 and MOLPRO implementations of open shell CCSD(T) are different. The MOLPRO implementation uses ROHF orbitals, while Gaussian 98 implementation

uses UHF orbitals. Again, for comparison to previous work, we also report some B3LYP calculations performed using Gaussian 98. The final energies were extrapolated to a complete basis set (CBS) limit using the following formula36

Eaug-cc-pVxZ ) ECBS +

B (nx + 1)4

where nT ) 3 and nQ ) 4 for the aug-cc-pVTZ and aug-ccpVQZ basis sets, respectively. One of the prime motivations for this study was to characterize the key features of the potential surface using multireference ab initio methods. This problem is a challenging one for multireference approaches as a full valence active space, with 15 electrons in 13 orbitals, (15e,13o), is not computationally tractable. Takayanagi22 reported a characterization of the surface using CASPT2 calculations with an (11e,11o) active space (this includes all valence orbitals except the nitrogen 2s orbitals) combined with a Shepard interpolation scheme for fitting the calculated points and deriving geometries and frequencies. In this work, we also use CASPT2 calculations to derive geometries and frequencies; however, we make use of analytic gradients37 rather than Shepard interpolation. The majority of the calculations was performed using (3e,3o), (5e,5o), and/or (7e,7o) active spaces, although some key parts of the surface were reexamined with (9e,9o) and (11e,11o) active spaces to verify convergence with respect to the size of the active space. Geometry optimizations and frequency calculations were performed using the aug-cc-pvtz basis set. Single point CASPT2, CAS+1+2+QC, and CCSD(T) calculations were then performed using the larger aug-cc-pvqz basis set. 2.2. Kinetics Methods. The number of available states for all of the transition states was determined at the energy E and total angular momentum J resolved level. These transition state theory (TST) calculations were based on rigid-rotor harmonicoscillator (RRHO) assumptions for the energy levels, except for the transition state connecting HNCN with H + NCN. For the latter transition state, the reverse reaction is a barrierless radical-radical reaction, for which RRHO assumptions are not expected to be appropriate. Thus, for that channel, the number of available states was instead obtained from variable reaction coordinate TST (VRC-TST).38-40 Variational effects were also considered for TS1 and TS3 via standard reaction path evaluations but were found to be insignificant. The canonical partition function for the reactants was similarly determined with RRHO assumptions. The rovibrational properties for the TST calculations were obtained from the CASPT2/aug-cc-pvtz calculations. The active spaces employed in these analyses are described in detail in the electronic structure results. The saddle point energies were obtained from the present best estimates, which are based on extrapolations of CCSD(T) or CI+1+2+QC calculations (as appropriate) with the aug-cc-pvtz and aug-cc-pvqz basis sets. The potential energy surface for the VRC-TST calculations in the NCN + H channel was obtained from an analytic fit to a grid of CI+1+2+QC/aug-cc-pvtz calculations. Essentially, identical results could be expected for a CASPT2 surface. These VRC-TST calculations are closely analogous to our prior calculations for the isoelectronic 3HCCCH + H system described in detail in ref 41. An empirical dynamical correction factor of 0.9 is also incorporated as discussed previously.42 At low temperatures, stabilization to HCNN dominates the kinetics observed in a number of experimental studies.43-46 This

524 J. Phys. Chem. A, Vol. 112, No. 3, 2008

Figure 1. Schematic representation of the stationary point energies (including zero point) on the CH + N2 doublet potential surface (see text for details).

stabilization is irrelevant to the question of prompt NO, and so we do not consider this aspect further. At temperatures of 1000 K and higher, stabilization in this HCNN well becomes insignificant. For such temperatures, sample master equation simulations indicate that the rate coefficient to produce H + NCN is essentially independent of pressure. The HNCN well, which is the deepest well, is still not deep enough to support any collisional stabilization at such temperatures. In essence, the small number of vibrational modes (six) correlates with rapid dissociation and ineffective collisional stabilization. Thus, the results presented next are for the collisionless limit obtained with the methodology described in ref 47. The latter approach employs an analytic solution to the two-dimensional master equation based on the inversion of the rate matrix within the steady-state approximation. It requires only the canonical partition function for the reactants and the transition state number of available states for each of the transition states. It is readily implemented at the E and J resolved level. 3. Results and Discussion 3.1. Electronic Structure Results. 3.1.1. OVerView. A schematic representation of the important stationary points is given in Figure 1. Here, we use the transition state numbering scheme from ref 25, and we simplify the figure by omitting TS4 and TS9, both of which are too high to be important in the thermal reaction. There are two distinct paths for the reaction CH + N2 f NCN + H. Path A, first proposed by Moskaleva and Lin17 and shown with a solid line in Figure 1, follows the route

CH + N2 f TS1 f Ring1 f TS2 f Ring2 f TS3 f HNCN f NCN + H Path B, proposed by Berman et al.,25 is shown with a dotted line in Figure 1 and follows the route

CH + N2 f TS7 f HNNC f TS8 f Ring3 f TS5 f HNCN f NCN + H The two paths converge at HNCN, which then dissociates to form the final products. There is a minor difference between path B as shown in Figure 1 and the results of Berman et al. Our results indicate that TS7 connects directly to the reactants, while Berman et al. has TS7 connected to Ring2. This may simply be a difference between the topology of the CCSD(T)

Harding et al.

Figure 2. Schematic representation of the stationary point energies (including zero point) on the CH + N2 quartet potential surface (see text for details).

TABLE 1: CASPT2/aug-cc-pvtz Frequencies (cm-1) and Rotational Constants (GHz) Nacta 2CH 4CH 1N 2 2Ring

1

2Ring

2

2Ring

3 2HCNN 2HNNC 2HNCN 2HCN

2 4HCNN 4HCN

2

3NCN 1HCN 2TS

1

2TS

2

2TS

3-C2V

2TS

3-Cs

2TS

5

2TS

7

2TS

8

4TS

q1

4TS

q2

3 3 4 3 5 5 7 5 5 11 7 5 6 4 7 5 11 11 5 5 5 7 7

A

B

C

438.0 43.0 41.5 42.1 606.2 570.2 632.9 317.1 108.4 95.6 50.8 42.7 185.9 334.2 53.8 45.9 56.5 114.5 75.5

59.2 29.7 28.1 30.1 12.1 12.0 11.1 11.1 14.4 13.1 11.8 44.1 21.6 29.6 11.3 11.1 19.3 25.4 20.2 10.0 11.7

18.3 16.8 18.0 11.9 11.7 10.9 10.7 12.7 11.5 15.3 17.9 10.7 10.7 14.2 16.4 15.2 9.2 10.2

ω1

ω2

ω3

ω4

ω5

ω6

2945 3134 2328 3186 3286 3537 3216 3417 3505 2739 3147 3042 1623 3490 3156 3272 2892 2726 3267 2353 3448 3097 3313

1546 1588 1482 1839 1931 1862 1493 1336 1266 1236 2109 1636 1512 1184 1451 1680 1518 1701 1957 1753

1102 1225 1064 1232 1368 1186 1063 1214 1222 420 654 962 1057 1036 1073 1202 978 1288 545 790

1050 902 1003 902 1111 1072 794 1035 1139 420 654 689 906 952 449 798 952 731 385 757

953 891 786 504 363 455 458 713 807

666 625 492 464 317 447 226 512 550

459 816 580 294 278 726 427 367 373

386i 1188i 243i 937i 584i 1088i 446i 648i 772i

a

Nact is the number of orbitals (and electrons) included in the active space.

and CASPT2 surfaces, but it is unlikely to have a significant impact on the kinetics. A third path exists that begins with a barrierless addition to form HCNN, but subsequent barriers on this path have been shown11 to be too high for it to compete with the first two paths. At higher temperatures, where the 4Σ state of CH becomes thermally populated, reaction can also occur on the quartet surface as shown in Figure 2. The rotational constants and vibrational frequencies for all of the stationary points shown in Figures 1 and 2 (and used in the kinetics calculations) are given in Table 1. One of the challenges in this study is the choice of active spaces used in the multireference calculations. Taking Ring1 as an illustrative example, the structure is believed to be of Cs symmetry. CASPT2 calculations with (3e,3o) or (9e,9o) active spaces yield a Cs structure. Calculations with (5e,5o) or (7e,7o) active spaces, however, yield a slightly distorted C1 structure owing to an imbalance between the orbitals included in the active space and those treated as closed shells. For Ring1, the differences in the structures and frequencies from the use of these four different active spaces are too small to be of significance in the kinetics



TABLE 2: Relative Energies (kcal/mol) and T1 Diagnostics of Stationary Points from aug-cc-pvqz Calculationsa CH + N2 4 CH + N2 2 Ring1 2 Ring2 2 Ring3 2 HCNN 2 HNNC 2 HNCN 2 HCN2 2 H + 3NCN 4 HCNN 4 HCN2 1 HCN + 4N 2 TS1 2 TS2 2 TS3-Cs 2 TS3-C2V

2

2

TS5 TS7 2 TS8 4 TSq1 4 TSq2 2

CASPT2b

CASPT2c

CAS+1+2+QC

CCSD(T)

T1

0.0(0.0) 13.9(13.1) -12.9(-11.8) -31.4(-30.4) 0.4(1.6) -36.4(-35.3) -28.6(-27.6) -70.1(-69.1) 2.7(3.1) 14.6(14.4) 4.4(4.7) -12.6(-12.9) -5.7(-7.1) 4.3(4.8) -11.1(-10.1) 2.8(3.2) 3.4(3.4)

0.0(0.0) 14.6(13.8)

0.0(0.0) 17.4(16.8)

-31.9(-30.9)

-27.2(-26.3)

0.0(0.0) 16.8(16.3) -12.2(-11.1) -27.6(-26.6) 2.8(4.0) -31.7(-30.5) -22.5(-26.4) -67.4(-66.3) 10.4(10.7) 21.7(21.7) 12.1(12.5) -4.9(-5.2) -1.2(-2.6) 8.9(9.5) -7.0(-5.9) 9.8(10.2) -11.2(-15.2)d 13.8(13.6)e 10.0(10.6) 18.2(18.9) 8.0(8.8) 32.1(31.8) 31.7(30.7)

0.012 0.016 0.014 0.015 0.014 0.040 0.025 0.020 0.023 0.013 0.029 0.021 0.012 0.029 0.048 0.034 0.063d 0.027e 0.041 0.027 0.024 0.023 0.027

7.7(8.3) 14.3(15.0) 5.0(5.7) 27.5(27.0) 24.6(23.6)

(-28.3) (-69.6) 2.1(2.5) 14.0(13.8)

9.3(9.4) 21.1(21.0)

-13.2(-13.6) -6.3(-7.7)

-4.9(-5.4) -0.5(-2.0)

(2.7) 2.8(2.8)

(9.0) 10.9(10.7)

(6.9) (13.7)

(10.0) (18.7)

a Numbers in parentheses are from aug-cc-pvtz calculations. The multireference calculations all employ an (11e,11o) active space. Geometries are all from caspt2/aug-cc-pvtz calculations. b Internally contracted. c Uncontracted. d On the basis of RHF configuration a. e On the basis of RHF configuration b.

calculations. For example, the largest difference in vibrational frequencies is less than 20 cm-1. Larger differences in frequencies and geometries with different active spaces are found for TS3 (see section 3.1.3.). The relative energies of the stationary points from (11e,11o)CASPT2, (11e,11o)-CAS+1+2+QC, and CCSD(T) calculations are given in Table 2. Also given in Table 2 are the T1 diagnostics48 from the CCSD(T) calculations. There is very good agreement between the CAS+1+2+QC and the CCSD(T) energies for all stationary points having a small T1 diagnostic ( 0 to motion toward Ring1.

further check on the calculations in this part of the potential surface, we followed the intrinsic reaction coordinate (IRC) from TS1 using B3LYP and then evaluated CASPT2, CAS+1+2+QC, and CCSD(T) energies using the geometries from the B3LYP IRC. The results are shown in Figure 4. Again, it is seen that on the reactant side of TS1, the CASPT2 energies are anomalous. We tentatively conclude that the new stationary points, TS0 and Int0, are artifacts of the CASPT2 calculations of unknown origin. Focusing now on TS1, the structures and frequencies for TS1 from a variety of electronic structure calculations are given in Table 4. All are found to be in good agreement with each other except for MP2, which gives an unrealistically high NN stretching frequency of 3211 cm-1. In the kinetics calculations, we use the structure and frequencies from the (7e,7o)-CASPT2/ aug-cc-pvtz calculation. 3.1.3. TS3. Qualitatively, the electronic structure of TS3 can be diagrammed as follows: having an essentially broken NN σ

bond and a three electron allylic resonance in the π system. Both the broken σ bond and the resonance provide challenges to the electronic structure theory. In terms of orbitals, the

The two configurations differ by a double excitation between the 6a1 and the 4b2 orbitals. These two configurations are needed to allow the NN σ bond to evolve smoothly from a highly overlapping pair of bonding orbitals to a weakly overlapping pair of radical orbitals. As one moves down the reaction path toward Ring2, configuration a becomes increasingly dominant, while motion in the opposite direction causes configuration b to become increasingly dominant. For Ring2, the coefficients of configurations a and b are 0.94 and -0.14, respectively, while at TS3, these coefficients become 0.52 and -0.77, respectively. The next most important configurations are the ones describing the allylic resonance (c and d)

1a122a123a124a125a121b222b223b226a121b112b111a21 (c) 1a122a123a124a125a121b222b223b224b221b112b111a21 (d) corresponding to a single excitation between the 1b1 and the 2b1 orbitals from configurations a and b. At Ring2, c has a coefficient of 0.15, while at TS3, c and d have coefficients of 0.13 and 0.15, respectively. We find that the geometry of TS3 depends dramatically on the electronic structure method used. To better understand as to why this is the case, we show, in Figure 5, some twodimensional contour plots of the PES in the vicinity of TS3. These plots were obtained with both CASPT2 and CAS+1+2+QC calculations and with both aug-cc-pvdz and aug-cc-pvtz basis sets employing a (5e,5o) active space in each case. Two points are worth noting: first, the CASPT2 and CAS+1+2+QC calculations yield potentials that look remarkably similar. Second, all calculations show a potential that is very flat near TS3. This flatness is what causes the sensitivity of the exact location of the saddle points to the details of the electronic structure method used. Two saddle points can be seen in each of the four PESs shown in Figure 5. In the aug-cc-pvdz calculations, the two saddle points are of Cs symmetry, lying slightly off the C2V axis, and are related to each other by symmetry. In the aug-cc-pvtz calculations, both saddle points lie on the C2V axis, one being the transition state for ring opening and (1,2)-hydrogen migration (TS3), and the second slightly lower saddle point being a transition state for (1,3)-hydrogen migration in HNCN. Because HNCN is of lower symmetry than Ring2, the reaction path for Ring2 f HNCN must contain a

TABLE 4: Geometries and Vibrational Frequencies for 2TS1 method

RCH (Å)

RCN (Å)

RNN (Å)

HCN angle (deg)

CNN angle (deg)

HCNN dihedral angle (deg)

ωi (cm-1)

(7e,7o)-CASPT2/aug-cc-pvdz (7e,7o)-CAS+1+2/aug-cc-pvdz (7e,7o)-CAS+1+2+QC/aug-cc-pvdz (7e,7o)-CASPT2/aug-cc-pvtz B3LYP/6-311G(d,p)a MP2/6-311G(d,p) CCSD(T)/aug-cc-pvdz CCSD(T)/aug-cc-pvtz b

1.104 1.103 1.109 1.109 1.103 1.104 1.115 1.096

1.485 1.486 1.524 1.465 1.572 1.679 1.560 1.535

1.222 1.212 1.209 1.209 1.162 1.156 1.197 1.173

122.3 121.4 119.7 123.4 115.1 111.9 116.5

86.1 88.3 85.8 86.4 86.2 80.5 86.3

31.8 22.7 34.0 26.5 30.7 38.0 31.3

3146, 1634, 978, 678, 496, 381i 3158, 1678, 977, 642, 448, 255i 3096, 1696, 989, 633,426, 370i 3156, 1636, 962, 689, 459, 386i 3000, 1900, 1037, 600, 301, 346i 3211, 3057, 1279, 704, 500, 549i 3032, 1745, 1002, 597, 388, 294i

a

Ref 18. b Ref 25.



Figure 5. Contour plots of the potential surface in the vicinity of TS3 as a function of the two HCN angles. The remaining degrees of freedom are optimized at each point using CASPT2. The contour increment is 0.5 kcal/mol. The blue contours represent energies below that of TS3, the red contours represent energies above that of TS3, and the black contours correspond to the energy of TS3. (a) CASPT2/aug-cc-pvtz, (b) CAS+1+2+QC/ aug-cc-pvtz, (c) CASPT2/aug-cc-pvdz, and (d) CAS+1+2+QC/aug-cc-pvdz.

TABLE 5: Geometries and Vibrational Frequencies for TS3 from Multireference Calculationsa method (3e,3o)-CASPT2/aug-cc-pvdz (5e,5o)-CASPT2/aug-cc-pvdz (5e,5o)-CAS+1+2+QC/aug-cc-pvdz (9e,9o)-CASPT2/aug-cc-pvdz (11e,11o)-CASPT2/aug-cc-pvdz (11e,11o)-CASPT2/cc-pvtz (from ref 22) (11e,11o)-CASPT2/cc-pvtz (this work) (5e,5o)-CASPT2/aug-cc-pvtz (5e,5o)-CAS+1+2+QC/aug-cc-pvtz (9e,9o)-CASPT2/aug-cc-pvtz 2

TS3-Cs HCN2 2 TS3-C2V 2

(11e,11o)-CASPT2/aug-cc-pvtz a

RCH (Å)

RCNa (Å)

RCNb (Å)

HCNa angle (deg)

HCNb angle (deg)

ωi (cm-1)

1.128 1.148 1.148 1.149 1.161 1.138 1.147 1.108 1.111 1.110 1.152 1.146 1.120

1.349 1.305 1.305 1.309 1.310 1.296 1.287 1.309 1.306 1.308 1.284 1.282 1.305

1.277 1.292 1.289 1.291 1.286 1.283 1.274 1.309 1.306 1.308 1.274 1.282 1.305

104.5 92.1 91.7 90.4 87.2 92.0 89.0 108.0 106.6 106.1 88.7 95.2 104.7

104.8 102.5 103.1 103.3 105.2 102.5 100.8 108.0 106.6 106.1 99.5 95.2 104.7

2904, 1474, 948, 850, 559, 554i 2804, 1480, 1064, 470, 406, 459i 2790, 1399, 1058, 472, 423, 511i 2803, 1407, 1045, 454, 423, 532i 2705, 1395, 1043, 445, 384, 896i 2508, 1504, 1092, 806, 269, 385i 2748, 1444, 1067, 459, 294, 902i 2985, 1424, 1062, 1029, 656, 361i 2990, 1123, 1058, 1000, 639, 334i 2968, 1176, 1039, 971, 606, 298i 2726, 1451, 1073, 449, 294, 937i 2739, 1493, 1063, 794, 458, 226 2892, 1184, 1036, 952, 580, 243i

All structures are planar.

TABLE 6: Geometries and Vibrational Frequencies for TS3 from Single Reference Calculationsa method

RCH (Å)

RCNa (Å)

RCNb (Å)

HCNa angle (deg)

HCNb angle (deg)

ωi (cm-1)

UB3LYP/6-311G** UHF/6-311G** UMP2/6-311G** CCSD(T)/aug-cc-pvtzb RCCSD(T)/aug-cc-pvdz

1.164 1.213 1.197 1.142 1.158

1.304 1.287 1.368 1.300 1.299

1.246 1.250 1.198 1.254 1.282

85.0 68.4 77.5 91.0 92.0

112.4 131.8 119.0 106.1 96.5

2502, 1494, 1060, 542, 509, 1212i 2326, 1472, 1105, 483, 478, 2495i 2230, 1914, 852, 572, 531, 1688i

a

2729, 1625, 1097, 493, 364, 797i

b

All structures are planar. Ref 25.

bifurcation. In the aug-cc-pvdz calculations, this bifurcation occurs between Ring2 and TS3, while in the aug-cc-pvtz calculations, it occurs between TS3 and HNCN. Geometries for TS3 calculated with a number of different methods and basis sets are listed in Tables 5 and 6. All of the

structures are planar, but for the reasons discussed previously, some methods yield saddle points having C2V symmetry, while others yield structures having Cs symmetry. Consider first the multireference calculations. In Table 5, results are shown for the (3e,3o), (5e,5o), (9e,9o), and (11e,11o) active spaces. The

528 J. Phys. Chem. A, Vol. 112, No. 3, 2008 active orbitals for the (3e,3o) calculation consist of 6a14b21a2 (in C2V symmetry). Note that this active space does not include configurations c or d describing the allylic resonance. When the C2V symmetry restriction is relaxed in the (3e,3o)-CASSCF calculation, the wavefunction breaks symmetry by converging to either one of the two localized resonance structures shown previously. The minimum active space required to prevent this artificial symmetry breaking is a (5e,5o) space consisting of the following orbitals: 6a11b12b14b21a2. Note that a (7e,7o) active space is also found to undergo symmetry breaking (because one of the CN σ bond pairs is included in the active space and one is not), while the (9e,9o) and (11e,11o) active spaces do not undergo symmetry breaking. This analysis suggests that the (5e,5o), (9e,9o), and (11e,11o) active spaces all give qualitatively reasonable results, while the (3e,3o) and (7e,7o) active spaces suffer from symmetry breaking artifacts. Also shown in Table 5 are comparisons of results obtained with three different basis sets, aug-cc-pvdz, cc-pvtz, and augcc-pvtz. The two smaller basis sets yield saddle point structures of Cs symmetry. The largest basis set, aug-cc-pvtz, is found to yield a transition state structure having C2V symmetry. With the largest basis set, we also find a dependence on the size of the active space. With the (5e,5o) and (9e,9o) active spaces, we find a single saddle point of C2V symmetry, while the (11e,11o) active space yields a shallow C2V minimum corresponding to a ring-opened intermediate, 2HCN2, a C2V saddle point between 2HCN2 and Ring2, and two Cs symmetry saddle points between 2HCN2 and HNCN. The relative energies of 2HCN2, TS3-Cs and TS3-C2V, in this calculation are 0.00, 0.04, and 1.01 kcal/mol, respectively. Consider now the TS3 structures obtained from single reference calculations (Table 6). B3LYP, UHF, and UMP2 are all found to give structures very distorted from C2V symmetry. The problem here is that the HF wavefunction yields a very poor description of the tri-radical character found in the multireference calculations. In fact, the two lowest RHF wavefunctions at the geometry of the (5e,5o)-CASPT2/aug-ccpvtz saddle point are a 2B2 and a 2A1 state having the orbital occupancies shown in d and e

Harding et al.

Figure 6. RHF and CCSD(T) energies and T1 diagnostics along a C2V MEP for ring opening of Ring2. The geometries are obtained from C2V constrained (5e,5o)-CASPT2/aug-cc-pvtz geometry optimizations. Solid lines are RHF energies, dashed lines are CCSD(T) energies, and dotted lines are the T1 diagnostics. The red curves correspond to configuration a and the blue to configuration b. The vertical black line denotes the position of TS3 from the (5e,5o)-CASPT2/aug-cc-pvtz calculation.

1a122a123a124a125a121b121b222b223b221a224b21 (d) 1a122a123a124a125a121b121b222b223b221a226a11 (e) with four π electrons rather than three. The CCSD(T) calculation, when constrained to have three π electrons, yields a saddle point structure only slightly distorted from C2V and in reasonable agreement with the structures from the multireference calculations. There is, however, an ambiguity in the CCSD(T) calculations at least when applied to C2V geometries in that the calculations can be based on either configuration a or configuration b. In Figure 6, we plot the RHF and CCSD(T) energies and T1 diagnostics corresponding to configurations a and b along a C2V constrained minimum energy path (MEP) for ring opening of Ring2. Both CCSD(T) calculations exhibit divergent behavior in the vicinity of TS3. In Figure 7, we plot the results of analogous UHF based calculations. The curves for configurations a and b are quite similar to those in Figure 6. Also shown in Figure 7 are calculations based on the lowest UHF solution at each point. In the vicinity of TS3, the lowest UHF solution is a highly spin contaminated A1 state (S2 ) 1.9). The UHF wavefunction is found to undergo symmetry breaking quite close to the Ring2 minimum. Given the extreme sensitivity of the TS3 structure to the electronic structure method used, none of the calculations shown

Figure 7. UHF and CCSD(T) energies along a C2V MEP for ring opening of Ring2. The geometries are obtained from C2V constrained (5e,5o)-CASPT2/aug-cc-pvtz geometry optimizations. Solid lines are UHF energies, and dashed lines are CCSD(T) energies. The red curves correspond to configuration a, the blue to configuration b, and the green to the lowest UHF solution. The vertical black line denotes the position of TS3 from the (5e,5o)-CASPT2/aug-cc-pvtz calculation.

in Tables 5 and 6 can be considered to be definitive regarding the question of a C2V or Cs saddle point. However, even if the actual saddle point is Cs, as in Figure 5c,d, variational effects in the transition state theory calculation will tend to shift the dynamical bottleneck off the saddle point, toward the C2V axis and toward Ring2, where the reaction valley is more constrained (higher frequencies and lower reaction path degeneracy). For this reason, we now focus on determining the energy of TS3C2V as accurately as possible. As noted in section 3.1.1., the energies of the stationary points (except for TS3) are calculated at the CCSD(T)/aug-cc-pvqz// CASPT2/aug-cc-pvtz level. For TS3, the above analysis and the large T1 diagnostic found in the CCSD(T)/aug-cc-pvqz//



TABLE 7: Energy Estimates from Multireference aug-cc-pvqz Calculations for TS3-C2W (kcal/mol)a method

Nact

CASPT2b CASPT2c CAS+1+2 CAS+1+2+QC CASPT2b CAS+1+2 CAS+1+2+QC

11 11 11 11 5 5 5

scheme 1 (direct)

scheme 2 (via Ring2)

scheme 3 (via quartet)

2.8(2.8) 3.4(3.4) 14.3(14.0) 10.9(10.7)

7.1(6.1) 7.2(6.2) 11.3(10.3) 10.5(9.3) 6.6(7.1) 11.3(11.7) 11.1(10.8)

11.1(11.4) 11.1(11.4) 11.9(12.1) 10.9(11.2) 10.3(10.6) 11.5(11.8) 10.8(10.6)

a Numbers in parentheses were obtained using the aug-cc-pvtz basis set. b Uncontracted. c Contracted.

CASPT2/aug-cc-pvtz calculation (0.063 for configuration a and 0.027 for configuration b) suggest that the CCSD(T) calculation may not be reliable for this saddle point. We note that several other stationary points listed in Table 2 have unusually large T1 diagnostics (TS2, 0.048; TS5, 0.041; and HCNN, 0.040), but the kinetics is not sensitive to these three stationary points. We also note that the T1 diagnostic, having been developed for CCSD calculations, is not a rigorous indication of problems with the CCSD(T) calculations. However, as noted above, the CCSD(T) energies show clear evidence of divergent behavior in the vicinity of TS3, and for this reason, we choose to rely on multireference methods for estimating the energy and rovibrational properties of this saddle point. Three schemes for calculating the energy of TS3-C2V relative to reactants were used: (i) direct multireference calculations of TS3-C2V relative to CH + N2, (ii) multireference calculation of TS3-C2V relative to Ring2 combined with CCSD(T) calculation of Ring2 relative to CH + N2, and (iii) multireference calculation of TS3-C2V relative to the quartet state, 4HCN2 (which lies close to TS3-C2V), combined with CCSD(T) calculation of the quartet relative to CH + N2. The results are summarized in Table 7. All calculations employ the (11e,11o)-CASPT2/aug-cc-pvtz C2V geometry for TS3-C2V given in Table 5. The single point calculations in Table 7 were all obtained using the aug-cc-pvqz basis set. The CASPT2 based estimates were found to vary from 2.8 to 11.1 kcal/mol depending on the scheme used and the size of the active space. The CAS+1+2+QC estimates show much smaller variations, ranging from 10.5 to 10.9 kcal/mol. In our kinetics calculations, we used a value of 10.9 kcal/mol (the average of the (11e,11o)-CAS+1+2+QC estimates extrapolated to the complete basis set limit). Inclusion of zero point corrections gives the result shown in Table 3 (12.9 kcal/ mol). It was not feasible to perform a (11e,11o)-CAS+1+2+QC/ aug-cc-pvqz calculation for TS3-Cs. Instead, we used the (11e,11o)-CASPT2/CBS calculation to fix the energy of TS3-Cs relative to TS3-C2V (-1.0 kcal/mol, without zero point and -1.9 kcal/mol, including zero point). However, the kinetics calculations indicate that TS3-Cs has no effect on the kinetics. These energy estimates (11.0 kcal/mol for TS3-Cs and 12.9 kcal/mol for TS3-C2V) are in good agreement with the earlier calculations of Moskaleva et al. (11.7 kcal/mol) and are in reasonable agreement with the results of Berman et al. (14.4 kcal/mol), although the geometry from the multireference calculations is very different from either of the earlier single reference calculations. Interestingly, our own CCSD(T) calculations (based on ROHF reference wavefunctions and the CASPT2 C2V geometry) yield energies of 12.2 or 17.7 kcal/mol depending on whether configuration a or b is used. The good agreement between the CCSD(T) calculation using configuration a and the multireference calculations must be considered as fortuitous

Figure 8. Contour plots of the H + NCN interaction potential. The NCN fragment is constrained to be collinear, and the plotting plane contains this molecular axis. The red contours denote positive energies, and the blue contours denote negative energies. Zero energy is defined to be the energy of the H + NCN asymptote. The contour increment is 1 kcal/mol. The blanked circle in the center hides parts of the surface that are not relevant to the association kinetics and were therefore not calculated. The left side (a) was calculated using (7e,7o)-CAS+1+2+QC/ aug-cc-pvtz and the right side (b) with (7e,7o)-CASPT2/aug-cc-pvtz.

because, as seen in Figure 6, the CCSD(T) energy for configuration a is clearly beginning to diverge for this geometry. 3.1.4. HNCN f NCN + H. The final step, the dissociation of HNCN to H + NCN, is the reverse of a barrierless radicalradical addition. The rate of this reaction will be determined by the long-range interaction potential between the H atom and the NCN fragment. A two-dimensional contour plot of this interaction potential is shown in Figure 8. As expected, the plot shows barrierless paths for addition to either nitrogen atom. Note that the left side of this plot is calculated using multireference CI, while the right side was calculated with second order perturbation theory. The fact that the two sides of the plot are nearly identical demonstrates that the two multireference methods are in agreement with each other. NCN is formed in the ground, 3Σg-, state. Consequently, the potential shown in Figure 8 is cylindrically symmetric, containing all the information needed for the variational transition state theory calculation. In Figure 9, we show one-dimensional potential curves for this NH bond cleavage reaction using both multireference and single reference, CCSD(T), methods. For the CCSD(T) calculations, we show potential curves calculated using both UHF orbitals, as implemented in Gaussian 98, and ROHF orbitals, as implemented in MOLPRO. The ROHF-RCCSD(T) potential curve exhibits a spurious barrier to the reverse reaction, while the UHF-CCSD(T) calculations are significantly less attractive than the multireference calculations in the kinetically sensitive region, RNH ) 2-3 Å. For example, at RNH ) 2.5 Å, the UHF/ CCSD(T), CASPT2, and CAS+1+2+QC interaction energies (relative to H + NCN) are -0.35, -1.2, and -1.1 kcal/mol, respectively. In this region, the UHF wavefunction is highly spin contaminated, having an S2 expectation value of 2.2 instead of the expected value for a pure doublet of 0.75. These results demonstrate, not unexpectedly, that CCSD(T) calculations are not appropriate for this part of the potential surface. See ref 49 for a recent review of the problems associated with using coupled-cluster theory for bond breaking reactions. 3.2. Kinetics Results. 3.2.1. Path A. The present theoretical predictions for path A are illustrated in Figure 10. These predictions are for the collisionless limit of the rate coefficient for the CH + N2 f H + NCN reaction. Included in the plot are the predictions obtained when considering each transition state individually, as well as the full calculation, which properly

530 J. Phys. Chem. A, Vol. 112, No. 3, 2008

Figure 9. One-dimensional potential curves for H + NCN f HNCN. The NCN radical is constrained to its equilibrium geometry. The HNC angle is fixed at 100°. All calculations employ the aug-cc-pvdz basis set. The solid lines are for the doublet state and the dashed line for the quartet state. The asymptotic energy for the ROHF-RCCSD(T) calculations is taken from the quartet state.

Figure 10. Plot of the predicted rate coefficient for CH + N2 f NCN + H via path A when considering different subsets of the transition states along the path.

incorporates the joint effect of each of the transition states. At low temperatures, the decomposition of HNCN to form H + NCN is clearly the rate-limiting step, while at higher temperature, TS3 provides the rate-limiting bottleneck. At intermediate temperatures (i.e., from about 500 to 1500 K), both transition states play an important role in the kinetics. The entrance channel transition state, TS1, is never the dominant bottleneck, although it does yield a minor reduction in the rate coefficient at the highest temperatures. Notably, the implementation of the CCSD(T) properties for TS3 yields a rate constant that is 2 times larger at 2000 K. 3.2.2. Path B. The predicted temperature dependence of the rate coefficient for reaction via path B is illustrated in Figure 11, together with the corresponding results for path A and for a calculation in which both path A and path B are considered jointly. For path B, TS7, the transition state for the initial addition to form HNNC, is the dominant bottleneck. The TS7 barrier energy of 19.5 kcal/mol is substantially greater than the

Harding et al.

Figure 11. Plot of the contributions from different pathways to the predicted rate coefficient for CH + N2 f NCN + H.

corresponding maximum saddle point energy of 12.9 kcal/mol for path A. Thus, the contribution from path B increases with increasing temperature, with the ratio of the path B to path A rate coefficients varying from 0.21, to 0.41, to 0.60 as the temperature increases from 1000, to 2000, and then to 3000 K. The two contributions are not directly additive, and so the joint rate coefficient is increased by the slightly reduced factors of 1.13, 1.34, and 1.53 for these same temperatures. 3.2.3. 4CH + N2. Another pathway that has not been considered previously involves the reaction on solely the quartet state. With this pathway, the spin crossing may be more efficient as it would arise through collision induced intersystem crossing of CH (2CH + M f 4CH + M). The quartet reaction proceeds via first the formation of a quartet HCNN complex, followed by NN fission to produce HCN + 4N. The entrance and exit barriers for this pathway are both quite large, and both present significant bottlenecks to the reaction. At the CCSD(T)/CBS level, they are 33.9 and 34.9 kcal/mol, respectively, relative to 2CH + N . As a result, the predicted contribution from the purely 2 quartet pathway is quite small, being less than 16% for temperatures below 3000 K. 3.2.4. Comparison with Experiment. There have been a number of experimental studies of the CH + N2 reaction at high temperatures. In particular, Lindackers et al.,50 Dean et al.,5 and, most recently, Hanson and co-workers24 have provided fairly direct studies of the kinetics for temperatures of 2000 K and higher. The last of these studies also examined the product distribution, concluding that H + NCN are the principal products, with a conservative lower bound to their branching ratio of 0.7. Meanwhile, Becker et al.51 directly studied the kinetics in the 850-890 K range. A number of more indirect studies have also been presented.4,52-54 Interestingly, a recent modeling study23 has suggested that the rate to produce NCN + H must be 1-2 orders of magnitude greater than predicted by Moskaleva and Lin17 to adequately predict flame species in low temperature flames (1000-1500 K). In Figure 12, the present theoretical predictions for the rate coefficient are compared with the various experimental observations as well as the prior theoretical predictions of Moskaleva and Lin.17 The present predictions are in remarkably good agreement with the very recent experimental measurements of Hanson and co-workers,24 which in turn are representative of the body of high temperature experimental data. Furthermore,


Figure 12. Plot of the experimental measurements and theoretical predictions for the CH + N2 f NCN + H rate coefficient. The black symbols denote direct experimental results, the colored lines and symbols denote modeling results, while the black solid and dashed lines denote the theoretical predictions.

near the key temperature of 1500 K, the present predictions are in quantitative agreement with the best fit from the recent modeling study of El Bakali and co-workers23 but do diverge somewhat toward lower temperatures. The experimental results of Becker et al.51 in the 850-900 K range are 2 orders magnitude greater than the present predictions. This discrepancy may be indicative of the difficulty of disentangling the addition and bimolecular reaction rate coefficients at these temperatures. The theoretical predictions of Moskaleva and Lin17 are a factor of 5 lower than the theoretical and experimental results at 2000 K. Furthermore, the discrepancy between the two theoretical predictions increases with decreasing temperature, reaching a factor of 30 at 1000 K. It is unclear as to why the two theoretical predictions are so discordant. The difference of 2.0 kcal/mol in the H + NCN reaction endothermicity provides some explanation, but this should correlate with only a factor of 2.7 difference at 1000 K. The present VRC-TST treatment of the H + NCN channel should be much more accurate than the approach employed by Moskaleva and Lin. Nevertheless, a factor of 10 error and such a strong temperature dependence to the discrepancy is somewhat surprising. Notably, the present multireference treatment of TS3 should actually yield a reduction in the rate coefficient because the barrier is higher and the entropy is lower. It is perhaps worthwhile to consider the dominant sources of uncertainty in the theoretical predictions. The uncertainties in the TS3 saddle point energy and in the H + NCN channel energy have the greatest affect on the kinetics since they are the two dominant bottlenecks. The effect of varying the TS3 and H + NCN energies by (2 kcal/mol is illustrated in Figure 13. These variations are more or less representative of the true uncertainties in the present theoretical estimates of these quantities. Note that a photodissocation study of Neumark and co-workers55 estimated the NCN heat of formation (at 0 K) as 111.4 ( 0.7 kcal/mol, which implies a CH + N2 f H + NCN reaction endothermicity of 21.2 ( 0.7 kcal/mol. The present theoretical estimate is 1.7 kcal/mol lower.


Figure 13. Plot of the sensitivity of the theoretical predictions for the CH + N2 f NCN + H rate coefficient to the TS3 and H + NCN energies. The blue, red, and green lines denote changes by (2 kcal/ mol of TS3, of H + NCN, or of both TS3 and H + NCN, respectively.

The TS3 variation changes the overall rate constant by a more or less constant factor of 1.25 over the 1000-3000 K range. In contrast, the variation in the H + NCN channel energy has a larger but more temperature-dependent effect; at 1000 K, it changes the rate by a factor of 1.7, while at 3000 K, it is only a 4% effect. This difference in behavior is related to the changing importance of the H + NCN and TS3 transition states. The joint variation in both energies implies uncertainties in the present rate predictions of factors of 2.4, 1.4, and 1.2 at 1000, 2000, and 3000 K, respectively. 4. Concluding Remarks The present ab initio chemical kinetics predictions for the CH + N2 f NCN + H rate coefficient are in quantitative agreement with the recent experimental measurements of Hanson and co-workers24 for temperatures of 1940 K and higher. Furthermore, these predictions provide strong support for the magnitude of the rate constant proposed in the recent modeling study of El Bakali et al., at least near the practically important temperature of 1500 K.23 The present predictions are wellreproduced over the 400-3000 K range by the modified Arrhenius form 6.80 × 10-16 T1.122 exp(-8819/T) cm3 molecule-1 s-1, with T in Kelvin. These predictions are expected to be accurate to within about a factor of 2 for the practically important temperature range of 1000-2000 K. Accurate measurements of the heat of formation for NCN would greatly reduce the uncertainties in these predictions. The application of multireference electronic structure methods was central to the resolution of the discrepancy between the recent experiments of Hanson and co-workers24 and the prior theoretical results of Moskaleva and Lin.17 The implementation of such multireference calculations within the VRC-TST approach was also a key component of the analysis, particularly at lower temperatures. At higher temperatures, TS3, the transition state for the formation of HNCN from its ring precursor, provides the dominant bottleneck. The HNNC pathway of Berman and co-workers25 makes a significant contribution to the kinetics at higher temperatures. In contrast, the contribution of a purely quartet pathway is negligible, at least for temperatures below 3000 K.

532 J. Phys. Chem. A, Vol. 112, No. 3, 2008 Acknowledgment. This work was supported by the Division of Chemical Sciences, Geosciences, and Biosciences, the Office of Basic Energy Sciences, the U.S. Department of Energy. The work at Argonne was supported under DOE Contract DE-AC0206CH11357. Sandia is a multiprogram laboratory operated by Sandia Corporation, a Lockheed Martin Company, for the United States Department of Energy under Contract DE-AC0494-AL85000. References and Notes (1) Fenimore, C. P. Proc. Combust. Inst. 1971, 13, 373. (2) Hayhurst, A. N.; Vince, I. M. Prog. Energy Combust. Sci. 1980, 6, 35. (3) Miller, J. A. Proc. Combust. Inst. 1996, 26, 461. (4) Miller, J. A.; Bowman, C. T. Prog. Energy Combust. Sci. 1989, 15, 287. (5) Dean, A. J.; Hanson, R. K.; Bowman, C. T. Proc. Combust. Inst. 1991, 23, 259. (6) Lindackers, D.; Burmeister, M.; Roth, P. Proc. Combust. Inst. 1991, 23, 251. (7) Glarborg, P.; Miller, J. A.; Kee, R. J. Combust. Flame 1986, 65, 177. (8) Bair, R. A. Annual Report; Theoretical Chemistry Group (Chemistry Division), Argonne National Laboratory: Argonne, IL, September 1985. (9) Manaa, R. M.; Yarkony, D. R. J. Chem. Phys. 1991, 95, 1808. (10) Manaa, R. M.; Yarkony, D. R. Chem. Phys. Lett. 1992, 188, 352. (11) Martin, J. M. L.; Taylor, P. R. Chem. Phys. Lett. 1993, 209, 143. (12) Walch, S. P. Chem. Phys. Lett. 1993, 208, 214. (13) Miller, J. A.; Walch, S. P. Int. J. Chem. Kinet. 1997, 29, 253. (14) Rodgers, A. S.; Smith, G. P. Chem. Phys. Lett. 1996, 253, 313. (15) Miller, J. A.; Pilling, M. J.; Troe, J. Proc. Combust. Inst. 2005, 30, 43. (16) Cui, L.; Morokuma, K.; Bowman, J. M.; Klippenstein, S. J. J. Chem. Phys. 1999, 110, 9469. (17) Moskaleva, L. V.; Lin, M. C. Proc. Combust. Inst. 2000, 28, 2393. (18) Moskaleva, L. V.; Xia, W. S.; Lin, M. C. Chem. Phys. Lett. 2000, 331, 269. (19) Smith, G. P. Chem. Phys. Lett. 2003, 367, 541. (20) Bise, R. T.; Hoops, A. A.; Neumark, D. M. J. Chem. Phys. 2001, 114, 9000. (21) Szpunar, D. E.; Faulhaber, A. E.; Kautzman, K. E.; Crider, P. E., II; Neumark, D. M. J. Chem. Phys. 2007, 126, 114311. (22) Takayanagi, T. Chem. Phys. Lett. 2003, 368, 393. (23) El Bakali, A.; Pillier, L.; Desgroux, P.; Lefort, B.; Gasnot, L.; Pauwels, J. F.; da Costa, I. Fuel 2006, 85, 896. (24) Vasudevan, V.; Hanson, R. K.; Bowman, C. T.; Golden, D. M.; Davidson, D. F. J. Phys. Chem. A 2007, in press. (25) Berman, M. R.; Tsuchiya, T.; Gregusova, A.; Perera, A.; Bartlett, R. J. J. Phys. Chem. A 2007, 111, 6894. (26) Dunning, T. H., Jr. J. Chem. Phys. 1989, 90, 1007. (27) Kendall, R. A.; Dunning, T. H., Jr.; Harrison, R. J. J. Chem. Phys. 1992, 96, 6796. (28) Woon, D. E.; Dunning, T. H., Jr. J. Chem. Phys. 1993, 98, 1358. (29) MOLPRO is a package of ab initio programs written by Werner, H.-J. and Knowles, P. J. with contributions from Almlof, J.; Amos, R. D.;

Harding et al. Berning, A.; Cooper, D. L.; Deegan, M. J. O.; Dobbyn, A. J.; Eckert, F.; Elbert, S. T.; Hampel, C.; Lindh, R.; Lloyd, A. W.; Meyer, W.; Nicklass, A.; Peterson, K.; Pitzer, R.; Stone, A. J.; Taylor, P. R.; Mura, M. E.; Pulay, P.; Schutz, M.; Stoll, H.; Thorsteinsson, T. The majority of calculations reported here was performed with version 2002.6. (30) Werner, H.-J.; Knowles, P. J. J. Chem. Phys. 1985, 82, 5053. Knowles, P. J.; Werner, H.-J. Chem. Phys. Lett. 1985, 115, 259. (31) Celani, P.; Werner, H.-J. J. Chem. Phys. 2000, 112, 5546. (32) Werner, H.-J.; Knowles, P. J. J. Chem. Phys. 1988, 89, 5803. (33) Knowles, P. J.; Werner, H.-J. Chem. Phys. Lett. 1988, 145, 514. (34) Knowles, P. J.; Hampel, C.; Werner, H.-J. J. Chem. Phys. 1993, 99, 5219. (35) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Zakrzewski, V. G.; Petersson, G. A.; Montgomery, J. A., Jr.; Stratmann, R. E.; Burant, J. C.; Dapprich, S.; Millam, J. M.; Daniels, A. D.; Kudin, K. N.; Strain, M. C.; Farkas, O.; Tomasi, J.; Barone, V.; Cossi, M.; Cammi, R.; Mennucci, B.; Pomelli, C.; Adamo, C.; Clifford, S.; Ochterski, J.; Petersson, G. A.; Ayala, P. Y.; Cui, Q.; Morokuma, K.; Malick, D. K.; Rabuck, A. D.; Raghavachari, K.; Foresman, J. B.; Cioslowski, J.; Ortiz, J. V.; Stefanov, B. B.; Liu, G.; Liashenko, A.; Piskorz, P.; Komaroni, I.; Gomperts, R.; Martin, R. L.; Fox, D. J.; Keith, T. A.; Al-Laham, M. A.; Peng, C. Y.; Nanayakkara, A.; Gonzalez, C.; Challacombe, M.; Gill, P. M. W.; Johnson, B. G.; Chen, W.; Wong, M. W.; Andreas, J. L.; Head-Gordon, M.; Replogle, E. S.; Pople, J. A. Gaussian 98; Gaussian, Inc.: Pittsburgh, PA, 1998. (36) Martin, J. M. L.; Uzan, O. Chem. Phys. Lett. 1998, 282, 16. (37) Celani, P.; Werner, H.-J. J. Chem. Phys. 2003, 119, 5044. (38) Klippenstein, S. J. J. Chem. Phys. 1992, 96, 367. (39) Georgievskii, Y.; Klippenstein, S. J. J. Chem. Phys. 2003, 118, 5442. (40) Georgievskii, Y.; Klippenstein, S. J. J. Phys. Chem. A 2003, 107, 9776. (41) Harding, L. B.; Klippenstein, S. J.; Georgievskii, Y. J. Phys. Chem. A 2007, 111, 3789. (42) Harding, L. B.; Georgievskii, Y.; Klippenstein, S. J. J. Phys. Chem. A 2005, 109, 4646. (43) Berman, M. R.; Lin, M. C. J. Phys. Chem. 1983, 87, 3993. (44) Brownsword, R. A.; Herbert, L. B.; Smith, I. W. M.; Stewart, D. W. A. J. Chem. Soc., Faraday Trans. 1996, 92, 1087. (45) Fulle, D.; Hippler, H. J. Chem. Phys. 1996, 105, 5423. (46) Le Picard, S. D.; Canosa, A. Geophys. Res. Lett. 1998, 25, 485. (47) Hahn, D. K.; Klippenstein, S. J.; Miller, J. A. Faraday Discuss. 2001, 119, 79. (48) Lee, T. J.; Taylor, P. R. Int. J. Quantum Chem. 1989, S23, 199. (49) Barlett, R. J.; Musial, M. ReV. Mod. Phys. 2007, 79, 291. (50) Lindackers, D.; Burmeister, M.; Roth, P. Proc. Combust. Inst. 1991, 23, 251. (51) Becker, K. H.; Engelhardt, B.; Geiger, H.; Kurtenbach, R.; Schrey, G.; Wiesen, P. Chem. Phys. Lett. 1992, 195, 322. (52) Blauwens, J.; Smets, B.; Peeters, J. Proc. Combus. Inst. 1977, 16, 1055. (53) Matsui, Y.; Yuuki, A. Jpn. J. Appl. Phys. 1985, 24, 598. (54) Matsui, Y.; Nomaguchi, N. Combus. Flame 1978, 32, 205. (55) Bise, R. T.; Choi, H.; Neumark, D. M. J. Chem. Phys. 1999, 111, 4923.

Revisited with Multireference Methods - American Chemical Society

Recommend Documents