Computational Tests of Models for Kinetic Parameters of Unimolecular

ARTICLE pubs.acs.org/JPCA

Computational Tests of Models for Kinetic Parameters of Unimolecular Reactions of Organophosphorus and Organosulfur Compounds David K. Hahn, Krishans S. RaghuVeer, and J. V. Ortiz* Department of Chemistry and Biochemistry, Auburn University, Auburn, Alabama 36849-5312, United States

bS Supporting Information ABSTRACT: A computational study of the kinetics of isomerization and elimination reactions of organophosphorus and organosulfur reactions is presented with a view to characterizing the predictive capabilities of widely applied techniques for processes that pertain to the destruction of chemical warfare agents. A set of 22 reactions has been studied, and the results have been compared to experimentally derived data. The BMK functional and the MG3S basis set have been used to compute minimum energy paths. Corrections have been added from CBS-QB3, CASSCF, and CASMP2 calculations. Thermal rate constants at experimental temperatures have been calculated with canonical variational transition state theory and small-curvature tunneling theory. The quality of these results may depend on recrossing of the variational transition state, the amount of radical or diradical character found in the minimum energy paths, or the accuracy of barrier heights.

’ INTRODUCTION Chemical warfare agents (CWAs) are organic compounds of phosphorus and sulfur. United Nations treaties require treaty signatory countries to destroy all CWAs in weapons and in bulk storage. One of the many ways of accomplishing this end is to employ thermal means such as pyrolysis or combustion. A clear understanding of these thermal processes requires knowledge of the thermochemical properties, bond energies, and kinetic parameters for unimolecular and bimolecular gas-phase reactions. Measurement of these properties of the actual CWAs in a laboratory is impractical because of safety and security issues associated with these compounds. In general, the literature contains thermochemical properties of phosphorus and sulfur compounds only to a limited extent and none is available for all of the CWAs themselves. The only example of a detailed study of experimental kinetics and reaction mechanisms for thermolysis of nerve agents is for the nerve agent sarin (GB), and no models exist for mustard (HD) and VX.1 Often, the industry has chosen safer phosphorus and sulfur compounds (simulants) that resemble the CWAs for actual laboratory and field measurements and extrapolate the results to derive the properties of the CWAs.2 However, there are drawbacks to such extrapolation procedures because variations in activation energies or absence of certain bonds in the simulants compared to a CWA, for example, can lead to entirely different reaction pathways for the CWA compared to their simulants. For example, tributyl phosphate or O,S-diethyl methylphosphonothioate have been used as VX simulants. In the first compound a PS bond is absent, and in the second compound the amino group found in VX is absent. Therefore, any reaction mechanisms derived for the thermolysis of these two simulants will have serious limitations in the prediction of mechanisms for VX. r 2011 American Chemical Society

The alternative to the simulant approach is to perform quantum chemical modeling (QCM) of the molecular structures and properties of CWAs. We have carried out this approach in several stages. The purpose of the modeling we have presented has been to determine which quantum and thermochemical models provide the best results compared to the best available experimental values of molecular geometries, vibrational frequencies, and thermochemical parameters such as heats of formation and bond dissociation energies. When the models are validated against experimental results those models and methods will be applied to CWAs to provide properties of CWAs with confidence. In the first stage of this effort we demonstrated that CBS-Q , CBS-QB3, G3X, G4, and W1 thermochemical models predicted the structures and thermochemical properties of small phosphorus and sulfur compounds with acceptable error margins.3 For further studies we have now identified CBS-QB3 and G4 models as best suited to predict geometries and thermochemical properties of larger phosphorus- and sulfur-containing compounds. In the second stage of this effort we chose a set of 12 compounds that can be used as simulants for CWAs,4 and six of these simulants were found in the experimental study by MacNaughton et al.2 These 12 compounds were also chosen because the literature contains measured structural and thermochemical properties.1,2,59 We applied the CBS-QB3 and G4 models to these 12 compounds to calculate the ground-state structures, vibrational frequencies, heats of formation, and bond dissociation energies of key bonds in the simulants. Heats of formation Received: July 5, 2011 Revised: November 1, 2011 Published: November 03, 2011 14143

dx.doi.org/10.1021/jp206344r | J. Phys. Chem. A 2011, 115, 14143–14152

The Journal of Physical Chemistry A

ARTICLE

calculated with either CBS-QB3 or G4 energies showed a mean absolute deviation from experiment of less than 1.4 and 1.2 kcal/mol, respectively. Calculated vibrational frequencies showed errors in the range of 0.31% compared to experimental values. Theoretical formation enthalpies were used to compute bond dissociation energies in all cases. These results for small and large molecules at once indicate that these thermochemical models can be successfully applied to predict ground-state properties of the large CWA molecules. This paper is the third stage of the QCM. In this paper, we present computational modeling of the kinetics of isomerization and elimination of groups from the parent compounds in which the central atom is phosphorus or sulfur. In some cases, we have also computed the kinetic parameters for elimination reactions from the first elimination products. Initially, the minimum energy structures and vibrational frequencies of the ground states of each molecule, transition states for elimination, and the product molecules were determined by density functional theory (DFT) calculations using the B3LYP functional and modest basis sets such as 6-31G(d,p), 6-311G++(2df, 2p), or cc-pVTZ. It is known that similar basis sets, the B3LYP functional, and Gn thermochemical models show large errors in the calculated reaction activation energies.10,11 In view of this, the minimum energy structures and transition states structures were reoptimized using the BMK density functional and the MG3S basis set.12 In the literature there is a plethora of density functionals that provide accurate Arrhenius parameters for unimolecular and bimolecular gas-phase reactions. The performance of many density functionals has been extensively reviewed by Truhlar and co-workers.10,1315 In the present work, kinetic parameters were calculated using the program POLYRATE.16 As in the previous papers, we chose of set of 18 compounds and reactions because experimental kinetics data are available in the literature, thus enabling us to compare calculated to experimental values. Three of these reactions came from the work of MacNaughton et al.2 We have conducted a computational investigation of the kinetics of thermal decomposition of phosphorus- and sulfur-containing molecules. Dual-level direct dynamics was applied to the following unimolecular reactions: CH3 CðSÞOCH3 f CH3 CðOÞSCH3

