Reaction of Phenyl Radicals with Acetylene: Quantum Chemical

Aug 23, 2003 - Reaction of Phenyl Radicals with Acetylene: Quantum Chemical Investigation of the Mechanism and Master Equation Analysis of the Kinetic...
0 downloads 0 Views 305KB Size
Download PDF
Published on Web 08/23/2003

Reaction of Phenyl Radicals with Acetylene: Quantum Chemical Investigation of the Mechanism and Master Equation Analysis of the Kinetics I. V. Tokmakov and M. C. Lin* Contribution from the Department of Chemistry, Emory UniVersity, Atlanta, Georgia 30322 Received February 18, 2003; E-mail: [email protected]

Abstract: The mechanism of the C6H5 + C2H2 reaction has been investigated by various quantum chemical methods. Electrophilic addition to the CC triple bond is found to be the only important mode of phenyl radical attack on acetylene. The initially formed chemically activated C6H5C2H2 adducts may follow several isomerization pathways in competition with collisional stabilization and H-elimination. Thermochemistry of various decomposition and isomerization channels is evaluated by the G2M method. For key intermediates, the following standard enthalpies of formation have been deduced from isodesmic reactions: 94.2 ( 2.0 kcal/mol (C6H5CHCH), 86.4 ( 2.0 kcal/mol (C6H5CCH2), and 95.5 ( 1.8 kcal/ mol (o-C6H4C2H3). The accuracy of theoretical predictions was examined through extensive comparisons with available experimental and theoretical data. The kinetics and product branching of the C6H5 + C2H2 reaction have been evaluated by weak collision master equation/Rice-Ramsperger-Kassel-Marcus (RRKM) analysis of the truncated kinetic model including only kinetically important transformations of the isomeric C8H7 radicals. Available experimental kinetic data can be quantitatively reproduced by calculation with a minor adjustment of the C6H5 addition barrier from 3.7 to 4.1 kcal/mol. Our predicted total rate constant, kR1 ) (1.29 × 1010)T0.834 exp(-2320/T) cm3 mol-1 s-1, is weakly dependent on P and corresponds to the phenylation process under combustion conditions (T > 1000 K).

I. Introduction

Over the past years, a growing concern about the quality of air in the industrial and urban zones has been reflected in tightening the pollution and emission controls, specifically targeting HC (hydrocarbon) air toxics, which constitute a class of the most widespread atmospheric pollutants. Polycyclic aromatic hydrocarbons (PAHs) account for the largest and very dangerous portion of the HC emissions because many of them are potent mutagens1,2 and carcinogens.3 And since they escape the combustion sources in a form of easily inhaled ultrafine particles (d < 0.1 µm), their presence in the atmosphere can cause acute and long-term respiratory effects. In addition, air deposition of PAH on soil and directly on plants is likely to be responsible for contamination of such food products as cereal and vegetables, which then become the main PAH sources in a human diet.4,5 A better understanding of the mechanisms of PAH and soot formation during HC combustion is essential for the development of more efficient combustion devices with minimal environmental impact. (1) Durant, J. L.; Busby, W. F.; Lafleur, A. L.; Penman, B. W.; Crespi, C. L. Mutat. Res. 1996, 371, 123-157. (2) Samanta, S. K.; Singh, O. V.; Jain, R. K. Trends Biotechnol. 2002, 20, 243-248. (3) Menzie, C. A.; Potocki, B. B.; Santodonato, J. EnViron. Sci. Technol. 1992, 26, 1278-84. (4) Oleszczuk, P. Pol. Arch. Ochr. Srodowiska 2002, 28, 107-118. (5) Enzminger, J. D.; Ahlert, R. C. EnViron. Technol. Lett. 1987, 8, 269-78. 10.1021/ja0301121 CCC: $25.00 © 2003 American Chemical Society

Chemical reactions responsible for the formation of PAH in hydrocarbon combustion and pyrolysis have been reviewed recently.6-8 At the initial stage, the first aromatic ring is formed from small aliphatic radicals and molecules either via radicalmolecule addition pathways (e.g., n-C4H3 + C2H2 f C6H5, n-C4H5 + C2H2 f C6H6 + H, n-C4H4 + C2H3 f C6H6 + H, etc.) or from radical-radical recombinations (e.g., C3H3 + C3H3 f C6H6/C6H5 + H, c-C5H5 + CH3 f C6H6 + H + H, etc.). The following stage is largely attributed to the HACA (Habstraction-acetylene-addition) mechanism,9-11 which represents the molecular growth of PAH as sequential additions of the C2H2 building blocks to aryl radicals (Ai) generated from the corresponding aromatic molecules (AiH) by H-abstraction: AiH + R f Ai + RH (most importantly, R ) H, OH). In the absence of reliable kinetic data for the reactions of large aromatics, present kinetic models typically utilize the rate constants determined for prototype reactions. A natural choice of the prototype reactions for the HACA mechanism is to (6) Frenklach, M. Phys. Chem. Chem. Phys. 2002, 4, 2028-2037. (7) Richter, H.; Howard, J. B. Prog. Energy Combust. Sci. 2000, 26, 565608. (8) Richter, H.; Howard, J. B. Phys. Chem. Chem. Phys. 2002, 4, 2038-2055. (9) Bockhorn, H.; Fetting, F.; Wenz, H. W. Ber. Bunsen-Ges. Phys. Chem. 1983, 87, 1067. (10) Frenklach, M.; Clary, D. W.; Gardiner, W. C.; Stein, S. E. Proc. Combust. Inst. 1984, 20, 887-901. (11) Frenklach, M.; Warnatz, J. Combust. Sci. Technol. 1987, 51, 265-83. J. AM. CHEM. SOC. 2003, 125, 11397-11408

9

11397

ARTICLES

Tokmakov and Lin

Table 1. Experimental Kinetic Data on C6H5 + C2H2 f Products kR1/cm3 mol-1 s-1

T/K

P

(4.0 × 1013) exp(-5082/T) (1.6 × 1013) exp(-4730/T) (2.2 × 1013) exp(-4730/T) e5 × 108

method

1000-1330 1.3-13 µbar VLPP/MS 1000-1330 1.3-13 µbar reevaluationa 1000-1330 1.3-13 µbar reevaluationb 297 laser photolysis/ laser absorption 8 e6 × 10 489 0.05 bar 11 (2.2 × 10 ) exp(-1560/T) 297-523 0.03 bar CRDS (1.0 × 1013) exp(-3850/T) 1130-1430 ∼3 bar shock tube/ UV absorption

16 18 17 18 19

a Taking into account both 1986 and 1988 sets of product yields measured in the experiments of Stein et al.16 and assuming kR2 ) 3.0 × 1012 cm3 mol-1 s-1. b Using the kR2 ) 5.7 × 1012 cm3 mol-1 s-1 recommended by Heckmann et al.19

consider reactions of the smallest aromatic radical (A1 ) phenyl):

C6H5 + C2H2 f products

(R1)

for C2H2 additions and C6H6 + R f C6H5 + RH for H-abstractions. Needless to say, accurate mechanism and kinetics for these prototype reactions are of great importance to the large-scale modeling of the PAH growth in combustion. Reactions of benzene with some typical combustion radicals (e.g., R ) H, CH3, OH) were the subjects of our previous investigations,12-15 where we had reviewed available experimental data and examined different theoretical approaches to the calculation of kinetic and thermodynamic parameters from first principles. In this work, our top choice methodologies will be employed to study the more complex mechanism and kinetics of reaction R1. But first we survey available experimental kinetic data16-19 (Table 1), which can be used to test the reliability of theoretical predictions. Stein and co-workers16 used the very-low-pressure pyrolysis (VLPP) of C6H5NO and C6H5SO2C2H3 to generate phenyl in a flow reactor connected to a quadrupole mass spectrometer. They determined the kR1 rate constant relative to the assumed rate constant, kR2 ) 3.2 × 1012 cm3 mol-1 s-1, for the recombination of the C6H5 radicals:

C6H5 + C6H5 f C6H5C6H5

(12) Mebel, A. M.; Lin, M. C.; Yu, T.; Morokuma, K. J. Phys. Chem. A 1997, 101, 3189. (13) Tokmakov, I. V.; Park, J.; Gheyas, S.; Lin, M. C. J. Phys. Chem. A 1999, 103, 3636. (14) Tokmakov, I. V.; Lin, M. C. Int. J. Chem. Kinet. 2001, 33, 633-653. (15) Tokmakov, I. V.; Lin, M. C. J. Phys. Chem. A 2002, 106, 11309-11326. (16) Fahr, A.; Mallard, W. G.; Stein, S. E. Proc. Int. Symp. Combust. 1986, 21, 825-831. (b) Fahr, A.; Stein, S. E. Proc. Int. Symp. Combust. 1988, 22, 1023-1029. (17) Preidel, M.; Zellner, R. Ber. Bunsen-Ges. Phys. Chem. 1989, 93, 1417. (18) Yu, T.; Lin, M. C.; Melius, C. F. Int. J. Chem. Kinet. 1994, 26, 10951104. (19) Heckmann, E.; Hippler, H.; Troe, J. Proc. Int. Symp. Combust. 1996, 26, 543-550. 9

