Kinetics of the Hydrogen Abstraction PAH + ∙OH → PAH Radical +

4 days ago - A Reaction Class Transition State Theory (RC-TST) augmented with Structure-Activity Rela-tionship (SAR) methodology is applied to predict...
0 downloads 0 Views 2MB Size
Subscriber access provided by McMaster University Library

A: Kinetics, Dynamics, Photochemistry, and Excited States

Kinetics of the Hydrogen Abstraction PAH + ·OH # PAH Radical + H2O Reaction Class: An Application of the Reaction Class Transition State The-ory (RC-TST) and Structure-Activity Relationship (SAR) Maciej Baradyn, and Artur Ratkiewicz J. Phys. Chem. A, Just Accepted Manuscript • DOI: 10.1021/acs.jpca.8b10988 • Publication Date (Web): 31 Dec 2018 Downloaded from http://pubs.acs.org on January 1, 2019

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 43 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Kinetics of the Hydrogen Abstraction PAH + ∙OH → PAH Radical + H2O Reaction Class: An Application of the Reaction Class Transition State Theory (RC-TST) and Structure-Activity Relationship (SAR)

Maciej Baradyn, Artur Ratkiewicz* Institute of Chemistry, University of Bialystok, ul. Ciolkowskiego 1K 15-245 Bialystok, Poland

*

Corresponding author. Email address: [email protected]

Tel: (+48) 857388250 Fax: (+48) 857388258

1 ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 2 of 43

Abstract A Reaction Class Transition State Theory (RC-TST) augmented with Structure-Activity Relationship (SAR) methodology is applied to predict high-pressure limit thermal rate constants for hydrogen abstraction by ∙OH radical from Polycyclic Aromatic Hydrocarbons (PAHs) reaction class in the temperature range of 300-3000K. The rate constants for the reference reaction of C6H6 + ∙OH → C6H5 + H2O, is calculated by the Canonical Variational Transition State Theory (CVT) with the Small Curvature Tunneling (SCT). Only the reaction energy is needed to predict RC-TST rates for other processes within the family, the parameters needed were obtained from M06-2X/cc-pVTZ data for a training set of 34 reactions. The systematic error of the resulting RC-TST rates is smaller than 50% in comparison with explicit rate calculations, which facilitates application of the proposed methodology to the Automated Reaction Mechanism Generators (ARMGs) schemes.

2 ACS Paragon Plus Environment

Page 3 of 43 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

1. Introduction It is well known, that Polycyclic Aromatic Hydrocarbons (PAH), mainly formed in combustion of fossil fuels, are ubiquitous in the atmosphere and significantly contribute to the formation of combustion generated particles, like soot1-4. They are also utilized in the petrochemical industry as important intermediates in coal conversion processes and byproducts in steam cracking units. Initializing the combustion process, the hydrogen atom metatheses by atoms (i.e., ∙O, ∙H,)4-10 and small radicals (i.e., ∙OH9-15, ∙CH39, 16-18, ∙HO29) are usually the most important channels for the fuel consumption, hence the complex models are highly sensitive to their rate constants12. Moreover, ∙OH radicals and multi-ring PAH compounds are among the key players in the competition between pyrolysis and oxidation processes13, thus the availability of accurate kinetic parameters is an essential requirement for reliable modeling of PAH growth and for proper capturing of the critical features of these systems in extended conditions. However, the detailed oxidation and pyrolysis kinetics of aromatics, reflecting the rich and intricate chemistry of the conjugated rings, are less well understood than those for aliphatic hydrocarbons. The kinetic data for such systems are generally sparse and only available for a narrow range of temperature, thus the thorough knowledge of appropriate rates is required to avoid even relatively small uncertainties in macroscopic simulations. There have been considerable amount of studies related to the title reaction class6, 10-13, 15, 19-28. The majority of these were concerned with the benzene + ∙OH system, due to its relative simplicity and practical and fundamental importance. The achievements in the kinetics of the ∙OH radical reaction with benzene acquired till 80’s, have been reviewed and evaluated by Atkinson29, 30. More recently, Tokmakov and Lin27 reported the electronic structure and kinetics calculations for possible channels of the benzene + ∙OH reaction. Their analysis was further extended by Chen et al.13, where reaction paths and elementary steps were presented with 3 ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 4 of 43

kinetic parameters evaluated with Quantum Rice-Ramsperger-Kassel (QRRK) theory31, 32 and master equation analysis. It is interesting to compare above calculations

with shock tube

experiments augmented with TST calculations of Seta et al.11, where rate constants for the reactions of ∙OH radicals with benzene and toluene were evaluated. Recently, Semenikhin et al.12 carried out a detailed analysis of the reactions of selected PAHs with radicals, including ∙OH. They utilized a Classical TST (CTST) theory and distinguished two reaction sites for PAH + ∙OH reaction, namely zigzag and armchair. Also, computational and experimental rate constants were reported for naphthalene23, 24, anthracene26 and phenanthrene14, 33. The reactions of ∙OH radicals with PAHs: PAH + ∙OH → ∙PAH + H2O (∙H-abstraction channel) PAH + ∙OH → PAH–OH

(∙OH addition channel)

are mutually competitive. The favored pathway depends mainly on temperature. In the low-T region, the ∙OH-addition channel, followed by the H atom abstraction from C6H6OH in secondary processes, dominates. However, for T > 500 K, reactant decay may be attributed nearly entirely to the abstraction channel11,

13, 23, 26, 27, 29

, thus this pathway is of crucial importance for proper

understanding of combustion of benzene. For the other members of the title reaction family it is known that, similarly like for benzene, the abstraction channel dominates under combustion conditions11, 23, 26, 34-36. However, only limited kinetic data are available only for the simplest members of the title reaction class. There are no experimental works concerned with PAHs consisting of more than three rings, which contribute significantly to the soot formation and the regeneration processes occurring in DPF (Diesel Particulate Filters) in modern cars. This raises a challenge of obtaining accurate kinetic parameters for the complex schemes. To meet this challenge, an automated method capable of generating reliable thermal rates of elementary steps 4 ACS Paragon Plus Environment

Page 5 of 43 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

is needed37. A strict theoretical framework, relaying on an extrapolation of accurate rate constants of a chosen (reference) reaction to the whole family, is provided by the Group Additivity (GA) for transition states38, 39 and Reaction Class Transition State Theory (RC-TST)40 methodologies. In this work, the RC-TST is applied to estimate the thermal rate constants of any processes belonging to the PAH + ∙OH → PAH radical + H2O reaction class. The basic assumption, extensively validated previously40, is that the averaged expressions, referred to as the RC-TST correlations, can be extended to all reactions within the family. These correlations are analytically derived from explicit DFT calculations for a set of selected processes within a reaction class, called further the Representative Training Set. After successful application of the RC-TST to investigate H abstractions from PAHs by H atoms14 and alkyl radicals (∙CH3 and ∙C2H5)18 it is expected that this methodology can be also effectively utilized for the title reaction family. To compose the Representative Training Set 34 reactions were selected, as shown in Table 1. The different reaction sites were distinguished based on their closest vicinity. The classification of the C6 abstraction sites, based on the molecular topology, is pictured and explained in the Figure S101 in the Supporting Info file. Specifically, the sites at the six-membered rings, called further the C6 sites, were divided into five types, two among which (namely α and β ) are considered as a main type and referred as “typical” further in this study. The other kinds of sites, namely zigzag ( armchair ( ) and mixed zigzag-armchair (

),

), are referred to as “atypical” or “others” hereafter,

see Figure 1 below and Figure S101 in the Supporting Info. In contrast to distinctions among C6 sites, all the H abstraction sites of 5-membered rings were found to be similar one to another, thus they were assumed as belonging to the same type, denoted as “C5”.

5 ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 6 of 43

b site

1

1

b

CH4

2

3

2

a

a site

Anthracene C10H14 - C14Anth

Naphthalene C10H8 - C10

Benzene C6H6

(a) Linear PAHs (acenes) 4 1 5

2

3

a

b

a

b

4

3

3

b

b

6

4

b

b Acenaphtalene C12H8 - C12

2

10 9

8

4 5

a

b

b

4

b

b

5

7

6

a

a

3

2

a

b

Pyrene C16H10 - C16Py

b

a

a

10

1

1

11

b

10

a

a

a

a

8

b

7

2

9

CH4

3

4 6

8

7

a

9

a

a

Tetraphene C18H12- C18Tp

Aceanthrylene C16H10 - C16Anth

8

b

a

2

12

9 6

a

Acephenanthrylene C16H10 - C16Phen

1

3

7

1

10

Phenantrene C14H10 - C14Phen

3

1

2

1

a

2

b

a

5

5

a

Cyclopenta[cd]pyrene C18H10 - C18

b 2

a

1

5

3

4

a

12

b

3

2

11

4

a

5

a

10

b

6

1

a

9

8

a

Dicyclopenta[cd,jk]pyrene C20H10 - C20

a

7

a

Cyclopenta[fg]benzo[ghi]perylene C24H12 - C24

(b) Nonlinear PAHs

Figure 1: List of species considered in this study with IUPAC nomenclature and short notation. The sites are labeled on the basis of the hydrogen abstraction sites classification (see text and Figure S101 in the Supporting Info).

6 ACS Paragon Plus Environment

Page 7 of 43 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Table 1: List of reactions forming the Representative Training Set, examined in this study. Notation of PAH reactants and products are explained in Figure 1.

No R1 (β) R2 (α) R3 (β) R4 (β) R5 (α) R6 (□) R7 (α) R8 (β) R9 (Δ) R10 (α) R11 (β) R12 (□) R13 (□) R14 (β) R15 (C5) R16 (C5) R17 (β) R18 (β) R19 (α) R20 (α) R21 (□) R22 (β) R23 (β) R24 (α) R25 (β) R26 (β) R27 (C5) R28 (C5) R29 (β) R30 (α) R31 (α) R32 (α)

Reaction ∙