ð1Þ

C3 H6 S f SCH2 þ C2 H4

ð2Þ

t C4 H9 SH f SH2 þ i C4 H8

ð3Þ

CH2 dCHCH2 SCH3 f SCH2 þ CH3 CHdCH2

ð4Þ

CH2 dCHCH2 SCH2 Cl f SCHCl þ CH3 CHdCH2

ð5Þ

CH3 CðSÞSC2 H5 f CH3 CðSÞSH þ C2 H4

ð6Þ

CH3 CðOÞSC2 H5 f CH3 CðSÞOH þ C2 H4

ð7Þ

CH3 SCðSÞOC2 H5 f CH3 SCðOÞSH þ C2 H4

ð8Þ

CH3 OCðSÞOC2 H5 f CH3 OCðOÞSH þ C2 H4

ð9Þ

CH3 OCðOÞSC2 H5 f CH3 OCðSÞOH þ C2 H4

ð10Þ

CH3 SðCH2 Þ4 Cl f tetrahydrothiophene þ CH3 Cl

ð11Þ

POðCH3 ÞðOC2 H5 Þ2 f POðOHÞðCH3 ÞðOC2 H5 Þ þ C2 H4

ð12Þ POðOC2 H5 Þ3 f POðOHÞðOC2 H5 Þ2 þ C2 H4

ð13Þ

PðOC2 H5 Þ3 f HPOðOC2 H5 Þ2 þ C2 H4

ð14Þ

POðC2 H5 ÞðOC2 H5 ÞðSC2 H5 Þ f PSðOHÞðC2 H5 ÞðOC2 H5 Þ þ C2 H4

ð15Þ PhPðCH3 ÞðCH2 CHCH2 Þ f PhPdCH2 þ CH3 CHCH2

ð16Þ

PhPðCH2 CHCH2 Þ2 f PhPdCHCHCH2 þ CH3 CHCH2

ð17Þ

POðCH3 ÞðOC2 H5 ÞðFÞ f PO2 ðCH3 Þ þ CH3 CH2 F

ð18Þ

ðCH3 Þ3 CPH2 f PH2 þ ðCH3 Þ3 C

ð19Þ

This test set consists of one isomerization (reaction 1), 17 eliminations (reactions 218), and one dissociation (reaction 19). All reactions except 12 and 13 involve bond making or bond breaking with a phosphorus or sulfur atom, and all reactions except 11 occur in one step. All reactions except 18 have experimental gasphase rate constant data available from literature sources. All studied reactions except perhaps 13 and 14 obey first-order kinetics under experimental conditions, indicating measurements taken in the high-pressure limit, suitable for comparison to variational transition state theory.