VOL. 125, NO. 37, 2003

Preidel and Zellner17 attempted to measure the phenyl radical kinetics by monitoring the continuous-wave laser absorption signal at 488 nm. This technique did not permit a reliable determination of the relatively slow C6H5 + C2H2 reaction rate. The first successful direct measurement of the total rate of reaction R1 was reported by Yu et al.,18 using the cavity ringdown spectrometry (CRDS) technique. They also gave a theoretical interpretation of their low-T experimental results and the high-T kinetic data of Stein et al.16 in terms of RRKM theory (Scheme 1), employing the energetic and molecular parameters either computed by the BAC-MP4 method (for stable molecules) or adjusted to reproduce the room T rate constant (for transition states). Wang and Frenklach20 considered a similar one-well RRKM model (Scheme 1) but used corrected semiempirical (AM1) energetic and molecular parameters. Apparently unaware of the results of Yu et al.,18 they fitted the transition state parameters to reproduce the less reliable room-temperature rate constant of Preidel and Zellner.17 The experimental studies described above provided important benchmark values of the total rate constant kR1, but to the best of our knowledge no attempt was made to characterize the intermediates or products of reaction R1 other than C6H5C2H. Some of the possible transformations of C6H5CHCH have been considered in two recent computational studies.21,22 Richter et al.21 performed density functional theory (DFT) calculations to construct the potential energy surface (PES) of reaction R1 and a QRRK analysis to derive P, T-dependent branching rate constants, whereas Moriarty et al.22 tested the capability of several quantum chemical methods (PM3, MP2, B3LYP, CASPT2, and G2MP2) to predict accurate energetics for the reversible isomerization between the C6H5CHCH and o-C6H4C2H3 radicals:

(R2)

More recently, the high T reactions of phenyl were studied by Heckmann et al.19 in reflected shock waves. They reported a higher value of kR2(T ) 1050-1450 K) ) 5.7 × 1012 cm3 mol-1 s-1. Upon reevaluation (with the higher value of kR2), the kR1 rate constants of Stein et al.16 show excellent agreement with the kinetic data of Heckmann et al.,19 despite some differences in the Arrhenius parameters listed in Table 1. Both Stein et al.16 and Heckmann et al.19 concluded that C6H5CCH + H (channel R1d in Scheme 1) are the only products of reaction R1 under their high T conditions.

11398 J. AM. CHEM. SOC.

Scheme 1

ref

In the latter study, unimolecular rate constants and their falloff behavior were determined from RRKM theory, but reaction R3 was not coupled with various isomerization and decomposition pathways and the effect of chemical activation was not included. Thus, the theoretical description of the mechanism and kinetics of reaction R1 remains incomplete. Our goal in the present study is to provide the following missing ingredients: (1) an extended PES including all kinetically important channels and calculated with chemical accuracy, and (2) a comprehensive RRKM-ME analysis23-25 of the evolution of the chemically activated C8H7 (20) Wang, H.; Frenklach, M. J. Phys. Chem. 1994, 98, 11465-89. (21) Richter, H.; Mazyar, O. A.; Sumathi, R.; Green, W. H.; Howard, J. B.; Bozzelli, J. W. J. Phys. Chem. A 2001, 105, 1561-1573. (22) Moriarty, N. W.; Brown, N. J.; Frenklach, M. J. Phys. Chem. A 1999, 103, 7127. (23) Gilbert, R. G.; Smith, S. C. Theory of Unimolecular and Recombination Reactions; Blackwell Scientific Publications: Oxford, U.K., 1990. (24) Forst, W. Theory of Unimolecular Reactions; Academic Press: New York, 1973.

Reaction of Phenyl Radicals with Acetylene

radicals and the resulting product branching for reaction R1. The RRKM-ME analysis here denotes a procedure consisting of the Rice-Ramsperger-Kassel-Marcus calculation of the microcanonical rate constants [k(E)] for all elementary reactions included in the kinetic model, coupling of these k(E) values with the collisional energy transfer rates by means of a timedependent one-dimensional (E-resolved) master equation (ME), and analysis of the evolution of reactive intermediates by solving the ME for each set of experimental conditions. II. Computational Methods A detailed search of the PES for the [C8H7] molecular system was performed with the B3LYP hybrid density functional26 at the 6-311++G(d,p) level.27 Tight convergence criteria were reinforced in both geometry and electronic wave function optimizations. Harmonic vibrational frequencies calculated at the same level of theory were used for zero-point energy (ZPE) corrections, characterization of the stationary points as minima or saddle points, and RRKM calculations of the microscopic rate constants [k(E)]. The calculated and available experimental vibrational data for acetylene and phenylacetylene,28 benzocyclobutadiene,29 pentalene,30 and phenyl radical31 are listed in the Supporting Information. From the present comparisons and our earlier examinations, the B3LYP harmonic frequencies of various hydrocarbons and their radicals are on average ∼2-3% higher than the experimental fundamentals. Scaling factors of the same magnitude have been proposed.32 We have used calculated frequencies without any adjustments, mainly because of their small deviations from the available experimental data and negligible effect of frequency scaling on the calculated thermodynamic functions and kinetic parameters. To obtain chemically accurate energetic parameters, higher level calculations were carried out on the lower level optimized structures. At the higher level limit, the (R)CCSD(T)/6-311+G(3df,2p) electronic energy was approximated within a framework of the G2M composite method,33 which employs a series of single-point calculations at the (R)CCSD(T),34 (P)MP4,35 and (R)MP236 levels of theory37 with various basis sets. In some cases, we have also performed analogous spinunrestricted calculations to approximate the (U)CCSD(T)/6-311+G(3df,2p) level of theory. All post-self-consistent field (SCF) calculations employed a frozen-core (FC) approximation. The G2M(xCC5, x ) R, U) schemes were chosen as the most accurate and still feasible for eight carbon atom open-shell systems: (25) Robinson, P. J.; Holbrook, K. A. Unimolecular Reactions; Wiley: New York, 1972. (26) Becke, A. D. J. Chem. Phys. 1993, 98, 5648. (b) Becke, A. D. Phys. ReV. A 1988, 38, 3098. (c) Lee, C.; Yang, W.; Parr, R. G. Phys. ReV. B 1988, 37, 785. (d) Stephens, P. J.; Devlin, F. J.; Chabalowski, C. F.; Frisch, M. J. J. Phys. Chem. 1994, 98, 11623. (27) Pople, J. A.; Head-Gordon, M.; Raghavachari, K. J. Chem. Phys. 1987, 87, 5968. (28) King, G. W.; So, S. P. J. Mol. Spectrosc. 1970, 36, 468-487. (b) Bacon, A. R.; Hollas, J. M.; Ridley, T. Can. J. Phys. 1984, 62, 1254. (29) Chapman, O. L.; Chang, C. C.; Rosenquist, N. R. J. Am. Chem. Soc. 1976, 98, 261. (30) Bally, T.; Chai, S.; Neuenschwander, M.; Zhu, Z. J. Am. Chem. Soc. 1997, 119, 1869-1875. (31) Lapinski, A.; Spanget-Larsen, J.; Langgard, M.; Waluk, J.; Radziszewski, J. G. J. Phys. Chem. A 2001, 105, 10520-10524. (32) Scott, A. P.; Radom, L. J. Phys. Chem. 1996, 100, 16502. (33) Mebel, A. M.; Morokuma, K.; Lin, M. C. J. Chem. Phys. 1995, 103, 7414. (34) Knowles, P. J.; Hampel, C.; Werner, H. J. J. Chem. Phys. 1993, 99, 521927 and references therein. (b) Knowles, P. J.; Hampel, C.; Werner, H. J. J. Chem. Phys. 2000, 112, 3106-7; (c) Watts, J. D.; Gauss, J.; Bartlett, R. J. J. Chem. Phys. 1993, 98, 8718. (35) Chen, W.; Schlegel, H. B. J. Chem. Phys. 1994, 101, 5957-68. (36) Knowles, P. J.; Andrews, J. S.; Amos, R. D.; Handy, N. C.; Pople, J. A. Chem. Phys. Lett. 1991, 186, 130. (37) Methods with prefixes R, P, and U use the same restricted formalism for closed-shell systems, but these prefixes specify different formalisms for open shell systems: RCCSD(T) here denotes a partially spin-adapted openshell coupled cluster theory (see ref 34a,b) with perturbation correction for triples defined in ref 34c [MOLPRO keyword RHF-RCCSD(T)]; PMP4 is an approximate spin-projected MP4(SDTQ) energy after annihilation of s + 1 to s + 4 spin states; RMP2 is a spin-restricted open-shell MP2 (Gaussian keyword ROMP2).

