Silane-initiated nucleation in chemically active plasmas: validation of density functionals, mechanisms, and pressure-dependent variational transition state calculations

Abstract

The growth of anionic silicon hydride clusters is a critically important process in nanodusty plasmas. In the current study, we focus on the formation of homologs of silylene (Si_n+1H_2n+2⁻, n = 3, 4) and silyl (Si_nH_2n+1⁻, n = 4, 5) anions via anion–neutral reaction pathways. Species like silyl or silylene anions and their related elementary reactions, which are involved in the formation of silicon hydride clusters, were not used in developing exchange–correlation (xc) density functionals (i.e., they were not included in the training set of semiempirical density functionals); therefore, we explored the accuracy of various widely used xc density functionals based on reaction energies and barrier heights. Among the 21 density functionals we tested, M06-2X has the best performance for a hybrid functional, and MN15-L has the best performance for a local functional. Thermal rate constants of the elementary reactions involved in the reaction mechanism are calculated using M06-2X and multistructural canonical variational transition state theory with the small-curvature tunneling approximation (MS-CVT/SCT). The pressure dependence of unimolecular isomerization reactions is treated with system-specific quantum RRK theory (SS-QRRK) and the Lindemann–Hinshelwood mechanism.

---

1. Introduction

The growth of nanoparticles in nanodusty plasmas is an active research field in plasma physics, chemistry, and engineering. Many physical and chemical processes are involved in the formation of nanoparticles in chemically active plasmas, including nucleation, isomerization, electron capture and ionization, and mass and momentum transport. Anion–neutral reactions are one of the major chemical reactions in silicon hydride clustering in silane-containing reactive plasmas.^1,2 Rate constants of anion–neutral reactions are used in building transport equations or boundary conditions in many simulation studies for investigating the distribution of the sizes of nanoparticles and their population in a plasma.^3–7 However, due to the complexity of the systems and the unusual conditions of the chemical reactions (i.e., plasma), accurate experimental measurements for these anion–neutral reactions are very difficult and hence are scarce. The majority of the reaction rates needed for modeling nanodusty plasmas are empirically estimated, and Agarwal and Girshick⁸ noted that “there is considerable uncertainty in the values of rate constants for anion–neutral reactions that are primarily responsible for clustering in this system, and thus any correction factor one might apply for the predicted nucleation rate would itself be highly uncertain”. Theoretical calculations of chemical structures, energetics, and thermal rate constants can play an important role in reducing the uncertainty involved in plasma modeling. Although fast empirical methods for estimating thermodynamic functions and rate constants exist,⁹ we do not use such approach in the current work; thermodynamic functions and thermal rate constants reported in this work are computed from ab initio calculations.

Silyl anion (Si_nH_2n+1⁻) reactions with silane and silylene anion (Si_n+1H_2n+2⁻) reactions with silane are the dominant anionic pathways for the formation of nanoparticles in plasmas. Higher homologs of silicon hydrides with branched chains are generated in these reactions, which proceed with the elimination of molecular hydrogen. The silicon hydrides formed via silyl anion–silane reactions are H₃SiH₂Si:⁻ → (H₃Si)₂HSi:⁻ → (H₃Si)₃Si:⁻ → (H₃Si)₃SiH₂Si:⁻, where “:” represents paired electrons on the terminal Si atom; the silicon hydrides formed via silylene anion–silane reactions are H₃SiHSi˙⁻ → H₃SiH₂SiHSi˙⁻ → (H₃Si)₂HSiHSi˙⁻ → (H₃Si)₃SiHSi˙⁻, where “˙” represents an unpaired electron on the terminal Si atom.

In the present study, we focus on clusters that contain no more than five silicon atoms. The chemical mechanism of the initial step of the clustering process, i.e., the formation of (H₃Si)₂HSi:⁻ and H₃SiH₂SiHSi˙⁻, has been investigated in a previous work.¹⁰ In the current work, since the studied kind of system is not represented in the training set of any density functional known to us, we carry out benchmark tests for various density functionals, and this provides an opportunity to investigate the transferability of various semiempirical exchange–correlation functionals. Chemical structures, energetics, and thermal rate constants in the mechanisms are computed using the best functional among the functionals we tested in this work.

Two unimolecular isomerization reactions are involved in the reaction mechanism we propose in this work, and the pressure dependences of their thermal rate constants are estimated using system-specific quantum Rice–Ramsperger–Kassel (SS-QRRK) theory combined with the Lindemann–Hinshelwood thermal activation mechanism. The SS-QRRK method was proposed recently¹¹ as a convenient way to use variational transition state theory to treat pressure dependences of chemical reaction rates, and it was applied to a chemical activation mechanism; the present article is the first application to a thermally activated unimolecular reaction. For a chemically activated unimolecular reaction, the low-pressure rate constants of the formation of the stabilized adduct are lower than the high-pressure-limit, and the rate constants of the further isomerization/dissociation of the formed adduct are larger than the high-pressure equilibrium rate constant; and in thermally activated unimolecular reactions, the low-pressure rate constants are smaller than the high-pressure equilibrium rate constant, and hence the pressure effect for this kind of reaction is called “falloff.”

2. Theoretical background

2.1. High-pressure-limit thermal rate constants

High-pressure-limit thermal rate constants are computed using multi-structural canonical variational transition state theory with the small-curvature tunneling approximation (MS-CVT/SCT) as follows:^12–15

k^MS-CVT/SCT = F_actκ^SCT₁Γ^CVT₁k^TST₁ (1)

where F_act is multi-structural torsional potential anharmonicity factor of activation^16,17 computed by MS-T(C) method;¹⁸κ^SCT₁ is a small-curvature tunneling transmission coefficient for the lowest-energy conformer of the transition state, Γ^CVT₁ is CVT variational transmission coefficient for the lowest-energy path, and k^TST₁ is the conventional transition state theory rate constant computed based on the lowest-energy path:where Q^‡ and Q^R are rigid-rotor–harmonic-oscillator partition functions for the transition state structure and reactant; Φ^R is the reactants partition function per unit volume. V^‡ is the barrier height, k_B is the Boltzmann constant, h is Planck's constant, and T is absolute temperature. Expressing the TST rate constant in the form of eqn (2), the rotational partition function does not contain rotational symmetry number; the rotational symmetry number is included in the reaction symmetry number σ_rxn as follows:¹⁹where σ^R and σ^‡ are the rotational symmetry number for the reactant and transition state structure, which are equal to the order of the rotational subgroup for polyatomic molecule; rotational symmetry number is 1 for heteronuclear diatomic molecule and 2 for homonuclear diatomic molecule. The reaction symmetry number we used here excludes the non-superimposable mirror-image conformers; the contributions from these mirror-image conformers are treated within MS-T method. However, the contributions from enantiomers should be included in the reaction symmetry number or the rotational partition function if the MS-T method is not applied.^19,20 The reaction symmetry numbers for forward and reverse reactions involved in the current study are tabulated in the ESI.†

2.2. System-specific quantum RRK theory with Lindemann–Hinshelwood theory

In this section, we consider the effect of pressure on the unimolecular reactions by using SS-QRRK theory¹¹ with the Lindemann–Hinshelwood mechanism.

For a unimolecular reaction, the following Lindemann–Hinshelwood thermal activation mechanism is considered:²¹