’ COMPUTATIONAL METHODS The stationary points and minimum energy paths associated with reactions 118 were computed at the BMK/MG3S level of theory. The BMK functional is a meta hybrid functional designed for application to chemical kinetics through the inclusion of the kinetic energy density and a large percentage of exact exchange (specifically, 42%).12 In addition, it is parametrized against a training set that includes 22 reaction barrier heights, two of which pertain to homolytic cleavage of a bond to phosphorus or sulfur, and 23 reaction energies.12 The MG3S basis set is 6-311+G(2df,2p) on hydrogen and first-row elements and a modified version of 6-311+G(3d2f) on phosphorus through chlorine obtained by the deletion of core polarization functions on non-hydrogen atoms and diffuse functions on H.17 To minimize the error in DFT energies from the numerical calculation of two-electron integrals, a very fine pruned grid (99 radial points, 590 angular points per radial point) was used for the BMK calculations. To minimize the error in calculated equilibrium geometries and vibrational frequencies, the cutoff value for determining convergence of optimized structures was 2.0 106 hartree/bohr for maximum force and 6.0 106 bohr for maximum displacement. All quantum chemical calculations were performed with Gaussian 09.18 For reactions 118, the minimum energy paths were calculated to s = (2.5 bohr from the saddle point in steps of 0.01 bohr using the Euler steepest-descents method.19 The generalized free energy was computed at the saddle point and every 0.05 bohr along the minimum energy path (MEP). Rate constants were determined from this data with canonical variational transition state theory20 (CVT) using partition functions evaluated within the harmonic oscillator and rigid rotor approximations. Internal 14144

dx.doi.org/10.1021/jp206344r |J. Phys. Chem. A 2011, 115, 14143–14152


ARTICLE

Figure 1. BMK/MG3S first-order saddle points of reactions 110.

rotation was treated as a harmonic oscillator in all cases. Vibrational partition functions, as well as zero-point energies, were evaluated with BMK/MG3S harmonic frequencies scaled by a factor of 0.9835 to remove systematic error.13 A multiplicative transmission coefficient that corrects the CVT rate constant for multidimensional tunneling and nonclassical reflection was computed at the centrifugal-dominant small-curvature tunneling

(SCT) level of theory.21 All kinetics calculations were performed with POLYRATE.16 Corrected rate constants were obtained with variational transition state theory with interpolated single-point energies (VTST-ISPE).16 This method corrects the energies along the MEP (but not the vibrational frequencies or moments of inertia) using higher level electronic structure calculations at various 14145



ARTICLE

Figure 2. BMK/MG3S first-order saddle points of reaction 11.

points on the MEP. CBS-QB322 single-point calculations were therefore performed on the BMK/MG3S-optimized reactants, saddle points, and products of reactions 110 and 1218. For reaction 11 , the MEP was corrected with a CASMP223 singlepoint calculation on the stationary points, using the MG3S basis set and an active space consisting of eight electrons distributed over eight orbitals (HOMO 3 to LUMO + 3). Pyrolysis of tert-butyl phosphine (C4H9PH2) is believed to proceed by a free radical mechanism because its rate and products are strongly affected by the input C4H9PH2 concentration.24 The first-order initiation step is believed to be dissociation of phosphino, based on the reaction (activation) temperature. In this study the transition state for dissociation of phosphino from tert-butyl phosphine, reaction 19, was taken as the structure obtained from a partial BMK/MG3S geometry optimization in which the tertiary carbonphosphorus bond length is held fixed and the remaining parameters are allowed to vary. The energy of this structure was corrected with a CASMP2(4,4)/MG3S singlepoint calculation, using an active space consisting of the two σ(PC) and σ *(PC) orbitals nearest the Fermi level. The rate constant for reaction 19 was then computed from conventional transition state theory using the BMK/MG3S geometry and scaled frequencies as before, and BMK/MG3S, CASSCF, or CASMP2 electronic energies, for the reactant and transition state. These calculations were repeated for several distances in the carbonphosphorus internuclear separation in order to find an approximate maximum in the generalized free energy, thereby obtaining a rough approximation to the variational value for the rate constant.

’ RESULTS The BMK/MG3S-optimized first-order saddle points of reactions 118 are shown as ball-and-stick models in Figures 1, 2, and 3. For saddle points with multiple conformations, the structure on the reaction path originating from the lowest energy conformation of the reactant is shown. All reactions in the test set except 2 and 10 have onefold reaction path degeneracy. Reaction 2 has twofold reaction path degeneracy since the reactant, thietane, has Cs symmetry in its ground-state equilibrium geometry (not shown), whereas the saddle point has C1 symmetry. The computed rate constant and A factor for this reaction are therefore weighted by a statistical factor of 2. Conversely, reaction 10 has C1 and Cs symmetry in the reactant and saddle point equilibrium geometry, respectively,

so the computed rate constant and A factor of this reaction is weighted by a factor of one-half. Classical barrier heights and reaction energies evaluated with BMK, CBS-QB3, CASSCF, or CASMP2 energies are shown in Table 1. The theoretical and experimental rate constants and Arrhenius parameters of each reaction are shown in Table 2. The rate laws are compared at 700 K, the approximate mean of the experimental temperatures. Comparisons for higher temperatures are shown in Table 3. Kinetic and thermochemical data may also be calculated using the BMK/MG3S electronic energies, equilibrium geometries, and unscaled harmonic frequencies of the reactants, saddle points, and products of reactions 118, which are included as Supporting Information. In what follows, separate consideration is given to certain reactions or groups of reactions in the test set, followed by an assessment of the accuracy of CVT/SCT and VTST-ISPE/SCT over the entire test set. Reaction 2. Elimination of ethylene from thietane has at least five pathways: three in which bonds are broken in a more or less concerted manner, and two in which bonds are broken in a highly sequential manner (but still occurring in a single step so that no intermediate is formed). Saddle point geometries on the concerted pathways and the high-energy “sequential” pathway are shown in Figure 4. At BMK/MG3S, the lowest pathway involves concerted bond breaking, with active CS and CC bond lengths of 2.47 and 1.82 Å, respectively, at the saddle point (structure 2b in Figure 4). The classical reaction barrier height of this saddle point is 69.7 kcal/mol. The lowest sequential pathway has active CS and CC bond lengths of 3.03 and 1.65 Å at the BMK/MG3S saddle point (structure 2 in Figure 1), which amounts to a 75.4 kcal/mol reaction barrier. However, upon application of the CBS-QB3 energy correction, the classical barrier height of structure 2 becomes 49.3 kcal/mol, while structures 2b, 2c, 2d, and 2e in Figure 4 amount to barriers of 70.0, 91.0, 74.2, and 66.5 kcal/mol, respectively. Structure 2 was therefore used to compute the kinetics of reaction 2, and the other pathways were assumed to be kinetically unimportant in the range of experimental temperatures (9801040 K). As Table 2 shows, the CVT/SCT rate constant for reaction 2 is k∞ = 4.94 109 s1 at 700 K. (Using pathway 2b in Figure 4, a value of k∞ = 1.53 107 s1 is obtained.) In contrast, the VTST-ISPE/SCT rate constant is k∞ = 0.15 s1, which is in good agreement with the experimental value of k∞ = 0.084 s1. Reaction 2 is also unusual among the present test subjects in that both Arrhenius parameters are fairly well-reproduced by theory. 14146



ARTICLE

Figure 3. BMK/MG3S first-order saddle points of reactions 1218.

The VTST-ISPE/SCT rate constants yield Ea = 48.2 kcal/mol and A = 1.7 1014 s1, whereas the experimental values are Ea = 48.3 kcal/mol and A = 1.0 1014 s1. This provides additional support for the correctness of the assumed pathway, which primarily involves CS bond breaking on the reactant side of the MEP, followed by CC bond breaking on the product side, as implied by saddle point structure 2 in Figure 1. Reaction 5. Elimination of propylene by chloromethyl allyl sulfide (CH2dCHCH2SCH2Cl) has two distinguishable pathways originating from conformations that differ in structure primarily by a 180° rotation in the CH2Cl group. The two conformations differ in BMK/MG3S total energy by 2.2 kcal/mol, so the more stable conformer comprises at least 84% of the Boltzmann population at experimental temperature (613673 K). The activation free energy at 298.15 K, ΔGq(298), is 40.8 and 41.0 kcal/mol, respectively, for the lower and higher energy conformer, while the conformational dependence of the rotational constants and vibrational frequencies is also small. The rate constant of reaction 5 was therefore calculated only for the

pathway originating from the lower energy conformer. Values of k∞ = 1.51 s1 and 46.3 s1 at 700 K are obtained for reaction 5 using CVT/SCT and VTST-ISPE/SCT, respectively, and the former result is in good agreement with the value of 0.989 s1 obtained from the observed rate law. Reaction 6. Elimination of ethylene by ethyl dithioacetate (CH3C(S)SC2H5) may have a four- or six-center transition state, corresponding to hydrogen transfer from the β-carbon of the ethoxy group to a mono- or dicoordinate sulfur atom, respectively, as both mechanisms yield the same products. The BMK/ MG3S energy barrier at the saddle point is 21.5 kcal/mol larger for the four-center pathway, so reaction 6 was assumed to have a six-center transition state, and the kinetics of the four-center pathway was neglected. Values of k∞ = 0.0454 s1 and 0.424 s1 at 700 K are obtained for reaction 6 from CVT/SCT and VTSTISPE/SCT, respectively, compared to the value of 0.0703 s1 from the observed rate law. Reactions 8, 9, and 10. The ethylene elimination reactions of the carbonates also have two pathways each, originating from 14147



ARTICLE

conformations that differ in structure by a 180° rotation around the SCSCH3 or SCOCH3 dihedral angle. This angle is cis in the BMK/MG3S global minimum conformation of the reactant of reactions 8 and 9 and trans in that of reaction 10. The difference in conformer total energies is 1.7, 3.4, and 4.5 kcal/mol for reactions 8, 9, and 10, so the global minimum conformation comprises over 78%, 92%, and 94% of the total Boltzmann population in the experimental temperature range of reactions 8, 9, and 10. Unfortunately, however, ΔGq(298) has a significant conformational dependence for these reactions; it is larger on the reaction path originating from the global minimum conformation— by 1.8, 1.8, and 3.8 kcal/mol for reactions 8, 9, and 10, respectively. The rate constants for these reactions were therefore computed for both the cis and trans configuration, and the thermal rate constant was obtained from the Boltzmann average. For these reactions the CVT/SCT rate constants are widely ranging in accuracy, being 6 and 62 times larger than experiment in the case of reactions 10 and 9, respectively. The corresponding VTST-ISPE/SCT rate constants are too large by a factor of 14 and 133. Reaction 11. Unimolecular decomposition of 1-chloro4-(methylthio)-butane involves a two-step mechanism wherein dissociation of the chlorine atom to form a diradical intermediate is followed by elimination of methyl chloride and ring closing to form tetrahydrothiophene.25 The saddle point geometry of the first and second step of this reaction is shown in Figure 2, and the energy of both barriers, the intermediate, and the overall reaction is shown in Table 1. At BMK/MG3S, the first step amounts to a 62.0 kcal/mol barrier, the intermediate lies at 18.8 kcal/mol, and

Table 1. Classical Barrier Heights and Reaction Energies (kcal mol1) ΔErxn

Vq rxn

BMK

corrected

a

correcteda

1

54.2

50.9

14.2

12.5

2

75.4

49.3

36.8

29.2

3 4

64.1 45.7

64.9 41.2

21.4 27.9

22.3 25.8

5

43.4

38.5

22.1

22.0

6

52.0

49.0

30.3

29.7

7

58.3

55.8

32.6

30.3

8

47.2

46.3

18.1

19.2

9

47.3

46.0

19.3

18.9

10

58.4

57.6

38.0

37.5

11

62.0 64.2

74.5 79.8

18.8 3.1

12.6 3

12

48.8

46.8

18.1

17.9

13

48.1

46.1

19.1

19.0

14

58.1

55.2

8.0

6.2

15

48.4

46.2

22.8

20.7

16

48.1

45.9

26.3

28.9

17

45.5

43.3

19.2

22.8

18

73.6

73.9

49.1

51.3

68.3

72.4

19 a

BMK

CBS-QB3 or CASMP2 (italicized).

Table 2. Kinetic Parameters for the Thermal Decomposition of Phosphorus and Sulfur Compounds Ea (kcal mol1) temp

CVT/ VTST-ISPE/

CVT/

VTST-ISPE/

rxn

(K)

SCT

SCT

exptl

SCT

SCT

exptl

SCT

SCT

1

678704

52.1

48.6

45.8

33.0

32.9

29.7

1.20 102

1.33 101

3.98 102c

0.3

3.4

9

1

p

CVT/ VTST-ISPE/

8.43 102d

6 108

1.8

exptl

2

9801040

72.5

48.2

48.3

33.0

32.7

32.2

4.94 10

1.51 10

3

9501230

60.7

61.0

55.0

34.4

34.3

30.6

9.28 105

6.93 105

1.34 104e

kcalcd/kexptla

kcalcd/kexptlb

0.7

0.5

4

649691

41.9

37.1

38.2

28.4

28.4

25.9

1.80 10

5.64 10

1.95 101f

0.9

28.9

5 6

613673 651716

39.1 47.3

34.3 44.1

34.4 45.2

28.5 30.9

28.5 30.9

24.7 29.8

1.51 100 4.54 102

4.63 101 4.24 101

9.87 101g 7.03 102h

1.5 0.6

47.0 6.0

7

763841

53.6

50.9

50.4

30.3

30.2

28.6

2.63 104

1.57 103

4.60 104i

0.6

3.4

1

3.53 101j

34.0

77.0

2.35 101j

61.8

133.0

2.86 104k

5.9

1

0

8

590620

40.6

39.6

39.7

31.7

31.8

27.5

1.20 10

2.71 10

9

570660

39.8

38.8

35.7

31.3

31.4

24.2

1.45 101

3.13 101

1

3

3

10

724775

52.6

51.4

50.7

31.4

31.4

28.3

1.67 10

3.89 10

11

583623

59.3

75.9

41.8

29.6

29.6

28.2

2.23 106

9.40 1011 1.54 101l

12

802907

42.9

40.6

45.3m

30.0

29.6

29.9m 4.24 101

1.54 100

7.17 102m,n

1

13

a

k(T = 700 K) (s1)

ln(As)

13.6

1 105

6.00 1010

5.9

21.5

775854 840940

42.6 43.0

40.5 40.9

47.4 44.8

30.6 30.8

30.6 30.8

32.2 29.9

9.83 10 9.43 101

4.39 10 4.19 100

1.58 101o 1.03 101p

6.2 9.2

27.7 40.8

700900

42.6

40.4

45.3m

30.6

30.5

31.0m 9.93 101

4.45 100

2.01 101m,q

5.0

22.2

1

0

570650

41.7

39.5

41.3

29.9

29.9

35.4

9.86 10

4.40 10

3.05 102r

0.003

0.01

14

620780

52.7

49.9

37.0

29.9

30.0

15.3

3.45 104

2.61 103

1.20 105r

28.7

217.1

15

673873

40.3

38.3

40.1

29.8

29.8

27.2

2.24 100

9.99 100

1.99 101r

11.3

50.2

7.31 102s

1.6

6.4

7.92 101s

3.8

14.2

1

0

1

16

729772

43.6

41.4

34.4

29.3

29.0

22.1

1.20 10

4.68 10

17

652694

37.4

36.1

34.2

28.0

28.3

24.3

2.98 100

1.12 101

18 19

600800 648723

69.0 98.1

69.5 34.9

63.0

31.7 42.8

31.8 36.7

40.3

1.74 108 1.22 108 8.75 1013 1.09 105

d

e

b

c

f

g

h

j

6.30 103t k

l

1 1010 m

2 107

CVT/SCT. VTST-ISPE/SCT. Ref 9. Ref 28. Ref 29. Ref 30. Ref 31. Ref 32. Ref 33. Ref 34. Ref 35. Ref 24. Estimate. Ref 7. o Ref 1. Ref 8. q Ref 6. r Ref 2. s Ref 36. t Ref 24. 14148

i

n



ARTICLE

Table 3. Rate Constants (s1) for the Thermal Decomposition of Phosphorus and Sulfur Compounds k(T = 1000 K) rxn

CVT/SCT

CVT/SCT

1.68 10

3.3

4.44 10

7.85 10

5.7

9.24 106

3 104

1.29 106

5.30 108

2 103

5

7

CVT/SCT

exptl

1.2

5.63 10

6 106

2.94 103

1

1

9.06 10

7.76 10

1.53 102

2.81 103

2

6

k(T = 2000 K) kcalcd/kexptl

kcalcd/kexptl

exptl

2

2

k(T = 1500 K)

6

8

kcalcd/kexptl

exptl 7

3

4.50 10

1.90 10

2.4

1.19 10

1.93 10

6.2

1.93 10

1.95 10

9.9

4

1.51 103

7.45 102

2.0

1.70 106

4.55 105

3.7

5.72 107

1.12 107

5.1

5

6.94 103

1.65 103

4.2

4.90 106

5.31 105

9.2

1.30 108

9.53 106

13.7

6

1.22 103

1.19 103

1.0

3.41 106

2.32 106

1.5

1.80 108

1.03 108

1.8

7

2.75 101

2.42 101

1.1

2.20 105

1.14 105

1.9

1.97 107

7.80 106

2.5

8 9

7.66 104 7.80 104

1.86 103 5.16 102

41.3 151.3

6.99 107 6.21 107

1.46 106 2.05 105

48.0 303.5

2.11 109 1.75 109

4.08 107 4.08 106

51.8 429.8

10

2.81 102

1.62 101

17.4

1.90 106

8.04 104

23.6

1.56 108

5.67 106

27.5

6 10

1.70 10

4

1.42 10

1 10

2.46 10

6

4.75 107

5 102

6

6

8

8

1

1

11

8.06 10

1.28 10

3 3

4

6

6

2

8

12

4.43 10

1.26 10

3.5

5.91 10

2.51 10

2.4

2.16 10

1.12 10

1.9

13

9.70 103

4.36 103

2.2

1.24 107

1.24 107

1.0

4.44 108

6.61 108

0.7

14

3.00 101

3.52 102

853.1

2.09 105

1.75 101

11943.0

1.74 107

3.89 102

44685.0

15

1.35 104

1.13 103

11.9

1.17 107

9.47 105

12.3

3.45 108

2.74 107

12.6

16 17

1.47 103 9.59 103

1.22 102 1.26 103

12.0 7.6

2.22 106 5.12 106

3.94 104 3.89 105

56.5 13.2

8.65 107 1.18 108

7.06 105 6.84 106

122.5 17.3

19

1.36 103

5.12 103

3 107

1.92 104

2.02 108

9 105

7.20 107

4.02 1010

2 103

3

the second step amounts to a 64.1 kcal/mol barrier relative to the reactant. The present implementation of CASMP2 changes the energy of these stationary points to 74.5, 12.6, and 79.8 kcal/mol relative to the reactant. The rate constant of reaction 11 was computed from the MEP of the second step, and the influence of the first step and intermediate on the kinetics was neglected. Table 2 shows that the CVT/SCT and VTST-ISPE/SCT rate constants for reaction 11 are too small by 5 and 10 orders of magnitude, respectively. The CVT/SCT and VTST-ISPE/SCT activation energies both greatly overestimate experiment. Reaction 12. Elimination of ethylene from dimethyl ethyl phosphonate (PO(CH3)(OC2H5)2) may have a four- or sixcenter transition state, corresponding to hydrogen transfer from the β-carbon of the ethoxy group to a mono- or dicoordinate oxygen atom. The six-center transition state was assumed, and the kinetics of the four-center transition state were neglected, based on the ∼20 kcal/mol difference in the size of the BMK/ MG3S energy barriers of these two mechanisms. The saddle point of the six-center transition state was found to have three conformations, reflecting different equilibrium configurations of the free ethyl substituent. The thermal rate constant for this reaction was therefore obtained from the Boltzmann average of the three microscopic rate constants, owing to a conformational dependence in ΔGq(298). Table 2 compares theoretical kinetic parameters to estimated values obtained by Glaude et al.6 If a comparison of the CVT/ SCT rate constant is made to the experimental results of Zegers and Fisher,7 one obtains kcalcd/kexptl = 10.3 at 800 K, but much poorer agreement is obtained when the comparison is made slightly outside the experimental temperature range, e.g., kcalcd/ kexptl = 44.7 at 700 K. Thus, the CVT/SCT Arrhenius parameters, A = 1.1 1013 s1 and Ea = 42.9 kcal/mol, repudiate the experimental values of A = 2.9 1016 s1 and Ea = 59.2 kcal/mol obtained by Zegers and Fisher. Glaude et al. assumed a value of A = 1.0 1013 s1 for reaction 12, based on the experimental A factor for ester decomposition, which is similar to reaction 12

in its rate and also proceeds through a six-center transition state. Using this value for A, along with the experimental rate constant at 800 K, a value of Ea = 45.3 kcal/mol is then obtained from the Arrhenius equation, in much better agreement with CVT/SCT. It follows from CVT/SCT theory and the estimated rate law of Glaude that kcalcd/kexptl = 5.9 for reaction 12 at 700 K. Reactions 13, 14, and 15. Reaction 13, the elimination of ethylene from triethyl phosphate (PO(OC2H5)3), was assumed to have a six-center transition state. This reaction, as well as the elimination of ethylene from triethyl phosphite (P(OC2H5)3) or the VX simulant O,S-diethyl ethylphosphonothioate (PO(C2H5)(OC2H5)(SC2H5)) (reactions 14 and 15, respectively), has at least four pathways, so the thermal rate constant for these reactions was approximated as the microscopic rate constant with the greatest statistical weight, without justification from relative conformer populations or ΔGq(298) values. In Table 2, the computed kinetic parameters for reaction 13 are compared to results from three experimental studies. MacNaughton et al., working in the 570650 K temperature range, observed quasi-first-order kinetics because of catalytic activity in a reaction byproduct.2 A very short exposure time in the reactor was therefore needed in order to obtain a stable reaction rate, and this compromised the reproducibility of their results. Zegers and Fisher observed the reaction over a temperature range of 706854 K,1 but a second process below 775 K was implied by the Arrhenius plot of their data. Consequently, their rate law was extracted from a short temperature range, which diminishes its precision. Tsang used shock tube measurements of the rate constant over a slightly greater temperature range of 840940 K to obtain a rate law.8 Glaude et al.6 estimated the rate law from an assumed value of Ea = 45.3 kcal/mol, based on the estimated Ea of reaction 12 and Zegers and Fisher’s experimental rate constant at 800 K. The results of Zegers, Tsang, and Glaude are supported by the CVT/SCT calculation, i.e., kcalcd/kexptl = 6, 9, and 5, respectively, at 700 K. The computed activation energy underpredicts the Ea 14149



ARTICLE

Figure 4. BMK/MG3S first-order saddle points of reaction 2.

values from these studies by 25 kcal/mol, and the computed A factor falls in between Tsang’s value and Glaude’s value. The experimental difficulties encountered by MacNaughton et al. are reflected in results that are much less consistent with preceding studies and the present calculations. From CVT/SCT and the MacNaughton study one obtains kcalcd/kexptl = 0.003 at 700 K. The MacNaughton and CVT/SCT activation energies agree to within 0.5 kcal/mol, but the corresponding A factors differ by a factor of 245, suggesting systematic error in the observed rate constant. MacNaughton’s study of reaction 14 was also affected by catalytic activity in the reaction byproducts, but to a lesser extent than reaction 13.2 The CVT/SCT rate constant for reaction 14 overestimates experiment such that kcalcd/kexptl = 28.7 at 700 K. The experimental A factor for reaction 14 has an unusually small value of A = 4.3 106 s1, compared to the CVT/SCT-derived value of A = 1.01 1013 s1. Reactions involving five-center transition states, such as reaction 14, typically have A factors between 1011 and 1015 s1.26 MacNaughton’s study of reaction 15 encountered no experimental difficulties,2 and in this case there is comparatively good agreement with CVT/SCT theory (kcalcd/kexptl = 11.3 at 700 K). The MacNaughton and CVT/SCT

activation energies of reaction 15 agree to within 0.2 kcal/mol, but the corresponding A factors differ by 1 order of magnitude. Reaction 17. Due to conformations of the free allyl group, elimination of propylene from diallylphenylphosphine has three pathways. The thermal rate constant for this reaction was obtained from the Boltzmann average of the three microscopic rate constants, owing to a conformational dependence in ΔGq(298). Values of k∞ = 2.98 s1 and 11.2 s1 at 700 K are obtained for reaction 17 with CVT/SCT and VTST-ISPE/SCT, respectively, compared to a value of 0.792 s1 from the experimental rate law. Reaction 18. Thermal decomposition of the sarin simulant O-ethyl methylphosphonofluoridate apparently has not been the subject of experimental study. Reaction 18, elimination of ethyl chloride from O-ethyl methylphosphonofluoridate, was therefore selected for computational study on the basis of its BMK/MG3S activation free energy in comparison to the following reactions:

14150

POðCH3 ÞðOC2 H5 ÞðFÞ f POðCH2 ÞðOC2 H5 Þ þ HF

ð18bÞ

POðCH3 ÞðOC2 H5 ÞðFÞ f cyc POCH2 CH2 ðOÞðCH3 Þ þ HF

ð18cÞ

POðCH3 ÞðOC2 H5 ÞðFÞ f POðOC2 H5 Þ þ CH3 F

ð18dÞ


The Journal of Physical Chemistry A Values of ΔGq(298) = 70.3, 74.0, 90.7, and 108.3 kcal/mol are obtained, respectively, for reactions 18, 18b, 18c, and 18d. The CVT/SCT Arrhenius parameters for reaction 18 are A = 5.8 1013 s1 and Ea = 69.0 kcal/mol. Reaction 19. No saddle point could be located for the dissociation of phosphino from tert-butyl phosphine at the BMK/ MG3S level of theory. A maximum in the CBS-QB3//BMK/ MG3S potential energy is found at a tertiary carbonphosphorus internuclear separation of 3.5 Å, but this result is unphysical in that the height of this barrier is smaller than the CBS-QB3// BMK/MG3S reaction energy. The CCSD(T) base energy of the CBS-QB3 method is known to fail when internuclear separations are sufficiently large, because of the omission of higher order clusters from the wave function or because of divergence by the perturbative expansion used to estimate the contribution from connected triples.27 The correction to the reactant and generalized transition state energies of reaction 19 was therefore performed with the CASMP2 method, and no saddle point or unphysical result is associated with PH2 dissociation at this level of theory. Unfortunately, however, the kinetics of reaction 19 is not reproduced with the present model. At 700 K, the generalized transition state for reaction 19 occurs at a phosphorustertiary carbon separation of RPC ∼ 5.2, 3.6, and 4.0 Å according to the BMK, CASSCF(4,4), and CASMP2(4,4) potential energies, respectively. The computed rate constant is in error such that kcalcd/kexptl ∼ 1010, 106, and 107 at 700 K, respectively, with the BMK, CASSCF(4,4), and CASMP2(4,4) potential energies. However, increasing the size of the active space from four electrons in four orbitals to 10 in 10 is found to significantly decrease the CASSCF and CASMP2 energy difference between the reactant and transition state. At RPC = 4.0 Å, the decrease is by about 9 kcal/mol. Using the CASMP2(10,10) energies, the rate constant computed at RPC = 4.0 Å 700 K therefore improves so that kcalcd/kexptl = 3 104 at 700 K. CASMP2 calculations with an even larger active space were not pursued in the present study due to the computational cost.

’ SUMMARY AND CONCLUSIONS Dual-level direct dynamics was applied to a set of 19 unimolecular reactions of phosphorus or sulfur molecules, 18 of which have observed rate laws. Information regarding the minimum energy paths of these reactions was obtained at the BMK/MG3S level, and corrections to the potential energy along these MEPs were taken from CBS-QB3, CASSCF, or CASMP2 calculations. Thermal rate coefficients were computed over experimental temperatures with canonical variational transition state theory or variational transition state theory with interpolated singlepoint energies, in conjunction with small-curvature tunneling theory. All reactions in the test set with a single pathway and little radical character in the transition state are well-reproduced by the CVT/SCT level of theory. More specifically, the CVT/SCT rate constants for reactions 1, 37, and 16 all fall within a factor of 3 of experiment at 700 K. Activation energies derived from these rate constants overestimate experiment by 25 kcal/mol, whereas the corresponding A factors are at least 1 order of magnitude larger than experiment in most cases. The comparatively good accuracy of CVT/SCT rate constants in these applications therefore results from a partial cancellation between the errors in these two parameters.

ARTICLE

CVT/SCT rate constants for reactions of molecules with more than one conformation are widely varying in accuracy. In cases such as reaction 17, the level of agreement with experiment is similar to that obtained for the single-conformation molecules, but in the case of reaction 9, the CVT/SCT rate constant is 60 times too large. The VTST-ISPE/SCT rate constants in Table 2 typically overshoot experiment and are less accurate than the CVT/SCT rate constants. The implementation of VTST-ISPE with CBSQB3 energies improves the accuracy of the rate constant only for reaction 2, which is unique among the present test subjects in that it has significant radical character in the transition state. The VTST-ISPE/SCT derived activation energies are usually more accurate than those of CVT/SCT and tend to overestimate experiment by 12 kcal/mol. (For reactions 12 and 13, Ea is more accurately computed with CVT/SCT, but these reactions only break bonds between first row elements.) Use of VTSTISPE for elimination reactions in which bonds to phosphorus or sulfur are broken therefore cannot be recommended. Dual-level methods that correct vibrational frequencies and moments of inertia along the MEP, as well as energies, are perhaps needed in order to improve on CVT/SCT rate constants.

’ ASSOCIATED CONTENT

bS

Supporting Information. Cartesian coordinates, total energies, and harmonic frequencies for reactants, products, and saddle points. This material is available free of charge via the Internet at http://pubs.acs.org.

’ AUTHOR INFORMATION Corresponding Author

*E-mail: [email protected].

’ ACKNOWLEDGMENT Computational resources were provided by the Alabama Supercomputer Center and the Arctic Region Supercomputing Center. This work was supported by the Defense Threat Reduction Agency. ’ REFERENCES (1) Zegers, E. J. P.; Fisher, E. M. Combust. Flame 1998, 115, 230–240. (2) MacNaughton, M. G.; Piccini, P. R.; Shannon, K. L.; Carpenter, M. A.; Crumlett, S.; Bohmann, J. A. Experimental Determination of Thermal Neutralization Parameters for Chemical Agents, Final Report— Thermal Decomposition of CWA; SwRI Project 13161; Southwest Research Institute: San Antonio, TX, 2010. (3) Hahn, D. K.; RaghuVeer, K. S.; Ortiz, J. V. J. Phys. Chem. A 2010, 114, 8142–8155. (4) Hahn, D. K.; RaghuVeer, K. S.; Ortiz, J. V. J. Phys. Chem. A 2011, 115, 8532–8539. (5) Pitzer, K. S.; Gwinn, W. D. J. Chem. Phys. 1942, 10, 428–440. (6) (a) Glaude, P. A.; Curran, H. J.; Pitz, W. J.; Westbrook, C. K. Proc. Combust. Inst. 2000, 28, 1749–1756. (b) Glaude, P. A.; Melius, C.; Pitz, W. J.; Westbrook, C. K. Detailed Chemical Reaction Mechanisms for Combustion and Inceneration of Organophosphorus and Fluoro-Organophosphorus Compounds. Proc. Combust. Inst. 2002, 29, 2469–2476. (7) Zegers, E. J. P.; Fisher, E. M. Combust. Sci. Technol. 1996, 116, 69–89. (8) Tsang, W. Combust. Sci. Technol. 1998, 138, 85–103. 14151


The Journal of Physical Chemistry A (9) Jeffers, P.; Dasch, C.; Bauer, S. H. Int. J. Chem. Kinet. 1973, 5, 545–551. (10) Zheng, J.; Zhao, Y.; Truhlar, D. G. J. Chem. Theory Comput. 2009, 5, 808–821. (11) Curtiss, L. A.; Raghavachari, K.; Redfern, P. C.; Rassolov, V.; Pople, J. A. J. Chem. Phys. 1998, 109, 7764–7776. (12) Boese, A. D.; Martin, J. M. L. J. Chem. Phys. 2004, 121, 3405–3416. (13) Zhao, Y.; Truhlar, D. G. Theor. Chem. Acc. 2008, 120, 215–241. (14) Lynch, B. J.; Truhlar, D. G. In Electron Correlation Methodology; Wilson, A. K., Peterson, K. A., Eds.; ACS Symposium Series; American Chemical Society: Washington, DC, 2007; Vol. 958, pp 153167. (15) Fernandez-Ramos, A.; Miller, J. A.; Klippenstein, S. J.; Truhlar, D. G. Chem. Rev. 2006, 106, 4518–4584. (16) Zheng, J.; Zhang, S.; Lynch, B. J.; Corchado, J. C.; Chuang, Y.-Y.; Fast, P. L.; Hu, W.-P.; Liu, Y.-P.; Lynch, G.-C.; Nguyen, K. A.; Jackels, C. F.; Ramos, A. F.; Ellingson, B. A.; Melissas, V. S.; Villa, J.; Rossi, I.; Coiti~no, E. L.; Pu, J.; Albu, T. V.; Steckler, R.; Garrett, B. C.; Isaacson, A. D.; Truhlar, D. G. POLYRATE 2010-A; University of Minnesota: Minneapolis, MN, 2010. (17) Lynch, B. J.; Zhao, Y.; Truhlar, D. G. J. Phys. Chem. A 2003, 107, 1384–1388. (18) 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, G. A.; Nakatsuji, H.; Caricato, M.; Li, X.; Hratchian, H. P.; Izmaylov, A. F.; Bloino, J.; Zheng, G.; Sonnenberg, J. L.; Hada, M.; Ehara, M.; Toyota, K.; Fukuda, R.; Hasegawa, J.; Ishida, M.; Nakajima, T.; Honda, Y.; Kitao, O.; Nakai, H.; Vreven, T.; Montgomery, J. A., Jr.; Peralta, J. E.; Ogliaro, F.; Bearpark, M.; Heyd, J. J.; Brothers, E.; Kudin, K. N.; Staroverov, V. N.; Kobayashi, R.; Normand, J.; Raghavachari, K.; Rendell, A.; Burant, J. C.; Iyengar, S. S.; Tomasi, J.; Cossi, M.; Rega, N.; Millam, J. M.; Klene, M.; Knox, J. E.; Cross, J. B.; Bakken, V.; Adamo, C.; Jaramillo, J.; Gomperts, R.; Stratmann, R. E.; Yazyev, O.; Austin, A. J.; Cammi, R.; Pomelli, C.; Ochterski, J. W.; Martin, R. L.; Morokuma, K.; Zakrzewski, V. G.; Voth, G. A.; Salvador, P.; Dannenberg, J. J.; Dapprich, € Foresman, J. B.; Ortiz, J. V.; Cioslowski, J.; S.; Daniels, A. D.; Farkas, O.; Fox, D. J. Gaussian 09, revision A.02; Gaussian Inc.: Wallingford, CT, 2009. (19) Melissas, V. S.; Truhlar, D. G.; Garrett, B. C. J. Chem. Phys. 1992, 96, 5758–5772. (20) Garrett, B. C.; Truhlar, D. G. Variational Transition State Theory. In Theory and Applications of Computational Chemistry: The First Forty Years; Dykstra, C. E., Frenking, G., Kim, K., Scuseria, G., Eds.; Elsevier: Amsterdam, The Netherlands, 2005; pp 6787. (21) Lu, D.-h.; Truong, T. N.; Melissas, V. S.; Lynch, G. C.; Liu, Y.-P.; Garrett, B. C.; Steckler, R.; Isaacson, A. D.; Rai, S. N.; Hancock, G.; Lauderdale, T. J.; Truhlar, D. G. Comput. Phys. Commun. 1992, 71, 235–262. (22) Montgomery, J. A., Jr.; Frisch, M. J.; Ochterski, J. W.; Petersson, G. A. J. Chem. Phys. 2000, 112, 6532–6542. (23) Hegarty, D.; Robb, M. A. Mol. Phys. 1979, 38, 1795–1812. (24) Li, S. H.; Larson, C. A.; Buchan, N. I.; Stringfellow, G. B. J. Electron. Mater. 1989, 18, 457–464. (25) Chuchani, G.; Martin, I.; Dominguez, R. M. Int. J. Chem. Kinet. 1987, 19, 683–690. (26) Gilbert, R. G.; Smith, S. C. Theory of Unimolecular and Recombination Reactions; Blackwell Scientific Publications: Oxford, U.K., 1990. (27) Piecuch, P.; Kowalski, K.; Fan, P.-D.; Pimienta, I. S. O. In Advanced Topics in Theoretical Chemical Physics; Maurani, J., Lefebvre, R., Br€andas, E., Eds.; Kluwer: Dordrecht, The Netherlands, 2003; pp 119206. (28) Bigley, D. B.; Gabbott, R. E. J. Chem. Soc., Perkin Trans. 2 1975, 317–320. (29) Tsang, W. J. J. Chem. Phys. 1964, 40, 1498–1505. (30) Martin, G.; Ropero, M.; Avila, R. Phosphorous, Sulfur Relat. Elem. 1982, 13, 213–220. (31) Martin, G.; Martinez, H.; Suhr, H.; Suhr, U. Int. J. Chem. Kinet. 1986, 18, 355–362. (32) Al-Awadi, N.; Bigley, D. B.; Gabbott, R. E. J. Chem. Soc., Perkin Trans. 2 1978, 1223–1224.

ARTICLE

(33) Oele, P. C.; Tinkelenberg, A.; Louw, R. Tetrahedron Lett. 1973, 23, 2375–2378. (34) Al-Awadi, N.; Bigley, D. B. J. Chem. Soc., Perkin Trans. 2 1982, 773–775. (35) Taylor, R. J. Chem. Soc., Perkin Trans. 2 1983, 291–296. (36) Martin, G.; Ocando-Mavarez, E.; Osorio, A.; Laya, M.; Canestrari, M. Heteroat. Chem. 1992, 3, 395–401.

14152


Computational Tests of Models for Kinetic Parameters of Unimolecular

Recommend Documents