C6H6 + OH → C6H5 + H2O C10H8 + OH ∙ → 2*C10H7∙ + H2O C12H8 + OH ∙ → 3*C12H7∙ + H2O C14-anthracene + OH ∙ → 1*C14-anthracene∙ + H2O C14-anthracene + OH ∙ → 2*C14-anthracene∙ + H2O C14-phenanthrene + OH ∙ → 1*C14 -phenanthrene∙ + H2O C14-phenanthrene + OH ∙ → 4*C14 -phenanthrene∙ + H2O C16-aceanthrylene + OH ∙ → 4*C16-aceanthrylene∙ + H2O C16-aceanthrylene + OH ∙ → 6*C16-aceanthrylene∙ + H2O C16-aceanthrylene + OH ∙ → 7*C16-aceanthrylene∙ + H2O C16-aceanthrylene + OH ∙ → 9*C16-aceanthrylene∙ + H2O C16-aceanthrylene + OH ∙ → 10*C16-aceanthrylene∙ + H2O C16-acephenanthrylene + OH ∙ → 1*C16-acephenanthrylene∙ + H2O C16-acephenanthrylene + OH ∙ → 2*C16-acephenanthrylene∙ + H2O C16-acephenanthrylene + OH ∙ → 4*C16-acephenanthrylene∙ + H2O C16-acephenanthrylene + OH ∙ → 5*C16-acephenanthrylene∙ + H2O C16-acephenanthrylene + OH ∙ → 8*C16-acephenanthrylene∙ + H2O C16-acephenanthrylene + OH ∙ → 9*C16-acephenanthrylene∙ + H2O C16-pyrene + OH ∙ → 2*C16-pyrene∙ + H2O C18H10 + OH ∙ → 1*C18H9 ∙ + H2O C18Tp + OH ∙ → 1*C18Tp + H2O C18Tp + OH ∙ → 4*C18Tp + H2O C18Tp + OH ∙ → 5*C18Tp + H2O C18Tp + OH ∙ → 6*C18Tp + H2O C18Tp + OH ∙ → 11*C18Tp + H2O C18Tp + OH ∙ → 12*C18Tp + H2O C20H10 + OH ∙ → 3*C20H9 + H2O C20H10 + OH ∙ → 4*C20H9 + H2O C24H12 + OH ∙ → 2*C24H9 + H2O C24H12 + OH ∙ → 3*C24H11 + H2O C24H12 + OH ∙ → 5*C24H11 + H2O C24H12+ OH ∙ → 6*C24H11 + H2O 7 ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

R33 (α) R34 (C5)

Page 8 of 43

C24H12+ OH ∙ → 8*C24H11 + H2O C24H12+ OH ∙ → 10*C24H11 + H2O

R1, the simplest process within the title family, is assumed to be the reference reaction. Of the 34 processes, eleven are H abstractions from , thirteen from  and five from “other” C6 positions. Five remaining H abstractions take place at the C5 sites.

2. Methodology 2.1. Reaction Class Transition State Theory (RC-TST) Since the principles of the RC-TST methodology were exhaustively described in the previous reports40, only main features, necessary to comprehend this study, are presented here. The assumption lies in the similarity of the reactive moieties of all reactions within a given family. The differences between their rate constants are considered as mainly caused by diverse interactions of the reactive center and substituents. The rate constant 𝑘(𝑇) of any reaction from a given class is related to the rate constant of a reference reaction kref(T) by a temperature-dependent coefficient f(T): 𝑘(𝑇) = 𝑓(𝑇) × 𝑘𝑟𝑒𝑓 (𝑇)

(1)

The rate constant of the reference reaction kref(T) may be taken directly from experiment or calculated at the high level of accuracy. Under the RC-TST framework, the f(T) is factorized into components, capturing differences (between reference process and reaction investigated) of quantities affecting the TST rate constants. 𝑓(𝑇) = 𝑓𝜎 × 𝑓𝑉 (𝑇) × 𝑓𝜅 (𝑇) × 𝑓𝑄 (𝑇) × 𝑓𝐻𝑅 (𝑇)

(2)

8 ACS Paragon Plus Environment

Page 9 of 43 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

where 𝑓𝜎 , 𝑓𝑉 (𝑇), 𝑓𝜅 (𝑇), 𝑓𝑄 (𝑇), and 𝑓𝐻𝑅 (𝑇) are the coefficients (factors) corresponding to: symmetry number, potential energy, tunneling, partition function and hindered rotations, respectively. Mathematically, all of the above factors are the ratios of the particular components of the TST41 formula for the arbitrary and reference reactions:

𝑓𝜎 =

𝜎

(3)

𝜎𝑟𝑒𝑓

fκ (𝑇) =

κ(T)

(4)

κref (𝑇)

𝑄‡ (𝑇) 𝑄‡ (𝑇) ( ‡ ) ( 𝑅 ) 𝑄𝑟𝑒𝑓 (𝑇) Φ (𝑇) fQ (𝑇) = = ‡ ‡ (𝑇) (𝑇) 𝑄𝑟𝑒𝑓 Φ𝑟𝑒𝑓 ( ‡ ) ( 𝑅 ) Φ (𝑇) Φ𝑟𝑒𝑓 (𝑇) (5)

𝑓𝑉 (𝑇) = 𝑒𝑥𝑝 (−

‡ Δ𝑉 ‡ − Δ𝑉𝑟𝑒𝑓

𝑘𝐵 𝑇

ΔΔ𝑉 ‡ ) = 𝑒𝑥𝑝 (− ) 𝑘𝐵 𝑇 (6)

‡ 𝑄𝐻𝑅 (𝑇) ‡ 𝑄𝐻𝑂 (𝑇)

‡ ‡ 𝑄𝐻𝑅 (𝑇) × 𝑄𝐻𝑂,𝑟𝑒𝑓 (𝑇) 𝑐(𝑇) 𝑓𝐻𝑅 (𝑇) = = ‡ = ‡ ‡ 𝑐𝑟𝑒𝑓 (𝑇) 𝑄𝐻𝑅,𝑟𝑒𝑓 (𝑇) 𝑄𝐻𝑂 (𝑇) × 𝑄𝐻𝑂,𝑟𝑒𝑓 (𝑇) ‡ 𝑄𝐻𝑂,𝑟𝑒𝑓 (𝑇)

(7) In the formulas above, σ is the degeneracy of the reaction path (reaction symmetry number), 𝛷R(T) and Q‡(T) symbolize the total partition function of the reactants and transition state, ‡

respectively, (T) is the transmission coefficient gauging the quantum tunneling, V is the 9 ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 10 of 43

classical reaction barrier height, c(T ) is the hindered internal rotation coefficient, equal to the ratio of the total partition function with explicit hindered rotations treatment QHR(T) to those calculated within the harmonic oscillator (HO) approximation QHO(T); T is the temperature in Kelvins and k B is the Boltzmann constant. In the RC-TST/SAR methodology the potential energy ‡

factor is calculated from Eq. (6), where the classical reaction barrier height V for the arbitrary reaction is obtained from the Linear Energy Relationship (LER), which utilizes structure−activity relationships (SAR) to approximate the above reaction class factors. The other factors, except for the exactly calculated symmetry factor, are averaged and then approximated with simple expressions, called further RC-TST correlations. The rate constants of any reactions in the title class can be predicted either from the SAR approach with the reaction energy and symmetry number or BHG (Barrier Height Grouping), where only the symmetry number is needed.

2.2. Computational Details All the electronic structure calculations were performed with the GAUSSIAN 09 program42. In the RC-TST, the relative reaction energies and barrier heights, which may be fast and relatively easily obtained on-the-fly, are more important than values calculated directly, even at the high theory level. It was shown previously40, that for some reaction classes, the reaction energies calculated even at as low as AM1 theory level, give satisfactory accuracy. For the vast majority of families considered in the past, including also these with PAHs2, 18, BH&HLYP/ccpVDZ was the method of choice for on–the–fly calculated reaction energies, the advantage of this functional over the most common B3LYP was also shown. However, we found this theory level as insufficient to properly capture the transition states properties of the members of the title reaction class. For that reason, in this study we utilized the M06-2X functional43, 44, designed 10 ACS Paragon Plus Environment

Page 11 of 43 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

primarily for the chemical kinetics. As shown further, the cc-pVTZ basis set works noticeable better than simpler cc-pVDZ for the family considered, thus we decided to use it here. Results for both stable forms and transition states were calculated for the lowest energy conformers. Normal mode analysis was made for stationary points to ensure their features, that is, the stable structure does not have imaginary vibrational frequencies, whereas transition states (TS) show one and only one imaginary frequency corresponding to the motion along the reaction coordinate. Geometries, energies, and frequencies (accessible at SI Tables S1-S160 and Figures S1-S80) were used to derive the RC-TST correlations and the Linear Energy Relationship (LER) between the reaction energies and barrier heights. To obtain the RC-TST correlations, the TST/Eckart rate constants for all processes from the Representative Training Set (cf. Table 1) were calculated with the TheRate45 and MSMC46 computer codes. In these calculations, overall rotations were treated classically, and vibrations were treated according to the harmonic oscillator (HO) approximation with the exception of the internal rotations of the –OH groups in TS, considered as the hindered rotations (HRs). Results were used to derive the RC-TST factors.

3. Results and discussion In this section, the calculations of the reference reaction rate constant are firstly reported. Then, derivation of the particular RC-TST factors from the electronic structure data for the Representative Training Set is described. Finally, three error analyses, estimating the accuracy of the reported RC-TST implementation, are presented.

3.1.

Reference reaction

11 ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 12 of 43

Since, according to eq. (1), the RC-TST rate constant for arbitrary process within the title class is proportional to the rate constant of the reference reaction, it is crucial to determine it with the maximum possible accuracy, which may be provided by experiment or high level calculations. In the previous RC-TST studies dealing with PAHs, the simplest reactions within a family, namely C6H6 + abstraction agent (∙H2, ∙CH317, ∙C2H517), were found to be appropriate (i.e. resulting in satisfactory agreement with available experimental data) reference processes. Their rate constants were accurately calculated from first principles and, for the limited T region, are known experimentally. As such, in this study we also decided to choose the simplest process within the title class, i.e. C6H6 + ∙OH  ∙C6H5 + H2O, as a reference reaction. 3.1.1. Potential energy surface Molecular geometry of the transition state of the C6H6 + ∙OH  ∙C6H5 + H2O process, optimized at the M06-2X/cc-pVTZ level of theory, is depicted in Figure 2. The TS was confirmed by both intrinsic reaction coordinate (IRC) calculations as well as analysis of the vibrational modes, which showed one and only one imaginary frequency, corresponding to the H transfer between C6H6 and ∙OH. Values calculated at the CBS-APNO level of theory and these by Chen et al.13 (B3LYP/6-31G(d,p)) are also presented for the sake of comparison. As it is seen, geometries predicted with both DFT functionals and the high level composite method are rather similar, except for the C-H and O-H distances, for which perceptible differences are observed. It is interesting to note, that the literature values from Chen et al.13 are closer to our CBS-APNO results rather than the M06-2X data. However, as it is shown below, the predicted M06-2X/ccpVTZ reaction energy and barrier are in a good accordance with those obtained at the high theory level, thus the M06-2X/cc-pVTZ was used to get the electronic structure data to derive the RC-TST factors. The reference reaction energies and zero point corrected barrier heights, 12 ACS Paragon Plus Environment