where A(T) represents the thermally equilibrated reactant at temperature T; M is the bath gas; A*(E) is the rovibrationally excited molecule with total rovibrational energy E; and P is the isomerization product. This mechanism includes the RRK assumption that energy in A* is rapidly statistically redistributed among modes subject only to the constraint of total energy E so that the reactivity of A* is simply a function of total energy. The rate constant of energization is k₁, and the rate constant of isomerization is k₂, and by the RRK assumption both of these rate constants depend on the total energy of A*; and k₁ also depends parametrically on temperature T. The rate constant of de-energization is k_c and is treated as temperature-dependent but energy-independent, which is the strong collision assumption. However the strong-collision assumption is mitigated in the present work by computing k_c as the product of the Lennard-Jones collision rate constant and a collision efficiency factor β_c, where the latter is computed by using Troe's modified collision model.^22,23

The pressure-dependent unimolecular reaction rate constant k_uni(T, p) for the above mechanism is:^24–26

where K(E,T) is the equilibrium constant of the first step, [M] is the concentration of bath gas, p is pressure, and E₀ is the threshold energy. Note that K(E,T) is mixed-ensemble equilibrium constant, representing the thermal equilibrium of species A* in a microcanonical ensemble at energy E with species A in a thermal ensemble with temperature T. Therefore,where ρ(E) is the density of states, and Q^A(T) is the rovibrational partition function of A.

In QRRK theory,^24,27 one models the states of A and A* as the discrete states of a system with s uncoupled harmonic oscillator modes, each with frequency [small omega, Greek, macron] in wave numbers (e.g., cm⁻¹). Then the integrals over E are replaced by sums over n, where n is the number of quanta excited at energy E (n = E/hc[small omega, Greek, macron], where the zero of energy is the potential energy of the equilibrium structure of A). With this model and with [M] given by the ideal gas law, eqn (4) can be written as:

where

m = E₀/hc[small omega, Greek, macron] (8)

and R is ideal gas constant, and c is the speed of light. Notice that in the limit of p → ∞, the high-pressure-limit unimolecular rate constant is recovered in eqn (4) and (6), and in the low-pressure limit, the unimolecular rate “constant” is no longer a constant (i.e., no longer independent of concentrations), but rather is proportional to [M] or p.

In SS-QRRK, we use eqn (6)–(8) with [small omega, Greek, macron] taken as the geometric mean of the s vibrational frequencies of A (where s is 3N − 6, and N is the number of atoms in A), and with E₀ and the microcanonical isomerization rate constant k₂(E) parameterized by using QRRK theory calibrated to match the MS-CVT/SCT canonical rate constant at high pressure. In QRRK theory, the microcanonical rate constant is a frequency factor A (with units of reciprocal time) times the probability that a molecule with n quanta of vibrational theory has at least m quanta in one of the modes (the reactive mode), where m is given in terms of the threshold energy by eqn (8). This yields²⁴

Substituting this into eqn (6), carrying out the sum, and taking the limit of p → ∞ yields the Arrhenius form:^19,24

k^QRRK_uni(T, p = ∞) = A exp(−mhc[small omega, Greek, macron]/RT) (10)

Thus, to parameterize QRRK theory, the parameter m is calculated from the MS-CVT/SCT Arrhenius activation energy E^MS-CVT/SCT_a(T):and the frequency factor A in QRRK theory is set equal the MS-CVT/SCT Arrhenius pre-exponential factor, i.e.,

A(T) = k^MS-CVT/SCT(T) exp[E^MS-CVT/SCT_a(T)/RT] (13)

Note that both E₀ and A depend on temperature in the parameterized rate expression. Therefore k₂(E) in eqn (6) becomes k₂(E,T) given by eqn (9) and (11)–(13). The sum in eqn (6) is evaluated with a step size of one quantum, and the factorials therefore all have non-integer arguments and are evaluated using gamma functions.

The de-energization rate constant is modeled as in our previous work,¹¹ using empirical Lennard-Jones parameters, the average energy transferred, and the energy dependence factor F_E of the density of states (which is the thermal population of unimolecular states above the threshold energy of the reactant normalized by a density of states factor at the threshold energy); these quantities are given in Section 3.4.

3. Computational details

3.1. Electronic structure calculations

Initial geometries of the reactants, products, and transition state structures are optimized with the M08-HX functional²⁸ and the MG3S basis²⁹ for benchmark studies; tight convergence criteria are used for both the SCF calculations and the geometry optimizations. For silicon atoms, the MG3S basis is equivalent to the 6-311+G(3d2f) basis, and for hydrogen atoms it is equivalent to the 311G(2p) basis.^30–32 All the density functional integrations are carried out with a grid of 99 radial shells around each atom and 974 angular points per shell.³³ All the electronic structure calculations are performed with a locally modified version of the Gaussian 09 software.^34,35

In order to test the accuracy of various popular density functionals for silyl or silylene anions reactions, a benchmark study was carried out based on classical barrier heights and classical energies of reactions (these are relative Born–Oppenheimer potential energies at stationary points and are exclusive of zero-point or thermal vibrational energy). Five reactions (as summarized in Table 1) were selected from our proposed reaction mechanisms for use in the benchmark study.

Table 1 Reactions used in benchmark study

	Chemical equation	Reaction type
R1	SiH3SiHSiH2^− + SiH4 → (SiH3)2SiHSiH2^− + H	Nucleophilic reaction
R2	(SiH3)2SiHSiH2^− + H → (SiH3)2SiHSiH^− + H2	Hydrogen (H) abstraction
R3	(SiH3)2SiH^− + SiH4 → (SiH3)2SiH2 + SiH3^−	Hydrogen (H) abstraction
R4	(SiH3)2SiH2 + SiH3^− → (SiH3)2SiHSiH2^− + H2	Hydrogen (H2) elimination
R5	(SiH3)2SiHSiH2^− → (SiH3)2SiSiH3^−	Intramolecular hydrogen migration

---

Single-point energy calculations are performed based on M08-HX/MG3S geometries using various density functionals (tabulated in Table 2) combined with the MG3S, jun-cc-pVTZ³⁶ and jul-cc-pVTZ³⁷ basis sets. Density functionals tested in the current work includes three local functionals (M11-L,³⁸ MN12-L,³⁹ and MN15-L⁴⁰) and 18 hybrid functionals (B3LYP,⁴¹ PBE0,⁴² TPSSh,⁴³ MGGA_MS2h,⁴⁴ MPW1K,⁴⁵ M05,⁴⁶ M06,⁴⁷ M05-2X,⁴⁸ M06-2X,⁴⁷ M08-HX,²⁸ M08-SO,²⁸ M11,⁴⁹ MN12-SX,⁵⁰ SOGGA11-X,⁵¹ B97-3,⁵² HSE06,⁵³ τHCTHhyb,⁵⁴ and ωB97X-D⁵⁵).

Table 2 Exchange–correlation functionals tested in the current work and their percentage of non-local Hartree–Fock exchange (% X)