ARTICLES

E[G2M(xCC5)] ) E[(x)CCSD(T)/6-311G(d,p)] + ∆E(+3df2p) (I-x) The present implementation of the G2M method is slightly different from the original version of Mebel et al.33 First, the 6-311G(d,p) basis set that was used originally with the B3LYP density functional for geometry optimization is now extended with diffuse functions (++). Although B3LYP geometric parameters only modestly depend on the basis set,33,38,39 in our experience, an inclusion of diffuse functions considerably improves the quality of vibrational frequencies for species with delocalized π-bonds and the accuracy of geometric parameters for transition states involving radicals. Second, the empirical higherlevel corrections (HLCs) defined by the number of R and β valence electrons are neglected. We note that individual HLCs have to be derived for each altered G2M scheme before it can be used for any nonisogyric reactions (e.g., to calculate atomization energies). However, all reactions considered in this study are isogyric (with a conserved number of electron pairs); hence, HLCs cancel in all relative energies. Third, in addition to the originally proposed (U)MP2 theory, we also used (R)MP2 calculations to evaluate the basis set extension term:

∆E(+3df2p) ) E[(y)MP2/6-311+G(3df,2p)] E[(y)MP2/6-311G(d,p)]

y ) U, R (II-y)

A less computationally demanding G2M(RCC6) version has been tested too:

E[G2M(RCC6)] ) E[(P)MP4/6-311G(d,p)] + ∆E(RCC) + ∆E(+3df,2p) (III) ∆E(RCC) ) E[(R)CCSD(T)/6-31G(d,p) - E[(P)MP4/6-31G(d,p)] (IV) Replacing (U)MP2 with (R)MP2 should help us recognize and cure possible deficiencies of (U)MP2 due to spin contamination. For similar reasons, the use of (R)MP2 in the basis set additivity approximation has been exercised previously in other model chemistries.40 The following notations will be used to distinguish between different G2M versions: G2M(xCC5,yMP2) where x, y ) R, U. Four possible models are defined by different combinations of eqs I-x and II-y. The G2M(RCC6,RMP2) model is defined by eqs III, IV, and II-R. We should comment that the G2M(UCC5,UMP2) model is similar to the G2(B3LYP/MP2/CC) model proposed by Bauschlicher and Partridge.39 For all DFT and most ab initio molecular orbital (MO) calculations the Gaussian 98 program package41 was used, with the exception of the (R)CCSD(T) calculations, which were done with MOLPRO 2000.42

III. Results and Discussion

1. Potential Energy Profile. Electrophilic addition to the CC triple bond is the only important mode of the phenyl radical (38) Baboul, A. G.; Curtiss, L. A.; Redfern, P. C.; Raghavachari, K. J. Chem. Phys. 1999, 110, 7650-7. (39) Bauschlicher, C. W.; Partridge, H. J. Chem. Phys. 1995, 103, 1788-91. (40) For example, G2MS(R); Froese, R. D. J.; Humbel, S.; Svensson, M.; Morokuma, K. J. Phys. Chem. A 1997, 101, 227-233. (41) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Zakrzewski, V. G.; 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.; Baboul, A. G.; Stefanov, B. B.; Liu, G.; Liashenko, A.; Piskorz, P.; Komaromi, I.; Gomperts, R.; Martin, R. L.; Fox, D. J.; Keith, T.; Al-Laham, M. A.; Peng, C. Y.; Nanayakkara, A.; Gonzalez, C.; Challacombe, M.; Gill, P. M. W.; Johnson, B.; Chen, W.; Wong, M. W.; Andres, J. L.; HeadGordon, M.; Replogle, E. S.; Pople, J. A. Gaussian 98, Revision A.7; Gaussian, Inc.: Pittsburgh, PA, 1998. (42) Amos, R. D.; Bernhardsson, A.; Berning, A.; Celani, P.; Cooper, D. L.; Deegan, M. J. O.; Dobbyn, A. J.; Eckert, F.; Hampel, C.; Hetzer, G.; Knowles, P. J.; Korona, T.; Lindh, R.; Lloyd, A. W.; McNicholas, S. J.; Manby, F. R.; Meyer, W.; Mura, M. E.; Nicklass, A.; Palmieri, P.; Pitzer, R.; Rauhut, G.; Schu¨tz, M.; Schumann, U.; Stoll, H.; Stone, A. J.; Tarroni, R.; Thorsteinsson, T.; Werner, H.-J. MOLPRO version 2000.1; University of Birmingham, Birmingham, U.K., 2001. J. AM. CHEM. SOC.

9

VOL. 125, NO. 37, 2003 11399

ARTICLES

Tokmakov and Lin

Figure 1. Potential energy diagram for the C6H5 + C2H2 reaction: phenylacetylene and benzocyclobutadiene branches. ZPE-corrected energies (kilocalories per mole) relative to C6H5 + C2H2 are calculated by the [G2M(RCC6,RMP2)] and G2M(RCC5,RMP2) methods.

Figure 2. Potential energy diagram for the C6H5 + C2H2 reaction: pentalene branch. ZPE-corrected energies (kilocalories per mole) relative to C6H5 + C2H2 are calculated by the [G2M(RCC6,RMP2)] and G2M(RCC5,RMP2) methods.

attack on acetylene. The H-abstraction mode is not competitive, because the C-H bond in C2H243 is ∼20 kcal/mol stronger than the C-H bond in C6H6.44 The addition of C6H5 to acetylene produces internally excited 2-phenylvinyl radicals (1), which have sufficient energy to undergo several isomerization and decomposition reactions. A complete network of feasible unimolecular transformations of 1 can be constructed by a systematic and explicit investigation of the [C8H7] PES. However, a purely combinatorial approach is prohibitively expensive due to a very large number of structural isomers for the [C8H7] molecular system and an even larger number of unimolecular transformations interconnecting these isomers. Therefore, we have to rely on chemical intuition in selecting the most probable reaction pathways to be considered in the PES search. Two basic strategies were used to ensure that all kinetically important channels are recovered. First, we investigated the most conceivable transformations of 1 and recovered those isomers that are connected to 1 via accessible barriers. The procedure was repeated for each isomer that was confirmed to be accessible. Second, we identified the most stable structures on the [C8H7] PES and tried to find the lowest energy channels to connect them to the structurally close isomers that had been confirmed to be accessible from the C6H5 + C2H2 entrance channel. In this manner, the global PES for the C6H5 + C2H2 reaction has been constructed, which can be roughly divided into three branches according to the topology of the carbon backbone: the phenylacetylene, benzocyclobutadiene, and pentalene branches. Each branch is named after a hydrocarbon that has the same carbon framework as the key intermediates of the corresponding branch. The global PES calculated at the highest levels of theory employed is schematically shown in Figures1 and 2. In the following, the elementary reactions included in each branch will be described in more detail. Unless specifically noted, the energetic parameters from Figures 1 and 2 will be used. A. Phenylacetylene Branch. The C6H5 + C2H2 reaction forms exclusively the 2-phenylvinyl radical (1) by the addition process. This compound has planar equilibrium geometry and