Page 13 of 43 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

calculated at different theory levels, are shown in Table 2. The compound model chemistries, such as CBS or G3, claimed to provide results with a chemical accuracy (≈1 kcal/mol with respect to experiment), previously received an extensive validation of their ability to model reactions involving PAHs 11, 13, 20, 21, 23, 27, 34, 35. The same is also true for the CCSD(T) method, which was also proven to work excellently for the hydrocarbon radicals40. For these reasons, the CBS-APNO, CBS-QB3 , G3B3 and CCSD(T) results are included in Table 2 and expected to be the most accurate. The reaction energies and ZPE corrected barrier heights are -3.4 and 4.4 kcal/mol for CBS-APNO, -4.0 and 4.0 for CBS-QB3 and -4.2 and 5.0 kcal/mol for G3B3, respectively. The results obtained with the CCSD(T)/cc-pVTZ are within the range of 1 kcal/mol to composite methods. Surprisingly, neither CCSD(T)/cc-pVDZ nor the CCSD(T)/CBS (extrapolation according to Truhlar’s scheme47) provide reasonable barrier. Among DFT methods, the results obtained with the M06-2X functional seem to be in the best accordance with both experiment and theoretical values from thermochemical recipes. Consequently, we used the M06-2X/cc-pVTZ theory level to calculate both the energies and hessians along the reaction path of the reference reaction and parameters needed to derive the RC-TST correlations. A complete failure of the B3LYP functional, which performs the most poorly showing very little (≈0.3-0.5 kcal/mol) 27, 48 or even no barrier (this study) for the reference reaction is also worth mentioning. This is in accordance with the previous studies, where B3LYP tends to underestimate reaction barriers (see, for example, ref.12, 49) . On the contrary, the BH&HLYP functional, which was previously widely used and extensively validated for different reactions classes29, seems to significantly overestimate the reference reaction barrier. These significant differences are due to the high spin contamination of the system under consideration. Since RC-TST relies on the relative barrier heights rather than absolute values, this spread does not remarkably perturb the

13 ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 14 of 43

final results, however. The great advantage of this methodology is that the relative reaction energies can be obtained at a fairly low theory level, such as the M06-2X /cc-pVTZ as applied here. It was shown previously50, 51, that the CCSD(T) method is fairly “resistant” to spin contamination, which is also confirmed by Table 2.

Figure 2: Geometry of the transition state of the C6H6 + ∙OH  ∙C6H5 + H2O reaction (bond lengths in Å and angles in degrees), optimized at the M06-2X/cc-pVTZ level. Values without brackets correspond to geometry calculated at the M06-2X/cc-pVTZ level of theory, values in parentheses are from the CBS-APNO optimization and numbers in square brackets are from work of Chen et al.4 at the B3LYP/6-31G(d,p) level of theory. Table 2: Calculated reaction energies (E) and zero point corrected barriers heights 𝑉𝑔𝑎 for the C6H6 + ∙OH  ∙C6H5 + H2O reaction (numbers are in kcal/mol). Since electronic data for the calculations of the rate constants are obtained at the M06-2X/cc-pVTZ level, energies corresponding to this level are provided with two decimal accuracy. Level of theory CBS-APNO CBS-QB3 G3B3

Classical reaction energy E

Classical barrier 𝑉𝐶

-3.4 -4.0 -4.2

7.7 6.1 8.3

Zero point corrected barrier 𝑉𝑔𝑎 4.4 4.0 5.0 14

ACS Paragon Plus Environment

Page 15 of 43 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

M06-2X/cc-pVTZ M06-2X/cc-pVDZ BMK/cc-pVTZ BH&HLYP/cc-pVTZ BH&HLYP/cc-pVDZ B3LYP/6-31G(d,p) QCISD(T)/cc-pVTZ//M06-2X/cc-pVTZ CCSD(T)/cc-pVDZ//M06-2X/cc-pVTZ CCSD(T)/cc-pVTZ//M06-2X/cc-pVTZ CCSD(T)/CBS//M06-2X/-cc-pVTZ experiment52 experiment53 experiment54 experiment55 B3LYP/(II)27 G3//B3LYP(II)27 G2M(CC,MP2)//B3LYP(II)27 PMP4(SDTQ)/6-311G(d,p)//B3LYP(II)27 G3(MP2)/B3LYP11 CBS-QB311 B3LYP/6-311+G(3df,2p)48 G3(MP2,CC)//B3LYP48 G3(MP2,CC)12

-5.57 -3.0 -4.1 0.5 3.2 -1.6 -2.0 3.5 -2.3 -4.9

-4.7 -4.2 -4.0 2.0 -2.9 -4.0 -5.5 -5.4 -5.2

5.19 4.3 5.6 11.3 10.3 2.1 6.6 9.9 6.5 4.5

3.23 2.4 3.3 9.1 8.2 -0.1 4.7 7.9 4.6 2.6 4.0 ± 3.0 3.3 5.1 4.5 0.3 5.3 6.6 8.9 7.77 3.51 0.42 5.5 4.3

Figure 3 illustrates the potential energies (relative to reactant values) along the MEP (Minimum Energy Path): classical potential VC and the vibrationally adiabatic ground-state potential 𝑉𝑔𝑎 , which is the sum of the two terms VC + ZPE. As mentioned before, hessians and ZPE corrections along the path were obtained at the M06-2X/cc-pVTZ level. As it is seen from Figure 3, the ZPE values typically lower the classical barrier by about 2 kcal/mol, which is about 25% of the classical energy VC. It is known, that the separation between both maxima reflects the magnitude of the variational effect in the calculated rate constants. A noticeable difference, as seen in Figure 3, may suggest that the variational effects may be of importance for the reference reaction. However, as pointed out below, the difference between TST and VTST rates are not significant in this case. 15 ACS Paragon Plus Environment

The Journal of Physical Chemistry

7

𝑉𝐶 (react=0) 𝑉𝑔𝑎 (react=0)

5

Energy (kcal/mol)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 16 of 43

3 1 -1 -3

-2

-1

0

1

2

3

Reaction coordinate s (amu0.5 bohr) -3 -5

Figure 3: Potential energy curves along the reaction coordinates of the C6H6 + ∙OH  ∙C6H5 + H2O reaction. 𝑉𝐶 represents the classical adiabatic ground-state potential, whereas 𝑉𝑔𝑎 is the vibrationally adiabatic ground-state potential. To acquire data needed to derive the RC-TST correlations for particular factors, the high pressure limit CTST (classical TST) rate constants were calculated for all processes from the Representative Training Set in the temperature region of 300–3000K, appropriate to combustion applications. 3.1.2. Rate constants Since there is no experimental data for the full T range, the thermal rate constant of the reference reaction was calculated with the POLYRATE 17b56 code. To achieve accuracy comparable with experiment, a canonical variational transition state theory (CVT), augmented with small curvature tunneling (SCT) correction methodology was chosen. The vibrations perpendicular to the reaction path were modelled by using internal curvilinear coordinates,

16 ACS Paragon Plus Environment

Page 17 of 43 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

defined by bond stretches, angle bends and torsions. Hessians were calculated at 200 (100 on each side) points along the MEP (Minimum Energy Path). The VTST calculations were performed for the reaction coordinate s range from -3 to 3. DFT energies of all 200 points within this range were scaled (by multiplying by a factor of ≈1.25) to match CCSD(T)/cc-pVTZ barrier. Energies of reactants and products were not scaled, thus the M06-2X/cc-pVTZ reaction energy was provided to the program. The rotation of –OH group in the transition state along the C-O bond is a low frequency torsional motion, thus it was treated explicitly as a hindered rotation by direct solving of the 1-D Schrödinger equation, as implemented in the MSMC code57. To account for reaction path degeneracy, a symmetry number of 6 was used58. The electronic partition function for ∙OH radical was corrected by considering the ground electronic state of OH radical59 with experimental spin-orbit coupling constant, A0 = 139.21 cm-1. The final CVT/SCT/HR rate constant, along with literature data, are plotted in Fig. 4 and fitted to an Arrhenius expression, as given below: 𝑘𝑟 (𝑇) = (1.93 × 10−21 )𝑇 3.19 𝑒𝑥𝑝(−733.5⁄𝑇)

(cm3 s-1 molecule-1)

(8)

The TST/HR rate constants were also plotted for comparison purposes. Since the CVT/HR rate is almost the same as the TST/HR, CVT/HR results are not included in Figure 4. The difference between the TST and CVT/SCT rate constants suggests substantial tunneling contribution in the low T region. Tunneling effects become less important when T > 1000 K. Unfortunately, it is not possible to firmly establish the error bars for the calculated rate constants. The only way to assess its accuracy is comparison with experimental data and/or other calculations. Experimental kinetic data from both direct measurements and fitting complex mechanism for the reference reaction are available in the NIST database60. Knispel et al.53, using resonance technique, derived reaction rate from fitting to a complex mechanism. Data reported 17 ACS Paragon Plus Environment

The Journal of Physical Chemistry

are in reasonable accordance with laser induced fluorescence experiment values reported later by Seta et al.11. As it is seen from Figure 4, our CVT/SCT data lie within error bars claimed by the experimentalists. Chen et al.13 reported high pressure limits of the Quantum Rice−Ramsperger−Kassel61 theory (QRRK) rate of the reference reaction. Their rates are generally noticeably higher than ours; the difference is minor for T>1000K, however. It is well known, that the 1D tunneling approaches often provide inaccurate description of the quantum tunneling, in most cases overestimating it. It may be thus expected, that the observed discrepancy is due to the more accurate tunneling treatment in our calculations. Generally, our CVT/SCT rates are in satisfactory agreement with literature data, including experiment. As such, their uncertainties are also believed to be comparable to these claimed by experimentalists and they will be used for estimation of the unknown rates within the title reaction family.

TST

-9.5

CVT/SCT Chen et al.

-10.5

log (kr) (cm3*molecule-1 - s-1)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 18 of 43

Seta et al. Knispel et al.

-11.5

-12.5

-13.5

-14.5

-15.5 0.0

0.5

1.0

1.5

2.0

2.5

3.0