Functional	% X	xcF	% X	xcF	% X
a The percentage of Hartree–Fock exchange increases from the first value listed for small interelectronic separation to 100% at large interelectronic separation. b The percentage of Hartree–Fock exchange decreases from 25% at small interelectronic separation to 0 at large interelectronic separation.
M11-L	0	ωB97X-D	22.2–100^a	SOGGA11-X	35.42
MN12-L	0	HSE06	25–0^b	MPW1K	42.8
MN15-L	0	MN12-SX	25–0^b	M11	42.8–100^a
MGGAMS2h	9	PBE0	25	M08-HX	52.23
TPSSh	10	B97-3	26.93	M06-2X	54
τHCTHhyb	15	M06	27	M05-2X	56
B3LYP	20	M05	28	M08-SO	56.79

---

The reference energy values are computed by CCSD(T)⁵⁶/CBS, where the complete basis set (CBS) limit is obtained by the following strategy:⁵⁷

E^CCSD(T)_CBS = E^MP2_CBS + (E^CCSD(T)_SBS − E^MP2_SBS) (14)

where SBS means the small basis set (which is aug-cc-pVTZ^37,58 in the current work), and the MP2⁵⁹/CBS energy is computed by respectively extrapolating the Hartree–Fock (HF) exchange energy and MP2 correlation energy as follows:^60–64where In the above equations, we use X = 4 (which is taken as aug-cc-pVQZ) and X − 1 = 3 (which is taken as aug-cc-pVTZ). Finally,

E^MP2_CBS = E^HF_CBS + E^corr_CBS (18)

3.2. Direct dynamics calculations

For the direct dynamics calculations, all the species are re-optimized by M06-2X/MG3S, which was selected as the level for direct dynamics calculations based on the benchmark tests (as will be discussed in Section 4.1).

Canonical variational transition state theory calculations with the small-curvature tunneling approximation were carried out in non-redundant internal coordinates^65,66 with a step size of 0.002 a₀ (note: 1 a₀ = 1 bohr = 0.5292 Å). The minimum energy paths (MEPs) are computed using Page–McIver algorithm⁶⁷ from −2.0 to +2.0 a₀. The RODS algorithm⁶⁸ was used to re-orient the generalized-transition-state-theory dividing surface. A scaling factor 0.970⁶⁹ was used to scale all the vibrational frequencies in the generalized normal mode calculations.

For bimolecular reactions with a negative barrier, the small-curvature tunneling transmission coefficient was computed using the ion–dipole complex as the initial state; the final bimolecular reaction rate constant is the product of the so-obtained tunneling transmission coefficient with the bimolecular reaction rate constant computed without tunneling. This is consistent with the fact that the tunneling calculation is performed for the high-pressure limit where the ion–dipole complex is fully thermalized (the issue of tunneling from the states of the precursor complex at lower energies than the bimolecular reactant ground state is discussed elsewhere^70,71). Multi-structural torsional anharmonicity (MS-T) rovibrational partition functions were computed using the MSTor program;⁷² VTST calculations were performed with the Polyrate⁷³ and Gaussrate⁷⁴ programs.

3.3. Torsional anharmonicity

The multi-structural torsional anharmonicity (MS-T) rovibrational partition functions are computed based on coupled effective torsional potentials. The local periodicities of –SiH₃ groups are set to be 3. The MS-T partition functions include the contributions from all the distinguishable conformational structures including non-superimposable mirror images.

3.4. Parameters used in de-energization

In the current study, Ar gas, which is commonly used in chemical vapor decomposition⁷⁵ (CVD) and in studying nucleation processes in plasmas,⁸ is selected as the bath gas used in estimation of the pressure-dependent rate constants of the unimolecular isomerization reactions. Lennard-Jones parameters ε/k_B and σ are taken for Ar as 120 K and 3.4 Å⁷⁶ and for Si₄H_n species as 254 K and 5.8 Å, as used in previous silane plasma modeling.⁷⁷

The energy transfer parameter 〈ΔE〉, which is the average energy transferred during both energization and de-energization processes and which is used for computing the collision efficiency factor, is chosen to be 740 cal mol⁻¹; this value has been used previously for modeling SiH₄ colliding with Ar.⁷⁸ The Si clusters we considered in the current work are larger than SiH₄, and therefore one might hypothesize that a larger energy transfer parameter 〈ΔE〉 should be used. To test this, we repeated the calculations for reaction RB, step 1 with a doubled energy transfer parameter of 1480 cal mol⁻¹, and we found that the obtained falloff curves (shown in Fig. S2 in ESI†) are not sensitive to this change. Doubling the energy transfer parameter leads to larger k(p) values (i.e., stronger collisions and smaller falloff effect) and at most a factor of 2 difference. (The maximum effect is at 1500 K, 0.001 bar.) The collision parameters and energy transfer parameters, in principle, could be determined from theoretical trajectory calculations;⁷⁹ we do not use such approach in the current work based on computational cost and the desire for a simple method that can be widely used in mechanism development.

The energy dependence factor F_E of the density of states can be directly computed using Troe's definition:²²

with density of states computed from MS-T partition function by inverse Laplace transform; alternatively, F_E can be computed by the empirical Whitten–Rabinovitch approximation using eqn (6) and (8) in Troe's work.⁸⁰F_E is used to compute the collision efficiency coefficient β_c using the following equation:

For falloff calculations on (SiH₃)₂SiHSiH⁻, we used F_E values computed from MS-T partition function in eqn (19) which yields 1.55, 1.84, and 2.85 at 298 K, 400 K, and 600 K, respectively, for comparison, F_E values computed using Whitten–Rabinovitch approximation at 298 K, 400 K, and 600 K are 1.43, 1.63, and 2.30, which agree with the values we used within 8, 11, and 19%, respectively. The F_E values computed from MS-T partition function for (SiH₃)₂SiHSiH₂⁻ are 1.60, 1.97, and 3.32 at 298 K, 400 K, and 600 K, respectively, and they are used in falloff calculations of RC step 3. The corresponding F_E values computed by the Whitten–Rabinovitch approximation at 298 K, 400 K, and 600 K are 1.45, 1.71, and 2.51, which agree with the values we used within 9, 13, and 24%, respectively.

4. Results and discussion

4.1. Benchmark of various density functionals

Although coupled cluster theory with single and double excitations and quasi-perturbative connected triple excitations, i.e., CCSD(T), is often viewed as the gold standard in quantum chemistry, it uses a single configuration state function as the reference wave function (i.e., it is a “single-reference” method), and therefore it might not be appropriate to use CCSD(T) as a benchmark for systems with strong multi-reference characters (such as bond dissociation⁸¹ and transition metal chemistry⁸²). To ascertain its expected reliability, we computed the T₁ diagnostic⁸³ values for all the species involved in the reactions that are used for benchmark study, since this diagnostic has been proposed as a measure of the suitability of singe-reference coupled cluster theory. The molecule with the largest T₁ diagnostic value is the transition state structure of reaction (R1) (Si₄H₁₀⁻), for which the value is 0.0252 computed by CCSD(T)/aug-cc-pVTZ; this value is much smaller than 0.045, which has been suggested^84–86 as a criterion for the applicability of single-reference methods to open-shell molecules, and hence we concluded that CCSD(T) can serve as a reference for testing other methods in this work.

The mean unsigned errors (MUEs) of various density functionals computed with respect to CCSD(T)/CBS for reactions (R1)–(R5) and the overall MUEs are shown in Fig. 1. MUEs are calculated based on zero-point-vibrational-energy-exclusive forward barrier height and energy of reaction. MUEs shown in Fig. 1 are computed using the MG3S basis; using larger basis sets such as jun-cc-pVTZ and jul-cc-pVTZ has a negligible effect on the MUEs. All the computed classical forward barrier heights, energies of reaction, and MUEs are tabulated in the ESI.†

