ARTICLE pubs.acs.org/JPCA
Theoretical Methodology for Prediction of Tropospheric Oxidation of Dimethyl Phosphonate and Dimethyl Methylphosphonate Marshall G. Cory,† DeCarlos E. Taylor,‡ Steven W. Bunte,‡ Keith Runge,§ Joseph L. Vasey,† and Douglas S. Burns*,† †
ENSCO, Inc., 4849 North Wickham Road, Melbourne, Florida 32940, United States United States Army Research Laboratory Weapons and Materials Research Directorate, RDRL-WMB-D, Aberdeen Proving Ground, Maryland 21005, United States § BWD Associates, LLC, 2901 NW 54th Avenue, Gainesville, Florida 32653-1819, United States ‡
ABSTRACT: Rate constants for the reactions of OH radicals with dimethyl phosphonate [DMHP, (CH3O)2P(O)H] and dimethyl methylphosphonate [DMMP, (CH3O)2P(O)CH3] have been calculated by ab initio structural methods and semiclassical dynamics modeling and compared with experimental measurements over the temperature range 250-350 K. The structure and energetics of reactants and transition structures are determined for all hydrogen atom abstraction pathways that initiate the atmospheric oxidation mechanism. Structures are obtained at the CCSD/6-31þþG** level of chemical theory, and the height of the activation barrier is determined by a variant of the G2MP2 method. A Transfer Hamiltonian is used to compute the minimum energy path in the neighborhood of the transition state (TS). This calculation provides information about the curvature of the potential energy surface in the neighborhood of the TS, as well as the internal forces that are needed by the semiclassical flux-flux autocorrelation function (SCFFAF) dynamics model used to compute the temperature-dependent reaction rate constants for the various possible abstraction pathways. The computed temperature-dependent rate curves frequently lie within the experimental error bars.
’ INTRODUCTION A primary pathway for the oxidation of volatile organic compounds (VOCs) that are emitted into the atmosphere is through hydrogen atom abstraction by OH radical.1-3 Organophosphorus (OP) compounds such as alkyl phosphates [(RO)3P(O)] and alkyl phosphonates [(RO)2P(O)R] where R = aryl or alkyl are widely used as pesticides, flame retardants, and plasticizers.4,5 These OP compounds and their precursors may be released into the atmosphere where they can undergo transport, dispersion, and chemical transformation. The ability to accurately and confidently compute the kinetics of potentially significant atmospheric reactions of toxic industrial compounds (TICs) is a critical component of environmental hazard assessment models. However, many of the critical kinetic data (e.g., rate constants for target reactions) do not exist and can be difficult and expensive to obtain experimentally. The approach described here is to use computational quantum chemistry (QC)/ quantum dynamics models to compute reaction rate constants, which can then be directly input into atmospheric chemistry modules within transport and dispersion (T&D) models. The addition of atmospheric chemistry via these reaction rate constants to environmental T&D models has resulted in the improved ability of the models to handle the fate and interaction of pollutants in the atmosphere, resulting in a more realistic representation of the contamination zone. The theoretical methods applied can predict the nature of the resulting byproducts. This information is valuable in programming advanced sensors to r 2011 American Chemical Society
warn of the release of a toxic compound even if the material has degraded by the time it reached the detection array. In addition, the results of these calculations can also be used to make a more scientifically defensible selection of simulants for challenging detectors, as well as assisting in the evaluation of both individual and collective protection systems. These calculations help provide a sound theoretical underpinning to the development of physical properties data, which is critical for choosing the proper simulant for use in synthesis, fate, and decontamination studies, as well as in the assessment of new potential threat agents and the selection of improved decontamination concepts. The calculations corroborate experimental data, provide an explanation for how products observed by experimentation are formed, and provide insight on the functional form of rate constant curves that are difficult to achieve in the lab.
’ THEORETICAL METHODS Electronic Structure. Ab initio electronic structure calculations were performed with the GAMESS,6 ACES II/III,7 and Gaussian038 suites of programs. The structure and energetics of reactants and transition structures are determined for all hydrogen atom abstraction pathways that initiate the atmospheric oxidation mechanism. Second-order Møller-Plesset perturbation theory (MP2)9 and Received: August 17, 2010 Revised: December 13, 2010 Published: February 22, 2011 1946
dx.doi.org/10.1021/jp107804m | J. Phys. Chem. A 2011, 115, 1946–1954
The Journal of Physical Chemistry A
ARTICLE
Pople’s 6-31þþG** basis set10,11 were used to determine which conformations of the parent DMHP and DMMP structures are energetically accessible at ambient temperature. The relevant lowest energy parent conformations were used as the starting points, and likely hydrogen abstraction transition-state (TS) structures are fully optimized by use of UMP212/6-311þþG** to convergence criteria better than 0.001 Å for bond lengths and 0.1 for angles. For each of the obtained TS structures, the activation energy (i.e., the energy of the TS minus the energy of the reactants) was determined. For those transition states that have an activation energy of e3 kcal/mol, frozen-core UMP2/ 6-311þþG** vibrational frequency calculations are performed. This vibrational frequency calculation serves a 2-fold purpose: (1) it verifies that the computed structure is indeed a first-order transition state (i.e., the signature of the Hessian contains one and only one negative eigenvalue) and (2) it generates a guess Hessian to be used in the subsequent transition-state search by single and double excitation coupled-cluster theory (CCSD)13,14 at the CCSD/6-31þþG** level of chemical theory. Energy Barriers. To achieve a high level of accuracy, we used an energy extrapolation that accounts for the well-understood and predictable deficiencies in the QC methods used. The reported energy barriers were computed by a variant of the G2MP2 protocol.15 This variant, hereafter labeled CCG2MP2, differs from that published in ref 15 in that all structures used in the extrapolation are relaxed at the CCSD/6-31þþG** level with the core orbitals frozen, and the zero-point energy of the CCSD optimized structures is obtained at the frozen-core MP2 level of theory, again by use of the 6-31þþG** basis set. The CCG2MP2 energy extrapolation on the CCSD converged structures consists of the following calculations: (a) (b) (c) (d)
single-point energy: QCISD(T)/6-311G** (QCI) single-point energy: MP2/6-311G** (SML) single-point energy: MP2/6-311þG(3df,2p) (LRG) zero-point energy from a vibrational frequency analysis: MP2/6-31þþG** (ZPE)
The extrapolated energy is then computed as15 EðCCG2MP2Þ ¼ EðQCIÞ þ EðLRGÞ - EðSMLÞ þ ZPE ð1Þ The activation barrier is computed as the difference in energy between the TS and parent structures as shown in eq 2: EA ¼ ETS - Ereactant
ð2Þ
Potential Energy Surface Determination. In addition to the activation barrier, other information characterizing the transition state is also needed. Specifically, the negative (repulsive) CCSD curvature of the minimum energy path at the transition state is also required. However, since a full CCSD vibrational frequency calculation cannot be completed in a timely manner, we compute the CCSD curvature of the repulsive mode at the transition state by fitting a harmonic potential to three CCSD single-point energies (one at the transition state and one on each “side” of the transition state along the reaction coordinate). In addition to the curvature of the repulsive mode at the transition state, the curvatures of the orthogonal modes, in the neighborhood of the transition state, are also required. Again, the molecular size (i.e., the total number of valence electrons and atomic basis functions) of the compounds under study precludes
the use of coupled-cluster theory to calculate this information; instead, we use a neglect of diatomic differential overlap (NDDO) Hamiltonian that is reparametrized to accurately reproduce coupled-cluster theory in order to compute the minimum-energy reaction path.16,17 In particular, we choose to use the geometry of the transition state, including particularly the bond lengths of the atoms directly involved in the hydrogen abstraction, for evaluating the success of the reparametrization. The strength of the repulsive mode calculated from the NDDO Hamiltonian, compared to that obtained from CCSD fitting, is a second major contributor to the determination of the reparametrization’s quality. In order to fit the coupled-cluster reference data, we use the AM1 variant of the NDDO Hamiltonian as our functional form. The one-electron contribution to the Fock matrix by use of the AM1 Hamiltonian for a charge distribution (μv) is given by X AA AA ð3Þ ¼ Uμv δμv ZB μA vA jsB sB Hμv B AB Hμv ¼
1 βμ þ βv Sμv 2
ð4Þ
where the superscripts A and B denote the atomic core on which the basis function is centered and Z is the atomic charge. In AM1, U (atomic orbital energy), β (“resonance” parameter), Slater orbital exponents (used to compute overlap integrals Sμv), and all nonzero one-center two-electron integrals are treated as parameters. The two-center two-electron integrals, which contribute to the oneelectron matrix elements above as well as the two-electron part of the Fock matrix, are computed in terms of the one-center twoelectron integrals on centers A and B by use of the point charge model of Dewar and Thiel.18 Finally, within AM1, the core repulsion between two atomic centers A and B is given by 0 0 EðA, BÞ¼ZA ZB ðsA sA jsB sB Þ 1 þ e-RA RAB þ e-RB RAB X 2 2 Ri, A e-bi, A ðci, A - RAB Þ þ Ri, B e-bi, B ðci, B - RAB Þ ð5Þ þ i
where R, a, b, and c are also taken as variable parameters. Using a simulated annealing algorithm, we vary the β, Slater orbital exponents, one-center two-electron repulsion integrals, R core exponents, and the a, b, and c parameters of the leading Gaussian core repulsion function for each atomic center such that the AM1 transition-state geometry and vibrational frequencies reproduce, to the extent possible within the given parameter space, the CCSD reference values. In this parametrization approach, the parameter sets for each atom in the molecule are individually determined. To clarify, each individual atom is assigned a set of AM1 parameters that are modified separately and allowed to converge to their own respective values. This is similar in spirit to the use of different basis sets on equivalent atom types within the same molecule, as is often done in ab initio modeling, in order to get a better description of the local geometric topology. It has been our experience that it is very difficult to fit equivalent bond types that have significantly different bond lengths (such as those that occur at transition states) with a single set of AM1 parameters for each atom type. The additional flexibility provided by using individual parameter sets for each atom has enabled us to fit the distorted geometries and varying bond lengths that are typically found in transition-state configurations for a wide class of molecules containing as many as 24 atoms. 1947
dx.doi.org/10.1021/jp107804m |J. Phys. Chem. A 2011, 115, 1946–1954
The Journal of Physical Chemistry A
ARTICLE
Table 1. Values of Bond Lengths and Imaginary Mode That Define the Reaction Coordinate for Hydrogen Abstraction from DMHP TSP1,2 and DMMP TS1 DMHP TSP1,2 coordinate
CCSD
fitted AM1
DMMP TS1 CCSD
fitted AM1
C-H bond (Å)
1.243
1.220
1.244
H-(OH) bond (Å)
1.274
1.287
1.264
1.272
-1915
-1965
-2008
-2025
repulsive mode (cm-1)
1.244
The AM1 parameters are varied until the objective function shown in eq 6 reaches a minimum where q corresponds to atomic coordinates and vibrational frequencies at the transition state. X ð6Þ wi qi, NDDO - qi, CCSD i
Each point is assigned a weight wi, and during the fit the points qi corresponding to coordinates defining the reaction path and repulsive eigenmode are weighted more heavily than the remaining internal coordinates and frequencies. Table 1 gives values of the internal coordinates defining the reaction path and the magnitude of the negative (repulsive) mode for the abstraction of hydrogen from the dimethyl methylphosphonate (DMMP) and dimethyl phosphonate (DMHP) TS1 systems. As shown in Table 1, the agreement between the ab initio and AM1 data is excellent; however, the AM1 calculation completes in less than 30 s. By use of the converged AM1 parameters, the minimumenergy path on each side of the saddle point is computed with a second-order Runge-Kutta integration algorithm. Dynamics Modeling. There are a number of available methods for the calculation of reaction rate constants. We choose to use a semiclassical version of the flux-flux autocorrelation function (FFAF) method as presented by Miller and co-workers.19 The approach is semiclassical in two senses. First, the activation barrier that is supplied to the Boltzmannized flux (vide infra) includes the effects of zero-point motion. Second, the value of the flux operator, a quantum mechanical operator, is generated with classical trajectories. As the FFAF method is, in principle, exact, we have chosen not to include possible correction factors, sometimes called “tunneling corrections”, in our calculations. Details of the dynamics modeling used, the semiclassical fluxflux autocorrelation function (SCFFAF) approach,20 have been published elsewhere and we include only a short summary of the approach’s salient features here. The thermal rate constant is related to the flux-flux autocorrelation function by eq 7: Z ¥ -1 Cff ðtÞ dt ð7Þ kðTÞ ¼ Qr ðTÞ 0
where
h i ^ ^ ^ ^ Cff ðtÞ ¼ tr e-βH=2 F^e-βH=2 eiHt=p F^e-iHt=p
ð8Þ
is the Boltzmannized flux-flux autocorrelation function with the ^ inverse temperature β = 1/(kBT), the Hamiltonian operator H, and the flux operator F^. This fully quantum mechanical expression can have a semiclassical realization that begins by casting the flux operator in the form of eq 9: aN pi p 1X B ^ ri sδðsÞ þ δðsÞ ri s Bi F ¼ 2 i ¼ 1 mi mi
ð9Þ
where s is the reaction coordinate whose zero is taken to be at the position of the transition state, and N is the number of atoms in the system. This expression is written in Cartesian coordinates so that the first N momentum entries, B p i, can be thought of as the x components for each atom, the second N as the y components, and the remainder as the z components. In any event, mi is the mass of the atom that corresponds to the ith momentum component. The semiclassical implementation of the equations for reaction rate constants, flux-flux autocorrelation function, and flux operator requires the generation of classical trajectories for the atomic degrees of freedom governed by the electronic potential energy surface. The efficient calculation of thousands to tens of thousands of trajectories is necessary for systems comprised of 15-30 atoms. Rather than attempting direct dynamics, we choose to generate the forces at equal spacing along the reaction coordinate, using the NDDO Hamiltonian. While the trajectories are not constrained to be along the reaction coordinate, the force at any point is obtained by Taylor expansion from points along the reaction coordinate. This Taylor expansion for the potential is expressed in terms of the Cartesian coordinates for each point along the reaction path, B q 0. For the current point along a given dynamic trajectory, B q , the potential is written:
X ∂V
qj - qj0 V ð qBÞ ¼ V ð qB0 Þ þ
∂qj
j q¼ B q0 B
1 X ∂2 V
þ qj - qj0 qk - qk0 ð10Þ
2 j, k ∂qj qk
q¼ B q0 B
where j and k are the points whose reaction coordinate are closest to the current value of the reaction coordinate, one with a slightly larger value and the other a slightly lower value of the reaction coordinate. The force is the negative gradient of the above expression. For the purpose of these semiclassical calculations the δ function is interpreted to be δðsÞ ¼
Θðs þ hÞ - Θðs - hÞ 2h
ð11Þ
Here h, the width of the dividing surface, is chosen such that the flux-flux autocorrelation function is stationary at early time and approaches zero at late times. In the current implementation of the SCFFAF method, the reaction coordinate is determined by diagonalizing the NDDO Hessian at the transition state and choosing the normal mode coordinates that correspond to the repulsive mode (imaginary frequency). Initial conditions for trajectories are chosen from a Gaussian distribution of momenta whose kinetic energy is scaled to sample energies up to the energy available at a given temperature range (e.g., near room temperature the kinetic energy is about 0.6 kcal/ mol or less). The total kinetic energy of any given trajectory is a randomly chosen percentage of the available thermal energy, so that kinetic energies between 0.0 and 0.6 kcal/mol are sampled for room-temperature calculations. The total energy of a given trajectory, the value to represent the Hamiltonian in the calculation of the flux-flux autocorrelation function, is the sum of the kinetic energy and the activation energy extrapolated from the CCSD geometry. The final reaction rate constants are calculated from between 1000 and 10 000 trajectories and are insensitive to the total number of trajectories. 1948
dx.doi.org/10.1021/jp107804m |J. Phys. Chem. A 2011, 115, 1946–1954
The Journal of Physical Chemistry A
ARTICLE
Figure 1. CCSD/6-31þþG** optimized structures for three DMHP parent structures.
Table 2. Comparison of DMHP and DMMP Parent Structures parent
ID
ECCSD/6-31þþG** (au)
ZPE (au)
E (ZPE-corrected) (au)
relative energy (kcal/mol)
Boltzmann weight (%)
DMHP
P1
-646.082 08
0.10334
-646.078 75
0.00
89.6
DMHP DMHP
P2 P3
-646.179 51 -646.178 73
0.10318 0.10310
-646.076 33 -646.075 93
1.52 1.95
7.0 3.4
DMMP
C1a
-685.394 76
0.13186
-685.262 90
0.00
100
Figure 2. MP2/6-31þþG** optimized structures for DMMP with C1 symmetry.
’ REVIEW OF EXPERIMENTAL METHODS Experimental data used to evaluate the theoretical methods and results discussed here have been previously reported.22,25,26 In general, temperature-dependent experiments are carried out in Teflon chambers at atmospheric pressure that are outfitted with a heating/cooling system to maintain the temperature to (2 K within the range of 280-350 K, and a Teflon-coated fan to mix the reactants during their introduction into the chamber. Rate constants for the reactions of OH radicals with DMHP and DMMP were measured by relative rate techniques in which the concentrations of the organophosphorus compound and a reference compound (whose OH radical reaction rate constant is reliably known) were measured in the presence of OH radicals. Reaction 1 :
OH þ target compound f degradation products
ð12Þ Reaction 2 :
OH þ reference compound f degradation products
ð13Þ If it is assumed that the target and reference compounds react only with OH radicals, then ! ! ½targett0 ½referencet0 k1 ln - Dt ¼ ln - Dt ð14Þ ½targett k2 ½referencet where [target]t0 and [reference]t0 are the concentrations of the target and reference compounds, respectively, at time t0; [target]t and [reference]t are the corresponding concentrations at time t; Dt is a factor to account for dilution caused by any
Figure 3. CCSD/6-31þþG** optimized structures for five DMHP TS’s.
additions to the chamber during the experiments (Dt = 0 for the OH radical reactions); and k1 and k2 are the rate constants for reactions 1 and 2, respectively. Specific details of the experimental methodology and results can be found in Atkinson and co-workers.22
’ RESULTS DMHP and DMMP Parent Structures. Three conformers of the DMHP parent were originally optimized at MP2/631þþG** and then refined by use of CCSD/6-31þþG** to identify the most stable structure(s). As shown in Figure 1, DMHP exists with different relative orientations of the O-CH3 groups. The ZPE-corrected CCSD/6-31þþG** energies of the three conformers lie close in energy (see Table 2) with P1 lying 1.52 and 1.95 kcal/mol lower in energy than P2 and P3, resulting in a Boltzmann-weighted relative number concentration of 90/ 7/3 at 300 K. On the basis of the relative energetics of the DMHP conformers, the P1 conformer was optimized at the CCSD/ 6-31þþG** level of chemical theory. Two conformers of the DMMP parent were optimized by use of MP2/6-31þþG** to identify the more stable structure (Figure 2). ZPE-corrected energies show the two C1 conformers are close in energy (EMP2,C1a = -685.386 87 au and EMP2,C1b = -685.384 08 au) with C1a lying 1.75 kcal/mol lower in energy, resulting in a Boltzmann-weighted relative number concentration of 95/5 at 300 K. The C1a conformer was then optimized by use of CCSD/6-31þþG**. 1949
dx.doi.org/10.1021/jp107804m |J. Phys. Chem. A 2011, 115, 1946–1954
The Journal of Physical Chemistry A
ARTICLE
Table 3. Relevant Characteristics for Important DMHP TS Structures TS
multiplicitya
EA (kcal/mol)
νrepulsive (cm-1)
TSP1,2
3
-0.95
-1915
TSP2,2
3
þ1.07
-2041
TSP1,1
1
þ0.13
TSP2,1
1
TSP3,1
1
rH-O(H) (Å) in TS
rX-H (Å) in TS
rX-H (Å) in parent
1.273
1.234
1.090
1.253
1.245
1.089
-1165
1.397
1.479
1.388
þ0.55
-1119
1.409
1.481
1.392
-0.060
-677
1.456
1.470
1.398
Abstraction from O-CH3
Abstraction from P-H
a
Multiplicity = number of rotationally equivalent hydrogen atoms at 300 K.
Rate Constants for the Reactions of OH Radicals with DMHP. Potential transition-state structures for abstraction of
hydrogen were identified for all seven of the hydrogen atoms contained in the three DMHP parent conformers. In each case, three orientations of the OH radical about each hydrogen atom (i.e., 3 7 = 21 geometry optimizations) were evaluated at MP2/6-31þþG** to identify the lowest energy TS structures. In many cases the resulting transition-state structures morphed into common configurations. Transition States. The five lowest energy TS’s were optimized at the CCSD/6-31þþG** level of theory by use of the ACES III QC package and are shown in Figure 3. These five TS’s include the lowest energy TS for each type of abstraction and all those within about 3 kcal/mol of the overall lowest energy TS. Two of the TS’s involve abstraction from the O-CH3 group (TSP1,2 and TSP2,2) and three involve abstraction of H from phosphorus (TSP1,1, TSP2,1, and TS3,1). Interestingly, at the TS the hydroxyl radical has an interaction with other oxygen atoms in the molecule (e.g., the PdO oxygen). Energetics. The computed barrier (EA) and repulsive mode (νrepulsive) for the five lowest energy TS’s are summarized in Table 3. TSP1,2 is the only pathway associated with abstraction from O-CH3 that leads to a negative barrier (indicating the presence of a van der Waals complex between the DMHP and the hydroxyl radical)21 and as such is expected to be a significant contributor to the overall kOH. The TSP3,1 associated with abstraction of the H from the P-H group is very slightly negative (-0.060 kcal/mol), but it is also associated with the least repulsive mode that was studied (-677 cm-1). Interestingly, the change in the C-H bond lengths associated with abstraction from an O-CH3 are greater (rC-H = 1.09 f 1.23-1.24 Å) than in the case of abstraction from P-H (where rP-H = 1.39 f 1.48 Å). Kinetics. Table 4 presents a summary of the temperaturedependent kOH kinetics (and error bars) measured by several groups including Atkinson and co-workers,22 Kleindeinst and Smith,24 and Andino.23 The experimental error bars for the temperature are (2 K. In addition, the theoretical rate curve for temperatures from 250 to 350 K for each of the four independent pathways is included, as is the total theoretical abstraction kOH. The data in Table 4 are plotted in Figure 4 and show the computed temperature-dependent rate curve for the DMHP system relative to experimental measurement. For DMHP, there are five transition states that contribute to the overall reaction rate constant. The total rate, shown in Figure 4 as a solid blue line, is simply the sum of the rate constants resulting from each individual transition state weighted by the number of equivalent hydrogen atoms (curves TS1 through TS3 in Figure 4) and agrees quite favorably with the previously reported experimental
data.6-10 Note that each of the experimental data points in Figure 4 is an individual measurement at a single temperature, whereas the theoretical calculations provide the full temperaturedependent rate constant. Calculations reveal that most of the decomposition is attributed to the abstraction from O-CH3 via TSP1,2 (∼88% of the total kOH,300). Consistent with the calculated barriers shown in Table 3, kOH,TSP2,2 (O-CH3) increases with increasing temperature. However, the overall negative temperature dependence is dominated by TSP1,2. Smith and Ravishankara21, Atkinson and co-workers22 propose that such reactions, as they have negative activation barriers, proceed through the formation of a van der Waals complex as in eq 15: OH þ ðCH3 OÞ2 PðOÞH T HO - ðCH3 OÞ2 PðOÞH vdw f H2 O þ CH3 OPðOÞðHÞOC• H2
ð15Þ
The decomposition then occurs from the van der Waals complex via breaking of the C-H bond to form H2O and the alkyl radical. Evidence of the effect of the van der Waals complex and electronic induction therefrom can be seen in the transitionstate geometries where the stationary point is encountered earlier as a function of the C-H (or P-H) bond distance. The quality of the theoretical model is demonstrated in the following: (1) The quantitative values are in good agreement with and, in many cases, within the experimental error bars from multiple laboratories. (2) The negative temperature dependence of the overall kinetics has been properly captured. Atkinson and coworkers22 reported an Arrhenius expression for DMHP with hydroxyl radical of kOH(T) = 1.01 10-12 e(474(159)/T cm3 3 molecule-1 3 s-1 for 278 < T < 351 K. Here the activation barrier fitted to the experimental data is about -0.94 kcal/mol, which is essentially the same as the activation barrier of the transition state (TSP1,2) that contributes most strongly to the overall theoretical rate constant, which is -0.955 kcal/mol. Clearly, including the contribution of the system dynamics is essential to appropriately describing this reaction in this temperature range. (3) The partitioning of the degradation products is in excellent agreement. Atkinson and co-workers22,26 estimated that ∼15% of the overall reaction occurred at the P-H bond. These calculations suggest (at 300 K) ∼8% of the overall decomposition occurs at the P-H bond. Rate Constants for the Reactions of OH Radicals with DMMP. Potential transition-state structures for abstraction of 1950
dx.doi.org/10.1021/jp107804m |J. Phys. Chem. A 2011, 115, 1946–1954
The Journal of Physical Chemistry A
ARTICLE
Table 4. kOH as a Function of Temperature from Experimental and Theoretical Calculations for DMHP theoretical kOH 1012 (cm3 3 molecule-1 3 s-1) experimental kOH 1012 (cm3 3 molecule-1 3 s-1)
P-H
TSP1,2
TSP2,2
TSP1,1
TSP2,1
TS3,1
total
250
7.90
0.20
0.28
0.081
0.17
8.64
260
7.23
0.21
0.28
0.082
0.16
7.96
270
6.63
0.22
0.27
0.083
0.15
7.35
280
6.09
0.23
0.26
0.083
0.15
6.81
290
5.62
0.23
0.25
0.083
0.14
6.32
temp (K)
278
ref 22
296
4.83 ( 0.25a
296 298
5.33 ( 0.43 5.01 ( 0.17b
300
ref 23
ref 24
5.44 ( 0.39
4.23 ( 0.16
4.94 ( 0.19
4.88 ( 0.49
298 4.76 ( 0.25
5.19
0.24
0.24
0.083
0.13
5.88
310
4.80
0.24
0.23
0.082
0.13
5.50
320
4.44
0.25
0.22
0.081
0.12
5.12
330
4.13
0.25
0.21
0.080
0.12
4.78
340 350
3.83 3.57
0.25 0.26
0.20 0.19
0.079 0.078
0.11 0.11
4.48 4.20
322
351 a
O-CH3
4.53 ( 0.17
3.84 ( 0.25
See ref 25. b See ref 26.
Figure 4. Experimental and theoretical rate constants for the reactions of OH radicals with DMHP.
hydrogen were identified for all nine of the hydrogen atoms contained in the lowest energy C1 DMMP parent. Transition States. The five lowest energy symmetry-unique TS’s were optimized at the CCSD/6-31þþG** level of theory and are shown in Figure 5. As with DMHP, these five TS’s include the lowest energy TS for each type of abstraction and all those within about 3 kcal/mol of the overall lowest energy TS. Four of the TS’s involve abstraction from the O-CH3 group (TS1-TS4) and one involves abstraction of H from the P-CH3 group (TS5). The hydroxyl radical in all of the obtained TS’s except TS3 have an interaction with other oxygen atoms in the DMMP (e.g., the PdO oxygen). Energetics. The computed barrier (EA) and repulsive mode (νrepulsive) values for the five lowest energy TS’s are summarized in Table 5.
Kinetics. Table 6 presents a summary of the temperaturedependent kOH kinetics (and error bars) measured by several groups including Atkinson and co-workers,22 Henley et al.,27 Kleindeinst and Smith,24 and Andino.23 The experimental error bars for the temperature are (2 K. In addition, the theoretical rate curve from 250 to 350 K for each of the five independent pathways is included, as is the total theoretical abstraction kOH. The data in Table 6 are plotted in Figure 6 and show the computed temperature-dependent rate curve for the DMMP system relative to experimental measurement. For DMMP, there are five transition states that contribute to the bulk of the overall reaction rate. The total rate shown in Figure 6 as a solid blue line is simply the sum of the rates resulting from each individual transition state weighted for the number of equivalent hydrogen atoms (curves TS1 through TS5 in Figure 6) and agrees quite favorably with the previously reported experimental data.22-24,27 Note that each of the experimental data points in Figure 6 is an individual measurement at a single temperature, whereas the calculated data provide the full temperature-dependent rate curve. Again, the computed rate curves show excellent agreement with the experimental results; however, the experimental data show a marked increase in the low-temperature region as compared to the theoretical prediction. Atkinson and coworkers22 attributed the higher apparent kOH data at low temperature to aerosol formation at lower temperatures. Calculations reveal that most of the decomposition is attributed to the abstraction from O-CH3 via TS2 (∼80% of the total kOH,300). Consistent with the calculated barriers shown in Table 5, kOH,TS1 and kOH,TS2 (O-CH3) decrease with increasing temperature and lead to the overall negative temperature dependence. Similar to DMHP, the negative barrier might suggest that these reactions proceed through the formation of a van der Waals 1951
dx.doi.org/10.1021/jp107804m |J. Phys. Chem. A 2011, 115, 1946–1954
The Journal of Physical Chemistry A
ARTICLE
Figure 5. CCSD/6-31þþG** optimized structures for five DMMP TS’s.
Table 5. Relevant Characteristics for Important DMMP TS Structures TS
multiplicitya
EA (kcal/mol)
νrepulsive (cm-1)
TS1
3
-0.24
-2008
TS2
3
-1.55
-2010
TS3 TS4
3 3
þ0.76 þ0.92
-1380 -1830
TS5
3
þ1.62
-2038
rH-O(H) (Å) in TS
rC-H (Å) in TS
rC-H (Å) in parent
1.265
1.245
1.090
1.265
1.245
1.091
1.354 1.292
1.199 1.226
1.090 1.090
1.245
1.258
1.090
Abstraction from O-CH3
Abstraction from P-CH3 a
Multiplicity = number of rotationally equivalent hydrogen atoms at 300 K.
O-CH3 site. These calculations suggest (at 300 K) ∼1% of the overall decomposition occurs at the P-CH3 bond.
complex21,22 as in eq 16: OH þ ðCH3 OÞ2 PðOÞCH3 T HO - ðCH3 OÞ2 PðOÞCH3 vdw f H2 O þ CH3 OPðOÞðCH3 ÞOC• H2
ð16Þ
The decomposition then occurs from the van der Waals complex via breaking of the C-H bond to form H2O and the alkyl radical. The quality of the theoretical model for DMMP is demonstrated in the following: (1) The quantitative values are in good agreement with and, in many cases, within the experimental error bars from multiple laboratories. (2) The negative temperature dependence of the overall kinetics has been properly captured. Atkinson and coworkers22 reported an Arrhenius expression for DMMP with hydroxyl radical of kOH(T) = 6.25 10-14 e(1538(112)/T cm3 3 molecule-1 3 s-1 for 283 < T < 348 K. Both the theoretical and experimental results indicate greater negative temperature dependence for DMMP compared to DMHP [i.e., (Ea/R)expt = 1538 vs 474]. For DMMP, the fitted activation barrier of -3.06 kcal/mol is quite different from the predominant transition-state activation barrier of -1.55 kcal/mol, so in this case including the contribution of the dynamics to the rate constant calculation is essential. (3) The partitioning of the degradation products is again in excellent agreement. Atkinson and co-workers22 estimated that most, if not all, of the overall reaction occurred at the
’ DISCUSSION Accurate calculation of thermal reaction rate constants for molecules of modest size, say a dozen or so atoms, is essential for the description of transport and dispersion of many environmentally significant chemicals. In this study, we have examined two organophosphorus compounds, DMHP and DMMP, and their interaction with hydroxyl radicals that constitutes the first step in the tropospheric oxidation pathway. While the current method is computationally intensive, requiring the use of high-level ab initio quantum chemistry to calculate accurate structures of molecules at the important stationary points along the minimum energy path that connects reactants to products, it does allow the semiclassical calculation of thermal rate constants that are unlikely to be calculable by fully quantum mechanical techniques. The computed energies of the stationary points are critically important as an error of (0.4 kcal/mol in the computed activation barrier height (the difference in the computed energies of the saddle point and reactant) can lead to a doubling or halving of the computed rate constant at 298 K. In view of this demand, an energy extrapolation was used that accounts for the wellunderstood and predictable deficiencies in the QC methods used. The only test of the computed activation barriers available to us is the calculated rate constant curves. On the basis of the reported thermal rate constants and their good agreement with 1952
dx.doi.org/10.1021/jp107804m |J. Phys. Chem. A 2011, 115, 1946–1954
The Journal of Physical Chemistry A
ARTICLE
Table 6. kOH as a Function of Temperature from Experimental and Theoretical Calculations for DMMP theoretical kOH 1011 (cm3 3 molecule-1 3 s-1) experimental kOH 1011 (cm3 3 molecule-1 3 s-1) temp (K)
ref 22
TS2
TS3
TS4
250
0.247
1.52
0.0167
260
0.229
1.33
0.0167
270
0.212
1.18
0.197
2.28 ( 0.15
278
1.83 ( 0.14
278
1.92 ( 0.21
280 283 288 296
1.04 ( 0.06
0.0233
0.0021
1.81
0.0235
0.0022
1.60
0.0165
0.0236
0.0023
1.43
1.04
0.0163
0.0236
0.0024
1.28
0.183
0.925
0.0161
0.0235
0.0025
1.15
0.170
0.826
0.0158
0.0232
0.0026
1.04
0.158
0.740
0.0154
0.0230
0.0027
0.939
1.03 ( 0.095c 0.632 ( 0.029b
1.12 ( 0.05
0.977 ( 0.096a 0.936 ( 0.117a
298 300 303 310
0.591 ( 0.122c
313
0.447 ( 0.097c
320
total
1.17 ( 0.17c
298 298
TS5
1.52 ( 0.06
290
0.776 ( 0.062
0.147
0.664
0.0151
0.0226
0.0027
0.852
330
0.137
0.599
0.0147
0.0222
0.0028
0.775
340
0.128
0.540
0.0143
0.0217
0.0028
0.707
0.119
0.489
0.0138
0.0212
0.0028
0.646
348
0.551 ( 0.031
350 a
P-CH3
TS1
278
other
O-CH3
See ref 23. b See ref 24. c See ref 27.
than the reactant energy. While this is indicative of a van der Waals complex, the agreement with experimental reaction rate constants suggests that this influence is appropriately accounted for in the SCFFAF approach. The discrepancy between the calculated and experimentally fitted activation barriers for the DMMP hydrogen abstraction reactions highlights the need for a full treatment of the dynamics of the system. Sensitivity of the reaction rate constants to temperature and the fact that the reactions slow as the temperature increases are important input information for atmospheric transport and dispersion modeling. The use of semiclassical dynamics coupled with high-level QC methods can allow for the examination of the tropospheric oxidation pathway for environmentally important molecules, as has been demonstrated. Figure 6. Experimental and theoretical rate constants for the reactions of OH radicals with DMMP.
experiment, our energy extrapolation appears to have performed well. The computed overall reaction rate constants [kOH(T)] for DMHP and DMMP compare well with experimental literature values over the measured temperature range.22,25 Furthermore, our predicted partitioning of the various possible pathways toward observed reaction products is in agreement with experimental observation. For each of the reactants, the dominant contribution to the total reaction rate constants arises from transition states whose activation energies are less
’ AUTHOR INFORMATION Corresponding Author
*E-mail:
[email protected]. Telephone: (321) 7757537.
’ ACKNOWLEDGMENT This work was supported by the United States Air Force. While this research has been funded by this agency, the results and content of this publication do not necessarily reflect the views and opinions of the funding agency. 1953
dx.doi.org/10.1021/jp107804m |J. Phys. Chem. A 2011, 115, 1946–1954
The Journal of Physical Chemistry A
’ REFERENCES (1) Atkinson, R.; Arey, J. Chem. Rev. 2003, 103, 4605–4638. (2) Atkinson, R. J. Phys. Chem. Ref. Data, Monograph 1 1989, 1–246. (3) Finlayson-Pitts, B. J.; Pitts, J. N., Jr. Chemistry of the Upper and Lower Atmosphere. Theory, Experiments, and Applications; Academic Press: New York, 2000. (4) The Pesticide Manual, 9th ed.; Worthing, C. R., Hance, R. J., Eds.; British Crop Protection Council: Farnham, U.K., 1991. (5) Toy, A. D.; Walsh, E. N. Phosphorus Compounds in Everyday Living; American Chemical Society: Washington, DC, 1987. (6) Schmidt, M. W.; Baldridge, K. K.; Boatz, J. A.; Elbert, S. T.; Gordon, M. S.; Jensen, J. H.; Koseki, S.; Matsunaga, N.; Nguyen, K. A.; Su, S. J.; Windus, T. L.; Dupuis, M.; Montgomery, J. A. J. Comput. Chem. 1993, 14, 1347–1363. (7) Stanton, J. F.; Gauss, J.; Watts, J. D.; Nooijen, M.; Oliphant, N.; Perera, S. A.; Szalay, P. G.; Lauderdale, W. J.; Kucharski, S. A.; Gwaltney, S. R.; Beck, S.; Balkova, A.; Bernholdt, D. E.; Baeck, K. K.; Rozyczko, P.; Sekino, H.; Hober, C.; , Bartlett, R. J. ACES II, a program product of the Quantum Theory Project, University of Florida. Integral packages included are VMOL (J. Alml€of and P. R. Taylor), VPROPS (P. Taylor), and ABACUS (T. Helgaker, H. J. Aa. Jensen, P. Jørgensen, J. Olsen, and P. R. Taylor). (8) 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 C.02; Gaussian Inc., Wallingford, CT, 2004. (9) Møller, C.; Plesset, M. S. Phys. Rev. 1934, 46, 618–622. (10) Hehre, W. J.; Ditchfield, R.; Pople, J. A. J. Chem. Phys. 1972, 56, 2257–2261. (11) Francl, M. M.; Petro, W. J.; Hehre, W. J.; Binkley, J. S.; Gordon, M. S.; DeFrees, D. J.; Pople, J. A. J. Chem. Phys. 1982, 77, 3654–3665. (12) Pople, J. A.; Nesbit, R. K. J. Chem. Phys. 1954, 22, 571–572. (13) Bartlett, R. J.; Purvis, G. D., III. Int. J. Quantum Chem. 1978, 14, 561–581. (14) Pople, J. A.; Krishnan, R.; Schlegel, H. B.; Binkley, J. S. Int. J. Quantum Chem. 1978, 14, 545–560. (15) Curtiss, L. A.; Raghavachari, K.; Pople, J. A. J. Chem. Phys. 1993, 98, 1293–1298. (16) (a) Gonzalez-Lafont, A.; Truong, T. N.; Truhlar, D. G. J. Phys. Chem. 1991, 95, 4618–4627. (b) Liu, Y.-P.; Lynch, G. C.; Truong, T. N.; Lu, D.-H.; Truhlar, D. G.; Garrett, B. C. J. Am. Chem. Soc. 1993, 115, 2408–2415. (c) Rossi, I.; Truhlar, D. G. Chem. Phys. Lett. 1995, 233, 231–236. (d) Corchado, J. C.; Truhlar, D. G. ACS Symp. Ser. 1998, 712, 106–127. (e) Chuang, Y.-Y.; Tadhakrishnan, M. L.; Fast, P. L.; Cramer, C. J.; Truhlar, D. G. J. Phys. Chem. A 1999, 103, 4893–4909. (17) Sekusak, S.; Cory, M. G.; Bartlett, R. J.; Sabljic, A. J. Phys. Chem. A 1999, 103, 11394–11405. (18) Dewar, M. J. S.; Thiel, W. Theor. Chim. Acta 1977, 46, 89–104. (19) Wang, H.; Thompson, W. H.; Miller, W. H. J. Chem. Phys. 1997, 107, 7194–7201. (20) Runge, K.; Cory, M. G.; Bartlett, R. J. J. Chem. Phys. 2001, 114, 5141–5148. (21) Smith, I. W. M.; Ravishankara, A. R. J. Phys. Chem. A 2002, 106, 4798–4807.
ARTICLE
(22) Aschmann, S. M.; Long, W. D.; Atkinson, R. J. Phys. Chem. A 2008, 112, 4793–4799. (23) Andino. Personal communication, 2004. (24) Kleindienst, T. E.; Smith, D. F. Chemical Degradation in the Atmosphere; Final report on the atmospheric chemistry of three important volatile chemical precursors to Applied Research Associates, Inc.; U.S. Air Force Contract F08635-93-C-0020; Armstrong Laboratory Environics Directorate, Tyndall AFB, FL, September 1996. (25) Aschmann, S. M.; Tuazon, E. C.; Atkinson, R. J. Phys. Chem. A 2005, 109, 11828–11836. (26) Martin, P.; Tuazon, E. C.; Atkinson, R.; Maughan, A. D. J. Phys. Chem. A 2002, 106, 1542–1550. (27) Henley, M. V.; Laine, P. L.; Forrester, Q. L.; Calidonna, S. E.; Weber, R. M. Kinetics of gas-phase reactions of DMMP and DIMP with atmospheric oxidants. Abstracts, 58th Southeast Regional Meeting of the American Chemical Society, Augusta, GA, November 2006.
1954
dx.doi.org/10.1021/jp107804m |J. Phys. Chem. A 2011, 115, 1946–1954