3.5

1000/T (1/K)

18 ACS Paragon Plus Environment

Page 19 of 43 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Figure 4: Arrhenius plot of the calculated and literature rate constants for the C6H6 + ∙OH  ∙C6H5 + H2O reaction in the temperature range of 300–3000K. Literature data are taken from references: 13 (Chen et al.), 11 (Seta et al.) and 53 (Knispel et al.).

3.2.

Reaction class parameters

This section explains the derivation of the RC-TST factors based on the electronic structure data obtained from the Representative Training Set at the M06-2X/cc-pVTZ level of theory. 3.2.1. Potential energy factor, fV This factor is calculated directly with Eq. 6, where Δ𝑉𝑟# and Δ𝑉𝑎# are the barriers of the reference reaction and reaction investigated, respectively. It was shown in the previous reports40, that within a given reaction family, a linear correlation between reaction energy and barrier height holds, which is comparable to the Evans–Polanyi linear free-energy relationship. Hence, a barrier for any process from the title family can be accurately predicted from its M06-2X/ccpVTZ reaction energy. The calculated energies Δ𝐸 and barriers Δ𝑉 # for all members of the Representative Training Set are given in Table S161 in the Supporting Info. The resulted LER’s are plotted in Figure 5(a-d) below:

19 ACS Paragon Plus Environment

The Journal of Physical Chemistry

5.7

4.6

(a)  sites

4.5

5.5

V‡ (kcal/mol)

V‡ (kcal/mol)

(b) sites

5.6

4.4 4.3 4.2 4.1 4.0

5.4 5.3 5.2

5.1

3.9

5.0

-5.8

-5.6

-5.4

-5.2

-5.0

-4.8

-5.7

-5.6

-5.5

E (kcal/mol)

-5.4

-5.3

-5.2

-5.1

-5.0

-4.5

E (kcal/mol) 4.14

6.90

4.12

(c) C5 sites

6.88

(d) Other sites

4.10 4.08

6.86

V‡ (kcal/mol)

V‡ (kcal/mol)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 20 of 43

6.84 6.82

4.06 4.04

4.02 4.00 3.98

6.80

3.96 6.78

3.94

-2.1

-2.0

-1.9

-1.8

-1.7

-1.6

-1.5

-7.5

E (kcal/mol)

-7.0

-6.5

-6.0

-5.5

E (kcal/mol)

Figure 5(a-d): Linear energy relationship (LER) between barrier heights (Δ𝑉 ‡ calculated at M062X/cc-pVTZ level of theory) and classical reaction energies (∆E, calculated at M06-2X/cc-pVTZ level of theory) for the H abstractions by the ∙OH radical from (a)  sites, (b)  sites, (c) C5 sites and (d) other sites. As it is seen from Figure 5, the nature of LER depends on the abstraction site. It is interesting to note, that barriers corresponding to -site abstractions are strongly correlated with those corresponding to abstractions, for which barriers are about 1 kcal/mol higher than barriers. Abstractions from C5 sites differ significantly from reactions occurring at the  and  positions, corresponding barriers and reaction energies are the highest from the whole family. The plot corresponding to reactions occurring at other sites is almost flat; however the barriers lower with increasing E, which is not the case for other subclasses. The linear fitting expression for the abstraction reactions at the  sites is:

20 ACS Paragon Plus Environment

Page 21 of 43 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Δ𝑉𝛼‡ = 0.747 × Δ𝐸 𝑀06−2𝑋 + 8.183

(9a)

where Δ𝐸 𝑀06−2𝑋 (kcal/mol) is the classical (no ZPE added) reaction energy calculated at the M062X/cc-pVTZ level. The relationships for the other kinds of sites, which do not follow the equation above, are given below: Δ𝑉𝛽‡ = 0.887 × Δ𝐸 𝑀06−2𝑋 + 10.118 for H abstraction at the  sites

(9b)

‡ Δ𝑉𝐶5 = 0.206 × Δ𝐸 𝑀06−2𝑋 + 7.218

(9c)

for H abstraction at the C5 sites

‡ Δ𝑉𝑜𝑡ℎ𝑒𝑟 = −0.0584 × Δ𝐸 𝑀06−2𝑋 + 3.687

for H abstraction at the other sites

(9d)

Table S161 in the Supporting Info shows the absolute deviations between the barriers calculated with eqs. 9(a-d) and these obtained directly from DFT. One may notice, that the differences are typically smaller than 0.1 kcal/mol, with exception of R15 (C16-acephenanthrylene + ∙OH → 4*C16-acephenanthrylene + H2O) with error of 0.78 kcal/mol. Since both reactants and transition state of R15 are quite similar to these for other H abstractions occurring at the 5membered rings, there is no straightforward explanation why this reaction shows such an unexpected deviation. Most probably, it may be considered just as a consequence of the empirical parameterization fits using a functional form. In spite of this aberration, the mean absolute deviation (MAD) of 0.06 kcal/mol is significantly smaller than the systematic error of the DFT calculations, which is definitely acceptable in the complex process modeling. The relevance of eqs. 9(a-d) is further confirmed by the small uncertainty factor f of 1.13, which means that every EDFT energy lies between ELER/f and ELER*f. It have to be pointed out that only the relative barriers are needed in RC-TST. To compute them, the M06-2X/cc-pVTZ classical barrier of the reference reaction R1 is needed. According to Table 2, this barrier is equal to 5.19 kcal/mol. As previously reported40, for a given reaction family the average barrier depends on the kind of active reaction 21 ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 22 of 43

site (for example H abstraction from the primary, secondary etc. carbon sites). Therefore, processes occurring at particular site may be considered as a reaction subclass. Such approximation is referred to as the Barrier Height Grouping (BHG). Since there is no site orders in the PAHs, abstraction sites are divided according to their physical properties (cf. Figure 1 and Table 1). For the title reaction class there are four groups, corresponding to four LER’s; namely the abstraction at ,  and other (“atypical”) C6 sites and at C5 sites with the average barrier height of 4.32, 5.29, 6.83 and 4.04 kcal/mol, respectively. For the BHG approximation, the maximum (0.29 kcal/mol) and medium (0.11kcal/mol) deviations are larger than those for the LER one. However, this approach has a significant advantage over the LER – only the reaction symmetry data are needed to estimate the unknown rate constants. This makes its actual application in the ARMG schemes, where symmetry data are necessary to account for the reaction multiplicity, much simpler then implementation of the LER method. In conclusion, the unknown barrier height of any process within the title family may be obtained with either LER or BHG approach. The estimated barriers are then used to calculate the potential energy factor with Eq. (6). The performance of both approximations is discussed further in this study. 3.2.2. Reaction Symmetry Number Factor, fσ This factor reflects the change of the number of symmetrically equivalent reaction paths between the reference process and the reaction investigated, and is the only one not dependent on temperature. As mentioned above, its value is readily available in the ARMGs62-64, thus there is no need for major modifications of their codes. For the title reaction class, the symmetry number of the reference reaction is equal to 6. Of course, fσ for the reference process is equal to 1. According to Eq. (3), its value for an arbitrary reaction is obtained by dividing its symmetry number by 6.

22 ACS Paragon Plus Environment

Page 23 of 43 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

These values for all processes from the Representative Training Set are available in Table S162 in the Supporting Info. 3.2.3. Tunneling Factor, f Being the ratio of the transmission coefficient  of the arbitrary reaction R to that of the reference reaction R1, the factor f quantifies the difference between tunneling effect in both processes. It was previously shown that, due to the cancellation of errors, these differences can be accurately captured with the one-dimensional Eckart method65. Results obtained for the Representative Training Set were then averaged and fitted to analytical expressions within the framework of the RC-TST/SAR methodology. Our investigation showed, that the averaged factor depends on the abstraction site. In particular, it was found that both “typical” and “other” C6 sites are characterized by similar imaginary frequencies, thus we decided to group  and “other” sites to one factor; error associated with this approximation does not exceed 2%, thus the total systematic error of the method is not significantly affected. A small (but observable) difference for sites is accounted for by deriving a separate equation. The significantly larger  coefficients for reactions occurring at C5 sites indicate that tunneling is more important here than for C6 abstractions; and leads to f significantly higher than unity for C5 sites. The calculated values of f at T = 300 K for the Representative Training Set (c.f. Table 1) with error analysis are listed in Table S162 in the Supporting Info . All the factors were fitted to the analytical expressions (with uncertainty factor f of less than 1.02) as given below and plotted in Figure 6: 𝑓𝜅(𝐶6 𝛼,𝑜𝑡ℎ𝑒𝑟𝑠) (𝑇) = 1 − 1.296 × exp(−𝑇⁄117.04)

(10a)

for abstraction by OH radicals from  and “others” C6 sites; 𝑓𝜅 (𝐶6 𝛽) (𝑇) = 1 + 0.225 × 𝑒𝑥𝑝(−𝑇⁄237.29)

(10b)

23 ACS Paragon Plus Environment

The Journal of Physical Chemistry

for abstraction by OH radicals from  C6 sites; 𝑓𝜅 (𝐶5) (𝑇) = 1.063 + 7.55 × 𝑒𝑥𝑝(−𝑇⁄149.19)

(10c)

for abstraction by OH radicals from C5 sites.

2.0 1.8 ALPHA

Tunneling factor f

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 24 of 43

1.6

BETA C5

1.4

OTHERS

1.2 1.0 0.8 0

500

1000

1500

2000

2500

3000

T(K)

Figure 6: Tunneling ratio factors fκ, as a function of temperature for particular subclasses of the title reaction class.

For an arbitrary reaction within a given subclass, the same expression may be reasonably assigned, with largest absolute deviation of 0.18 for reaction R16 and the largest percentage deviation smaller than 15%. For T=300K, the MAD (mean absolute deviation) is less than 5% in comparison with Eckart coefficients. It is well known, that the tunneling contribution to the rate constants decrease with temperature, thus the errors resulting from the RC-TST correlations also decrease to reach MAD less than 2% for T=500K. Results from Fig. 6 and Table S162 also indicate, that for real combustions systems tunneling contributions for reactions occurring at C6 sites may be safely neglected. However, such approximation may cause substantial errors for abstractions from C5 sites (underestimation). 24 ACS Paragon Plus Environment

Page 25 of 43 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