Fig. 1 Mean unsigned errors (MUEs, in kcal mol⁻¹) for reactions (R1)–(R5) and the overall MUEs of various density functionals. MUEs shown in this figure (computed with the MG3S basis) are based on classical forward barrier and energy of reaction; CCSD(T)/CBS values are used as references. The functionals are in order of increasing overall MUE.

Among the 25 tested density functionals, the global-hybrid meta-GGA functional M06-2X performs the best for all the reactions; the MUEs for reactions (R1)–(R5) and the overall MUE are 0.34, 0.35, 0.22, 0.76, 0.95 and 0.52 kcal mol⁻¹, respectively. They are all better than the value of 1 kcal mol⁻¹, which is usually quoted as “chemical accuracy”. M06-2X, M05-2X, M08-SO, and MPW1K are the four best performing functionals and all have overall MUEs within the chemical accuracy criterion. Among the three tested local density functionals (M11-L, MN12-L and MN15-L), MN15-L possesses the smallest overall MUE, which is 1.25 kcal mol⁻¹; it is encouraging that this surpasses the accuracy of many popular hybrid density functionals (such as B3LYP and ωB97X-D) and the recently developed MGGA_MS2h. This is encouraging for the very new MN15-L functional since there is no silyl or silylene-related species in its training set and since local functionals are usually less accurate than hybrid functionals for reaction energies and barrier heights.

Note that the density functional calculations converge more rapidly with respect to basis set than CCSD(T), and so CBS extrapolations were not needed. A great advantage of density functional theory in this respect is that the reaction path calculations and the required Hessians along the reaction path are affordable even with the reasonably well-converged MG3S basis set.

4.2. Proposed reaction mechanisms

In this work, we consider the silicon hydride clusters formed from silylene or silyl anions reactions with silane that involve no more than 5 silicon atoms. To be more specific, the composite reactions we studied are:

(RA): SiH₃SiHSiH₂⁻ + SiH₄ → (SiH₃)₂SiHSiH⁻ + H₂

(RB): (SiH₃)₂SiHSiH⁻ + SiH₄ → (SiH₃)₃SiSiH⁻ + H₂

(RC): (SiH₃)₂SiH⁻ + SiH₄ → (SiH₃)₃Si⁻ + H₂

(RD): (SiH₃)₃Si⁻ + SiH₄ → (SiH₃)₃SiSiH₂⁻ + H₂

In the above reactions, the ground-state spin multiplicities for all the silylene anions (i.e., SiH₃SiHSiH₂⁻, (SiH₃)₂SiHSiH⁻ and (SiH₃)₃SiSiH⁻) are doublet and for silyl anions (i.e., (SiH₃)₂SiH⁻, (SiH₃)₃Si⁻ and (SiH₃)₃SiSiH₂⁻) are singlet. The growth of (SiH₃)₃SiSiH⁻ starts from SiH₃SiHSiH₂⁻, which is produced during the initial polymerization reaction of SiH₄ + Si₂H₄⁻; (SiH₃)₂SiH⁻, which leads to the formation of (SiH₃)₃SiSiH₂⁻, is generated by the reaction SiH₄ + Si₂H₅⁻.

Potential energy diagrams (relative Born–Oppenheimer potential energy E in kcal mol⁻¹ with respect to reactants) for reaction mechanisms of reactions RA–RD are shown in Fig. 2–5. The elementary steps involved and their classical energies of reaction (ΔE) are listed in Table 3. Energetic values are computed based on geometries optimized at M06-2X/MG3S and single point energies calculated with the same method. All the elementary reactions are bimolecular reactions except for the first step of RB and the third step of RC, which are unimolecular isomerization reactions (intramolecular hydrogen transfer). We considered ion–dipole complex in four of the elementary bimolecular reactions, in which the barrier is negative: RA step 2, RB step 3, RC step 1, and RD step 1; the ion–dipole complex is used in computing tunneling transmission coefficient of the bimolecular reactions.

Fig. 2 Potential energy diagram for reaction mechanism of reaction RA; relative energies (E in kcal mol⁻¹, with respect to reactants SiH₃SiHSiH₂⁻ + SiH₄) are computed at M06-2X/MG3S level.

Fig. 3 Potential energy diagram for reaction mechanism of reaction RB; relative energies (E in kcal mol⁻¹, with respect to reactants (SiH₃)₂SiHSiH⁻ + SiH₄) are computed at M06-2X/MG3S level.

Fig. 4 Potential energy diagram for reaction mechanism of reaction RC; relative energies (E in kcal mol⁻¹, with respect to reactants (SiH₃)₂SiH⁻ + SiH₄) are computed at M06-2X/MG3S level.

Fig. 5 Potential energy diagram for reaction mechanism of reaction RD; relative energies (E in kcal mol⁻¹, with respect to reactants (SiH₃)₃Si⁻ + SiH₄) are computed at M06-2X/MG3S level.

Table 3 Elementary steps in the proposed reaction mechanisms for reactions RA–RD, and their classical energies of reactions (kcal mol⁻¹) at M06-2X/MG3S level

Reaction	Step	Chemical equation	ΔE (kcal mol^−1)
A	1	SiH3SiHSiH2^− + SiH4 → (SiH3)2SiHSiH2^− + H	27.1
A	2	(SiH3)2SiHSiH2^− + H → (SiH3)2SiHSiH^− + H2	−31.3
B	1	(SiH3)2SiHSiH^− → (SiH3)2SiH2Si^−	−3.0
B	2	(SiH3)2SiH2Si^− + SiH4 → (SiH3)3SiSiH2^− + H	31.6
B	3	(SiH3)3SiSiH2^− + H → (SiH3)3SiSiH^− + H2	−31.2
C	1	(SiH3)2SiH^− + SiH4 → (SiH3)2SiH2 + SiH3^−	24.7
C	2	(SiH3)2SiH2 + SiH3^− → (SiH3)2SiHSiH2^− + H2	−20.3
C	3	(SiH3)2SiHSiH2^− → (SiH3)2SiSiH3^−	−12.6
D	1	(SiH3)3Si^− + SiH4 → (SiH3)3SiH + SiH3^−	34.5
D	2	(SiH3)3SiH + SiH3^− → (SiH3)3SiSiH2^− + H2	−24.6

---

4.3. Multi-structural and torsional potential anharmonicity

Multi-structural and torsional potential anharmonicity factors of activation for the forward reactions are tabulated in Table 4, and for the reverse reactions are in ESI.† MS-T standard-state reaction enthalpy ( in kcal mol⁻¹) and reaction Gibbs free energy ( in kcal mol⁻¹) for all the reactions at 298 K and 1000 K are shown in Table 6, and their values at various temperatures are tabulated in ESI.† The numbers of distinguishable conformers found in the conformational searches for the species that have multiple structures are given in ESI.†

Table 4 MS-T factors for activation at various temperatures for forward reactions computed at M06-2X/MG3S level

T/K	RA		RB			RC			RD
	Step 1	Step 2	Step 1	Step 2	Step 3	Step 1	Step 2	Step 3	Step 1	Step 2