can have its vinylic CR-H and Cβ-H bonds either in cis (Z) or in trans (E) orientation. The 1(Z) form is more stable by only ∼1.0 kcal/mol. The lowest energy path between the E and Z forms of 1 goes through TS(E-Z) with a quasilinear CR-Cβ-H fragment. Its energy relative to 1(Z) is ∼ 4.0 kcal/mol, whereas the initially formed 2-phenylvinyl radicals are expected to have internal energy in excess of 40 kcal/mol. Hence, they will be dynamically equilibrated between the E and Z forms. In addition to the Cβ-H inversion described above, the molecular structure of 1 is flexible with respect to internal rotation about the C1CR bond. The barrier for internal rotation in 1 is ∼2.7 kcal/ mol. The corresponding torsional potential investigated at several high levels of theory is discussed separately in the Appendix. Among several possible rearrangements of 1, the most facile is the [1,4] H-shift, which can be regarded as an intramolecular H-abstraction from an ortho position of the aromatic ring by the vinylic radical center localized on the terminal Cβ atom. This channel yields 2-vinylphenyl radicals (3) via TS3 (see Figure 1). Similar to 1, the torsional motion around the C1-CR bond in 3 is hindered by ∼3.0 kcal/mol. The corresponding torsional potentials for 1 and 3 are discussed together in the Appendix. An important conclusion of the conformational analyses for 1 and 3 is that their most stable conformations are planar and they have the unpaired electron directly pointing toward the migrating H-atom; i.e., the mutual orientation of the reaction centers is favorable for the [1,4] H-shift. Consequently, the Cs molecular symmetry is preserved during this reaction, and the whole π-electronic system remains conjugated in TS3. All these factors facilitate the rearrangement between 1 and 3. Indeed, the barrier for this reaction (TS3) has a relatively small value (∼28 kcal/mol) for an H-shift involving the cleavage of a vinylic C-H bond. The unpaired electron localized on the terminal Cβ atom of 1 can attack the aromatic π-electronic cloud at ipso and ortho positions, producing intermediates 4 and 5, respectively. We have also examined the meta and para positions of the aromatic ring as possible targets of this intramolecular radical attack. However, these latter channels appear to be strongly unfavorable on the basis of their high endothermicity (>30 kcal/mol relative to the C6H5 + C2H2, according to DFT estimates). The low stability of the product molecules is easy to appreciate, because

(43) Berkowitz, J.; Ellison, G. B.; Gutman, D. J. Phys. Chem. 1994, 98, 2744. (44) Davico, G. E.; Bierbaum, V. M.; DePuy, C. H.; Ellison, G. B.; Squires, R. R. J. Am. Chem. Soc. 1995, 117, 2590-9. 11400 J. AM. CHEM. SOC.

9

VOL. 125, NO. 37, 2003

Reaction of Phenyl Radicals with Acetylene Scheme 2

they are very strained bicyclic radicals, not stabilized by conjugation. Therefore, only channels leading to 4 and 5 have been characterized at higher levels of theory and included in the mechanism of reaction R1. Spiro[2.5]octatrienyl radical (4) is destabilized by a very strained cyclopropenoid ring, which can be easily opened by breaking either one of its single C-C bonds. Thus, 4 is a shallow minimum (less than 3 kcal/mol deep) on the reaction profile for a relatively facile phenyl migration in 1. The latter process has a barrier of 28.5 kcal/mol (TS4), which is similar in magnitude to the energy of TS3 discussed above. 2-Phenylvinyl radical has to overcome a slightly higher barrier of 33.5 kcal/ mol (TS5) to produce a relatively stable intermediate 5, consisting of the fused four- and six-member rings. Both 4 and 5 are stabilized by delocalization of the π-electronic density in their C6 rings, but the unsaturated C4 ring in 5 is less strained than the C3 ring in 4. Furthermore, the π-orbitals of the C4 ring partially overlap with the C6 π-orbitals in 5, whereas the π-orbital of the C3 ring in 4 is orthogonal to the π-orbitals of the adjacent C6 ring. Thus, lower strain and higher degree of π-electron delocalization explain the higher stability of 5 compared to 4. Both isomers are easily accessible from the reactants, and as such they should be included in the mechanism of reaction R1. The [1,2] H-shift via TS2 transforms 1 into the more stable 1-phenylvinyl radical (2) where the unpaired electron is delocalized over the whole π-electronic system (Scheme 2). It has C2V symmetry with terminal CH2 moiety lying in the plane perpendicular to the aromatic ring. Because the C1-CR bond in 2 has a strong double character, internal rotation about this bond is very difficult and for this reason the harmonic oscillator model is appropriate to calculate the contribution from this motion to thermochemical functions. From both 1 and 2, the H-elimination pathways [via TS6(R) and TS6(β), respectively] lead to C6H5C2H. Although the reverse reaction clearly favors an attack by H-atoms at the terminal Cβ atom, this channel is not easily accessible from 1 due to the relatively high barrier of 45.0 kcal/mol (TS2) separating 1 and 2. An alternative pathway to isomer 2 is the [1,3] H-shift in 3 via TS7, but it involves an even higher barrier. Therefore, H-elimination via TS6(R) is expected to be the major phenylacetylene-producing branch of the C6H5 + C2H2 reaction. Nevertheless, at high T, a fraction of the chemically activated intermediates 1 will have sufficient energy to overcome TS2, which means that 1-phenylvinyl radical (2) and its consequent transformations need to be considered in the construction of the global PES. In fact, only one of the secondary reactions of 2 is kinetically important, namely, the H-elimination via TS6(β); other unimolecular transformations of 2, including the [1,2] and [1,3] H-shifts via TS2 and TS7, are not competitive due to high barriers. Remarkably, these latter transition states lie much higher than the previously discussed TS3. The possible reasons are that the [1,2] and [1,3] H-shifts in 2, unlike the [1,4] H-shift

ARTICLES

via TS3, involve very strained transition states and break the molecular symmetry. In addition, the CR-Cβ π-bond needs to be broken and then reconstructed in the orthogonal plane when going from 2 to either 1 or 3. The associated energetic expenses raise the energies of TS2 and TS7 by more than 17 kcal/mol relative to TS3. B. Benzocyclobutadiene and Pentalene Branches. In pursuit of further isomerization pathways, we have studied some plausible rearrangements of 3 and 5. The aromatic radical center in 3 can attack the CdC bond attached to the aromatic ring at the ortho position with respect to this radical center. Energetically, this intramolecular radical addition via TS8 is very close to the cyclization of 1 via TS5, where the attack occurs in the opposite direction (from the side chain onto the ring). Both reactions yield isomers (5 and 8) that have a new bicyclic framework of C atoms classified as [4.2.0]. In the following, we will separate all species having the [4.2.0] bicyclic framework and transformations between them into the benzocyclobutadiene branch or simply the [4.2.0] branch. The intramolecular cyclizations via TS5 and TS8 connect this branch to the phenylacetylene branch discussed above. Within the [4.2.0] branch, all isomers (5-8) can be obtained by adding the H-atom to benzocyclobutadiene at different positions or by sequential [1,2] H-shifts via TS9-TS11. Our calculations, however, predict that TS9-TS11 lie more than 15 kcal/mol above C6H5 + C2H2, which renders these transition states, as well as intermediates 6 and 7, inaccessible. The H-elimination channels from 5 and 8 were not included in our final kinetic model either, because of their high endothermicity. As mentioned earlier, an explicit characterization of the [C8H7] PES was beyond the scope of this study. Nevertheless, we have gone to some extra length to locate the most stable [C8H7] isomers. The global minimum on the [C8H7] PES corresponds to isomer 13, which has a [3.3.0] bicyclic structure consisting of two fused five-member rings (see Figure 2). This radical, along with 11 and 12, are products of the H atom addition to pentalene at three nonequivalent positions. All three isomers can interchange by relatively facile [1,2] H-shifts. We will refer to isomers 11-13 and their transformations as the pentalene branch or simply the [3.3.0] branch (Figure 2). Since a considerable rearrangement of the carbon skeleton is required in order to link this branch to C6H5 + C2H2, the question about its accessibility could not be answered solely by chemical intuition. Thus, we have further examined the [C8H7] PES with a focus on finding the lowest energy channel that connects isomers 11-13 to the structures already included in the reaction mechanism. The lowest energy channel is likely to be the one that involves a minimal number of rearrangements and bond-breaking reactions. Isomers 5 and 11 have the same framework of C-H and C-C bonds, except for the central C-C bond. They are the structurally closest pair of isomers that belong to different branches, one of which is the [3.3.0] branch. Therefore, we restricted our PES search to isomerization channels between 5 and 11. Two different pathways have been found (see Figure 2). Breaking the central C-C bond in either 5 or 11 gives the same eight-member ring structure 9, thus providing one obvious connection channel. The same outcome can be obtained in a concerted manner via a tricyclic intermediate 10. The C-C bonds of the central three-member ring in 10 are very weak. J. AM. CHEM. SOC.

9

VOL. 125, NO. 37, 2003 11401

ARTICLES

Tokmakov and Lin Table 2. Thermochemical Parameters of Selected Molecules and Radicals Relevant to This Study

Scheme 3