3.2.4. Partition Function Factor, fQ The total partition function includes the translational, rotational, electronic, vibrational and internally rotational components. This separation is also reflected in the partition function factor fQ.. Among its components, two first (translational and rotational) do not depend on temperature; the temperature dependent part originates mainly from different coupling of the reactive moiety with the substituents, and it is related to vibrational and internally rotational parts40. The physical nature of the partition function is difficult to interpret. It is known from the previous studies49, that fQ decreases with increasing interaction strength between reactive moiety and its substituents. For the reactions with aliphatic species, its values for tertiary reaction sites is smaller than for secondary sites, fQ for secondary sites are smaller than primary ones etc. In this study aromatic species are considered, thus these findings do not simply apply to our results. However, some similarities may be found – as the reference reaction occurs with a single ring, and all others with multiple ring species, the stronger coupling with the increasing number of rings may be also expected. For that reason, for arbitrary process within a title class, the coupling between substituents is stronger than it is for the reference reaction and, consequently, the fQ value is expected to be smaller than 1. Indeed, as it is seen in Figure 7, it is true for fQ factor for  and “atypical” sites, for which averaged factors are significantly lower than unity. It may be concluded, that the physical nature of these sites are far different to the C6H6 sites. This becomes even more of a case for “atypical” abstractions sites, which is due to the significant difference in their nature in comparison with benzene.  sites are the most similar to these in benzene, thus appropriate values of fQ are expected to be close to unity. Indeed, its value is about 1.1 in the whole temperature range. For C5 sites, the partition function factor is larger than unity, which may suggest smaller coupling between these sites and the abstraction agent than what occurs at benzene. 25 ACS Paragon Plus Environment

The Journal of Physical Chemistry

1.8 1.6

Partition function factor fQ

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 26 of 43

1.4 1.2 Alpha

1.0

Beta 0.8

C5 Others

0.6 0.4 0.2 0.0 0

500

1000

1500

2000

2500

3000

T [K]

Figure 7: Partition function factors fQ as a function of temperature for different H abstraction sites. Since the hindered rotation (HR) modes are handled separately, their contributions were not included in the calculations of partition function. The temperature dependence of the averaged partition function factors for particular subclasses, are plotted in Figure 7 and fitted as below: 𝑓𝑄 (𝐶6,𝛼) (𝑇) = 0.56 − 0.2 × 𝑒𝑥𝑝(−𝑇⁄308.8)

(11a)

for abstraction by OH radical at C6, sites; 𝑓𝑄 (𝐶6 𝛽,) (𝑇) = 1.1

(11b)

for abstraction by OH radical at C6, sites; 𝑓𝑄 (𝐶5,) (𝑇) = 1.42 − 0.55 × 𝑒𝑥𝑝(−𝑇⁄241.86)

(11c)

for abstraction by OH radical at C5 sites and 𝑓𝑄 (𝐶6, 𝑜𝑡ℎ𝑒𝑟𝑠) (𝑇) = 0.33 − 0.15 × 𝑒𝑥𝑝(−𝑇⁄312.5)

(11d)

26 ACS Paragon Plus Environment

Page 27 of 43 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

for abstraction by OH radical at C6 “others” (atypical), i.e. zigzag, armchair and mixed zigzag– armchair sites. The small uncertainty factor of 1.02, associated with using approximate formulas 11(a-d), further confirms the applicability of the RC-TST/SAR correlations. 3.2.5. Hindered Rotation (HR) Factor, fHR As all the RC-TST factors, the hindered rotations factor fHR(T) captures the changeability of a certain quantity (here it is the contribution of low frequency motions to rate constants) from the reference reaction to another arbitrary process within a given family. As pointed out in the Section 3.1.2, the low frequency modes are already treated separately in the rate constant of the reference reaction. The hindered rotation factor is than a measure of the influence of other rings in the PAH structure (they may be considered as substituents) to the contribution of the hindered rotors to the total rate constants. To capture such influence we utilized the approach proposed by East et al. and Kilpatrick et al.66, as implemented in the MSMC code57. The rotation of the –OH group along the C-OH axis in the transition states was explicitly considered as a low frequency hindered rotor. Our investigation showed, that potential energy curves for this motion for benzene (corresponding to reference reaction R1) differs from that for the multi-ring structures (any other reaction within the title family), thus the hindered rotation factor is expected to differ from unity. As it is seen from the results (see Figures S81-S100 in the Supporting Info for rotational potential and the HR/HO values), this expectation was confirmed by calculations. We also found, that the fHR does not depend on the kind of the active site. As such, for simplicity, only one factor (averaged over reactions R2-R34) is used for the whole training set and its value is plotted in Figure 8. The hindered rotations factor is smaller than one in the whole temperature range (which means that low frequency motions are less important for reactions involving multi-ring species than for benzene

27 ACS Paragon Plus Environment

The Journal of Physical Chemistry

+ ∙OH process) and increases with the temperature. For the temperature range of 300–3000K, fHR was fitted to the following expression, with r2 > 0.99 and uncertainty factor of less than 1.01: 𝑓𝐻𝑅 (𝑇) = −1.6 × 10−8 × 𝑇 2 + 8.75 × 10−5 × 𝑇 + 0.745

(12)

0.90 Hindered rotations factor fHR

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 28 of 43

0.85

0.80

0.75

0.70 0

500

1000

1500 T [K]

2000

2500

3000

Figure 8: Hindered rotations factor fHR as a function of temperature.

3.3.

Rate Constants Prediction

In previous sections, analytical expressions (RC-TST correlations) for all factors needed for the application of the RC-TST/SAR to calculate rate constant of an arbitrary reaction within the title family were established. To calculate the rate constant of an arbitrary reaction within the PAH + ∙OH → PAH radical + H2O family, the following steps have to be taken: (1) utilize Eq. (6) to calculate the potential energy factor fV, with the barrier of the reference process equal to 5.15 kcal/mol. The barrier of reaction investigated can be calculated either with LER (Eqs. 9(a-d)) or BHG approaches; (2) use Eq. (3) to calculate the symmetry number factor of the reaction investigated or see Table S162; (3) evaluate the tunneling factor using Eqs. (10a-c); (4) use Eqs. 11(a-d) to compute the partition function factor fQ; (5) apply Eq. 12 to calculate the hindered

28 ACS Paragon Plus Environment

Page 29 of 43 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

rotations factor fHR ; and (6) according to Eq. (8) obtain the final rate constants by taking a product of the rate constant of the reference reaction given by Eq. (8) and the five factors calculated above. All the RC-TST/SAR parameters needed are summarized in the Table 3 below. If, instead of Eqs. 9(a-d), BGH methodology is applied, an averaged barrier for the particular subclass should be used, namely 4.32, 5.29, 6.83 or 4.04 kcal/mol for the abstractions at , , other (“atypical”) C6 sites and C5 sites, respectively. The resulting rate is then denoted as RC-TST/BHG. Table 3 summarizes the RC-TST parameters for this reaction class. If the BHG barrier heights and average values for other factors are used, the rate constants are denoted by RC-TST/SAR-BHG. For simplicity, the RC-TST/SAR-BHG correlations (per reaction site, units are cm3/molecule∙s), valid for any reaction belonging to the title class, are presented below. Within this framework, the rate constants can be estimated without any further calculations as: 𝑘(𝑇) = 2.1 × 10−22 × 𝑇 3.4 × 𝑒𝑥𝑝(−411.2⁄𝑇)

(13a) for abstraction from C6 sites

𝑘(𝑇) = 1.1 × 10−22 × 𝑇 3.46 × 𝑒𝑥𝑝(−415.9⁄𝑇)

(13b) for abstraction from C6 sites

𝑘(𝑇) = 1.7 × 10−22 × 𝑇 3.45 × 𝑒𝑥𝑝(−795.4⁄𝑇)

(13c) for abstraction fromC5 sites

𝑘(𝑇) = 2.7 × 10−22 × 𝑇 3.35 × 𝑒𝑥𝑝(−420.9⁄𝑇)

(13d)

for

abstraction

fromothers

(“atypical”) C6 sites. The rates calculated with the above equations as a function of temperature are plotted in Figure S102 in the Supporting Info file. The accuracy of the BHG approach is obviously inferior to that of the LER approximation. The ease of ready to use rate expressions for any reaction within the class may compensate the lower accuracy, however. It should be noted, that large values of the tunneling coefficients at low temperature may lead to significant deviations in the fitted rate coefficients. Hence, the rates fitted at high (500K ≤ T ≤ 3000K) temperature range would be more 29 ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 30 of 43

accurate for combustion applications. The RC-TST parameters for this temperature region are provided in Tables S163 (RC-TST/SAR-BHG) and S164 (RC-TST/SAR-LER) in the Supporting Info file.

Table 3: Parameters and formulations of the RC-TST/SAR-LER method for the PAH + •OH → PAH radical + H2O reaction family (C6H6 + •OH → •C6H5 + H2O is the reference reaction). 𝑘𝑅𝐶−𝑇𝑆𝑇/𝑆𝐴𝑅 (𝑇) = 𝑓𝜎 × 𝑓𝜅 (𝑇) × 𝑓𝑄 (𝑇) × 𝑓𝑉 (𝑇) × 𝑓𝐻𝑅 (𝑇) × 𝑘𝑟𝑒𝑓 (cm3/molecule∙s)

fs f(T)

calculated explicitly from the symmetry of reactions (see Table 3) 𝑓𝜅(𝐶6 𝛼,𝑜𝑡ℎ𝑒𝑟𝑠) (𝑇) = 1 − 1.296 × exp(−𝑇⁄117.04)

𝑓𝜅(𝐶6,𝛽) (𝑇) = 1 + 0.225 × 𝑒𝑥𝑝(−𝑇⁄237.29) 𝑓𝜅(𝐶5) (𝑇) = 1.063 + 7.55 × 𝑒𝑥𝑝(−𝑇⁄149.19)

fQ(T)