298	1.41	0.57	0.93	1.91	1.68	0.52	1.64	0.75	0.70	0.52
300	1.41	0.57	0.93	1.90	1.68	0.52	1.64	0.75	0.70	0.52
400	1.50	0.61	0.92	1.70	1.69	0.47	1.65	0.82	0.66	0.57
500	1.55	0.66	0.96	1.52	1.72	0.42	1.68	0.91	0.62	0.61
600	1.56	0.70	1.02	1.36	1.74	0.38	1.72	1.02	0.58	0.66
700	1.55	0.75	1.09	1.24	1.77	0.35	1.76	1.13	0.55	0.70
800	1.53	0.80	1.17	1.13	1.80	0.32	1.82	1.25	0.52	0.74
900	1.50	0.85	1.25	1.05	1.82	0.29	1.87	1.37	0.49	0.78
1000	1.47	0.90	1.33	0.98	1.84	0.27	1.92	1.48	0.46	0.82
1500	1.31	1.13	1.71	0.74	1.95	0.20	2.17	2.03	0.36	1.00

---

The errors for computing the thermal rate constants introduced by ignoring the multi-structural and torsional anharmonicity effects vary between different reactions and temperatures. For reactions RA step 1, RB step 3 and RC step 2, the errors are only slightly temperature dependent; the averaged errors (over 298–1500 K temperature range) for thermal rate constants computed without including MS-T effects for these reactions are respectively 32%, 43%, and 44%, with a standard deviation of 4%, 3%, and 5% between various temperatures. For reaction RB step 1, the MS-T effect can be ignored for temperatures below 800 K, at which temperatures the errors are all smaller than 15%; at 1500 K, one would have an error of 42% in thermal rate constant if the MS-T factor were not included. For reaction RC step 1 at 1500 K, the multi-structural CVT/SCT rate constant is a factor of 5 smaller than the single-structural CVT/SCT rate constant; the ratio of the single-structural harmonic oscillator rovibrational partition function (SS-HO) to the MS-T rovibrational partition function is 0.138 for the TS and is 0.696 for the reactant (SiH₃)₂SiH⁻, while the ratio of MS-HO partition function to SS-HO partition function is 2 for the TS and is 1 for (SiH₃)₂SiH⁻, which indicates stronger couplings between torsional modes in the TS than in (SiH₃)₂SiH⁻.

4.4. High-pressure-limit thermal rate constants

Calculated MS-CVT/SCT rate constants for all the forward reactions are listed in Table 5, and for reverse reactions in the ESI.† The following equations are used to fit the thermal rate constants:^70,87where A, n, E and T₀ are fitting parameters, T is temperature in Kelvin, and R is the ideal gas constant (1.9872 × 10⁻³ kcal mol⁻¹ K⁻¹). Fitting parameters for forward reactions are shown in Table 6, and those for the reverse reactions are in the ESI.†

Table 5 MS-CVT/SCT rate constants for forward reactions computed at M06-2X/MG3S level at various temperatures. For bimolecular reactions, units of rate constants are cm³ molecule⁻¹ s⁻¹; for unimolecular reactions, units are s⁻¹

T/K	RA		RB			RC			RD
	Step 1	Step 2	Step 1	Step 2	Step 3	Step 1	Step 2	Step 3	Step 1	Step 2
298	9.08 × 10^−32	3.71 × 10^−10	1.35 × 10^−7	2.23 × 10^−34	3.30 × 10^−10	6.55 × 10^−28	3.09 × 10^−21	3.57 × 10^−6	7.30 × 10^−41	2.01 × 10^−18
300	1.22 × 10^−31	3.55 × 10^−10	1.78 × 10^−7	3.14 × 10^−34	3.19 × 10^−10	8.27 × 10^−28	3.35 × 10^−21	4.64 × 10^−6	9.97 × 10^−41	2.09 × 10^−18
400	8.25 × 10^−27	8.27 × 10^−11	6.33 × 10^−3	1.26 × 10^−28	8.49 × 10^−11	6.27 × 10^−24	1.14 × 10^−19	9.35 × 10^−2	1.58 × 10^−35	1.20 × 10^−17
500	8.20 × 10^−24	4.09 × 10^−11	3.82 × 10^+0	3.32 × 10^−25	4.51 × 10^−11	1.59 × 10^−21	1.45 × 10^−18	3.89 × 10^+1	2.64 × 10^−32	5.00 × 10^−17
600	9.32 × 10^−22	2.76 × 10^−11	2.86 × 10^+2	7.15 × 10^−23	3.12 × 10^−11	6.90 × 10^−20	9.62 × 10^−18	2.19 × 10^+3	4.28 × 10^−30	1.57 × 10^−16
700	2.97 × 10^−20	2.18 × 10^−11	6.47 × 10^+3	3.60 × 10^−21	2.48 × 10^−11	1.08 × 10^−18	4.20 × 10^−17	4.15 × 10^+4	1.79 × 10^−28	4.03 × 10^−16
800	4.26 × 10^−19	1.89 × 10^−11	6.88 × 10^+4	7.17 × 10^−20	2.15 × 10^−11	8.80 × 10^−18	1.39 × 10^−16	3.85 × 10^+5	3.18 × 10^−27	8.91 × 10^−16
900	3.52 × 10^−18	1.74 × 10^−11	4.39 × 10^+5	7.66 × 10^−19	1.97 × 10^−11	4.61 × 10^−17	3.73 × 10^−16	2.21 × 10^+6	3.13 × 10^−26	1.76 × 10^−15
1000	1.98 × 10^−17	1.66 × 10^−11	1.96 × 10^+6	5.25 × 10^−18	1.84 × 10^−11	1.78 × 10^−16	8.69 × 10^−16	9.04 × 10^+6	2.04 × 10^−25	3.20 × 10^−15
1500	4.54 × 10^−15	1.71 × 10^−11	1.88 × 10^+8	2.15 × 10^−15	1.74 × 10^−11	1.16 × 10^−14	1.60 × 10^−14	7.18 × 10^+8	7.78 × 10^−23	2.83 × 10^−14

---

Table 6 Fitting parameters for MS-CVT/SCT rate constants for forward reactions computed by M06-2X/MG3S and MS-T reaction enthalpy and reaction Gibbs free energy at 298 K and 1000 K^a

Step	RA		RB			RC			RD
	Step 1	Step 2	Step 1	Step 2	Step 3	Step 1	Step 2	Step 3	Step 1	Step 2
a The standard-state pressure is one bar. For bimolecular (bimol) reactions, the units of parameter A are cm^3 molecule^−1 s^−1; for unimolecular (unimol) reactions, the unit is s^−1. The parameters T0 and E are in units of K and kcal mol^−1 respectively, and n is unitless; enthalpy and free energy are in units of kcal mol^−1. Reactions that are endothermic (endo) at 0 K are fit using eqn (21), and those that are exothermic (exo) at 0 K are fit using eqn (21).
	Endo	Exo	Exo	Endo	Exo	Endo	Exo	Exo	Endo	Exo