One of them is shared in 10 by the adjacent three- and fourmember rings, the other by the adjacent three- and five-member rings. Radical 11 is produced from 10 by an almost effortless breaking of the first C-C bond, whereas the other C-C bond has a dissociation barrier (TS13) of ∼11 kcal/mol that separates 10 from 5. Isomerization channels via 9 and 10 complete our global PES for reaction R1. The lowest energy path that connects the [4.2.0] and [3.3.0] branches goes over TS13, intermediate 10, and TS14. After ZPE corrections are included, the relative energy of TS14 drops slightly below that of 10, which makes TS13 the bottleneck of this channel with an effective barrier of ∼6.0 kcal/ mol relative to the C6H5 + C2H2. We did not explore any further isomerization/decomposition pathways of 11-13, because the pentalene branch did not have a significant fraction in the product distribution of reaction R1 (vide infra). The mechanism that includes the lowest energy paths of reaction R1 can be expressed by Scheme 3. It is based on a truncated version of the PES for reaction R1 featuring only seven product channels and six intermediates. We have excluded several high-energy pathways from the [4.2.0] branch. The [3.3.0] branch is not expected to be kinetically significant either; therefore, we have not performed a detailed analysis of the fate of radicals 11-13 but considered the isomerization of 5 to 11 as an irreversible channel to estimate the contribution of the [3.3.0] branch (channel R1g). Scheme 3 was used as a basis for our rate constant calculations discussed later in this article. 2. Isodesmic Enthalpies of Formation for Radicals 1-3. To assess the quality of theoretical predictions, we derived accurate enthalpies of formation for radicals 1-3 through the following isodesmic reactions:45

C6H5CHCH + C2H4 f C6H5CHCH2 + C2H3

(R4)

C6H5CCH2 + C2H4 f C6H5CH2 + H2CCCH2 (R5) C6H4CHCH2 + C6H6 f C6H5CHCH2 + C6H5 (R6) The 0 K enthalpies of reactions R4-R6 calculated at several correlated levels of theory are given in Table 3. Although these hypothetical reactions contain similar types of chemical bonds in the reactants and products, post-unrestricted Hartree-Fock (UHF) correlated methods do not take full advantage of error cancellation in the evaluation of the enthalpies of reactions R4R6. For example, (U)MP2 theory predicts an unreasonably high exothermicity for reactions R4-R6, and this deficiency is not (45) Hehre, W. J.; Radom, L.; Schleyer, P. v. R.; Pople, J. A. Ab initio Molecular Orbital Theory; John Wiley & Sons: New York, 1986. 11402 J. AM. CHEM. SOC.

9

VOL. 125, NO. 37, 2003

a Evaluated in the present work based on the B3LYP/6-311++G(d,p) molecular parameters. b Isodesmic enthalpies of formation from eqs V-VII. Enthalpies in brackets are calculated from ∆fH°0(C6H5), ∆fH°0 (C2H2), and G2M relative energies shown in Figures 1 and 2.

completely cured even at the higher levels of theory, such as (P)MP4 and (U)CCSD(T) with the 6-311G(d,p) basis set. Slow convergence of the post-UHF correlated methods is a consequence of an unbalanced spin contamination of the UHF wave functions46 for open-shell species on the right- and left-hand sides of reactions R4-R6.

Reaction of Phenyl Radicals with Acetylene

ARTICLES

Table 3. Enthalpies of Isodesmic Reactions R4-R6 Calculated by Various Methods computational methodsa

∆H°R4(0 K)

∆H°R5(0 K)

∆H°R6(0 K)

B3LYP/6-311++G(d,p) (U)MP2/6-311G(d,p) (U)MP2/6-311+G(3df,2p) (R)MP2/6-311G(d,p) (R)MP2/6-311+G(3df,2p) (P)MP4/6-31G(d,p) (P)MP4/6-311G(d,p) (U)CCSD(T)/6-311G(d,p) (R)CCSD(T)/6-31G(d,p) (R)CCSD(T)/6-311G(d,p) G2M(UCC5,UMP2) G2M(RCC6,RMP2) G2M(RCC5,RMP2)

-0.46 -27.43 -28.10 -0.68 -0.81 -6.51 -6.67 -3.30 -0.29 -0.43 -3.96 -0.58 -0.55

-4.20 -8.47 -9.16 -2.26 -2.88 -3.42 -3.42 -3.08 -2.60 -2.66 -3.77 -3.22 -3.28

0.68 -9.78 -10.27 0.35 0.18 -2.36 -2.37 -0.96 0.52 0.47 -1.45 0.34 0.31

a All reaction enthalpies include ZPE corrections calculated at the B3LYP/6-311++G(d,p) level of theory and are given in kilocalories per mole.

A better approach to evaluate the enthalpies of reactions R4R6 is to use a restricted open-shell formalism. Indeed, both (R)MP2 and (R)CCSD(T) calculations give reasonable predictions of ∆H°Rn(0 K) (n ) 4, 5, 6), which are also close to the B3LYP estimates. At the highest G2M(RCC5,RMP2) level of theory employed, we obtain ∆H°R4(0 K) ) -0.6 ( 1.0 kcal/mol, ∆H°R5(0 K) ) -3.3 ( 1.0 kcal/mol, and ∆H°R6(0 K) ) 0.3 ( 1.0 kcal/mol, where an assumed theoretical uncertainty of (1.0 kcal/mol has been included. The isodesmic enthalpies of formation for radicals 1-3 are then defined as follows:

∆fH°0(1) ) ∆fH°0(C6H5CHCH2) + ∆fH°0(C2H3) ∆fH°0(C2H4) - ∆H°R4(0 K) ) ∆fH°0(C6H5CHCH2) + D°0(C2H3-H) ∆fH°0(H) - ∆H°R4(0 K) (V) ∆fH°0(2) ) ∆fH°0(C6H5CH2) + ∆fH°0(C3H4) ∆fH°0(C2H4) - ∆H°R5(0 K) (VI) ∆fH°0(3) ) ∆fH°0(C6H5CHCH2) + ∆fH°0(C6H5) ∆fH°0(C6H6) - ∆H°R6(0 K) ) ∆fH°0(C6H5CHCH2) + D°0(C6H5-H) ∆fH°0(H) - ∆H°R6(0 K) (VII) where the theoretically evaluated enthalpies of reactions R4R6 need to be combined with the experimental enthalpies of formation of styrene,47 allene,48 ethylene,47 and benzyl radical49 from Table 2 and the D°0(C2H3-H) and D°0(C6H5-H) bond dissociation energies. The latter quantities have been reevaluated very recently by Ervin and DeTuri:50 D°0(C6H5-H) ) 111.33 ( 0.45 kcal/mol and D°0(C2H3-H) ) 109.15 ( 0.65 kcal/mol. They used the negative ion thermochemistry cycle to derive (46) Nobes, R. H.; Pople, J. A.; Radom, L.; Handy, N. C.; Knowles, P. J. Chem. Phys. Lett. 1987, 138, 481. (b) Wong, M. W.; Radom, L. J. Phys. Chem. 1995, 99, 8582. (47) NIST Chemistry WebBook; NIST Standard Reference Database Number 69; Linstrom, P. J., Mallard, W. G., Eds.; National Institute of Standards and Technology: Gaithersburg, MD, July 2001 (http://webbook.nist.gov). (48) Computational Chemistry Comparison and Benchmark Database; NIST Standard Reference Database Number 101; Johnson, R. D., III, Ed.; National Institute of Standards and Technology: Gaithersburg, MD, September 2002 (http://srdata.nist.gov/cccbdb/). (49) Ellison, G. B.; Davico, G. E.; Bierbaum, V. M.; DePuy, C. H. Int. J. Mass Spectrom. Ion Proc. 1996, 156, 109-131 and references therein. (50) Ervin, K. M.; DeTuri, V. F. J. Phys. Chem. A 2002, 106, 9947-9956.

these bond dissociation energies from the corresponding gasphase acidities that were previously determined for both ethene51 and benzene44 relative to the deprotonation enthalpy of ammonia. The reevaluated values are ∼0.5 kcal/mol lower than the originally reported ones because of the revised anchor acidity of ammonia. Substitution of all auxiliary thermochemical data into eqs V-VII leads to the isodesmic enthalpies of formation for radicals 1, 2, and 3 listed in Table 2. The rather conservative error ranges of these estimates are sums of the theoretical and experimental uncertainties of all thermodynamic parameters used in the isodesmic reaction analyses. Our value of ∆fH°298(1) is similar to that deduced by Wang and Frenklach20 from semiempirical calculations. Richter et al.21 calculated isodesmic enthalpies of formation for 1 and 3 using the BLYP density functional with the cc-pvdz basis set. Their value of ∆fH°298(3) virtually coincides with our estimate (see Table 2), but their value of ∆fH°298(1) is outside the error limits of our best estimate with a deviation of 2.4 kcal/mol. The origin of this discrepancy is difficult to track, because no details of the isodesmic reaction analysis were given by Richter et al.21 Experimental and isodesmic enthalpies of formation of various species from Table 2 have been used to calculate the benchmark reaction enthalpies, which can be compared to the values predicted by different theoretical methods. The latter are given in the Supporting Information. Written in the notation from Scheme 3, the 0 K enthalpies of various channels of reaction R1 are ∆H°R1a(0 K) ) -39.6 ( 2.8 kcal/mol, ∆H°R1b(0 K) ) -47.5 ( 2.8 kcal/mol, ∆H°R1c(0 K) ) -38.3 ( 1.6 kcal/mol, and ∆H°R1d(0 K) ) -9.7 ( 1.2 kcal/mol. The G2M values obtained by all three versions considered in this study are very consistent and reproduce the benchmark values within 1.2 kcal/mol. Taking into account the good agreement between different G2M versions, we used only the least expensive G2M(RCC6,RMP2) method to calculate the energetics of some secondary isomerization channels. The results are summarized in Figures 1 and 2. Using relative energies from these figures and experimental enthalpies of formation of C2H2 and C6H5, we estimated the enthalpies of formation for several intermediates and products. They are included in Table 2. The most recent estimates of the standard enthalpies of formation for benzocyclobutadiene52 and pentalene53 agree reasonably well with our predictions. Earlier reference values54-56 of ∆fH°298([4.2.0]C8H6) appear to be somewhat too high. The enthalpies of several reaction channels have been estimated by Richter et al.21 Their ∆H°R1a(298 K) ) -38.7 kcal/ mol, ∆H°R1c(298 K) ) -39.8 kcal/mol, ∆H°R1d(298 K) ) -9.7 kcal/mol, ∆H°R1e(298 K) ) -22.9 kcal/mol, and ∆H°R1f(298 K) ) -45.1 kcal/mol agree well with our best values ∆H°R1a(298 K) ) -40.6 ( 2.8 kcal/mol, ∆H°R1c(298 K) ) -39.3 ( 1.6 kcal/mol, ∆H°R1d(298 K) ) -9.4 ( 1.2 kcal/mol, ∆H°R1e(298 K) ∼ -25.8 kcal/mol, and ∆H°R1f(298 K) ∼ -46.9 kcal/ (51) Ervin, K. M.; Gronert, S; Barlow, S. E.; Gilles, M. K.; Harrison, A. G.; Bierbaum, V. M.; DePuy, C. H.; Lineberger, W. C.; Ellison, G. B. J. Am. Chem. Soc. 1990, 112, 5750. (52) Broadus, K. M.; Kass, S. R. J. Am. Chem. Soc. 2000, 122, 10697-10703. (53) Rogers, D. W.; McLafferty, F. J. J. Phys. Chem. A 2000, 104, 93569361. (54) Schulman, J. M.; Disch, R. L. J. Am. Chem. Soc. 1993, 115, 11153-7. (55) Rogers, D. W.; McLafferty, F. J.; Podosenin, A. V. J. Phys. Chem. A 1996, 100, 17148-17151. (56) Glukhovtsev, M. N.; Laiter, S.; Pross, A. J. Phys. Chem. A 1995, 99, 68286831. J. AM. CHEM. SOC.

9

VOL. 125, NO. 37, 2003 11403

ARTICLES

Tokmakov and Lin

Table 4. Transition-State Theory Rate Constantsa for Elementary Reactions Included in Scheme 3 reaction

log A

n

Ea/kcal mol-1

C6H5 + C2H2 f 1 1 f C6H5 + C2H2 1f2 2f1 1f3 3f1 1f4 4f1 1f5 5f1 1 f C6H5C2H + H 2 f C6H5C2H + H 3f8 8f3 5 f [3.3.0] bicycles

6.43 14.13 12.69 13.82 10.31 9.66 10.85 13.27 11.21 12.66 11.58 13.09 11.00 12.95 12.14

2.05 0.34 0.45 0.14 0.70 0.81 0.67 0.05 0.43 0.18 0.82 0.55 0.43 0.26 0.28

3.72 45.71 45.74 53.12 27.50 26.27 28.56 3.11 33.17 17.81 38.91 42.58 30.86 37.78 30.64

a

∂gi(E, t) )ω ∂t

