J. Phys. Chem. A 2010, 114, 10717–10725
10717
Isomerization and Decomposition of a Model Nerve Agent: A Computational Analysis of the Reaction Energetics and Kinetics of Dimethyl Ethylphosphonate Debasish Mandal, Bhaskar Mondal, and Abhijit K. Das* Department of Spectroscopy, Indian association for the CultiVation of Science, JadaVpur, Kolkata 700 032, India ReceiVed: July 7, 2010; ReVised Manuscript ReceiVed: August 24, 2010
The gas-phase isomerization and decomposition reactions of dimethyl ethylphosphonate (DMEP) are investigated using the CBS-QB3 method followed by the calculation of rate constant for all reaction pathways using Rice-Ramsperger-Kassel-Marcus (RRKM) theory. Three conformational isomers C1, C2, and C3 are identified for DMEP having stability order C1 (0.00) < C2 (0.11) < C3 (1.87). Each conformer can isomerize via H-transfer reaction and can decompose via CH3OH, CH2, and H2 elimination at higher temperatures. The conformers C1, C2, and C3 can isomerize to IM1, IM2, and IM3, respectively, via different pathways. Decomposition requires much higher activation energy and therefore higher temperatures than the corresponding isomerization. For instance, the most stable conformer C1 isomerizes to IM1 twice as fast as decomposing to P1 + CH2 at 1000 K whereas the least stable conformer C3 isomerizes to IM3 104 times faster than decomposing to P5 + CH3OH. Only one decomposition channel is identified for C1 and two different decomposition channels are identified for C2 as well as for C3. The decompositions of C2 and C3 to P2 + CH3OH and P5 + CH3OH, respectively, are predicted to be more favorable thermodynamically as well as kinetically over the other decomposition channels within the temperature range 1000-3000 K. For the lack of experimental data, we have calculated the low as well as high-pressure limit rate constants for the decomposition reactions of DMEP. In addition, consistent and reliable enthalpies of formation at 298.15 K (∆fH°298.15) have been computed for all the species involved in the isomerization and decomposition reaction of DMEP. The results obtained for DMEP for the reaction mechanism and energetics are compared with that for DMHP and DMMP. 1. Introduction Decomposition of organophosphorus compounds (OPCs) in flames is the underlying chemical process for the disposal of toxic and hazardous chemical wastes, including the disposal of chemical warfare agents (CWAs) and solvent extraction processes.1-9 The organophosphorus materials are extremely toxic nerve and blistering agents, and therefore, these compounds are very dangerous to handle in the laboratory. Less toxic OPCs such as dimethyl hydrogenphosphonate (DMHP), dimethyl methylphosphonate (DMMP), and dimethyl ethylphosphonate (DMEP) are used to simulate the chemical properties of CWAs.10-12 The organophosphates are also used as plasticizers, flame retardants, and lubricants.13,14 A large number of experimental and theoretical investigations have been carried out on these compounds to study their conformational changes, thermochemical properties, gas-phase pyrolysis, and pyrolysis in H2/O2 flame. High-level thermochemical analysis of OPCs, like DMMP, can be found in the literature.15,16 Matrix isolation infrared and theoretical studies were performed to analyze the conformers of DMHP.17 Conformations of DMMP in liquid and gas phases and also in argon, nitrogen, and carbon monoxide matrices were investigated by Barnes et al.18 Recently. Yang et al. performed an excellent theoretical study to investigate the gas-phase unimolecular decomposition of DMMP.19 Kinetic models have also been developed for DMMP to study their decomposition in H2/O2 flame.20 The catalytic decomposition of chemical warfare agents using metal and metal oxides is also an important area of research.21-23 * Corresponding author. E-mail:
[email protected].
DMEP is a typical organophosphorus toxic liquid like DMHP and DMMP. Despite this fact, only the atmospheric aspect of DMEP regarding its reaction with the OH radical and NO3 radical was investigated by Aschmann et al.24 and the mass spectrometric fragmentation reactions of protonated DMEP was elucidated by Bell et al.25 In these early works only the reaction of DMEP with radicals and its fragmentation have been investigated, but no attention has been paid to study its unimolecular decomposition in flame temperature, isomerization, and thermochemistry. As far our knowledge goes, we could not find such a study on DMEP in existing literature. In this article, we have presented a high-level computational study to build up the potential energy surface (PES) for the gas-phase unimolecular decomposition and isomerization pathways of DMEP including the thermochemical calculations by high-level CBS-QB3 composite quantum chemical method on the species involved in the reaction pathways. We have modeled the unimolecular isomerization and decomposition pathways of DMEP by analyzing the existing reaction pathways of other CWAs like DMHP, DMMP, VX, GB, etc. because all these CWAs are examined to be similar in their physical and chemical properties. Following the fragmentations and reactions of protonated DMEP in an ion trap mass spectrometer,25 we can assume the destruction pathways of CWAs via CH3OH elimination and the isomerization pathways via 1,3 and 1,4-H transfer. The destruction of CWAs via elimination of CH2, CO, and H2 was reported in literature.19 To have a clear understanding of the unimolecular reaction pathways of DMEP, the kinetics study is also carried out using RRKM master equation and thermodynamic consideration.
10.1021/jp106270d 2010 American Chemical Society Published on Web 09/10/2010
10718
J. Phys. Chem. A, Vol. 114, No. 39, 2010
2. Computational Methods All ab initio and DFT molecular orbital theory calculations have been performed using Gaussian 03 suite of quantum chemistry program.26 For electronic structure calculations we have used density functional B3LYP (DFT-B3LYP)27-29 method in conjunction with 6-311++G(2d,2p) basis set and MP230-33 method in conjunction with the 6-311++G(d,p) basis set for all the elements involved. All calculations have been performed on the low-lying singlet species with restricted wave function. The B3LYP method consists of Becke’s three-parameter hybrid exchange functional combined with Lee-Yang-Parr correlation functional. The MP2 is the second order Møller-Plesset perturbation theory which includes double excitations. The frozencore (FC) approximation is employed during MP2 calculations. The 6-311++G(2d,2p) basis set is a valance triple-ξ quality with double polarization and double diffuse functions on all atoms and 6-311++G(d,p) involves single polarization and single diffuse functions on all atoms. The connecting first order saddle points, that are the transition states between the equilibrium geometries, are obtained using synchronous transitguided quasi-Newton (STQN) method. The normal-mode analysis has been carried out at the same level of theories for equilibrium as well as transition state geometries and characterized as minima (no imaginary frequency) or as a transition state (one imaginary frequency). A parallel intrinsic reaction coordinate calculation (IRC) has been performed with all transition states to confirm whether these transition states connect the right minima or not.34 The high-accuracy energy prediction multilevel method CBS-QB335 is used to obtain more reliable energies which can be utilized for the calculations of thermodynamic parameters. The CBS-QB3 method uses B3LYP/CBSB7 geometries and vibartional frequencies with appropriate scaling for accurate single-point energy calculations. The frequencies used in CBS-QB3 method are scaled by a factor of 0.99, and these frequencies are further used for thermodynamic calculation. We have used the CBS-QB3 method as implemented in Gaussian 03 without any modification. A one-dimensional singlet potential energy surface (PES) is constructed with the CBS-QB3 relative energies calculated with respect to the lowest energy conformer C1. For the analysis of the PES and reaction energetics, CBS-QB3 calculated energies in kcal/mol are used throughout this work. The enthalpies of formation for all the species involved in the title isomerization and decomposition reaction at 298.15 K (∆fH°298.15) are calculated using CBS-QB3 electronic energies. The atomization scheme36 has been employed for the calculation of ∆fH°298.15, and for this purpose we have used the literature values of ∆fH°298.15 for C [171.29 kcal/mol], H [52.10 kcal/ mol], O [59.51 kcal/mol], and P [75.62 kcal/mol].37 The free energy changes for reaction (∆rG) are calculated as a function of temperature using MP2/6-311++G(d,p) electronic energies and the statistical mechanical principles with the contribution from translational, rotational, and vibrational motion. The vibrational and translational frequencies are modeled using the rigid rotor harmonic oscillator (RRHO) approximation with the frequencies calculated at the MP2/6-311++G(d,p) level. Low-frequency vibrations corresponding to the rotation of the methyl groups are omitted from the RRHO analysis and treated as hindered rotors. The potential energy scans about the rotors are carried out at the PM3 semiempirical level and the rotation barrier height is calculated on the basis of the maximum value of the potential energy curve. As the estimated barrier height is higher than 1.4RT at a given temperature, it is
Mandal et al. concluded that the rotation is hindered. The contribution from internal rotation is estimated using the MOLTRAN program.38 The unimolecular decomposition and isomerization rate constants are calculated as a function of temperature and pressure using RRKM theory with a time-independent solution of the master equation.39 Collision energy transfer is treated with exponential down model with 〈∆Edown〉 ) 500 cm-1 where argon is used as a bath gas. During the calculation of the rate constants the effect of quantum mechanical tunneling is taken into account using Eckart’s tunneling correction.40 We have calculated the barrier width (in amu0.5 Å) along the reaction coordinate by fitting the IRC curves with one-dimensional Eckart’s potential.
V(x) )
(
eu B A+ u 1+e 1 + eu
)
where u ) 2πx/l, A ) E1 - E-1, and B ) E11/2 + E-11/2. The constants E1 and E-1 represent the barrier heights relative to the reactants and products, respectively. Taking into account the experimental temperature and pressure for the title reaction, we have calculated low- and high-pressure limit rate constants for different reaction schemes at a particular temperature window. The limiting high-pressure rate constants, k∞, are calculated using canonical transition state theory (CTST) and the low-pressure-limiting expression of the RRKM master equation is used to calculate the limiting low-pressure rate constants, k0. The calculated low-pressure and high-pressure limit rate constants are fitted to a modified form of Arrhenius expression
( )
ki ) ATn exp -
Ea RT
and the Arrhenius parameters A, n, and Ea are calculated for all reaction channels. All the kinetic computations are performed with ChemRate code with MP2/6-311++G(d,p) geometries, frequencies,and CBS-QB3 enthalpies of formation.41 ChemRate determines moments of internal rotors based upon molecular structure and connectivity, and these are subsequently employed in evaluating the contribution of the internal rotor to the partition function of the molecule. 3. Results and Discussion The optimized geometries of all the reactants, transition states (TS), and products are displayed in Figure 1 with the geometrical parameters calculated at the B3LYP/6-311++G(2d,2p) and MP2/6-311++G(d,p) levels of theory. The singlet potential energy surface is constructed using the relative energies calculated from the CBS-QB3 level of theory and presented in Figure 2. From Figure 2 it is observed that there are five decomposition channels from the direct decomposition of DMEP (C1, C2, and C3 conformers) involving the products P1 + CH2, P2 + CH3OH, P3 + H2, P5 + CH3OH, and P6 + CH3OH. Moreover, a secondary decomposition channel from P3 is also possible, leading to the dissociation products P4 + CO. We have also computed three different isomerization channels to IM1, IM2, and IM3 from the conformers C1, C2, and C3, respectively. Among the five decomposition channels from the primary decomposition of DMEP, there are three channels associated with methanol CH3OH elimination, one channel associated with hydrogen H2 elimination, and the rest associated with methylene CH2 elimination.
Energetics and Kinetics of Dimethyl Ethylphosphonate
J. Phys. Chem. A, Vol. 114, No. 39, 2010 10719
Figure 1. Optimized geometries with geometrical parameters for the species involved in the unimolecular reaction of DMEP. (B3LYP/6311++G(2d,2p) and MP2/6-311++G(d,p) values are in normal and bold faces, respectively.)
TABLE 1: Relative Energies (kcal/mol) Calculated at Different Theoretical Levels for the Species Involved in the Isomerization and Decomposition Reaction Pathways of DMEP
Figure 2. CBS-QB3 potential energy surface for unimolecular reactions of DMEP.
3.1. Conformers of DMEP. We have located three different conformers of DMEP having very low-energy separation. This observation is very consistent with the existing literature for the important analogues of DMEP, i.e., DMHP and DMMP.
species
B3LYP/ 6-311++G(2d,2p)
MP2/ 6-311++G(d,p)
CBS-QB3
C1 C2 C3 IM1 IM2 IM 3 P1 + CH2 P2 + CH3OH P3 + H2 P4 + H2 + CO P5 + CH3OH P6 + CH3OH TS1 TS2 TS3 TS4 TS5 TS6 TS7 TS8 TS9 + H2
0.00 0.21 2.34 49.12 63.48 44.29 80.89 56.60 81.22 41.44 51.30 56.12 65.71 85.83 71.44 67.29 89.47 63.48 84.97 94.08 95.92
0.00 0.04 1.79 45.18 59.88 44.64 80.72 62.49 84.83 35.06 53.14 58.02 65.89 90.09 70.20 66.79 95.78 63.91 83.25 91.91 99.77
0.00 0.11 1.87 48.72 62.26 45.82 86.48 68.79 83.39 51.55 59.92 64.48 63.95 89.27 71.55 72.99 93.41 65.01 90.19 99.23 105.63
10720
J. Phys. Chem. A, Vol. 114, No. 39, 2010
TABLE 2: Standard Enthalpies at 298.15 K (H°298.15), Atomization Enthalpies (Hat), and Standard Enthalpies of Formation at 298.15 K (∆fH°298.15) Using the CBS-QB3 Level for the Species Involved in the Isomerization and Decomposition Reactions of DMEP species
H°298.15 (hartree)
Hat (kcal/mol)
∆fH°298.15 (kcal/mol)
experimental (kcal/mol)
C1 C2 C3 IM1 IM2 IM3 P1 P2 P3 P4 P5 P6 TS1 TS2 TS3 TS4 TS5 TS6 TS7 TS8 TS9 CH3OH CO H2 CH2
-725.118802 -725.118606 -725.115979 -725.041879 -725.020388 -725.045317 -685.912804 -609.472221 -723.820446 -610.691051 -609.487556 -609.480107 -725.017772 -724.976175 -725.005532 -725.00249 -724.969918 -725.015389 -724.975694 -724.961209 -723.785083 -115.53568 -113.178703 -1.162778 -39.06603
-1721.92 -1721.80 -1720.15 -1673.65 -1660.16 -1675.80 -1453.30 -1164.40 -1531.51 -1304.91 -1174.02 -1169.35 -1658.52 -1632.42 -1650.84 -1648.93 -1628.50 -1657.02 -1632.11 -1623.03 -1509.32 -487.92 -257.53 -105.33 -180.78
-207.85 -207.73 -206.08 -159.58 -146.10 -161.74 -215.00 -90.02 -121.75 -126.23 -99.65 -94.97 -144.46 -118.35 -136.77 -134.87 -114.43 -142.96 -118.05 -108.96 -159.17 -48.23 -26.46 -1.03 94.98
-48.99 ( 2.40a -26.41 ( 0.04b 0.00 92.34c
a
Reference 44. b Reference 45. c Reference 46.
Ewig and Van Wazer42 located three stable conformers of dimethyl hydrogen phosphonate (DMHP). Barnes et al.18 studied the conformers of dimethylmethyl phosphonate (DMMP) in liquid and gas phases and also in argon, nitrogen, and carbon monoxide matrices. They identified the existence of three conformers for DMMP and concluded that depending upon the polarity of the matrix used the conformer distribution varied during the deposition process. Yang et al. have also located three main conformer of DMMP at ab initio and density functional level of theories.19 For conformational analysis, the geometry optimization has been started from different geometries of DMEP involving different spatial orientation of the methoxy group (-OCH3). All the conformers are generated by the internal rotation of P-O and C-C bonds with very low rotational barriers. Conformer 1 (C1) has the lowest energy among the all three conformers and follows the energy order C1 (0.00) < C2 (0.11) < C3 (1.87) calculated at the CBS-QB3 level. The standard enthalpies of formation at 298.15 K (∆fH°298.15) values for C1, C2, and C3 are very close to each other and the differences are within 1 kcal/mol (refer to Table 2). The calculated dipole moments of C1, C2, and C3 are 2.356, 2.334, and 4.7871 D, respectively. The three conformers differ from each other by the orientation of the ethyl group and can be interconvertable through internal rotation of the C-C and O-C bond. The conformer 3 (C3) has an orientation of one methoxy group, which is different from that of C1 and C2, and therefore contains a higher value of dipole moment compared to those for C1 and C2. C3 can be converted to C1 by an internal rotation of O-C bond and to C2 by C-C and O-C rotations. 3.2. Reaction Pathway Map for the Isomerization and Dissociation Reactions of DMEP. 3.2.1. Isomerization Pathways. In this section we have computed different isomerization pathways of C1, C2, and C3 producing isomer 1 (IM1), isomer
Mandal et al. 2 (IM2), and isomer 3 (IM3). The isomerization pathways from C1, C2, and C3 can be classified in two categories, 1,3 H-transfer and, 1,4 H-transfer through cyclization. The isomers IM1, IM2, and IM3 lie above 48.71, 62.14, and 43.95 kcal/mol from their precursors C1, C2, and C3, respectively. These isomers are very different from each other with respect to their electronic structures and energies. The CBS-QB3 calculated enthalpies of formation for IM1, IM2, and IM3 are -159.58, -146.10, and -161.73 kcal/mol, respectively. The isomerization pathways from C1, C2, and C3 are depicted in Figure 2 with the isomerization TSs. All isomerization channels are calculated to be highly endothermic but occur almost exclusively at relatively low-temperature region. Isomerization of C1 to IM1. The formation of IM1 from C1 occurs via 1, 4 H-transfer from one of the -OCH3 group to the O atom of PdO producing a OsH bond. The IM1 contains a three-membered backbone like -PsOsCH2. A very high activation barrier 63.94 kcal/mol is encountered via a five membered transition state TS1. In TS1 the transfer of a H atom maintains a distance 1.798 and 1.053 Å from the corresponding C and O atoms, respectively. The low-pressure isomerization rate for this isomerization is calculated to be twice that of the dissociation pathway from C1 at 1000 K. The low-pressure rate expression for isomerization reaction from C1 is k0 ) 7.31 × 10-2T-2.192exp(-16.62 kcal/mol-1 /RT) cm3 mol-1 s-1 (T ) 1000 - 3000K)
Isomerization of C2 to IM2. IM2 is formed form C2 via a 1, 4 H-transfer reaction from the terminal carbon of the -C2H5 group to the O atom of PdO producing a O-H bond. A three membered ring -P-CH2-CH2, different from the former case, builds the structure of IM2 resulted from the H-transfer reaction. A five-membered transition state TS3 involving activation barrier 71.44 kcal/mol is responsible for this step of conversion. In TS3 bond length of the breaking C-H bond is 1.859 Å and that of the forming O-H bond is 1.039 Å. This isomerization is almost four times faster than the decomposition reactions from C2 at 1000 K and the low-pressure isomerization rate is calculated to be k0 ) 2.28 × 103T-3.393 exp(-22.28 kcal/mol-1 /RT) cm3 mol-1 s-1 (T ) 1000 - 3000 K)
Isomerization of C3 to IM3. A simple 1,3 H-shift from the -CH2 group to the O atom of PdO leads to the formation of IM3, which involves a four-membered transition state TS6. The moving H atom in TS6 is 1.612 and 1.216 Å away from C and O atoms, respectively. The activation barrier for this step of conversion is also very high, 63.14 kcal/mol, which is similar to the former cases. The calculated low-pressure limit rate constant for this isomerization k0 ) 3.18 × 10-2T-2.281 × exp(-16.13 kcal/mol-1 /RT) cm3 mol-1 s-1 (T ) 1000-3000 K)
is 104 times larger than the corresponding decomposition rates of C3 at 1000 K. 3.2.2. Decomposition Pathways. Decomposition from Conformer C1. The conformer C1 is the most stable among the three conformers of DMEP. We have identified only one unimolecular decomposition channel from C1 through the elimination of methylene resulting in the product P1. The
Energetics and Kinetics of Dimethyl Ethylphosphonate TABLE 3: Temperature Dependent Free Energy Changes (∆rG, kcal/mol) for All the Unimolecular Decomposition Reactions from DMEP reaction channel decomposition from C1
decomposition from C2
decomposition from C3
T (K)
channel 1
channel 2
1500 1800 2000 2500 1500 1800 2000 2500 1500 1800 2000 2500
13.21 -0.32 -9.18 -30.59 14.34 4.43 -2.22 -19.03 9.17 0.8 -4.79 -18.52
37.34 29.73 21.54 21.54 17.77 9.70 -0.75 -14.10
optimized geometries and potential energy diagram for this decomposition channel are shown in Figures 1 and 2, respectively. In this decomposition reaction, one H atom is cleaved from -OCH3 and added to the O atom of the same group; besides, a simultaneous O-C bond cleavage occurs through a three-membered transition state TS2. The IRC from TS2 directly gives P2 + CH2 and C1 as two well-defined minima. We could not find any intermediate structure between TS2 and P2 + CH2. This observation is different from that found in the case of the conformer II of DMMP.19 In the H-transfer transition state TS2, the moving H atom is situated 1.552 and 1.027 Å away from C and O atoms, respectively. As expected, the transition vector for the imaginary frequency 1002i cm-1of the TS2 is governed by the movement of the transferring H atom. The dissociation products are 86.48 kcal/mol above C1. This dissociation processes is associated with a very high activation barrier 89.27 kcal/mol (CBS-QB3 value). Temperature dependent free energy change (∆rG) for this process shows that this reaction becomes exergonic at higher temperatures (g1800 K) and therefore favored thermodynamically at higher temperatures. Table 3 summarizes the ∆rG values at different temperatures. Decomposition from Conformer C2. Two dissociation channels are identified from the thermal decomposition of C2, which is the next thermodynamically stable conformer after C1. First unimolecular decomposition from C2 yields P2 via CH3OH elimination and the second one is associated with the formation of P3 via H2 elimination. The first channel occurs through a simultaneous breaking of C-H and P-O bonds and then the formation of an O-H bond resulted in the elimination of CH3OH. It passes through a four-membered nonplanar transition state TS4, where the P-O bond distance is elongated from 1.62 to 1.94 Å and the moving H atom maintains distances 1.645 and 1.122 Å from C and O atoms, respectively. The imaginary frequency 1119i cm-1 is associated with the transition vector dominated by the movement of the transferring H atom. The activation barrier for this decomposition pathway is 72.99 kcal/ mol. It is observed from Table 3 that this reaction becomes exergonic at higher temperatures (g2000 K) and therefore also favored thermodynamically at higher temperatures. Another decomposition channel from C2 is associated with the elimination of H2 through the formation of P3, which involves a simultaneous rupture of two C-H bonds from one -OCH3 group. It occurs via the transition state TS5 with barrier height 93.40 kcal/mol, which is 20.41 kcal/mol larger than the barrier for TS4. The TS5 consists of a three-membered C-H-H moiety where the two eliminating H atoms are 0.815 Å away from each other. The movement of the eliminating H atom,
J. Phys. Chem. A, Vol. 114, No. 39, 2010 10721 which is 1.454 Å away from C atom, is responsible for the transition vector associated with the imaginary frequency 901i cm-1. Interestingly, this decomposition channel, unlike the former one, is endergonic even at higher temperatures (∼2500 K) (refer to Table 3) and therefore least important. Moreover, the product P3 formed in this way is higher in energy by 14.60 kcal/mol than P2. A secondary decomposition channel from P3 is also identified and predicted to occur readily because it has a low activation barrier of 22.25 kcal/mol and high exothermicity of 31.38 kcal/ mol even at 0 K. It involves CsH and PsO bond cleavages and a consecutive H atom transfer to the O atom of PdO. This step occurs via a five-membered planar transition state TS9 in which the PsO distance changes from 1.71 to 2.20 Å as we move from P3 to TS9. Decomposition from Conformer C3. Like the decomposition pathways from C2, we have also identified two different decomposition channels from C3. These two reaction pathways are associated with the elimination of CH3OH. In the first unimolecular decomposition channel from C3, the elimination of CH3OH occurs through three-membered cyclization resulting in the formation of P5. This channel involves the transfer of one H atom from an -OCH3 group and simultaneous P-OCH3 bond cleavage of other -OCH3 with P-O-CH2 ring closure resulting in P5 formation. This step of decomposition occurs via a five-membered nonplanar transition state TS7. The breaking P-O bond increases from 1.61 to 1.99 Å as we move from C3 to TS7. Moreover, the forming P-C distance decreases from 2.61 to 2.09 Å from C3 to TS7, which clearly indicates the formation of a three-membered ring via elimination of CH3OH. The TS7 is associated with activation barrier 90.18 kcal/mol. The movement of the transferring H atom, which is 1.557 and 1.116 Å away from C and O atoms, respectively, is responsible for transition vector associated with the imaginary frequency 1183i cm-1. The calculated ∆rG values, displayed in Table 3, indicate that this reaction is also thermodynamically favorable at higher temperatures (g2000 K). In the second unimolecular decomposition of C3, the elimination of CH3OH occurs through three-membered cyclization resulting in the formation of P6. The product P6 is less stable than P5 by 4.56 kcal/mol. This channel involves the transfer of one H atom from the terminal carbon of the ethyl group to the O atom of the methoxy group and the corresponding cleavage of P-O and simultaneous ring closure P-CH2-CH2 resulting in P6 formation. This step of decomposition occurs via a fivemembered transition state TS8. The breaking P-O bond increases from 1.61 to 2.08 Å from C3 to TS8. The movement of the transferring H atom, which is 1.617 and 1.086 Å away from C and O atoms, respectively, is responsible for transition vector associated with the imaginary frequency 887i cm-1. The activation barrier associated with TS8 for this second decomposition path is 99.23 kcal/mol, which is 9.05 kcal/mol higher than that associated with TS7. The ∆rG values are calculated to be negative at higher temperatures (g2000 K) (refer to Table 3). 3.3. Thermochemistry. The important thermochemical parameter like standard enthalpies of formation at 298.15 K (∆fH°298.15) is calculated using CBS-QB3 energies for all the species involved in the unimolecular reaction of DMEP and are tabulated in Table 2. The high-accuracy energy prediction model CBS-QB3 is tested well with organophosphorus compounds and reported to produce very accurate results. As there is lack of experimental ∆fH°298.15 data for most of the species involved in this study, the reliability of the method is verified
10722
J. Phys. Chem. A, Vol. 114, No. 39, 2010
by comparing the calculated enthalpy of formation values for CH3OH, CO, and CH2 with existing literature43-45 and a good agreement is observed (refer to Table 2). Furthermore, the accuracy of the employed method is also checked with the DMMP molecule, which is an important analogue of DMEP. We have calculated the ∆fH°298.15 value -204.37 kcal/mol for DMMP, and this value agrees very well with the BAC-MP4 value -198.00 kcal/mol of Melius et al.15 From Table 2 we found that all the calculated ∆fH298.15 values are negative and very close to each other. Also the stability order of C1, C2, and C3 is clearly reflected from the ∆fH298.15 values. The calculated standard enthalpies of formation in Table 2 are used for the kinetics calculations. 3.4. Reaction Kinetics. From the above discussion on the decomposition mechanism of different DMEP conformers, one can conclude that the channel with lowest activation barrier is more important regarding product formation than that with higher activation barrier, and therefore, the lowest energy barrier channel is kinetically most important. So we have skipped the detailed discussions on the decomposition channels with higher activation barriers. The rate constants for the decomposition of DMEP are calculated with the microcannonical RRKM theory by solving one-dimensional master equation (using Nesbet method of iteration) coupling both pathways from same conformer. From the calculation of free energy changes for all the reaction channels presented in Table 3, it is found that noncatalytic thermal decompositions of the title compound are thermodynamically not feasible at temperatures up to 1800 K. So, drastic experimental conditions are required for the decomposition. Again, the decomposition of organophosphorus compounds in H2/O2 flames clearly indicates that it is a case of lowpressure and high-temperature decomposition. It is important to mention here that several experimental studies on the decomposition of DMMP in H2/O2 flames can be found in the literature.7,20,43 Due to the lack of experimental data regarding the actual condition for decomposition reactions of DMEP, we choose the temperature and pressure limit from thermodynamic consideration. For the calculation of rate constants we have selected the ranges of temperature and pressure are 1000-3000 K and 10-4-102 atm, respectively. Argon is used as a bath gas with σ ) 3.542 Å and ε ) 39.95 K. The tunneling barrier width is calculated by the nonlinear curve fitting of the IRC curve with Eckart’s potential and the calculated barrier width (in amu0.5 Å) is used to estimate the tunneling coefficient. Rate Constants for the Decomposition of C1. The Arrhenius fitted high-pressure (k∞1 ) and low-pressure (k01) limit rate expressions for C1 decomposition are as follows:
k∞1 ) 4.46 × 109T0.76 exp(-90.47 kcal/mol-1 /RT) s-1
Mandal et al.
Figure 3. Variation of the Arrhenius fitted rate constants with temperature at different pressures for C1 f P1 + CH2.
TABLE 4: Arrhenius Parameters at Different Pressures for Selected Channels reaction C1 f P1 + CH2
C2 f P2 + CH3OH
C3 f P5 + CH3OH
a
pressure (atm) a
Pf0 10-4 10-3 10-2 10-1 100 101 102 Pf∞ P f 0a 10-4 10-3 10-2 10-1 100 101 102 Pf∞ P f 0a 10-4 10-3 10-2 10-1 100 101 102 Pf∞
A (s-1) 5.01 × 10 3.54 × 1099 8.12 × 1083 8.51 × 1062 5.24 × 1040 6.16 × 1022 1.31 × 1013 1.84 × 1010 4.46 × 109 2.34 × 1004 1.38 × 1089 5.75 × 1086 6.02 × 1092 3.32 × 1092 9.99 × 1083 1.77 × 1068 4.67 × 1047 7.41 × 109 7.24 × 1018 3.80 × 1056 1.25 × 1035 5.88 × 1018 1.04 × 1011 1.06 × 109 5.62 × 108 5.24 × 108 6.76 × 108 18
Ea (kcal mol-1)
n
41.40 153.63 145.10 131.02 114.71 100.91 93.26 90.96 90.47 23.56 107.14 110.67 121.44 126.95 125.22 116.07 101.82 72.77 41.20 125.01 109.11 96.45 90.30 88.70 88.48 88.46 88.44
-7.53 -24.63 -20.05 -14.08 -7.84 -2.85 -0.17 0.62 0.76 -3.84 -22.66 -21.58 -22.83 -22.38 -19.77 -15.13 -9.29 1.24 -7.60 -13.02 -7.00 -2.45 -0.30 0.24 0.32 0.33 0.33
A has a unit of cm3 mol-1 s-1.
k01 ) 5.01 × 1018T-7.53 exp(-41.40 kcal/mol-1 /RT) cm3 mol-1 s-1
The apparent rate constants for this decomposition reaction have been calculated within the temperature and pressure limit stated previously. The temperature variation of the Arrhenius fitted apparent rate constants at different pressures are displayed in Figure 3. From Figure 3 we see that the high-pressure limit (HPL) rate constant lines almost coincide with the line at 1 atm and at lower temperatures but are distinguishable at higher temperatures. The fitted kinetic parameters at different pressures are included in Table 4. The values in Table 4 show that the decomposition reaction is highly temperature dependent at low pressures.
Rate Constants for Decomposition of C2. From the earlier discussions of the reaction mechanism we find that C2 decomposes into two product channels; one is P2 + CH3OH (2a) and the other is P3 + H2 (2b) with activation barriers 72.99 and 93.40 kcal/mol, respectively. According to the values of the activation barrier, the product channel P2 + CH3OH is clearly dominating and as a result the branching ratio of the remaining channel is very low for the entire range of temperature and pressure. The Arrhenius fitted high-pressure (k2∞) and lowpressure (k20) limit rate expressions for the two possible decomposition channels from C2 are as follows:
Energetics and Kinetics of Dimethyl Ethylphosphonate
J. Phys. Chem. A, Vol. 114, No. 39, 2010 10723
k∞2a ) 7.41 × 109T1.24 exp(-72.77 kcal/mol-1 /RT) s-1 k∞2b ) 2.81 × 109T1.55 exp(-93.63 kcal/mol-1 /RT) s-1 0 ) 2.34 × 104T-3.84 exp(-23.56 kcal/mol-1 /RT) cm3 mol-1 s-1 k2a
k02b ) 2.80 × 1023T-8.89 exp(-47.26 kcal/mol-1 /RT) cm3 mol-1 s-1
The Arrhenius parameters at different pressures are also calculated and presented in Table 4 for channel 2a. The branching ratios are calculated using rate constants, and their linear variation with temperature for channels 2a and 2b are shown in Figure 4a. Figure 4a indicates that the channel with lower activation barrier is exclusive and the branching ratio negligibly decreases with temperature. To illustrate the effect of pressure on reaction rate, we have plotted in Figure 4b the Arrhenius fitted apparent rate constants with temperature at different pressures for channel 2a. From Figure 4b we can conclude that pressure has a significant effect on reaction rates for channel 2a. The high-pressure line is separated from the 102 atm line at higher temperatures but merges to the 102 atm line at relatively lower temperatures. The Arrhenius fitted kinetic parameters for reaction reaction 2a are tabulated in Table 4. Rate Constants for Decomposition of C3. Like C2, we have also identified two different decomposition channels for C3; one resulted in P5 + CH3OH (3a), and the other resulted in P6 + CH3OH (3b). Channels 3a and 3b are associated with very high activation barriers 90.18 and 99.23 kcal/mol, respectively. The activation barrier for the channel producing P5 + CH3OH is lower than that for the other channel, and the product P5 is also thermodynamically more stable than P6. This makes the P5 + CH3OH producing channel thermodynamically as well as kinetically more favorable, and therefore, only this channel is taken into account for detailed kinetic consideration. The calculated high- and low-pressure limit rate constants for C3 decomposition are as follows:
k∞3a ) 6.76 × 108T0.333 exp(-88.44 kcal/mol-1 /RT) s-1 ∞ k3b ) 3.46 × 108T0.37 exp(-97.04 kcal/mol-1 /RT) s-1
0 k3a ) 7.24 × 1018T-7.60 exp(-41.20 kcal/mol-1 /RT) cm3 mol-1 s-1
k03b ) 4.89 × 1029T-10.46 exp(-54.70 kcal/mol-1 /RT) cm3 mol-1 s-1
Based on the rate constants, the branching ratios are calculated and plotted in Figure 5a as a linear function of temperature for both channels 3a and 3b. Figure 5a shows the exclusive formation of P5 + CH3OH over the applied temperature range. The pressure effect on the rate constant for the exclusive channel 3a is demonstrated in Figure 5b. The Arrhenius fitted apparent rate constants are plotted in Figure 5b as a function of temperature at different pressures. As seen from Figure 5b, the high-pressure rate line is distinguishable from the other pressure lines and also has a constant separation from others within the total temperature range. This is different from the pressure dependency of the decomposition nature for the conformers C1 and C2. The Arrhenius parameters for reaction channel 3a at different pressures are given in Table 4. The pressure-dependent
Figure 4. (a) Branching ratio for channels 2a and 2b. (b) Variation of the Arrhenius fitted rate constants with temperature at different pressures for channel 2a.
rate parameters (A, n, Ea) for the construction of kinetic models have been determined for important reaction pathways in the DMEP decomposition at P ) 0.0001-100 atm. It should be noted that some of the values of pre-exponential factors (A) and power of T (n) reported in Table 4 are quite high. These parameters are fits to the calculated rate constants at a given pressure over a temperature range as simply the best fits, as far as the program can perform, to the data over the temperature range. They have little physical significance; they are only functional to describe the rate constants over the temperature range. Concerning the values of the resulting rate constants, they are meaningful in the context of DMEP decomposition. 4. Conclusions The structures, isomerization, and decomposition mechanism of dimethyl ethylphsophonate DMEP have been studied using high-level B3LYP/6-311++G(2d,2p), MP2/6-311++G(d,p), and CBS-QB3 methods. Thermochemical and kinetic parameters of the isomerization and decomposition of DMEP have been determined using the CBS-QB3 method and RRKM master equation, respectively. The conformers and isomerization and decomposition products are analyzed thermochemically by calculating their standard enthalpies of formation at 298.15 K (∆fH°298.15). The reaction kinetics is modeled over the temperature and pressure ranges 1000-3000 K and 10-4-102 atm, respectively. Our results show that DMEP exhibits three
10724
J. Phys. Chem. A, Vol. 114, No. 39, 2010
Mandal et al. S1 and S2 and the harmonic frequencies for all species are tabulated in Tables S3 and S4. This material is available free of charge via the Internet at http://pubs.acs.org. References and Notes
Figure 5. (a) Branching ratio for channels 3a and 3b. (b) Variation of the Arrhenius fitted rate constants with temperature at different pressures for channel 3a.
conformational isomers that can isomerize and decompose differently. For all the conformers, isomerization is found to be the lower energy processes compared to the dissociations. On the basis of the calculated low- and high-pressure limit rate constants for unimolecular decompositions, it is observed that all the conformers dissociate through only one exclusive pathway. Calculation of temperature dependent free energy changes for unimolecuar decompositions (∆rG) leads to the conclusion that all the reactions become exergonic only at temperatures g1800 K. Therefore, the experimental fact of very high-temperature decomposition of organophosphorus compounds in H2/O2 flames is quantitatively verified in our work. Moreover, the gas-phase isomerization and decomposition mechanism, explored in this work, can provide a clear understanding of the unimolecular reactions of DMEP in the absence of a catalyst. Acknowledgment. D.M. and B.M. are very grateful to the Council of Scientific and Industrial Research (CSIR), Government of India, for providing research fellowships. We are also thankful to the Annesur Rahaman Center for High Performance Computing, IACS, for providing us with computational time. Thanks are due to the reviewer for his constructive comments to improve the quality of the manuscript. We thank Dr. Mokrushin, NIST, U.S.A., for his constant support to execute programs with ChemRate code. Supporting Information Available: Cartesian coordinates for optimized stationary points on PES are tabulated in Tables
(1) Dempsey, C. R.; Oppelt, E. T. Air Waste 1993, 43, 25. (2) U. S. Army’s AlternatiVe Demilitarization Technology; Report for Congress, Department of the Army, Program Manager for Chemical Demilitarization, April 11, 1994. (3) AlternatiVe Technologies for the destruction of chemical Agents and munitions; National Research Council, National Academy Press: Washington, DC, 1993. (4) Demilitarization and Disposal of U.S Chemical Warfare Agent and Munitions; Conference on Disarmament, Ad Hoc Committee on Chemical Weapons, CD/CW/WP.265, Dec 11, 1989. (5) Johnston Atoll Chemical Agent Disposal System. Green CoVer ReView Draft, Final Second Supplemental Environmental Impact Statement for the storage and Ultimate Disposal of the European Chemical Munitions Stockpile, 1990. (6) Korobeinichev, O. P.; Ilyin, S. B.; Boshova, T. A.; Shvartsberg, V. M.; Chernov, A. A. Combust. Flame 2000, 121, 593. (7) Korobeinichev, O. P.; Bolshova, T. A.; Shvartsberg, V. M.; Chernov, A. A. Combust. Flame 2001, 125, 744. (8) Macdonald, M. A.; Gouldin, F. C.; Fisher, E. M. Combust. Flame 2001, 124, 668. (9) Nogueira, M. F.; Fisher, E. M. Combust. Flame 2003, 132, 352. (10) Cuisset, A.; Mouret, G.; Pirali, O.; Ray, P.; Caizer, F.; Nouali, H.; Demaison, J. J. Phys. Chem. B 2008, 112, 12516. (11) Templeton, M. K.; Weinberg, W. H. J. Am. Chem. Soc. 1985, 107, 97. (12) Henderson, M. A.; Jin, T.; White, J. M. J. Phys. Chem. 1986, 90, 4607. (13) Troy, A. D.; Walsh, E. N. Phosphorus compound In EVeryday LiVing; American Chemical Society: Washington, DC, 1987. (14) Hance, C. R., Ed. The Pesticide Manual, 9th ed.; Worthing, British Crop Protection Council: Surrey, U.K., 1991. (15) Melius, C. http://z.ca.sandia.gov/∼melius/, 1997. (16) Sullivan, P. A.; Sumathi, R.; Green, W. H.; Tester, J. W. Phys. Chem. Chem. Phys. 2004, 6, 4296. (17) Sundararajan, K.; Sankaran, K. J. Phys. Chem. A 2008, 112, 5917. (18) Barnes, A. J.; Lomax, S.; Van Der Veken, B. J. J. Mol. Struct. 1993, 99, 137. (19) Yang, L.; Shroll, R. M.; Zhang, J.; Lourderaj, U.; Hase, W. L. J. Phys. Chem. A 2009, 113, 13762. (20) Werner, J. H.; Cool, T. A. Combust. Flame 1999, 117, 78. (21) Henderson, M. A.; White, J. M. J. Am. Chem. Soc. 1988, 110, 6939. (22) Moss, J. A.; Szczepankiewicz, S. H.; Park, E.; Hoffmann, M. R. J. Phys. Chem. B 2005, 109, 19779. (23) Templeton, M. K.; Weinberg, W. H. J. Am. Chem. Soc. 1985, 107, 774. (24) Aschmann, S. M.; Tuazon, E. C.; Atkinson, R. J. Phys. Chem. A 2005, 109, 11828. (25) Bell, A. J.; Ferrante, F.; Hall, S. E.; Mikhalilov, V.; Mitchell, D.; Timperley, C. M.; Watts, P.; Williams, N. Int. J. Mass. Spectrom. 2008, 269, 46. (26) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Montgomery, J. A., Jr.; Vreven, T.; Kudin, K. N.; Burant, J. C.; Millam, J. M.; Iyengar, S. S.; Tomasi, J.; Barone, V.; Mennucci, B.; Cossi, M.; Scalmani, G.; Rega, N.; Petersson, G. A.; Nakatsuji, H.; Hada, M.; Ehara, M.; Toyota, K.; Fukuda, R.; Hasegawa, J.; Ishida, M.; Nakajima, T.; Honda, Y.; Kitao, O.; Nakai, H.; Klene, M.; Li, X.; Knox, J. E.; Hratchian, H. P.; 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.; Ayala, P. Y.; Morokuma, K.; Voth, G. A.; Salvador, P.; Dannenberg, J. J.; Zakrzewski, V. G.; Dapprich, S.; Daniels, A. D.; Strain, M. C.; Farkas, O.; Malick, D. K.; Rabuck, A. D.; Raghavachari, K.; Foresman, J. B.; Ortiz, J. V.; Cui, Q.; Baboul, A. G.; Clifford, S.; Cioslowski, J.; Stefanov, B. B.; Liu, G.; Liashenko, A.; Piskorz, P.; Komaromi, I.; Martin, R. L.; Fox, D. J.; Keith, T.; Al-Laham, M. A.; Peng, C. Y.; Nanayakkara, A.; Challacombe, M.; Gill, P. M. W.; Johnson, B.; Chen, W.; Wong, M. W.; Gonzalez, C.; Pople, J. A. Gaussian 03, Revision E.01; Gaussian, Inc.: Pittsburgh, PA, 2004. (27) Becke, A. D. J. Chem. Phys. 1993, 98, 5648. (28) Lee, C. T.; Yang, W. T.; Parr, R. G. Phys. ReV. B 1988, 37, 785. (29) Miehlich, B.; Savin, A.; Stoll, H.; Preuss, H. Chem. Phys. Lett. 1989, 157, 200. (30) Møller, C.; Plesset, M. S. Phys. ReV. 1934, 46, 618. (31) Head-Gordon, M.; Pople, J. A.; Frisch, M. J. Chem. Phys. Lett. 1988, 153, 503. (32) Frisch, M. J.; Head-Gordon, M.; Pople, J. A. Chem. Phys. Lett. 1990, 166, 275.
Energetics and Kinetics of Dimethyl Ethylphosphonate (33) Frisch, M. J.; Head-Gordon, M.; Pople, J. A. Chem. Phys. Lett. 1990, 166, 281. (34) (a) Gonzales, C.; Schlegel, H. B. J. Chem. Phys. 1989, 90, 2154. (b) Gonzales, C.; Schlegel, H. B. J. Phys. Chem. 1990, 94, 5523. (35) Montgomery, J. A.; Frisch, M. J., Jr.; Ochterski, J. W.; Petersson, G. A. J. Chem. Phys. 1999, 110, 2822. (36) Nicolaides, A.; Rauk, A.; Glukhovtsev, M. N.; Radom, L. J. Phys. Chem. 1996, 100, 17460. (37) Lias, S. G.; Bartmess, J. E.; Liebman, J. F.; Holmes, J. L.; Levin, R. D.; Mallard, W. G. J. Phys. Chem. Ref. Data 1988, 17 (1), . (38) Ignatov, S. K. Moltran v.2.5 - Program for molecular visualization and thermodynamic calculations; University of Nizhny Novgorod, 2004; http://ichem.unn.ru/Moltran. (39) (a) Klippenstein, S. J. J. Chem. Phys. 1992, 96, 367. (b) Klippenstein, S. J.; Marcus, R. A. J. Chem. Phys. 1987, 87, 3410. (c) Wardlaw,
J. Phys. Chem. A, Vol. 114, No. 39, 2010 10725 D. M.; Marcus, R. A. Chem. Phys. Lett. 1984, 110, 230. (d) Wardlaw, D. M.; Marcus, R. A. J. Chem. Phys. 1985, 83, 3462. (40) Eckart, C. Phys. ReV. 1930, 35, 1303. (41) Mokrushin, V.; Bedanov, V.; Tsang, W.; Zachariah, M.; Knyazev, V. ChemRate, Version 1.5.8; National Institute of Standards and Testing: Gaithersburg, MD, 2009. (42) Van Wazer, J. R.; Ewig, C. S. J. Am. Chem. Soc. 1986, 108, 4354. (43) NIST Chemistry WebBook, http://webbook.nist.gov/. (44) Cox, J. D.; Wagman, D. D.; Medvedev, V. A. CODATA Key Values for Thermodynamics; Hemisphere Publishing Corp.: New York, 1984, 1. (45) Chase, M. W., Jr. NIST-JANAF Themochemical Tables, 4th ed. J. Phys. Chem. Ref. Data, Monograph 1998, 9, 1–1951. (46) Korobeinichev, O. P.; Shvartsberg, V. M.; Il’in, S. B. Combust., Explos., Shock WaVes 1997, 33, 270.
JP106270D