Molecularity	Bimol	Bimol	Unimol	Bimol	Bimol	Bimol	Bimol	Unimol	Bimol	Bimol
ln A	−33.986	−29.040	21.176	−33.713	−28.1662	−31.596	−38.010	24.442	−52.742	−38.301
n	4.7020	1.8197	2.7059	4.7496	1.3979	3.0955	4.8766	1.9478	5.4789	4.6543
T 0	110.23	50.00	107.81	99.92	42.27	114.19	186.33	10.00	104.93	255.65
E	18.406	−3.673	18.594	21.577	−3.2529	15.207	5.856	21.252	19.512	2.932
	24.4	−29.4	−2.9	28.5	−29.7	24.0	−21.2	−11.9	33.4	−25.2
	26.9	−30.3	−2.4	33.3	−33.2	24.3	−19.3	−12.2	33.7	−20.9
	25.7	−30.5	−2.1	29.6	−30.6	23.5	−21.0	−10.9	32.2	−24.8
	27.0	−32.7	−0.5	29.4	−31.5	21.9	−22.6	−8.2	29.1	−26.6

---

Fig. 6 shows the computed small-curvature tunneling (SCT) transmission coefficients at various temperatures for the third step in reaction RC. At low temperature, the thermal rate constant can be increased by an order of magnitude due to quantum tunneling. The SCT tunneling transmission coefficient decays rapidly as temperature increases; κ^SCT decreases from 11.1 at 298 K to 2.0 at 500 K, and to 1.1 at 1500 K.

Fig. 6 Small-curvature tunneling (SCT) transmission coefficient of reaction RC step 2 at various temperatures.

4.5. High-pressure-limit activation energy

The Arrhenius activation energies are defined byand are obtained by putting eqn (21) into eqn (22), which yields

Computed activation energies for both forward and reverse reactions are shown in ESI.† The temperature dependence of activation energies of forward reactions is depicted in Fig. 7.

Fig. 7 High-pressure-limit activation energies for forward reactions at various temperatures.

The temperature dependence of activation energy computed based on MS-CVT/SCT rate constants could be quite different from the one computed from conventional transition state theory (TST), because of the effects of tunneling, recrossing, and multi-structural torsional anharmonicity. For instance, in step 2 of reaction RC, the MS-CVT/SCT activation energy is 7.4 kcal mol⁻¹ at 298 K, and it increases by 12.5 kcal mol⁻¹ from 298 K to 1500 K; the conventional single-structural TST (without tunneling) activation energy at 298 K is 10.9 kcal mol⁻¹, which is 3.5 kcal mol⁻¹ higher than MS-CVT/SCT activation energy at 298 K. We also plot the SS-TST, SS-TST/SCT, MS-CVT, and MS-CVT/SCT rate constant for RC step 2 at various temperatures in Fig. 8. The SS-TST/SCT natural logarithm rate constant curve decreases slower than the straight-line-shape of SS-TST curve as 1000 K/T increases; at high temperature, tunneling is negligible and SS-TST overlaps with SS-TST/SCT curve, while at low temperature, SS-TST curve is significantly lower than SS-TST/SCT curve. The MS-CVT/SCT curve differs negligibly from the SS-CVT/SCT curve due to the small MS-T effect.

Fig. 8 The computed high-pressure-limit SS-TST, SS-TST/SCT, MS-CVT, and MS-CVT/SCT bimolecular thermal rate constants (cm³ molecule⁻¹ s⁻¹) for reaction RC step 2 at various temperatures (K).

4.6. Falloff effects for unimolecular isomerization reactions

Step 1 of reaction RB and step 3 of reaction RC are unimolecular isomerization reactions. The predicted falloff curves for these two reactions at various temperatures (in K) and pressures (in bar) are shown in Fig. 9. Falloff curves are plotted as log₁₀[k(p)/k(p = ∞)] versus pressure, where k(p) is the thermal unimolecular rate constant at pressure p, and k(p = ∞) is the high-pressure-limit rate constant computed by MS-CVT/SCT theory. The falloff curves for pressures from 10⁻⁵ to 1 bar are shown in ESI.†

Fig. 9 Predicted falloff curves for thermal rate constants of unimolecular isomerization reactions at various temperatures (K) and pressures (bar). Solid lines are for reaction RB step 1; dashed lines are for reaction RC step 3.

At low and middle temperatures (T < 600 K), falloff effects for these two unimolecular reactions are negligible. At 600 K, the rate constant of RB step 1 at 1000 bar is 2.83 × 10² s⁻¹, and at 0.01 bar it is 2.66 × 10² s⁻¹; the rate constant of RC step 3 at 600 K and 1000 bar is 2.19 × 10³ s⁻¹, while at the same temperature at 0.01 bar it is 1.90 × 10³ s⁻¹. For RB step 1, the collision efficiency coefficient β_c at 298 K, 600 K, and 1000 K are 0.34, 0.14 and 0.03, respectively; and for RC step 3, collision efficiency coefficients at 298 K, 600 K, and 1000 K are 0.33, 0.12 and 0.02, respectively.

At high temperature, falloff effects become more significant. At 1500 K, the rate constant of RB step 1 is 1.45 × 10⁸ s⁻¹ at 1000 bar and is 2.66 × 10⁶ s⁻¹ at 1.0 bar, so the 1 bar result is a factor of 0.014 smaller than the high-pressure-limit. The rate constant of RC step 3 at 1500 K is 2.82 × 10⁸ s⁻¹ at 1000 bar and is 1.25 × 10⁶ s⁻¹ at 1 bar, so the high-pressure-limit rate constant is 574 times larger than the one at 1.0 bar.

4.7. Pressure-dependent activation energy

The activation energies for the unimolecular isomerization reactions RB step 1 and RC step 3 depend on both temperature and pressure. We show the pressure-dependent activation energies at 600 K, 800 K, and 1000 K in Fig. 10. At the highest pressures shown (10³ bar), the activation energies are at the high-pressure limit. At 600 K and around 0.1 bar, where falloff effects on the rate become notable, the activation energies also start falling as the pressures decreases. At 800 K and 1000 K, the activation energy is significantly pressure dependent; the activation energy decreases almost linearly with respect to the pressure. At the transition pressure p_1/2,²⁵ which is defined as the pressure at which the unimolecular rate constant is half of the high-pressure-limit value, the ratio of E_a(p = p_1/2) to E_a(p = ∞) is 0.6 for RB step 1 (p_1/2 = 0.01 bar) and 0.6 for RC step 3 (p_1/2 = 0.03 bar) at 800 K, and it is 0.5 for RB step 1 (p_1/2 = 0.3 bar) and 0.4 for RC step 3 (p_1/2 = 1 bar) at 1000 K.

Fig. 10 Predicted pressure-dependent activation energies (kcal mol⁻¹) for unimolecular isomerization reactions RB step 1 and RC step 3 at 600 K, 800 K, and 1000 K. Solid lines are for reaction RB step 1; dashed lines are for reaction RC step 3.

5. Summary

In the current work, we tested various exchange–correlation density functionals for an important system in nanodusty plasmas, in particular silylene and silyl anions reacting with silane molecules. Among the functionals we tested in this work, M06-2X is the most successful hybrid functional and MN15-L is the most successful local functional. Reaction mechanisms for the growth of silicon hydride clusters have been proposed, and the thermal rate constants of the elementary reactions involved in the reaction mechanisms were computed using multi-structural canonical variational transition state theory with the small-curvature tunneling approximation. Two unimolecular isomerization reactions are involved in the reaction mechanism, and their pressure dependent thermal rate constants were estimated based on system-specific quantum RRK (SS-QRRK) theory combined with Lindemann–Hinshelwood theory.