Fitted to the modified Arrhenius form k ) ATn exp(-Ea/RT).

mol calculated from the standard enthalpies of formation of the reactants and products (Table 2). Our best estimates of the barrier for reaction R3, ER3(0 K) ∼ 27-28 kcal/mol, are close to the G2MP2 estimates of this parameter (28-29 kcal/mol) obtained by Moriarty et al.22 However, our calculations predict a similar stability for radicals 1 and 3 [∆H°R3(0 K) ) 1.2 kcal/ mol], unlike the G2MP2 method, which favors 1 by ∼6 kcal/ mol. This discrepancy is in part due to the choice of nonplanar conformations of 1 and 3 as the most stable in the study of Moriarty et al.22 Torsional potentials for 1 and 3 computed in this work at higher levels of theory disagree with that choice (see Appendix). To summarize, a comparison of our most reliable theoretical energetic parameters with the available benchmark values confirms the good quality of the former. We expect a similar quality of the G2M predictions for other energetic parameters shown in Figures 1 and 2. 3. Rate Constant Calculations. Statistical theory rate constant calculations for elementary reactions depicted in Scheme 3 were performed with the ChemRate program available from NIST.57 Molecular parameters listed in the Supporting Information were employed for the partition function, sum, and density of states computations followed by TST calculations of elementary rate constants (see Table 4) and RRKM calculations of microscopic rate constants:

k(E) )

m* N*(E - E0) m hF(E)

(VIII)

where m*/m is the path degeneracy due to optical isomers, N*(E - E0) is the sum of states of the transition state relative to the reaction threshold energy E0, and F(E) is the density of states of the active intermediate. The calculations of F(E) and N*(E - E0) were carried out with a 10 cm-1 increment via the modified Beyer-Swinehart algorithm.58 Several intermediates and transition states in Scheme 3 have among their internal degrees of freedom torsional motions hindered by small barriers. The hindered rotor treatment59 was applied to evaluate their contributions to the statistical functions. (57) Mokrushin, V.; Bedanov, V.; Tsang, W.; Zachariah, M.; Knyazev, V. ChemRate Version 1.19; National Institute of Standards and Technology: Gaithersburg, MD, 2002. (58) Astholz, D. C.; Troe, J.; Wieters, W. J. Chem. Phys. 1979, 70, 5107. (59) Knyazev, V. D. J. Phys. Chem. A 1998, 102, 3916. 11404 J. AM. CHEM. SOC.

9

VOL. 125, NO. 37, 2003

To calculate the effective bimolecular rate constants and product distributions, we have to analyze on a microcanonical level the interplay of chemical activation, isomerization, and decomposition channels for the present multiple quantum well system and also properly account for the energy transfer effects. The latter are very important because reactive intermediates produced by reaction R1, while having initial distributions that peak at energies above several decomposition and isomerization thresholds, are potentially long-lived radicals that can be stabilized in relatively deep wells. A rigorous way of predicting the kinetics of such systems is to solve the time-dependent master equation (ME), which denotes a set of coupled integrodifferential equations of motion for populations of specific energy levels of the reactive intermediates:

∫E∞ Pi(E, E′)gi(E′, t) dE′ - ωgi(E, t) 0i

ki(E)gi(E) + r(E, t) (IX) where gi(E, t) is the population of energy level E in well i at time t, ω is the collision frequency, E0i is the ground-state energy of well i, Pi(E, E′) is the transition probability for a molecule in well i with energy E′ to go on collision to another state in the same well with energy E, ki(E) is the total rate constant of decay via all isomerization and decomposition channels open from well i at energy E, and r(E, t) is the rate of formation of species i with energy E from the chemical activation and isomerization channels. In the present analysis, the single chemical activation channel provides a steady supply of reactive intermediates 1q from C6H5 and C2H2. Both reactants have Boltzmann distribution functions. Since we are only interested in the initial product branching of reaction R1, an “infinite sink” approximation is used for the R1d product channel, i.e., the concentrations of H and C6H5C2H are so small that the reverse H + C6H5C2H reaction can be neglected. In fact, this condition is satisfied even in the advanced stages of the C6H5 + C2H2 reaction in the study of Heckmann et al.,19 because H-atoms produced via R1d react more readily with phenyl to form benzene rather than adding to C6H5C2H. We assume that all energy transfer acts are induced by weak molecular collisions and the energy transfer probabilities are given by the standard “exponential-down” model23 with an empirical value of 400 cm-1 for 〈∆E〉down (average energy loss per collision of the active component with a bath gas molecule). The sensitivity of the calculated rate constants to this parameter has been examined. The frequency of collisions was derived from the Lennard-Jones (L-J) parameters of Ar [σ(Ar) ) 3.54 Å, ε/kB(Ar) ) 93.3 K]60 and C8H7 [σ(C8H7) ) 5.70 Å, ε/kB(C8H7) ) 546 K]. The latter values are obtained from an empirical relationship between the L-J parameters and molecular weight established for a series of aromatic hydrocarbons,61 and they are very similar to the L-J parameters for phenylacetylene and styrene estimated from their boiling points.61 The ME was solved in a matrix form (for an array of discrete states Ej, each with width δE, and energy-dependent functions represented by vectors) with a method based on the Householder’s tridiagonalization algorithm.62 The energy bin size δE ) 100 cm-1 was (60) Reid, R. C.; Prausnitz, J. M.; Sherwood, T. K. The Properties of Gases and Liquids, 3rd ed.; McGraw-Hill: New York, 1977. (61) Wang, H.; Frenklach, M. Combust. Flame 1994, 96, 163. (62) Wilkinson, J. H.; Reinsch, C. Linear Algebra; Springer: New York, 1971.

