Chemical Physics Letters 425 (2006) 201–204 www.elsevier.com/locate/cplett
Theoretical study of the thermochemistry and kinetics of the addition of silyl radical to ethylene Owen J. Clarkin b
a,1
, Gino A. DiLabio
b,*
a Department of Chemistry, Carleton University, 1125 Colonel By Drive, Ottawa, Ont., Canada K1S 5B6 National Institute for Nanotechnology, National Research Council, Molecular Scale Devices Group, W6-010 ECERF, 9107-116th Street, Edmonton, Alta., Canada T6G 2V4
Received 1 February 2006; in final form 10 May 2006 Available online 13 May 2006
Abstract The thermochemistry and kinetics of the H3Si + C2H4 ! H3Si-CH2CH2 reaction were studied using UQCISD(T)/aug-cc-pVTZ// UQCISD/aug-cc-pVDZ. Our results show that the Arrhenius activation energy, Ea, for the reaction is 3.3 kcal/mol, substantially lower than the 6 kcal/mol value determined experimentally by Loh et al. [S.K. Loh, D.B. Beach, J.M. Jasinski, Chem. Phys. Lett. 169 (1990) 55]. Loh et al.’s high Ea was the result of an estimate of the experimental A-factor that was two orders of magnitude too high. Our calculated rate constant of 1.41 · 106 M1 s1 is in very good agreement with the upper bound of 1.81 · 106 M1 s1 determined by Loh et al. A study of the thermochemistry and kinetics of the addition reaction using 15 density functional theory methods with 6-31G(d) basis sets indicates that BHandHLYP predicts results in best agreement with high-level theory. Published by Elsevier B.V.
1. Introduction The theoretical modeling of silicon-centered radicals is able to yield insights into the reactivity of incompletely terminated silicon surfaces toward organic molecules. Our long-standing interests in this area relate to the ability of vinyl (and aldehyde) substituted molecules to form linear nanostructures on silicon surfaces via a self-directed growth process [1,2]. The first step in the self-directed growth mechanism is the addition of the vinyl (aldehyde) moiety to a silicon surface dangling bond. The simplest model for the first step of the self-directed growth of alkenes on silicon is the addition of H3Si to ethylene according to H3 Si + C2 H4 ! H3 Si–CH2 C H2
*
ð1Þ
Corresponding author. Fax: +1 780 492 8632. E-mail address:
[email protected] (G.A. DiLabio). 1 Present address: Department of Chemistry, Queen’s University, 90 Bader Lane, Kingston, Ont., Canada K7L 3N6. 0009-2614/$ - see front matter. Published by Elsevier B.V. doi:10.1016/j.cplett.2006.05.030
Despite the simplicity of this chemical system, there are inconsistent reports on the fundamental quantities associated with the reaction. For example, an upper bound for the gas-phase reaction rate constant of 1.81 · 106 M1 s1 [3] was measured by Loh et al. These authors assumed a pre-exponential factor (A-factor) of 6 · 1010 M1 s1 to obtain an Arrhenius activation energy, Ea, of 6 kcal/mol for the addition. However, the pre-exponential factor used by Loh et al. is ca. two orders of magnitude higher than the value of 2 · 108 M1 s1 for the analogous methyl radical addition reaction [4]. The only theoretical study of (1), performed by Bottoni [5], reported an Ea of 6.46 kcal/mol (energy barrier of 4.31 kcal/mol), in surprisingly good agreement with the Loh et al. result. Arthur and Potzinger [6] noted some other inconsistencies with some of Bottoni’s predictions for reactions of H3Si + C3H6 and Me3Si + C2H4. It should be noted that the theoretical methods used by Bottoni were considered to be of high quality at the time of his study. In this work, we re-examine the thermochemistry and kinetics of (1) using high-level correlated ab initio techniques. Density functional theory (DFT) methods are
202
O.J. Clarkin, G.A. DiLabio / Chemical Physics Letters 425 (2006) 201–204
currently the only practical approaches that can be used for modeling reactions on silicon surfaces. We therefore benchmark some common DFTs for their ability to predict the key quantities for (1).
3. Results and discussion
2. Computational methods
The key calculated data for (1) are given in Table 1. The UQCISD/6-31G(d) method predicts a TS SiAC bond ˚ shorter than those found with length that is ca. 0.05 A UQCISD using 6-311+G(d, p) and aug-cc-pVDZ. These differences are consistent with the findings of Chuang et al [8]. The use of the two larger basis sets yield TS bond separations that are in agreement. Arrhenius activation energies and classical reaction barriers at the UQCISD level span a fairly broad range of ca. 2.5 kcal/mol, as do the reaction enthalpies. Computing single-point energies with UQCISD(T)/augcc-pVTZ with the three different UQCISD geometries gives barrier and enthalpy results that are within 0.2 kcal/mol. The slight geometry differences at the lower level of theory do not have a large effect on the computed single-point energetics. We have previously made similar conclusions with respect to this point in relation to bond dissociation enthalpy calculations [10]. Our highest-level calculations (UQCISD(T)/aug-cc-pVTZ//UQCISD/aug-cc-pVDZ) give ˚ Ea = 3.3 kcal/mol, and DH = 20.6 kcal/ RzSi–C ¼ 2:63 A, mol for the reaction. Bottoni reported [5] a QCISD(T)/6-311G(d, p)//MP2/631G(d) Ea = 6.46 kcal/mol which is close to our UQCISD/6-311+G(d, p) value and ca. 3.2 kcal/mol larger than our best calculated value, see Table 1. We repeated Bottoni’s calculations and found that the errors in his results arise from two sources: Firstly, the use of MP2/ 6-31G(d) for geometry optimizations yields a RzSi–C ¼ ˚ which compares reasonably well to our values. 2:58 A, However, the silyl moiety is oriented such that the orbital containing the unpaired electron is not pointing toward the p-orbital of the ethylene, see Fig. 1, but is off by ca. 11 compared to our QCISD(T)/6-311+G(d, p) value. The poorly predicted MP2 TS geometry is partly or completely due to the large degree of spin contamination (ÆS2æ = 0.94) in the TS wavefunction. This is an indication that the results obtained from the MP2 method are unreliable [8]. However, the poorly predicted TS geometry contributes only ca. 0.3 kcal/mol to the error in Ea. More importantly, the Ea determined using UQCISD(T)/6311G(d, p)//UQCISD/6-31G(d) is higher than our best value by 2.3 kcal/mol. The latter indicates that the basis sets used by Bottoni for his single-point calculations, although substantial at the time of his study, are too small to obtain an accurate Ea. As a check of our methods, we calculated the Ea for the analogous H3C + C2H4 addition reaction, for which there exists substantial experimental and theoretical data [4]. As was determined for the silyl reaction, UQCISD calculations with small basis sets yield Ea values that are too high by up to 4 kcal/mol (viz. Ea (UQCISD/6-31G(d)) = 11.3 kcal/mol). Our highest-level approach gives a
All calculations were performed using the GAUSSIAN-98 suite of programs [7]. For the highest-level calculations, geometry optimization, and frequency calculations were carried out using unrestricted singles and doubles quadratic configuration interaction (UQCISD) with basis sets ranging in size from 6-31G(d) to aug-cc-pVDZ. Transition states (TS) were verified by ensuring that the structures have one imaginary frequency and that the mode connects the reactant and product. Single-point energies were calculated at the QCISD(T)/aug-cc-pVTZ level of theory. Recent work by Chuang et al. on hydrogen atom transfer reactions [8], and the results of some of our benchmarking efforts (see Table 1), suggest that our chosen methods should provide good results. DFT calculations were performed using various methods with 6-31G(d) basis sets. This basis set is selected because the modeling of addition reactions on silicon surfaces necessitates the use of small basis sets. Test calculations indicate that the predicted energetics are not necessarily converged with respect to basis set with DFT/ 6-31G(d). A temperature of 298.15 K was used for all calculations. The Ea values were determined using the standard equation for a general bimolecular reaction leading to a TS, i.e., A + B ! [A B]à [9]: Ea ¼ DE0 þ DEtherm þ RT where DE0 is the zero-point corrected reaction energy barrier and DEtherm is the difference in the thermal contributions to the energies of the reactants and TS. A-factors are obtained assuming that rate constants obtained from Arrhenius theory (Arr) and conventional transition state theory (CTST) are equivalent: k Arr ¼ A-factor eðEa =RT Þ ¼ k CTST ¼
k B T QTS ðDE0 =RT Þ e : h QA QB
Table 1 ˚ Arrhenius activation energy (Ea/ Transition state bond length (RzSi–C =A), kcal/mol), classical barrier height (Ec/kcal/mol) and reaction enthalpy (DH/kcal/mol) for (1) at various levels of theory Method/basis set
RzSi–C
Ea (EC)
DH
UQCISD/6-31G(d) UQCISD/6-311+G(d, p) UQCISD/aug-cc-pVDZ UQCISD(T)/aug-cc-pVTZ// UQCISD/6-31G(d) UQCISD(T)/aug-cc-pVTZ// UQCISD/6-311+G(d, p) UQCISD(T)/aug-cc-pVTZ// UQCISD/aug-cc-pVDZ
2.59 2.64 2.63 2.59
6.9 6.2 4.5 3.4
(5.2) (4.5) (2.8) (1.6)
17.5 20.1 18.5 20.4
2.64
3.5 (1.7)
20.4
2.63
3.3 (1.6)
20.6
3.1. Highly-correlated wavefunction treatment of reaction (1)
O.J. Clarkin, G.A. DiLabio / Chemical Physics Letters 425 (2006) 201–204
203
˚ , hH –Si–C ¼ 129:0 , hSi–C–C = 101.2 and the (b) MP2/6Fig. 1. Optimized TS geometry for (1) at the: (a) QCISD(T)/6-311+G(d, p) level: RSi–C = 2.64 A ˚ , hH –Si–C ¼ 140:1 , hSi–C–C = 96.6. 31G(d) level: RSi–C = 2.58 A
Ea = 8.25 kcal/mol, in fairly good agreement with the experimental value of 7.4 kcal/mol given in Ref. [4]. The rate of (1) can be computed via transition state theory (TST) using the partition functions (Q) for the reaction participants [9]. We calculated Q values using UQCISD/ aug-cc-pVDZ and their ratios yield an A-factor 3.81 · 108 M1 s1. Combined with our best calculated value for Ea, the rate constant for (1) is computed to be 1.41 · 106 M1 s1, in very good agreement with the upper bound of 1.81 · 106 M1 s1 estimated by Loh et al. [3]. Loh et al. derived Ea = 6 kcal/mol for reaction (1) by assuming an A-factor of 6 · 1010 M1 s1. Using our calculated A-factor for (1), which is in excellent agreement with that of the H3C + C2H4 addition reaction [4], along with Loh et al.’s measured rate constant yields a more consistent experimental upper bound to Ea of 3.2 kcal/mol. 3.2. DFT Treatments of reaction (1)
Table 2 ˚ A-factor (/108 M1 s1), Arrhenius Transition state bond length (RzSi–C =A), activation energy (Ea/kcal/mol), rate constant (k/106 M1 s1), and reaction energy (DH/kcal/mol) for (1) using DFT/6-31G(d) DFT
RzSi–C
A-factor
Ea
G96LYP B1LYP BLYP B3LYP BHandHLYP BHandH B3P86 B3PW91 MPW1PW91 MPW1K PBE1PBE PBE PW91 BB1K B1B95 Referencec
2.80 2.72 2.79 2.73 2.69 2.87 2.82 2.79 2.80 2.76 2.82
13.2 6.3 10.7 6.9 5.1 13.1a 8.3 7.2 6.7 5.5 7.0
4.2 3.0 2.1 2.5 3.3
a
Transition state bond lengths, A-factors, Eas, rate constants (k), and DHs for the H3Si + C2H4 addition reaction using 15 DFT methods with 6-31G(d) basis sets are presented in Table 2. These small basis sets were chosen because they are the largest practical ones for our modeling of organics on extended silicon surfaces, which typically consist of several hundred atoms. The DFT methods tested, not surprisingly, produce results of variable quality. All of the DFT methods predict RzSi–C values that are larger than the high-level ab initio reference value, with the exception of PBE and PW91 which predict no saddle point at all along the reaction path. The deviations in DFT RzSi–C values with respect to the ref˚ (BHandH) to +0.06 A ˚ erence value range from +0.24 A (BHandHLYP). The A-factors, which depend parametrically on molecular geometry and predicted vibration frequencies, are also wide ranging with the largest value being ca. 3.5 times (G96LYP) larger than the reference value. B1B95, BB1K, and BHandHLYP predict A-factors that are within 1.7 times that of the reference value. The predicted Ea values span a range of ca. 3.3 kcal/mol, with
b c
a
1.2 2.3 1.5 1.7 0.9
b
b
b
b
b
b
2.70 2.74 2.63
4.0 5.1 3.8
2.5 2.2 3.3
k
DH
1.1 4.1 33.5 9.5 2.0 a
106.8 15.6 55.9 30.8 151.8 b b
6.2 12.7 1.4
13.2 18.4 14.8 18.4 22.4 30.3 22.2 20.8 22.8 25.5 23.6 20.4 20.9 22.1 19.8 20.6
Negative barrier height with respect to reactants. No saddle point found on minimum energy pathway. Best estimates from QCISD(T) calculations. See Table 1.
the BHandHLYP value being in best agreement with the reference value. The reactions rate constants are exponentially dependent on Ea and linearly dependent on A-factors. Collectively, the DFT rate constants are within about two orders of magnitude of the reference value, kref. The G96LYP k is in best agreement with the reference value for kref due to a fortuitous cancellation of errors. The contribution to the rate constant by the large G96LYP A-factor is compensated by the overly large barrier height. Only BB1K, B1LYP, and BHandHLYP predict both A-factors and Ea values, and therefore k values, that are in good agreement with the reference values. Overall, BHandHLYP displays the best performance for the kinetic parameters associated with reaction (1). Similarly good results were found for the H3C + C2H4 addition reaction (data not shown) using BHandHLYP/6-31G(d). These findings are consistent with those of Durant [11]. The B3PW91, PBE,
204
O.J. Clarkin, G.A. DiLabio / Chemical Physics Letters 425 (2006) 201–204
PW91, and B1B95 methods predict DH(1) to within ca. 1 kcal/mol of the reference value. BHandHLYP predicts a DH(1) value that is too low by 1.8 kcal/mol.
data. The predicted Ea value is in perfect agreement with the reference data while the DH is too low by 1.8 kcal/mol. References
4. Summary The thermochemistry and kinetics of the silyl radical + ethylene addition reaction were re-examined using high-level correlated ab initio techniques. Our results show that the Arrhenius activation energy is 3.3 kcal/mol, substantially lower than the 6 kcal/mol value previously determined experimentally [3]. The erroneous experimental Ea was the result of an estimate of the experimental A-factor in Ref. [3] that was two orders of magnitude too high. However, the calculated and experimental rate constants are in excellent agreement. The thermochemistry and kinetics of the addition reaction were also studied using 15 DFT methods with 6-31G(d) basis sets. The results indicate that BHandHLYP performs the best for the addition reaction. The A-factor and reaction rate constant are predicted to within a factor of 1.5 of the high-level ab initio reference
[1] See, for example P.G. Piva, G.A. DiLabio, J.L. Pitters, J. Zikovsky, M. Rezeq, S. Dogel, W.A. Hofer, R.A. Wolkow, Nature 435 (2005) 658. [2] J.L. Pitters, I. Dogel, G.A. DiLabio, R.A. Wolkow, J. Phys. Chem. B 110 (2005) 2159, and references therein. [3] S.K. Loh, D.B. Beach, J.M. Jasinski, Chem. Phys. Lett. 169 (1990) 55. [4] For an excellent review, see: H. Fischer, L. Radom, Angew. Chem., Int. Ed. 40 (2001) 1340. [5] A. Bottoni, J. Phys. Chem. A 101 (1997) 4402. [6] N.L. Arthur, P. Potzinger, Organometallics 21 (2002) 2874. [7] M.J. Frisch et al., GAUSSIAN 98, Revisions A.6, A.7, A.11, A.13, Gaussian Inc., Pittsburg, PA, 1998. [8] Y.-Y. Chuang, E.L. Coitin˜o, D.G. Truhlar, J. Phys. Chem. A 104 (2000) 446. [9] K.J. Laidler, Chemical Kinetics, third edn., Harper and Row, 1987. [10] E.R. Johnson, O.J. Clarkin, G.A. DiLabio, J. Phys. Chem. A 107 (2003) 9953, and references therein. [11] J.L. Durant, Chem. Phys. Lett. 256 (1996) 595.