This work provides guidance for choosing density functionals for studying anion–neutral reactions in the silane-based reactive plasma. The methodology for computing thermal rate constants, whose values are rarely available experimentally, is also useful in estimating input kinetic data in plasma modeling for engineering applications. The extension of the SS-QRRK method to unimolecular isomerizations is also of more general use; it may be applied, for example, to atmospheric chemistry and combustion chemistry.

Acknowledgements

The authors appreciate helpful discussion with Steven L. Girshick and Mark J. Kushner and valuable contributions from Prasenjit Seal. This work is supported by the National Science Foundation under award no. CHE11-24752.

References

1. U. V. Bhandarkar, M. T. Swihart, S. L. Girshick and U. R. Kortshagen, J. Phys. D: Appl. Phys., 2000, 33, 2731.
2. U. V. Bhandarkar, U. R. Kortshagen and S. L. Girshick, J. Phys. D: Appl. Phys., 2003, 36, 1399.
3. S. J. Warthesen and S. L. Girshick, Plasma Chem. Plasma Process., 2007, 27, 292.
4. P. Agarwal and S. L. Girshick, Plasma Sources Sci. Technol., 2012, 21, 055023.
5. A. Gallagher, Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top., 2000, 62, 2690.
6. L. Ravi and S. L. Girshick, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2009, 79, 026408.
7. S. Nunomura, I. Yoshida and M. Kondo, Appl. Phys. Lett., 2009, 94, 071502.
8. P. Agarwal and S. L. Girshick, Plasma Chem. Plasma Process., 2014, 34, 489.
9. R. V. de Vijver, N. M. Vandewiele, P. L. Bhoorasingh, B. L. Slakman, F. S. Khanshan, H.-H. Carstensen, M.-F. Reyniers, G. B. Marin, R. H. West and K. M. Van Geem, Int. J. Chem. Kinet., 2015, 47, 199.
10. J. L. Bao, P. Seal and D. G. Truhlar, Phys. Chem. Chem. Phys., 2015, 17, 15928.
11. J. L. Bao, J. Zheng and D. G. Truhlar, J. Am. Chem. Soc., 2016, 138, 2690.
12. T. Yu, J. Zheng and D. G. Truhlar, Chem. Sci., 2011, 2, 2199.
13. B. C. Garrett and D. G. Truhlar, J. Chem. Phys., 1979, 70, 1593.
14. B. C. Garrett and D. G. Truhlar, Acc. Chem. Res., 1980, 13, 440.
15. D. G. Truhlar, A. D. Isaacson and B. C. Garrett, in Theory of Chemical Reaction Dynamics, ed. M. Baer, CRC Press, Boca Raton, 1985, pp. 65–137.
16. J. L. Bao, R. Meana-Pañeda and D. G. Truhlar, Chem. Sci., 2015, 6, 5866.
17. J. L. Bao, P. Sripa and D. G. Truhlar, Phys. Chem. Chem. Phys., 2016, 18, 1032.
18. J. Zheng and D. G. Truhlar, J. Chem. Theory Comput., 2013, 9, 1356.
19. A. Fernández-Ramos, B. A. Ellingson, R. Meana-Pañeda, J. MC Marques and D. G. Truhlar, Theor. Chem. Acc., 2007, 118, 813.
20. N. M. Vandewiele, R. V. de Vijver, K. M. Van Geem, M.-F. Reyniers and G. B. Marin, J. Comput. Chem., 2015, 36, 181.
21. (a) H.-H. Carstensen and A. M. Dean, Compr. Chem. Kinet., 2007, 42, 101; (b) J. W. Bozzelli, A. Y. Chang and A. M. Dean, Int. J. Chem. Kinet., 1997, 29, 161; (c) A. Y. Chang, J. W. Bozzelli and A. M. Dean, Z. Phys. Chem., 2000, 214, 1533.
22. J. Troe, J. Chem. Phys., 1977, 66, 4745.
23. J. Troe, J. Chem. Phys., 1977, 66, 4758.
24. K. J. Laidler, Chemical Kinetics, McGraw-Hill, New York, 2nd edn, 1965, p. 145.
25. K. A. Holbrook, M. J. Pilling and S. H. Robertson, Unimolecular Reactions, John Wiley & Sons, New York, 2nd edn, 1996, ch. 2.
26. T. Baer and W. L. Hase, Unimolecular Reaction Dynamics, Oxford University Press, New York, 1996, p. 5.
27. A. M. Dean, J. Phys. Chem., 1985, 89, 4600.
28. Y. Zhao and D. G. Truhlar, J. Chem. Theory Comput., 2008, 4, 1849.
29. B. J. Lynch, Y. Zhao and D. G. Truhlar, J. Phys. Chem. A, 2003, 107, 1384.
30. R. Krishnan, J. S. Binkley, R. Seeger and J. Pople, J. Chem. Phys., 1980, 72, 650.
31. T. Clark, J. Chandrasekhar, G. W. Spitznagel and P. Schleyer, J. Comput. Chem., 1983, 4, 294.
32. M. J. Frisch, J. A. Pople and J. S. Binkley, J. Chem. Phys., 1984, 80, 3265.
33. V. I. Lebedev and L. Skorokhodov, Russ. Acad. Sci. Dokl. Math., 1992, 45, 587.
34. M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Robb, J. R. Cheeseman, G. Scalmani, V. Barone, B. Mennucci, G. A. Petersson, H. Nakatsuji, M. Caricato, X. Li, H. P. Hratchian, A. F. Izmaylov, J. Bloino, G. Zheng, J. L. Sonnenberg, M. Hada, M. Ehara, K. Toyota, R. Fukuda, J. Hasegawa, M. Ishida, T. Nakajima, Y. Honda, O. Kitao, H. Nakai, T. Vreven, J. A. Montgomery, Jr., J. E. Peralta, F. Ogliaro, M. Bearpark, J. J. Heyd, E. Brothers, K. N. Kudin, V. N. Staroverov, R. Kobayashi, J. Normand, K. Raghavachari, A. Rendell, J. C. Burant, S. S. Iyengar, J. Tomasi, M. Cossi, N. Rega, J. M. Millam, M. Klene, J. E. Knox, J. B. Cross, V. Bakken, C. Adamo, J. Jaramillo, R. Gomperts, R. E. Stratmann, O. Yazyev, A. J. Austin, R. Cammi, C. Pomelli, J. Ochterski, R. L. Martin, K. Morokuma, V. G. Zakrzewski, G. A. Voth, P. Salvador, J. J. Dannenberg, S. Dapprich, A. D. Daniels, O. Farkas, J. B. Foresman, J. V. Ortiz, J. Cioslowski and D. J. Fox, Gaussian 09 (Revision D. 01), Gaussian, Inc., Wallingford, CT, 2009.
35. Y. Zhao, R. Peverati, K. Yang, H. Xiao, H. Yu and D. G. Truhlar, MN-GFM version 6.5 computer program module, University of Minnesota, Minneapolis, 2015.
36. E. Papajak and D. G. Truhlar, J. Chem. Phys., 2012, 137, 064110.
37. E. Papajak and D. G. Truhlar, J. Chem. Theory Comput., 2010, 6, 597.
38. R. Peverati and D. G. Truhlar, J. Phys. Chem. Lett., 2012, 3, 117.
39. R. Peverati and D. G. Truhlar, Phys. Chem. Chem. Phys., 2012, 14, 13171.
40. H. S. Yu, H. Xiao and D. G. Truhlar, J. Chem. Theory Comput., 2016, 12, 1280.
41. P. J. Stephens, F. J. Devlin, C. F. Chabalowski and M. J. Frisch, J. Phys. Chem., 1994, 98, 11623.
42. C. Adamo and V. Barone, J. Chem. Phys., 1999, 110, 6158.
43. V. N. Staroverov, G. E. Scuseria, J. Tao and J. P. Perdew, J. Chem. Phys., 2003, 119, 12129.
44. J. Sun, R. Haunschild, B. Xiao, I. W. Bulik, G. E. Scuseria and J. P. Perdew, J. Chem. Phys., 2013, 138, 044113.
45. B. J. Lynch, P. L. Fast, M. Harris and D. G. Truhlar, J. Phys. Chem. A, 2000, 104, 4811.
46. Y. Zhao, N. E. Schultz and D. G. Truhlar, J. Chem. Phys., 2005, 123, 161103.
47. Y. Zhao and D. G. Truhlar, Theor. Chem. Acc., 2008, 120, 215.
48. Y. Zhao, N. E. Schultz and D. G. Truhlar, J. Chem. Theory Comput., 2006, 2, 364.
49. R. Peverati and D. G. Truhlar, J. Phys. Chem. Lett., 2011, 2, 2810.
50. R. Peverati and D. G. Truhlar, Phys. Chem. Chem. Phys., 2012, 14, 16187.
51. R. Peverati and D. G. Truhlar, J. Chem. Phys., 2011, 135, 191102.
52. T. W. Keal and D. J. Tozer, J. Chem. Phys., 2005, 123, 121103.
53. A. V. Krukau, O. A. Vydrov, A. F. Izmaylov and G. E. Scuseria, J. Chem. Phys., 2006, 125, 224106.
54. A. D. Boese and N. C. Handy, J. Chem. Phys., 2002, 116, 9559.
55. J. D. Chai and M. Head-Gordon, Phys. Chem. Chem. Phys., 2008, 10, 6615.
56. J. A. Pople, M. Head-Gordon and K. Raghavachari, J. Chem. Phys., 1987, 87, 5968.
57. Y. Zhao and D. G. Truhlar, J. Phys. Chem. A, 2005, 109, 6624.
58. T. H. Dunning, J. Chem. Phys., 1989, 90, 1007.
59. M. Head-Gordon, J. A. Pople and M. J. Frisch, Chem. Phys. Lett., 1988, 153, 503.
60. D. G. Truhlar, Chem. Phys. Lett., 1998, 294, 45.
61. A. Karton and J. M. L. Martin, Theor. Chem. Acc., 2006, 115, 330.
62. A. Halkier, T. Helgaker, P. Jørgensen, W. Klopper, H. Koch, O. Jeppe and A. K. Wilson, Chem. Phys. Lett., 1998, 286, 243.
63. Y. Zhao, H. T. Ng, R. Peverati and D. G. Truhlar, J. Chem. Theory Comput., 2012, 8, 2824.
64. J. L. Bao, H. S. Yu, K. Duanmu, M. A. Makeev, X. Xu and D. G. Truhlar, ACS Catal., 2015, 5, 2070.
65. K. A. Nguyen, C. F. Jackels and D. G. Truhlar, J. Chem. Phys., 1996, 104, 6491.
66. C. F. Jackels, Z. Gu and D. G. Truhlar, J. Chem. Phys., 1995, 102, 3188.
67. T. Joseph, R. Steckler and D. G. Truhlar, J. Chem. Phys., 1987, 87, 7036.
68. J. Villà and D. G. Truhlar, Theor. Chem. Acc., 1997, 97, 317.
69. I. M. Alecu, J. Zheng, Y. Zhao and D. G. Truhlar, J. Chem. Theory Comput., 2010, 6, 2872.
70. J. Zheng and D. G. Truhlar, Faraday Discuss., 2012, 157, 59.
71. J. Zheng and D. G. Truhlar, J. Chem. Theory Comput., 2013, 9, 2875.
72. J. Zheng, S. L. Mielke, K. L. Clarkson and D. G. Truhlar, MSTor computer program, version 2011, University of Minnesota, Minneapolis, MN, 2011.
73. J. Zheng, S. Zhang, B. J. Lynch, J. C. Corchado, Y.-Y. Chuang, P. L. Fast, W.-P. Hu, Y.-P. Liu, G. C. Lynch, K. A. Nguyen, C. F. Jackels, A. F. Ramos, B. A. Ellingson, V. S. Melissas, J. Vill_a, I. Rossi, E. L. Coitino, J. Pu, T. V. Albu, R. Steckler, B. C. Garrett, A. D. Isaacson and D. G. Truhlar, POLYRATE, version 2010-A, University of Minnesota, Minneapolis, 2010.
74. J. Zheng, S. Zhang, J. C. Corchado, Y. Y. Chuang, E. L. Coitino, B. A. Ellingson and D. G. Truhlar, GAUSSRATE, version 2009-A, University of Minnesota, Minneapolis, 2010.
75. Y. Shi, Acc. Chem. Res., 2015, 48, 163.
76. A. Rahman, Phys. Rev., 1964, 136, A405.
77. M. J. Kushner, J. Appl. Phys., 1988, 63, 2532.
78. A. Dollet, S. de Persisb and F. Teyssandier, Phys. Chem. Chem. Phys., 2004, 6, 1203.
79. (a) A. W. Jasper, J. A. Miller and S. J. Klippenstein, J. Phys. Chem. A, 2013, 117, 12243; (b) A. W. Jasper, K. M. Pelzer, J. A. Miller, E. Kamarchik, L. B. Harding and S. J. Klippenstein, Science, 2014, 346, 1212.
80. J. Troe, J. Phys. Chem., 1979, 83, 114.
81. K. R. Yang, A. Jalan, W. H. Green and D. G. Truhlar, J. Chem. Theory Comput., 2013, 9, 418.
82. X. Xu, W. Zhang, M. Tang and D. G. Truhlar, J. Chem. Theory Comput., 2015, 11, 2036.
83. T. J. Lee and P. R. Taylor, Int. J. Quantum Chem., Symp., 1989, 23, 199.
84. J. C. Rienstra-Kiracofe, W. D. Allen and H. F. Schaefer, J. Phys. Chem. A, 2000, 104, 9823.
85. M. Martinez-Avila and J. Peiro-Garcia, J. Phys. Chem. Solids, Lett. Sect., 2003, 370, 313.
86. N. Lambert, N. Kaltsoyannis, S. D. Price, J. Zabka and Z. Herman, J. Phys. Chem. A, 2006, 110, 2898.
87. J. Zheng and D. G. Truhlar, Phys. Chem. Chem. Phys., 2010, 12, 7782.

---

Footnote

† Electronic supplementary information (ESI) available: Classical forward barriers and energy of reactions; mean unsigned errors (MUEs) computed by various model chemistries; reaction symmetry numbers; MS-T factors of activation for reverse reactions; MS-CVT/SCT rate constants for reverse reactions; fitting parameters for MS-CVT/SCT rate constants of reverse reactions; activation energies. See DOI: 10.1039/c6cp00816j

---