Reaction of Phenyl Radicals with Acetylene

Figure 3. Experimental and calculated total rate constants for the C6H5 + C2H2 reaction. The recommended effective total rate constant (solid curve) is based on the fitted value of E1(0 K) ) 4.1 kcal/mol and the G2M value of E2(0 K) ) -2.2 kcal/mol. Dashed line represents the HPL of kR1, not corrected for the decomposition of 1 back to the reactants. Dotted curves in plot A show the sensitivity of kR1 to a variation of E1 from 3.7 kcal/mol (G2M value) to 4.4 kcal/mol (B3LYP value). Dotted curves in plot B show the kR1(0.1 Torr) and kR1(1 atm) calculated from the fitted value of E1(0 K) ) 4.1 kcal/mol and the B3LYP value of E2(0 K) ) 1.2 kcal/mol. Experimental data: (9) from ref 18.

used in the ME computations; the matrix size was up to 3840 × 3840 to ensure the convergence at high T. More details about the implementation of the time-dependent weak collision ME/ RRKM analysis in ChemRate are available in a series of publications by Tsang and co-workers.63 We used theoretical reaction barriers and enthalpies computed by the G2M(RCC5,RMP2) method (Figures 1 and 2) with only one minimal adjustment: E1(0 K), the 0 K entrance barrier (TS1), was raised from 3.7 to 4.1 kcal/mol. This minor empirical correction allowed us to quantitatively (within a factor of 2) account for the available experimental kinetic data illustrated in Figures3 and 4. The effective total rate constant for reaction R1 is weakly dependent on pressure and can be expressed as kR1 ) (1.29 × 1010)T0.834 exp(-2320/T) cm3 mol-1 s-1. Its effective nature is due to fact that a fraction of reactive intermediates may decompose back to the reactants. This process becomes increasingly important at high T, causing a deviation of the effective rate constant from the TST predictions (see Figure 3A). (63) Bedanov, V. M.; Tsang, W.; Zachariah, M. R. J. Phys. Chem. 1995, 99, 11452. (b) Tsang, W.; Bedanov, V.; Zachariah, M. R. J. Phys. Chem. 1996, 100, 4011. (c) Knyazev, V. D.; Tsang, W. J. Phys. Chem. A 2000, 104, 10747-10765. (d) Knyazev, V. D.; Tsang, W. J. Phys. Chem. A 1999, 103, 3944-3954.

ARTICLES

Figure 4. Branching rate constants for reaction R1. Plot A: (s) kR1d(1 atm); (‚‚‚) kR1d(0.1 Torr) and kR1d(100 atm). Experimental data (see Table 1): (O) from ref 19; (×) reevaluated kR1 from ref 16. Plot B: (s) kR1d(1 atm) and kR1s(1 atm) ) kR1a + kR1b + kR1c + kR1e + kR1f. Fitted rate constants (in cubic centimeters per mole per second): kR1s(1 atm) ) (4.80 × 1044)T-9.9 exp(-8870/T), kR1g(1 atm) ) (4.85 × 1030)T-5.5 exp(-13 630/T), kR1d(0.1 Torr) ) (2.93 × 1012)T0.18 exp(-3280/T), kR1d(1 atm) ) (2.66 × 1032)T-5.3 exp(-11 970/T), kR1d(100 atm) ) (5.10 × 1044)T-8.5 exp(-19 690/T). Previous theoretical estimates of kR1d at P ) 1 atm are taken from ref 18 ([), ref 20 (9), and ref 21 (4).

In addition to the total rate constant, the present analysis provides an insight into the mechanism and product branching of reaction R1 under different simulated experimental conditions. For practical applications, the most important kinetic parameter is the branching fraction of the H-elimination channel R1d. The T-dependence of the kR1d branching rate constant at different pressures is illustrated in Figure 4A. Provided that P is high enough for rapid collisional energy transfer and the T is low enough to prevent reactivation, reactive intermediates can be trapped in one of the relatively deep wells. For this reason, kR1d strongly depends on P at low T, and channel R1d contributes negligibly to the product distribution of reaction R1 at high P and low T. As shown in Figure 4A, kR1d becomes essentially independent of P at T > 1000 K, which marks the threshold of thermal stability of the reactive intermediates. At high T, even though reactive intermediates may suffer numerous collisions with bath gas, they are not stabilized because a large fraction of these collisions are activating. As a result, reaction R1 proceeds in a steady-state regime, where the newly formed C8H7 adducts quickly decompose either to C6H5C2H + H or back to the reactants, and kR1 ∼ kR1d, independent of P. The major fraction of C6H5C2H + H is produced directly from 1; J. AM. CHEM. SOC.

9

VOL. 125, NO. 37, 2003 11405

ARTICLES

Tokmakov and Lin

the remainder (up to 13% at 2000 K) is a contribution of the decomposition flux through TS6(β). Our calculated kR1d rate constant at P ) 1 atm is compared to the values recommended in the previous studies in Figure 4B. We notice that our predicted values of kR1d at low T are systematically higher than all previous recommendations. This is in part due to different energetics, e.g., a lower relative energy of TS6(R), E2(0 K), employed in our model and in part due to different methodological choices. Unlike the simplified steadystate analysis, the weak collision ME/RRKM modeling carried out in the present study allows us to accurately take into account the reversible nature of the isomerizations between different reactive intermediates and their thermal decomposition. The difference between the kR1 and kR1d rate constants is almost entirely due to the production of stabilized C8H7 isomers at low T (see Figure 4B). Their composition is determined by a competition of the isomerization pathways with decomposition and stabilization. Common trends can be derived from Figure 5, where the composition of the stabilized radicals is shown as a function of reaction time. At high P and low T, the stabilized products are composed almost exclusively of isomer 1, because collisional stabilization of the chemically activated 1q occurs faster than isomerization or decomposition. However, as T increases, the product distribution becomes more diverse. At the same time, the individual branching rate constants become time-dependent, and for this reason they are not always welldefined. At low P and low T, reactive intermediates have a better chance to isomerize before being trapped in one of the wells. Three wells (isomers 1, 3, and 8) effectively compete for reactive intermediates under these conditions. Phenyl migration via intermediate 4 takes place in the same regimes as the isomerization reactions discussed above. This channel does not produce any new chemical species and therefore does not bear any kinetic consequences for reaction R1. However, this process is ubiquitous to all molecules containing a radical site in the β-position to an aryl group. Hence, phenyl migration in 1 is an important prototype reaction. Our calculations indicate that pentalene branch R1g is essentially inaccessible from C6H5 + C2H2. The fraction of reactive intermediates reaching the [3.3.0] bicyclic structures increases with T (Figure 4B) but remains negligible even at 2000 K. To complete the study, we tested the sensitivity of our predicted total and branching rate constants to small variations in the most critical energetic parameters. The effect of replacing the fitted E1(0 K) ) 4.1 kcal/mol with either the G2M(RCC5,RMP2) value of 3.7 kcal/mol or the B3LYP/6-311++G(d,p) barrier of 4.4 kcal/mol is illustrated in Figure 3A. The effective total rate constant kR1 changes by a factor of 2 at room temperature and becomes less and less sensitive to the magnitude of the entrance barrier as the T increases. The kR1 values at high T are sensitive to the height of the H-elimination barrier. We have designated the 0 K energy of TS6(R) relative to the C6H5 + C2H2 reactants as E2(0 K). The theoretical estimates of this parameter strongly depend on the level of theory employed. Although we believe the G2M estimates to be the most accurate, we have tested the sensitivity of kR1 with respect to raising E2(0 K) to 1.2 kcal/mol (B3LYP value). The effect on kR1 is 2-fold (see Figure 3B): (1) kR1 drops by a factor of 2 at 1000 ( 100 K and to a lesser extent at other 11406 J. AM. CHEM. SOC.

9

VOL. 125, NO. 37, 2003

T, and (2) kR1 becomes more P-dependent in the T ) 5001000 K range. This behavior can be anticipated, since a higher E2 barrier reduces the fraction of the chemically activated intermediates that undergo H-elimination, while increasing the alternative channels’ yields (mainly decomposition back to the reactants, but also isomerizations and deactivation). Finally, we have tested our results with respect to a variation of the assumed value of 〈∆E〉down ) 400 cm-1 by (100 cm-1. Since the effective total rate constant is essentially independent of P under a wide range of experimental conditions, it is essentially unaffected (less than 5%) by the variations in 〈∆E〉down. The branching rate constants, on the other hand, are more sensitive, especially in the low-T region where they are more P-dependent. For example, kR1d(1 atm) becomes 30% faster at 700 K and 50% faster at 500 K if 〈∆E〉down ) 300 cm-1; the same rate constant becomes 20% slower at 700 K and 25% slower at 500 K if 〈∆E〉down ) 500 cm-1. Unfortunately, no experimental P-dependent kinetic data are available for reaction R1; therefore, it remains unclear which value of 〈∆E〉down is more appropriate. IV. Conclusion

