Xiao
Shan
*^{a},
Mark R.
Sambrook
^{b} and
David C.
Clary
*^{a}
^{a}Physical and Theoretical Chemistry Laboratory, Department of Chemistry, University of Oxford, South Parks Road, Oxford, OX1 3QZ, UK. E-mail: xiao.shan@chem.ox.ac.uk; david.clary@chem.ox.ac.uk
^{b}CBR Division, Defence Science and Technology Laboratory (Dstl), Porton Down, Salisbury, SP4 0JQ, UK
First published on 17th December 2019
It is very difficult to perform experiments on the physical parameters for the thermal decomposition of chemical nerve agents such as VX and computations, therefore, are useful. The reaction dynamics of the gas-phase pericyclic hydrogen transfer of the nerve agent VX is studied computationally. The geometries of the stationary structures are calculated at M06-2X/jul-cc-pVTZ level of theory. Single point energy calculations are carried out at the CBS/QB3 level to correct the energy barriers. Canonical reaction rate constants are calculated as a function of temperature. The one-dimensional semiclassical transition state theory is used to analyse the quantum tunneling effects. A reduced-dimensional hindered rotor model is proposed, tested, and applied to calculate the vibrational partition functions. It is found that the ester (O-side) and thioester (S-side) side chains of VX undergo pericyclic H-transfer reactions that result in decomposition of the molecule. The S-side reaction is favoured both kinetically and thermodynamically and dominates the pyrolysis over the temperature range from 600 K to 1000 K. It is predicted that VX completely decomposes in 2 s at temperatures above 750 K.
In order to overcome these challenges, we first employ density functional theory (DFT) for the quantum chemistry calculations. In our previous studies,^{21,22} we tested the performance of various combinations of DFT functionals and basis sets on calculating the quantum tunneling contribution to reaction rate constants of H-transfer reactions of VX simulants. In the current study, we use the functional with the best performance, M06-2X,^{24} and the basis set jul-cc-pVTZ^{25} in the calculation. The quantum tunneling contribution to the reaction rate constants is calculated using the semiclassical transition state theory (SCTST) developed by Miller and co-workers.^{26–31} The SCTST calculates the overall quantum tunneling probability by summing over a series of TS vibrational configuration-dependent tunneling probabilities. The shape of each individual configuration-dependent potential energy barrier is determined by the vibrational frequencies and anharmonic constants of the TS.^{30,32–35} Obtaining the anharmonic constants^{30} and number/density of states^{32} of the TS are the two main computational costs in a full-dimensional (FD) SCTST calculation. These costs increase dramatically as the number of atoms in a system increases.^{36–38} In addition, because of the size of the molecule, a great number of vibrations do not directly contribute to the reaction, and the coupling anharmonic constants between these modes and the reaction mode are very small. In practice, the numerical computations of these anharmonic constants can introduce unnecessary error in a FD calculation. In the current study, we employ the one-dimensional (1D) SCTST,^{21,36,39–41} which is a much less expensive and yet sufficiently accurate approximation to the FD SCTST. We have used the 1D SCTST in a combined experimental and theoretical study of a pericyclic H-transfer reaction of sarin,^{21} similar to the title reactions, and found excellent agreement between the two.
An additional computational challenge one faces when working with relatively large molecules is the treatment of the partition functions of the internal rotational degrees of freedom (DoFs).^{42–46} These DoFs are often referred to as the hindered rotors (HRs), because they behave similarly to a normal vibrational mode at low temperature, while at high temperature they should be treated as a rotation. It has been shown that one could solve the Schrödinger equation for one^{43,47–49} to two^{50,51} dimensional HRs provided that the relative PES is known. Recently, this approach has been extended by Truhlar and co-workers to deal with the multi-dimensional (MD) HRs cases,^{52} where the HR modes are treated as pair-wise coupled. All of the above methods are applicable at a great computational cost especially for systems as large as VX. Alternatively, a relatively more cost efficient method proposed by Truhlar and co-workers,^{43–46,53} the multi-structure all-structure (MS-T) method, can be easily adapted into the SCTST rate constant calculation.^{21,22}
The predicted full two-step mechanism of VX thermal decomposition through pericyclic H transfer reaction has been shown in ref. 22. The current study is focused on the first step reactions; a schematic illustration of the mechanism is shown in Fig. 1. The S- and O-side reactions yield O-ethyl methylphosphonothioic acid (EMPTA) and EA 2192, respectively. EA 2192 is also a possible product of VX hydrolysis,^{8–12} however, EMPTA can only be formed via pericyclic H transfer of VX. Although the decomposition of VX via pericyclic H transfer is somewhat similar to the hydrolysis, especially for the products, it should be noted that the two types of reactions are fundamentally different in the position of the bond cleavage. In hydrolysis, the PO–X (X = O or S) bond is broken forming an alcohol or thioalcohol leaving group, while in pericyclic H transfer, the POX–CH_{2} (X = O or S) bond is broken forming an alkene. EA 2192 formed in the hydrolysis of VX can have either S or R chirality, while the pericyclic H transfer can only product one specific enantiomer. R- and S-VX can yield only S- and R-EA 2192, respectively.
The remainder of the paper is organised in the following way. In Section 2, we briefly discuss the theoretical derivation of 1D SCTST and the implementation of the reduced-dimensional hindered rotor method. The details of computational and experimental methods employed in the current study are given in Section 3. Results and discussions of our measurements and calculations of the IR spectra, as well as calculated stationary point geometries, reaction energetics, reaction rate constants, and the survival percentage of VX are given in Section 4. Our key findings in this study are summarised in Section 5.
(1) |
(2) |
The vibrational perturbation theory (VPTn) to the second order is commonly used for this expression.^{27–29,32,35,54,56,57} Recently, Stanton et al.^{58} showed that it is also possible to apply the VPT4 in a 1D analysis. The VPT4 requires up to 6th order derivatives of the potential energy surface with respect to the vibrational modes, however it does not improve the VPT2 results by a great deal^{58} especially in a high temperature situation such as the pericyclic H transfer. In addition, for systems as large as VX, the computational cost must be taken into consideration. We therefore use the expression of θ(E_{v}) based on the VPT2 calculation, it is given by
(3) |
Ω = Im(ħω_{R}), |
It is known that this VPT2 expression of θ(E_{v}) is inaccurate for reactant or product asymptotes known as the deep tunneling (DT) region.^{62} In the DT correction proposed by Wagner et al.^{62}eqn (3) is replaced by a piece-wise continuous asymmetric Eckart barrier model. It improves the results at temperatures typically lower than 500 K.^{21,36,39–41,54} The additional computational cost of constructing the DT barrier is only the reverse barrier height.
(4) |
For a system that consists a total of J torsional conformers, the conformational rovibrational partition function, Q_{con-rovib}(T), is given by
(5) |
To carry out the MS-T treatment, the information of the Hessian matrix at each of the torsional conformers of a molecule is required. The number of calculations scales up very quickly as the number of HR modes in a molecule increases. For molecules of the size of VX, there are a total of 12 HR modes in VX. Four of them are CH_{3} group rotation and can be approximated as having a C_{3} symmetry, we then have 8 hindered rotors each have 2–3 conformer geometries to analyse, which translates to between 128 (= 2^{8}) and 6561 (= 3^{8}) geometries. In our previous study of decomposition of GB,^{21} we employed a single structure approximation to the MS-T treatment. In the following test study, this method is referred to as the single-structure (SS)-HR. In the SS-HR, instead of calculating all of the torsional conformers, only the lowest energy conformer is considered in the treatment. In practice, it is impossible to find the lowest energy conformer without computing the energies for all conformers. However, the multi-structure reference-structure (MS-RS) method^{45} can usually find a low energy conformer close to the lowest energy one.
In the current study, adapting the reduced dimensionality framework,^{63–65} we propose a method that is somewhat in between SS-HR and MS-T. Before we going into details of this method, it is helpful to discuss the locality of a vibrational mode.^{41,66} For a particular vibrational mode, Q, it is given by
(6) |
(7) |
In general, higher ζ value suggest greater single atom movement in the vibrational mode. Taking a simulant molecule, O,S-diethyl methylphosphonothiolate, as an example, there are two possible pericyclic H-transfer sites from the CH_{3} group on either the ester or the thioester side chain. We show the calculation results for the atomic locality of the reaction mode for the two TSs in Fig. 2. It can be seen that in either cases, the movement of the atoms in this vibrational mode is dominated by the H atom transferring from the attached C atom to the O atom.
Based on this observation, we propose a computationally less expensive HR treatment, the reduced-dimensional HR (RD-HR) method. In this method, only the HR modes in the reactant that change the relative position of the transferring H atom are treated explicitly as in the MS-T method. These set of HR modes are referred to as the active HR modes. At each of the conformers resulting from an active HR mode rotation, the rest of HR modes, (the spectator HR modes) are treated the same as in a simple SS-HR method. In Fig. 3, we illustrated the RD-HR treatment to the O-side TS for the reaction of O,S-diethyl methylphosphonothiolate. The active and spectator HR modes are marked as the red and blue arrows, respectively. For this reaction, using the RD-HR method, the rotations related to the thioester group as well as the CH_{3} group attached to the P-atom are not active. In this way, the number of conformers for the reactant and TS of the reaction is reduced from 16 to 2 and from 6 to 1, respectively. The same can be applied to the S-side reaction yielding 6 and 1 conformers for the reactant and TS, respectively.
Recently, Wu et al.^{67} proposed and tested a multi-structural 2-dimensional torsion (MS-2DT) method to compute accurate partition functions at a smaller computational cost than the MS-T method on C_{6}–C_{8} alkanes. Although it is possible to extend the application to VX, it is still more expensive than our method. For a 5 three-fold rotatable torsions, the MS-2DT method requires 90 initial conformers,^{67} while the RD-HR method based on the MS-RS method would require scanning through only 15 conformer. It should be noted that the goal of our RD-HR method is to sufficiently accurately and efficiently estimating the rate constant of reactions with presence of multiple hindered rotors affecting the mechanism, it should not be used to reproduce accurate MS-T partition functions.
We show in Fig. 4 the comparison between the RD-HR method and MS-T results for both reactions presented in Fig. 2. We also included the results calculated using the SS-HR method and purely harmonic oscillator (HO) models using the reactant and TS conformers of the lowest energy. Details of their structures and energies can be found in ref. 22. Note that in the logarithmic scale, a deviation of 1.0 is 10 times over- or underestimation of the rate constant. It can be seen that the RD-HR has a good agreement to the MS-T results (note that the MS-T results are simply the x-axes of the plots) for temperature typically lower than 1000 K, yielding a rate constant that is within ∼25% when comparing to the MS-T. On the other hand, both the SS-HR and the HO results are much worse comparing to the RD-HR ones. Therefore, in the following study of VX reactions, only the RD-HR method is used.
Fig. 4 Comparison of rate constants for (a) O-side and (b) S-side pericyclic H transfer reaction for O,S-diethyl methylphosphonothiolate calculated using various approximated hindered rotor models. The rate constants for the MS-T results, as well as the geometries, relative energies, vibrational frequencies and anharmonic constants of reactant and TS's can be found in ref. 22. |
S-side TS | O-side TS | R | |
---|---|---|---|
a The values are taken as the average of the respective atomic distances. | |||
O–P | 1.5997 (1.5996) | 1.5733 | 1.6001 |
S–P | 1.9964 (1.9980) | 2.1010 | 2.1065 |
PO | 1.5349 (1.5265) | 1.5379 | 1.4810 |
O–H | 1.2742 (1.2885) | 1.2789 | >4^{a} |
H–C | 1.3427 (1.3317) | 1.3522 | 1.0912^{a} |
O–H–C | 176.04 (173.46) | 169.61 | — |
In Table 2, we report the adiabatic barrier heights for the forward and reverse reactions. The reactions for both sides have significant forward barrier heights over 176 kJ mol^{−1}, approximately 10 kJ mol^{−1} higher than the barrier for the sarin reaction.^{21} Similar to sarin and VX simulants,^{21,22} both reactions are endothermic. The S-side reaction is favoured both kinetically and thermodynamically by having a lower forward barrier height by ∼14 kJ mol^{−1} and a smaller endothermicity by ∼20 kJ mol^{−1}. It should be noted that, unlike in the sarin case,^{21} the two stereoisomers of the S-side TS are very different in energy due to steric hindrance. The corresponding reactions should be treated independently. In addition, due to the existence of the chiral C atom in the TS, the two H atoms attached to it should be treated as in different chemical environments in the reactant partition function calculation especially for treating the HR modes. We also report the relative RD-HR reactant conformer energies in Table 2. We can see that the S-side reaction has more conformers than the O-side, although the two sides have the same number of active HR modes. The rotation of CH_{3} on the O-side ester can be treated as a simple three-fold symmetric top, and hence only the conformers of the other two active HR modes are needed in the calculation. On the other hand, all the possible conformers of the equivalent rotation on the S-side have to be treated individually. Note that this does not increase the number of reactant conformers for the S-side reaction by a great deal due to the steric hindrance of the (i-Pr)_{2}NCH_{2}CH_{2}– group.
S-side | O-side | |
---|---|---|
ΔV_{f} | 176.24 (186.19) | 190.64 |
ΔV_{r} | 127.95 (137.90) | 112.52 |
conf1 | 0.00 | 0.00 |
conf2 | 0.3403 | 0.2281 |
conf3 | 10.23 | 3.954 |
conf4 | 13.36 | 13.09 |
conf5 | 2.799 | 14.13 |
conf6 | 17.34 | — |
conf7 | 11.50 | — |
Fig. 6 Comparison between calculated IR spectrum of VX and experimental results, ©Crown Copyright (2019), Dstl. |
The values of the third and fourth order derivatives of the reaction mode calculated numerically with Richardson extrapolation up to the 2nd order are reported in Table 3. We can see that for the two third order derivatives, the extrapolation does not converge. In addition, large deviation can be found for the 0th order extrapolation results, this is because the barriers for both reactions are significantly asymmetric. As shown in Table 2, for instance, the forward barrier of the O-side reaction is approximately twice as large as the reverse barrier. As the displacement geometry moving further away from the TS, the calculated energies used in the numerical differentiation would have a more realistic representation of the shape of the overall PES. In other words, the Richardson extrapolation intended to improve the value calculated at the smallest displacements would no longer provide a sufficiently accurate result. Therefore, in the current study, we fitted the single point energies calculated at a total of eight displacement geometries with a fourth order polynomial to obtain the higher order derivatives. The results are reported in Table 3. We can see that these higher order derivatives from polynomial fitting for either reaction sides are quite different from the highest order Richardson extrapolated results. The residual sum of squares (RSS) of the fittings are also shown in Table 3, we can see the fittings are very good.
S-side TS | ||||
---|---|---|---|---|
ΔQ | Richardson extrapolation order | Polynomial fitting | ||
0 | 1 | 2 | ||
f _{3,R} × 10^{7} | ||||
0.03 | 5.737 | 6.244 | 6.388 | |
(−7.147) | (−7.260) | (−7.323) | ||
0.06 | 4.216 | 4.084 | — | 4.590 |
(−6.807) | (−6.320) | — | (−8.320) | |
0.12 | 4.613 | — | — | |
(−8.268) | — | — | ||
f _{4,R} × 10^{7} | ||||
0.03 | 3.528 | 3.479 | 3.477 | |
(0.6789) | (0.1280) | (0.02173) | ||
0.06 | 3.674 | 3.508 | — | 4.180 |
(2.332) | (1.722) | — | (4.238) | |
0.12 | 4.173 | — | — | |
(4.160) | — | — |
O-side TS | ||||
---|---|---|---|---|
ΔQ | Richardson extrapolation order | Polynomial fitting | ||
0 | 1 | 2 | ||
f _{3,R} × 10^{7} | ||||
0.03 | 6.838 | 7.202 | 7.277 | |
0.06 | 5.746 | 6.078 | — | 4.826 |
0.12 | 4.750 | — | — | |
f _{4,R} × 10^{7} | ||||
0.03 | 2.417 | 2.354 | 2.349 | |
0.06 | 2.608 | 2.427 | — | 3.156 |
0.12 | 3.149 | — | — |
The calculated vibrational frequencies (ω_{R}) and anharmonic constants (x_{RR}) for the reaction mode vibration and the anharmonic correction (G_{0}) to the PES are presented in Table 4. We show results from both the Richardson extrapolation and fitted polynomial data. It can be seen that for both the S- and O-side reactions, compared to the barrier heights reported in Table 2, the effect of G_{0} on reaction rate constant is negligible. The error in the calculation of G_{0} from using in the Richardson extrapolated higher order derivatives is therefore not significant. However, the errors in x_{RR} are more noticeable. For both reactions, the x_{RR} from the extrapolated derivatives are ∼20% smaller in absolute value comparing to the ones from the fitted derivatives, resulting in a narrower barrier and an overestimated tunneling probability. Therefore, for the remainder of the paper, our results will be based on the higher order derivatives obtained from polynomial fitting.
S-side TS | ||
---|---|---|
ω _{R} | 1457.1i | |
(1452.4i) | ||
x _{RR} | Extrapolation | Fitting |
−113.02 | −132.57 | |
(−7.0731) | (−139.10) | |
G _{0} | Extrapolation | Fitting |
−0.3285 | −0.3915 | |
(−8.408 × 10^{−4}) | (−0.4002) |
O-side TS | ||
---|---|---|
ω _{R} | 1418.1i | |
x _{RR} | Extrapolation | Fitting |
−84.12 | −106.77 | |
G _{0} | Extrapolation | Fitting |
−0.2377 | −0.3130 |
Since both the S-side and O-side reactions are unimolecular, it is easy to calculate the percentage survival of VX (VX %survival) as a function of temperature and reaction time using the following equation:
(8) |
We plot in Fig. 8(a) the total VX %survival as a function of temperature, where we have k(T) = k_{S-side}(T) + k_{O-side}(T). The results of a total of four different reaction times are presented. From observing each individual curve, we can see that the longer the reaction the lower temperature required for complete decomposition. Alternatively, we can also measure the decomposition completeness at a particular temperature after different reaction times. For instance, at 725 K, there are ∼30% and ∼8.8% of VX left after 1 s and 2 s reaction time, respectively. The inset of Fig. 8(a) shows a comparison between the 1D SCTST (solid black curve) and TST (dashed red curve) results after 2 s reaction time. The 1D SCTST estimated rate constants are ∼1.5 times greater than the TST ones in the temperature range shown, and the quantum tunneling effects still have contributions to the reaction. Note that similar results were found for the sarin pericyclic H-transfer reaction,^{21} where the 1D SCTST results had excellent agreement to the experimental data.
Fig. 8 Percentage survival of VX at various temperatures and reaction times for the (a) total reaction after various reaction times and comparison between the TST and 1D SCTST in the inset and (b) comparison between the S-side and O-side reaction and their contribution to the total reaction. The %survival of GB after 2 s calculated using the rate constants reported previously.^{21} |
In Fig. 8(b), we show the VX %survival after 2 s reaction calculated using the total reaction rate constant. The GB %survival estimated based on the rate constant reported in ref. 21 are also shown for comparison. To achieve the same percentage of decomposition VX in general requires a temperature that is ∼ 50 to ∼60 K higher than GB. This agrees with the finding that the smallest barrier found in the VX reactions, 176.24 kJ mol^{−1}, is approximately 10 kJ mol^{−1} higher than the barrier of the GB reaction, which is 166 kJ mol^{−1}.^{21}
In addition, we show the VX %survival after 2 s of the S-side and O-side reaction in Fig. 8(b). Comparing to the total %survival, the reaction is dominated by the S-side reaction, which yields O-ethyl methylphosphonothiolate (EMPTA). Although EMPTA is a VX precursor,^{72} its toxicity is considered to be low. Conversely, cleavage of the ester side-chain furnishes the breakdown product EA 2192, which is considered to be nearly as toxic as VX and to be relatively stable to, for example, hydrolysis.^{1,3} It should also be noted that hydrolysis of VX can proceed via either O-side or S-side chain cleavage, with the ratio of the products formed dependent on the solution pH. For example, Munro et al.^{9} described hydrolysis at pH < 6 and >10 as dominated by P–S cleavage, while pH 7–10 by P–O cleavage. The thermal degradation, presented in the current study, on the other hand shows that the S-side reaction is much more strongly favoured than the O-side, in agreement with results presented by Riches and co-workers^{73} on the treatment of VX-contaminated soil using accelerant-based fires. In addition, it is possible that under the pyrolysis condition both EMPTA and EA 2192 can further breakdown via pericyclic H-transfer reaction to form methylphosphonothioic S-acid and methylphosphonothioic acid, respectively,^{22} while it is difficult to perform further hydrolysis of EA 2192.^{3}
Both the ester (O-side) and thioester (S-side) side chains can carry out pericyclic H-transfer reactions that results in breaking down of the VX molecule. The S-side reaction was favoured both kinetically and thermodynamically and dominated the reaction mechanism in the temperature range of pyrolysis (∼600 K to ∼1000 K). We examined the %survival of the molecules as a function of time and temperature based on the 1D SCTST rate constant. In the gas-phase high pressure limit, at temperature above 750 K, VX would completely decompose after 2 s. It should be noted that the P-atom containing products of all the reactions can carry out further pericyclic H-transfer reactions and decompose to smaller less toxic compounds. The mechanisms and reaction rate constants of some of these decompositions have been reported previously.^{22}
Footnote |
† Electronic supplementary information (ESI) available. See DOI: 10.1039/c9cp05109k |
This journal is © the Owner Societies 2020 |