for C6 and C6 “others” ( C6 sites

,

and

for C6 sites

for C5 sites

𝑓𝑄(𝐶6,𝛼) (𝑇) = 0.56 − 0.2 × 𝑒𝑥𝑝(−𝑇⁄308.8)

for C6 sites

𝑓𝑄(𝐶6,𝛽) (𝑇) = 1.1

for C6 sites

𝑓𝑄(𝐶5) (𝑇) = 1.42 − 0.55 × 𝑒𝑥𝑝(−𝑇⁄241.86)

for C5 sites

𝑓𝑄(𝐶6,𝑜𝑡ℎ𝑒𝑟𝑠) (𝑇) = 0.33 − 0.15 × 𝑒𝑥𝑝(−𝑇⁄312.5)

for “others” (

fHR(T)

fHR (T) = −1.6 × 10−8 × T 2 + 8.75 × 10−5 × T + 0.745

for all sites

fV(T)

exp  ( V ‡  V ‡ref ) / ( kBT )

,

and

) C6 sites

(zero-point energy correction is not included)

Δ𝑉𝛼‡ = 0.748 × Δ𝐸 𝑀06−2𝑋 + 8.183 (kcal/mol) Δ𝑉𝛽‡ = 0.887 × Δ𝐸 𝑀06−2𝑋 + 10.118 (kcal/mol) ‡ Δ𝑉𝐶5 = 0.206 × Δ𝐸 𝑀06−2𝑋 + 7.218 (kcal/mol) ‡ Δ𝑉𝑜𝑡ℎ𝑒𝑟 = −0.0584 × Δ𝐸 𝑀06−2𝑋 + 3.687 (kcal/mol)

for C6 sites for C6 sites for C5 sites for “others” (

,

and

≠ ∆𝑉𝑟𝑒𝑓 = 5.19 kcal/mola

30 ACS Paragon Plus Environment

) C6 sites

Page 31 of 43 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

kref (T ) (cm3/molec ule∙s) a

𝑘𝑟𝑒𝑓 (𝑇) = (1.93 × 10−21 )𝑇 3.19 𝑒𝑥𝑝(−733.5⁄𝑇) (cm3 s-1 molecule-1)

Barrier height for the reference reaction R1 calculated at the M06-2X /cc-pVTZ level of theory.

3.4.

Error Analyses

The first error analysis compares RC-TST results with available literature data. As mentioned in the Introduction, only very limited experimental data exist for the title reaction class. Actually, the only experimental values available are for naphthalene34, anthracene26 and phenanthrene33. These studies report on the total rate of the PAH + OH reaction which, as mentioned above, is the sum of H abstraction and ∙OH addition channels. However, since the addition channel becomes negligible for T>500K29, only experimental data for T>700K are considered here. Furthermore, the total rate constants reported by experimentalists cannot be directly compared to the rates reported in this study, which are individual rate constants for single abstraction sites. Instead, the total RC-TST rates, being the sum of individual abstraction rates per site, are presented here. The above mentioned reactions include individual channels corresponding to ,  and zigzag abstraction sites. We are not aware of any experimental data for abstractions from C6 armchair and C5 sites. Recently, Semenikhin et al. 12 presented the CTST kinetic data for the H abstractions from PAHs chrysene and benzo[a]pyrene, consisting of four and five aromatic rings, respectively. The electronic structure data were calculated at the G3(MP2,CC) theory level. In contrast to the methodology proposed here, they considered only two kinds of H abstraction sites, namely zigzag and armchair. Whereas the “armchair” was defined in the same way as in Fig. 1, the zigzag sites are actually all the C6 sites other than “armchair”, i.e. , , zigzag and mixed zigzag-armchair. 31 ACS Paragon Plus Environment

The Journal of Physical Chemistry

The same definition of the armchair sites facilitates direct comparison with the RC-TST results for “atypical” sites. Computational results are also available for naphthalene24, 48. Figure 9(a-d) shows Arrhenius plots of the calculated and measured rate constants for reactions (9a) C10H8 + ·OH → ∙C10H7 + H2O34; (9b) C14H10 (anthracene) + ·OH → ∙C14H9 + H2O26; (9c) C14H10 (phenanthrene) + ·OH → ∙C14H9 + H2O33 and (9d) benzo(a)pyrene (BaP) + ·OH → ∙BaP (armchair) + H2O12. Error bars claimed by the experimentalist are provided for their data. As one may see from Figures 9(a-d), the rate constants obtained by the RC-TST/SAR approach with both LER and BHG approximations are in reasonable agreement with the literature results. However, there is significant discrepancy between our data and these reported by Shiroudi et al.24.

-9

-9

b) Anthracene + ∙OH → Anthracene*(total) + H2O

a) Naphtalene + ∙OH → Naphtalene*(total) + H2O log{k(T){cm3molecule-1s-1}

RC-TST/SAR RC-TST/BHG

-11

Lorenz 83

-12

Shiroudi 2014 total

RC-TST/SAR

-10

log{k(T){cm3molecule-1s-1}

-10

RC-TST/BHG

-11

Ananthula 2006

-12 -13

-13

-14

-14

-9

0.5

1.0

1.5 2.0 1000/T [1/K]

2.5

3.0

c) Phenantrene + ∙OH→ Phenantrene*(total)+ H2O RC-TST/SAR

-10

RC-TST/BHG

-11

0.0

3.5

Ananthula 2007

-12

-10

0.5

1.0

1.5 2.0 1000/T [1/K]

2.5

3.0

3.5

d) BaP+ ∙OH → BaP* (armchair)+ H2O RC-TST/SAR

-11