The mechanism of the phenyl radical reaction with C2H2 has been investigated quantum-chemically at several theoretical levels. The energetic parameters predicted at the highest theoretical level [G2M(RCC5)] agree with available benchmark values within 1.2 kcal/mol, which is well within the desired chemical accuracy limits. The mechanism of reaction R1 and thermochemistry of individual species established in this study considerably improve and extend the results of previous theoretical investigations. The rate constants for the C6H5 + C2H2 reaction have been derived from the ME analysis of the comprehensive kinetic model including the most important product branches. The available experimental kinetic data can be reproduced within its scatter, if 3.7 < E1 < 4.4 kcal/mol and -2.2 < E2 < 1.2 kcal/mol, i.e., within the range of the G2M and B3LYP predictions of the relative energies of TS1 and TS6(R), respectively. Under combustion conditions (T > 1000 K), the exclusive products of reaction R1 are phenylacetylene and H-atoms. At low T, the reactive intermediates can be deactivated primarily in well 1 but also in wells 3 and 8. Stabilized isomeric C8H7 radicals can serve as active agents in the mass growth reactions with C2H2, C2H4, and other light unsaturated hydrocarbons and radicals, thus contributing to the PAH formation. The most important isomerizations are 1 f 3, 3 f 8, and C6H5 migration via a short-lived intermediate 4. The pentalene branch containing the most stable C8H7 isomers remains kinetically unimportant virtually under any experimental conditions. Nevertheless, our extensive PES lays the groundwork for detailed studies of the mechanism and product distribution for other reactions (such as H + pentalene, C5H5 + C3H2, etc.) that are more likely to access intermediates 11-13. The PES for reaction R1 is also a good starting point for a study of the mechanism and product distribution of the C6H5C2H + H reaction; however, it needs to be appended with three H-abstraction and four H-addition channels, which will be deferred to a later publication. Acknowledgment. We are grateful for the support of this work from the Division of Chemical Sciences, Office of Basic

Reaction of Phenyl Radicals with Acetylene

ARTICLES

Figure 5. Time-dependent composition of the intermediates produced by the C6H5 + C2H2 reaction at selected P and T.

Energy Sciences, Department of Energy, through Contract DEFGO2-97ER14784. Also, we are thankful to the Cherry L. Emerson Center of Emory University for the use of its resources, which is in part supported by a National Science Foundation grant (CHE-0079627) and an IBM Shared University Research Award.

Appendix: Conformational Analysis of 2-Phenylvinyl (1) and 2-Vinylphenyl (3) Radicals

To establish reliable thermochemical parameters for intermediates 1 and 3, their internal rotational profiles have been explored in detail. As a benchmark, we have also studied internal J. AM. CHEM. SOC.

9

VOL. 125, NO. 37, 2003 11407

ARTICLES

Tokmakov and Lin

data. This potential features a hindering barrier for internal rotation of 3.0 kcal/mol, which separates two quasi-planar minima (barrier to planarity is less than 0.01 kcal/mol). The B3LYP and MP2 potentials are qualitatively similar to the benchmark CCSD(T) potential, the main deviations being a slightly higher hindering barrier at the B3LYP/6-311++G(d,p) level and a strongly twisted geometry of the lowest energy conformation at the MP2/6-31G(d,p) level. The disagreement between MP2 and other computational methods becomes much larger when they are applied to the open-shell species 1 and 3. As shown in Figure 6, B3LYP, PMP4, and (R)CCSD(T) predictions are qualitatively similar to each other and to the benchmark potential of styrene, in the sense that global maxima are located at τ ∼ 90° and global minima are at τ ) 0°. Unexpectedly, UMP2 theory predicts the global minima at τ ) 90°.65 UMP2 calculations fail to produce physically meaningful torsional profiles for 1 and 3, because of the adverse effect of spin contamination, which is higher for the UHF wave functions of planar conformations of 1 and 3 and minimal at τ ) 90°. Our best estimates of the hindering barriers for internal rotation in 1 and 3 were obtained from the (R)CCSD(T)/6-311G(d,p) calculations (shown as crosses in Figure 6). After ZPE corrections are included, V0(1) ) 2.7 kcal/mol and V0(3) ) 3.0 kcal/mol. For both radicals, their torsional potentials have two minima (τ ∼ 0° and τ ∼ 180°), which are equivalent for 1 but different for 3. The τ ∼ 180° conformation of 3 is ∼0.7 kcal/ mol less stable than the τ ∼ 0° conformation because of the repulsive interactions of the vinyl group with the o-hydrogen of the aromatic ring. These interactions are maximized when vinyl is in cis orientation to the o-hydrogen (τ ) 180°). The (P)MP4 and (R)CCSD(T) single-point calculations indicate that the second minimum of 3 may have a slightly twisted geometry. However, the barrier to planarity appears to be very small, as in the case of styrene. Figure 6. Molecular structures (τ ) 0 conformations) and internal rotational profiles for cis-C6H5CHCH (1), o-C6H4C2H3 (3), and C6H5CHCH2 calculated at different levels of theory (ZPE correction is not included): (9) B3LYP/ 6-311++G(d,p); (b) UMP2/6-31G(d,p); (0) UMP2/6-311G(d,p); (]) UMP2/6-311+G(3df,2p); (4) PMP4/6-31G(d,p); (3) PMP4/6-311G(d,p); (×) RCCSD(T)/6-311G(d,p). Geometries were optimized at the B3LYP/ 6-311++G(d,p) and UMP2/6-31G(d,p) levels; all other methods employed the B3LYP geometries. Dashed line represents the CCSD(T)/extrap.//CCD/ cc-pVDZ profile obtained in ref 64.

rotation in styrene, which is structurally very similar to 1 and 3. The results are shown in Figure 6. The molecular structure of styrene has been reviewed recently by Sancho-Garcia and Perez-Jimenez.64 Their final form of the classical torsional potential (dashed line in Figure 6) was calculated at the CCSD(T) level extrapolated to the complete basis set; it was also corroborated by the available spectroscopic (64) Sancho-Garcia, J. C.; Perez-Jimenez, A. J. J. Phys. B: At. Mol. Opt. Phys. 2002, 35, 1509-1523 and references therein.

11408 J. AM. CHEM. SOC.

9

VOL. 125, NO. 37, 2003

Supporting Information Available: Tables S1-S5 contain the molecular parameters for all species and transition states calculated in this study; Tables S6-S10 contain detailed energetics of all stationary points at various theoretical levels (print). This material is available free of charge via the Internet at http://pubs.acs.org. JA0301121 (65) In all post-UHF calculations we used a quadratically convergent SCF procedure (Bacskay, G. B. Chem. Phys. 1981, 61, 385-404) with tight convergence criteria [e.g., UMP2(SCF ) qc, tight)]. The structures of 1 and 3 optimized by the default UMP2(SCF ) DIIS) algorithm have been described in detail by Moriarty et al.22 They argued that the most stable conformations of 1 and 3 are located at some intermediate values of 0 < τ < 90°, but they did not calculate complete torsional potentials. We believe that neither UMP2(SCF ) qc, tight) nor UMP2(SCF ) DIIS) can afford reliable torsional profiles for 1 and 3, because both methods are adversely affected by spin contamination. Furthermore, given the fact that the optimized geometry of styrene is essentially indistinguishable from the planar conformation at higher levels of theory, there is no reason to expect a large deviation from planarity for the global minima of either 1 or 3, where the energy of π-conjugation is very similar to that in styrene, but the steric repulsions are weakened by removing one H-atom.