log{k(T){cm3molecule-1s-1}

0.0

log{k(T){cm3molecule-1s-1}

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 32 of 43

RC-TST/BHG

-12

Semenikhin 2017

-13

-13

-14

-14

-15 0.0

0.5

1.0

1.5 2.0 1000/T [1/K]

2.5

3.0

3.5

0.0

0.5

1.0

1.5 2.0 1000/T [1/K]

2.5

3.0

3.5

32 ACS Paragon Plus Environment

Page 33 of 43

Figure 9(a-d): Arrhenius plots of the calculated total rate constants using RC-TST/SAR and RCTST/BHG approaches for four reactions, along with the available literature values: Lorenz 83 – ref.34, Shiroudi 2014 –ref.24, Ananthula 2006 – ref.26, Ananthula 2007 – ref.26 and Semenikhin 2017 – ref.12.

A more systematic way of assessment of the performance of the method is comparison of the results computed directly from formulas (2-7) with approximations: RC-TST/SAR/LER (Table 3) and RC-TST/SAR/BHG (eqs. 13(a-d)).We performed such analysis for 33 reactions from the Representative Training Set (with exception of the reference reaction R1). Results are plotted in Figure S103(a-b) in the Supporting Info, detailed discussion of these results is also provided in the SI file.

Total factor (M06-2X) Total factor (BHG) Tunneling factor Partition function factor Potential energy factor (M06-2X) Potential energy factor (BHG) Hindered rotations factor

0.40 0.35 0.30 0.25

Error

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

0.20 0.15 0.10 0.05 0.00 0

500

1000

1500

2000

2500

3000

T [K]

33 ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 34 of 43

Figure 10: Mean absolute errors (%) of the total RC-TST factors f(T) (Eq. (2)) obtained with LER and BHG approaches and its components as a function of the temperature.

The last error analysis is on the systematic errors introduced by particular components of the total RC-TST factor (cf. Eq. 2). The deviations between values computed with formulas from Table 3 and these calculated directly from the M06-2X/cc-pVTZ data were determined for every reaction from the Representative Training Set for the temperature range of 300-3000K and then averaged over the whole family, with exception of the reference reaction R1. Results as functions of the temperature are plotted in Figure 10. The error of the total factor is affected by the errors of its components, namely potential energy factor fV(T), tunneling factor f(), partition function factor fQ(T) and hindered rotations factor fHR(T). As it is seen from Fig. 10, the error comes mainly from the potential energy and partition function components. The former one, also mainly responsible for the temperature dependence of the total error, decreases with temperature and makes this dependence perceptible only for T500K. As discussed previously, error introduced by applying the BHG approximation is higher than that associated with LER (28% vs. 18% at 300K, respectively). Of all the factors, the tunneling one f() introduces the smallest (less than 5%) error, decreasing with temperature. The error related to the HR coefficient is less than 10% for the whole T range, whereas the one introduced by the partition function factor is the largest of individual component errors, even bigger than these for total LER and (for T ≈ 1500 K) also total BHG factor. This can be explained by a mutual error cancelation. Overall, the LER approach is significantly more accurate than the BHG in the whole temperature range, with error remaining almost constant for T>800K.

3.5.

Rate based rules 34 ACS Paragon Plus Environment

Page 35 of 43 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

The RC-TST correlations, summarized in Table 3, enable one to draw wider conclusions regarding the dependence of the rate coefficient on the molecular size for a homologous reaction series and a comparison among different sites. Figure 11 shows a plot of the average RC-TST/SAR rate constants (T=300K) for H abstractions from different sites as a function of the number of rings in the PAH molecule. It is seen, that drawing simple conclusions regarding the dependence of the rate constants on the molecule’s size is not straightforward. There is a clear tendency for  and C5 sites – for  sites the averaged rate decreases with the number of rings, whereas for C5 it remains almost constant. There is no observable tendency for “atypical” sites. However, since the chemical properties of such sites are very diverse, their systematization is problematic. Therefore, the performance of the RC-TST methodology for these kinds of H abstractions is noticeably worse than for ‘‘typical’’ ones; their description should be considered rather as semi-quantitative only. The relations between averaged rates corresponding to particular sites are also not simple – rates for C5 abstractions are definitely the lowest, they are about 4-5 times slower than abstractions from

 sites. Abstractions from  and “other” sites show the highest rates. It is noteworthy that in spite of the lack of clear tendencies in the averaged rates, there is a significant discrepancy between particular single rates being averaged. This shows, that the simple rate based rules67, widely reported in the literature, although effectively utilized in modeling a broad range of chemistries67, may generally not be relevant to the challenge of PAHs, which is due to the complicated chemistry of the conjugated multi-ring structures.

35 ACS Paragon Plus Environment

The Journal of Physical Chemistry

2.5E-14

other sites kaverage(300K) (cm3molecule-1s-1)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 36 of 43

2.0E-14



1.5E-14

1.0E-14



C5

5.0E-15

0.0E+00 0

1

2

3

4

5

6

7

8

No. of rings in PAH Figure 11: Averaged RC-TST/SAR rate constants (T=300K) of the H abstractions from PAH by ∙OH radical from particular sites as a function of the number of rings in the PAH molecule. The average is taken over both the representative training set (cf. Table 1) and reactions not present in the training set, but considered in Figure 9 above.

4. Conclusions In this report, the kinetic RC-TST/SAR methodology was applied for the reliable on-thefly prediction of thermal rate constants of the PAH + ∙OH →PAH radical + H2O reaction family, rate based analysis on the basis of the RC-TST results was also performed. PAH abstraction sites were divided into four main groups. There is no clear tendency regarding the dependence of the rate constants versus PAH molecule size. However, such trend is noticeable about the dependence 36 ACS Paragon Plus Environment

Page 37 of 43 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

on the kind of the abstraction site, namely kC5 < k < k in average. The “atypical” sites (i.e. other than ,  or C5) were found to be too different in their chemical nature to exhibit any observable tendency. The proposed methodology can be effectively utilized in the Automated Reaction Mechanism Generator (ARMG) schemes for construction of the comprehensive mechanisms of the pyrolysis, partial oxidation and combustion of the PAHs–including fuels; it is also useful in the modeling of the process of soot formation and growth., The error analyses showed satisfactory agreement of the results in comparison with available experiments. The systematic error of the method was found to be within a factor of two when compared to the explicit high level rate calculations. Generally, it was determined that the RC-TST/SAR(LER), where only the symmetry (number of symmetrically equivalent reaction paths) data and reaction energy are necessary, and RC-TST/SAR(BHG), where only symmetry data are needed, are an encouraging and fairly comprehensive frameworks to predict the rate constant of any reaction within a title family. The worse performance of the BHG methodology can be offset by the handiness of ready to use rate expressions for any family member.

5. Associated Content * S Supporting Information

Tables: The optimized geometries and frequencies of all species calculated at the M06-2X/ccpVTZ level of theory for the representative set. Classical reaction energies and absolute deviations between the calculated barrier heights from LER and BHG approaches at 0K. Tunneling and symmetry factors for the reactions from the Representative Training Set. Parameters and formulations of RC-TST/SAR method for the title reaction class in the high temperature range 500-3000K.

37 ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 38 of 43

Figures: Optimized geometries of all species calculated at the M06-2X/cc-pVTZ level of theory for the representative set. The hindrance potential and HO/HR values for selected processes. BHG rates plots as a function of temperature. Some error statistics. RC-TST/SAR rate constants as a function of the number of rings.

6. Author Information Corresponding Author: Artur Ratkiewicz, e-mail address: [email protected].

7. Conflicts of interest There are no conflicts of interest to declare.

8. Acknowledgments The authors would like to thank the Computational Center of the University of Bialystok (Grant GO-008) for providing access to the supercomputer resources and the GAUSSIAN 09 program.

9. References 1. Watson, M. D.; Fechtenkötter, A.; Müllen, K. Big Is Beautiful−“Aromaticity” Revisited from the Viewpoint of Macromolecular and Supramolecular Benzene Chemistry. Chem. Rev. 2001, 101 (5), 12671300. 2. Violi, A.; Truong, T. N.; Sarofim, A. F. Kinetics of Hydrogen Abstraction Reactions from Polycyclic Aromatic Hydrocarbons by H Atoms. J. Phys. Chem. A 2004, 108 (22), 4846-4852. 3. Allen, J. O.; Dookeran, N. M.; Smith, K. A.; Sarofim, A. F.; Taghizadeh, K.; Lafleur, A. L. Measurement of Polycyclic Aromatic Hydrocarbons Associated with Size-Segregated Atmospheric Aerosols in Massachusetts. Env. Sci. Tech. 1996, 30 (3), 1023-1031. 4. Hemelsoet, K.; Van Speybroeck, V.; Marin, G. B.; De Proft, F.; Geerlings, P.; Waroquier, M. Reactivity Indices for Radical Reactions Involving Polyaromatics. J. Phys. Chem. A 2004, 108 (35), 72817290. 5. Saeys, M.; Reyniers, M.-F.; Marin, G. B.; Van Speybroeck, V.; Waroquier, M. Ab Initio Calculations for Hydrocarbons: Enthalpy of Formation, Transition State Geometry, and Activation Energy for Radical Reactions. J. Phys. Chem. A 2003, 107 (43), 9147-9159. 6. Hemelsoet, K.; Moran, D.; Van Speybroeck, V.; Waroquier, M.; Radom, L. An assessment of theoretical procedures for predicting the thermochemistry and kinetics of hydrogen abstraction by methyl radical from benzene. J. Phys. Chem. A 2006, 110 (28), 8942-51.

38 ACS Paragon Plus Environment

Page 39 of 43 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

7. Saeys, M.; Reyniers, M. F.; Van Speybroeck, V.; Waroquier, M.; Marin, G. B. Ab Initio Group Contribution Method for Activation Energies of Hydrogen Abstraction Reactions. ChemPhysChem 2006, 7 (1), 188-99. 8. Hou, D.; You, X. Reaction kinetics of hydrogen abstraction from polycyclic aromatic hydrocarbons by H atoms. Phys. Chem. Chem. Phys. 2017, 19 (45), 30772-30780. 9. Battin-Leclerc, F.; Curran, H.; Faravelli, T.; Glaude, P., A. Specificities Related to Detailed Kinetic Models for the Combustion of Oxygenated Fuels Components. Chapter 4 in Cleaner Combustion: Developing Detailed Chemical Kinetic Models; Battin-Leclerc, F.; Simmie, J. M.; Edward Blurock, E., Eds.; Springer-Verlag: London 2013; pp 93-109. 10. Frenklach, M.; Liu, Z.; Singh, R. I.; Galimova, G. R.; Azyazov, V. N.; Mebel, A. M. Detailed, sterically-resolved modeling of soot oxidation: Role of O atoms, interplay with particle nanostructure, and emergence of inner particle burning. Combust. Flame 2018, 188, 284-306. 11. Seta, T.; Nakajima, M.; Miyoshi, A. High-Temperature Reactions of OH Radicals with Benzene and Toluene. J. Phys. Chem. A 2006, 110 (15), 5081-5090. 12. Semenikhin, A. S.; Savchenkova, A. S.; Chechet, I. V.; Matveev, S. G.; Liu, Z.; Frenklach, M.; Mebel, A. M. Rate constants for H abstraction from benzo(a)pyrene and chrysene: a theoretical study. Phys. Chem. Chem. Phys. 2017, 19 (37), 25401-25413. 13. Chen, C.-C.; Bozzelli, J. W.; Farrell, J. T. Thermochemical Properties, Pathway, and Kinetic Analysis on the Reactions of Benzene with OH:  An Elementary Reaction Mechanism. J. Phys. Chem. A 2004, 108 (21), 4632-4652. 14. Edwards, D. E.; Zubarev, D. Y.; Lester, W. A.; Frenklach, M. Pathways to Soot Oxidation: Reaction of OH with Phenanthrene Radicals. J. Phys. Chem. A 2014, 118 (37), 8606-8613. 15. Gao, L. G.; Zheng, J.; Fernández-Ramos, A.; Truhlar, D. G.; Xu, X. Kinetics of the Methanol Reaction with OH at Interstellar, Atmospheric, and Combustion Temperatures. J. Am. Chem. Soc. 2018, 140 (8), 2906-2918. 16. Hemelsoet, K.; Van Speybroeck, V.; Moran, D.; Marin, G. B.; Radom, L.; Waroquier, M. Thermochemistry and kinetics of hydrogen abstraction by methyl radical from polycyclic aromatic hydrocarbons. J. Phys. Chem. A 2006, 110 (50), 13624-31. 17. Li, C. Y.; Agarwal, J.; Wu, C. H.; Allen, W. D.; Schaefer, H. F. Intricate Internal Rotation Surface and Fundamental Infrared Transitions of the n-Propyl Radical. J. Phys. Chem. B 2015, 119 (3), 728-735. 18. Mai, T. V. T.; Ratkiewicz, A.; Le, A.; Duong, M. v.; Truong, T. N.; Huynh, L. K. On-the-fly kinetics of hydrogen abstraction from polycyclic aromatic hydrocarbons by methyl/ethyl radicals. Phys. Chem. Chem. Phys. 2018, 20 (36), 23578-23592. 19. Li, S.-H.; Guo, J.-J.; Li, R.; Wang, F.; Li, X.-Y. Theoretical Prediction of Rate Constants for Hydrogen Abstraction by OH, H, O, CH3, and HO2 Radicals from Toluene. J. Phys. Chem. A 2016, 120 (20), 34243432. 20. Pelucchi, M.; Cavallotti, C.; Faravelli, T.; Klippenstein, S. J. H-Abstraction reactions by OH, HO2, O, O2 and benzyl radical addition to O2 and their implications for kinetic modelling of toluene oxidation. Phys.Chem. Chem. Phys. 2018, 20 (16), 10607-10627. 21. Yao, Q.; Cao, X.-M.; Zong, W.-G.; Sun, X.-H.; Li, Z.-R.; Li, X.-Y. Potential Energy Surface for Large Barrierless Reaction Systems: Application to the Kinetic Calculations of the Dissociation of Alkanes and the Reverse Recombination Reactions. J. Phys. Chem. A 2018. 22. Buras, Z. J.; Chu, T.-C.; Jamal, A.; Yee, N. W.; Middaugh, J. E.; Green, W. H. Phenyl radical + propene: a prototypical reaction surface for aromatic-catalyzed 1,2-hydrogen-migration and subsequent resonance-stabilized radical formation. Phys. Chem. Chem. Phys. 2018. 23. Shiroudi, A.; Deleuze, M. S.; Canneaux, S. Theoretical Study of the Oxidation Mechanisms of Naphthalene Initiated by Hydroxyl Radicals: The OH-Addition Pathway. J. Phys. Chem. A 2014, 118 (26), 4593-4610. 39 ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 40 of 43

24. Shiroudi, A.; Deleuze, M. S. Theoretical Study of the Oxidation Mechanisms of Naphthalene Initiated by Hydroxyl Radicals: The H Abstraction Pathway. J. Phys. Chem. A 2014, 118 (20), 3625-3636. 25. Sirjean, B.; Fournet, R. Theoretical study of the reaction 2,5-dimethylfuran+H→products. Proc. Combust. Inst. 2013, 34 (1), 241-249. 26. Ananthula, R.; Yamada, T.; Taylor, P. H. Kinetics of OH Radical Reaction with Anthracene and Anthracene-d10. J. Phys. Chem. A 2006, 110 (10), 3559-3566. 27. Tokmakov, I. V.; Lin, M. C. Kinetics and Mechanism of the OH + C6H6 Reaction:  A Detailed Analysis with First-Principles Calculations. J. Phys. Chem. A 2002, 106 (46), 11309-11326. 28. Yao, Q.; Cao, X.-M.; Zong, W.-G.; Sun, X.-H.; Li, Z.-R.; Li, X.-Y. Potential Energy Surface for Large Barrierless Reaction Systems: Application to the Kinetic Calculations of the Dissociation of Alkanes and the Reverse Recombination Reactions. J. Phys. Chem. A 2018, 122 (21), 4869-4881. 29. Atkinson, R. Kinetics and mechanisms of the gas-phase reactions of the hydroxyl radical with organic compounds under atmospheric conditions. Chem. Rev. 1986, 86 (1), 69-201. 30. Atkinson, R. J. Kinetics and mechanisms of the gas-phase reactions of the hydroxyl radical with organic compounds. Phys. Chem. Ref. Data, Monogr. 1989, 1, 1-246. 31. Lay, T. H.; Bozzelli, J. W.; Seinfeld, J. H. Atmospheric Photochemical Oxidation of Benzene:  Benzene + OH and the Benzene−OH Adduct (Hydroxyl-2,4-cyclohexadienyl) + O2. J, Phys, Chem. 1996, 100 (16), 6543-6554. 32. Dean, A. M. Predictions of pressure and temperature effects upon radical addition and recombination reactions. J. Phys. Chem. 1985, 89 (21), 4600-4608. 33. Ananthula, R.; Yamada, T.; Taylor, P. H. Kinetics of OH radical reaction with phenanthrene: New absolute rate measurements and comparison with other PAHs. Int. J. Chem. Kin. 2007, 39 (11), 629-637. 34. Lorenz, K.; Zellner, R. Kinetics of the Reactions of OH-Radicals with Benzene, Benzene-d6 and Naphthalene. Ber. Bunsen-Ges. Phys. Chem. 1983, 87, 629−636. 35. Atkinson, R.; Aschmann, S. M. Kinetics of the reactions of naphthalene, 2-methylnaphthalene, and 2,3-dimethylnaphthalene with OH radicals and with O3 at 295 ± 1 K. Int. J. Chem. Kin. 1986, 18 (5), 569-573. 36. Uc, V. H.; Alvarez-Idaboy, J. R.; Galano, A.; García-Cruz, I.; Vivier-Bunge, A. Theoretical Determination of the Rate Constant for OH Hydrogen Abstraction from Toluene. J. Phys. Chem. A 2006, 110 (33), 10155-10162. 37. Shi, S. Advances in modeling hydrocarbon cracking kinetic predictions by quantum chemical theory: A review. Int. J. Energy Res. 2018, 42 (10), 3164-3181. 38. Sumathi, R.; Green Jr., W. H. A priori rate constants for kinetic modeling. Theor. Chem. Acc. 2002, 108 (4), 187-213. 39. Van de Vijver, R.; Sabbe, M. K.; Reyniers, M.-F.; Van Geem, K. M.; Marin, G. B. Ab initio derived group additivity model for intramolecular hydrogen abstraction reactions. Phys. Chem. Chem. Phys. 2018, 20 (16), 10877-10894. 40. Ratkiewicz, A.; Huynh, L. K.; Truong, T. N. Performance of First-Principles-Based Reaction Class Transition State Theory. J. Phys. Chem. B 2016, 120 (8), 1871-1884. 41. Bao, J. L.; Truhlar, D. G. Variational transition state theory: theoretical framework and recent developments. Chem. Soc. Rev. 2017, 46 (24), 7548-7596. 42. Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Scalmani, G.; Barone, V.; Mennucci, B.; Petersson et al., J. Gaussian 09, Revision A.1, Gaussian, Inc.: Wallingford CT, 2009. 43. Zhao, Y.; Truhlar, D. G. Density Functionals with Broad Applicability in Chemistry. Acc. Chem. Res. 2008, 41 (2), 157-167. 44. Zhao, Y.; Truhlar, D. G. The M06 suite of density functionals for main group thermochemistry, thermochemical kinetics, noncovalent interactions, excited states, and transition elements: two new 40 ACS Paragon Plus Environment

Page 41 of 43 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

functionals and systematic testing of four M06-class functionals and 12 other functionals. Theor. Chem. Acc. 2008, 120 (1), 215-241. 45. Duncan, W. T.; Bell, R. L.; Truong, T. N. TheRate: Program for ab initio direct dynamics calculations of thermal and vibrational-state-selected rate constants. J. Comp. Chem. 1998, 19 (9), 10391052. 46. Duong, M. V.; Nguyen, H. T.; Truong, N.; Le Thong, N.-M.; Huynh, L. K. Multi-Species MultiChannel (MSMC): An Ab Initio-based Parallel Thermodynamic and Kinetic Code for Complex Chemical Systems. Int. J. Chem. Kin. 2015, 47 (9), 564-575. 47. Truhlar, D. G. Basis-set extrapolation. Chem. Phys. Lett. 1998, 294 (1), 45-48. 48. Kislov, V. V.; Islamova, N. I.; Kolker, A. M.; Lin, S. H.; Mebel, A. M. Hydrogen Abstraction Acetylene Addition and Diels−Alder Mechanisms of PAH Formation:  A Detailed Study Using First Principles Calculations. J. Chem. Theor. Comp. 2005, 1 (5), 908-924. 49. Huynh, L. K.; Panasewicz, S.; Ratkiewicz, A.; Truong, T. N. Ab Initio Study on the Kinetics of Hydrogen Abstraction for the H + Alkene → H2 + Alkenyl Reaction Class. J. Phys. Chem. A 2007, 111 (11), 2156-2165. 50. Villà, J.; Corchado, J. C.; González-Lafont, A.; Lluch, J. M.; Truhlar, D. G. Variational TransitionState Theory with Optimized Orientation of the Dividing Surface and Semiclassical Tunneling Calculations for Deuterium and Muonium Kinetic Isotope Effects in the Free Radical Association Reaction H + C2H4 → C2H5. J. Phys. Chem. A 1999, 103 (26), 5061-5074. 51. Ratkiewicz, A. Kinetics of the C–C bond beta scission reactions in alkyl radicals. Phys. Chem. Chem. Phys. 2011, 13 (33), 15037-15046. 52. Perry, R. A.; Atkinson, R.; Pitts, J. N. Kinetics and mechanism of the gas phase reaction of hydroxyl radicals with aromatic hydrocarbons over the temperature range 296-473 K. J. Phys. Chem. 1977, 81 (4), 296-304. 53. Knispel, P.; Koch, R.; Siese, M.; Zetzsch, C. Adduct formation of OH radicals with benzene, toluene, and phenol and consecutive reactions of the adducts with NOx and O2. Ber. Bunsen-Ges. Phys. Chem. 1990, 94 (1375-1379). 54. Lin, S. C.; Kuo, T. C.; Lee, Y. P. Detailed rate coefficients and the enthalpy change of the equilibrium reaction OH+C6H6↔MHOC6H6 over the temperature range 345–385 K. J. Chem. Phys. 1994, 101 (3), 2098-2105. 55. Tully, F. P.; Ravishankara, A. R.; Thompson, R. L.; Nicovich, J. M.; Shah, R. C.; Kreutter, N. M.; Wine, P. H. Kinetics of the reactions of hydroxyl radical with benzene and toluene. J. Phys. Chem. 1981, 85 (15), 2262-2269. 56. Zheng, J.; Bao, J. L.; Meana-Pañeda, R.; Zhang, S.; Lynch, B. L.; Corchado, J. C.; Chuang, Y.-Y.; Fast, P. L.; Hu, W.-P.; Liu et al.. Polyrate 17-B: Computer Program for the Calculation of Chemical Reaction Rates for Polyatomics2017. 57. Le, T. H. M.; Do, S. T.; Huynh, L. K. Algorithm for auto-generation of hindered internal rotation parameters for complex chemical systems. Comput. Theor. Chem. 2017, 1100, 61-69. 58. Fernández-Ramos, A.; Ellingson, B. A.; Meana-Pañeda, R.; Marques, J. M. C.; Truhlar, D. G. Symmetry numbers and chemical reaction rates. Theor. Chem. Acc. 2007, 118 (4), 813-826. 59. Huber, K. P.; Herzberg, G. Molecular Spectra and Molecular Structures, IV. Constants of Diatomic Molecules; Van Nostrand Reinhold: New York 1979. 60. Manion, J. A.; Huie, R. E.; Levin, R. D.; Burgess Jr., D. R.; Orkin, V. L.; Tsang, W.; McGivern, W. S.; Hudgens, J. W.; Knyazev, V. D.; Atkinson et. al.. NIST Chemical Kinetics Database, NIST Standard Reference Database 17, Version 7.0 (Web Version), Release 1.6.8, Data version 2015.12, National Institute of Standards and Technology, Gaithersburg, Maryland, 20899-8320. Web address: http://kinetics.nist.gov/.

41 ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 42 of 43

61. Westmoreland, P. R.; Howard, J. B.; Longwell, J. P.; Dean, A. M. Prediction of rate constants for combustion and pyrolysis reactions by bimolecular QRRK. AIChE J. 1986, 32 (12), 1971-1979. 62. Van de Vijver, R.; Vandewiele, N. M.; Bhoorasingh, P. L.; Slakman, B. L.; Seyedzadeh Khanshan, F.; Carstensen, H.-H.; Reyniers, M.-F.; Marin, G. B.; West, R. H.; Van Geem, K. M. Automatic Mechanism and Kinetic Model Generation for Gas- and Solution-Phase Processes: A Perspective on Best Practices, Recent Advances, and Future Challenges. Int. J. Chem. Kinet. 2015, 47 (4), 199-231. 63. Kim, Y.; Kim, W. Y. Universal Structure Conversion Method for Organic Molecules: From Atomic Connectivity to Three-Dimensional Geometry. Bull. Kor. Chem. Soc. 2015, 36 (7), 1769-1777. 64. Rangarajan, S.; Kaminski, T.; Van Wyk, E.; Bhan, A.; Daoutidis, P. Language-Oriented Rule-Based Reaction Network Generation and Analysis: Algorithms of RING. Comp. Chem. Eng. 2014, 64, 124-137. 65. Truong, T. N.; Duncan, W.; Tirtowidjojo, M. A reaction class approach for modeling gas phase reaction rates. Phys. Chem. Chem. Phys. 1999, 1 (6), 1061-1065. 66. Kilpatrick, J. E.; Pitzer, K. S. Energy Levels and Thermodynamic Functions for Molecules with Internal Rotation. III. Compound Rotation. J. Chem. Phys. 1949, 17 (11), 1064. 67. Carstensen, H.-H.; Dean, A. M. Rate Constant Rules for the Automated Generation of Gas-Phase Reaction Mechanisms. J. Phys. Chem. A 2009, 113 (2), 367-380.

42 ACS Paragon Plus Environment

Page 43 of 43 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47

The Journal of Physical Chemistry

ACS Paragon Plus Environment