Calculations on the unimolecular decomposition of the nerve agent VX

Abstract

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.

---

1 Introduction

The V-series of organophosphorus (OP) nerve agents are highly toxic to organisms including humans due to inhibition of the acetycholinesterase enzyme.^1–3 Various V-series agents were developed into chemical warfare agents (CWAs).^1,3 It is, therefore, important to study and mitigate against such synthetic toxic molecules. Numerous researches have been devoted to investigate different methods to degrade V-series agents including, among others, surface sorption,^4–7 hydrolysis,^8–12 perhydrolysis,¹³ novel catalysts,^14–18 binding to metal ions,¹⁹ and incineration. In recent studies, the thermally activated pericyclic H-transfer reaction in gas-phase via a 6-membered-ring transition state (TS) has been investigated for the G-series OP nerve agent sarin (GB),^20,21 and for VX simulants.²² Note that this mechanism is only possible for gas-phase reactions; in the solution state this is not a feasible degradation pathway.^1,3 In combination with earlier studies,²³ the results for unimolecular decomposition showed that the 6-membered-ring TS is energetically favoured with a second most stable TS somewhere over 100 kJ mol⁻¹ higher on the potential energy surface (PES). The structure of VX allows such calculations to be carried out on both the ester and thioester side chains. In our previous work,²² we studied the reaction mechanism of VX using a series of simulant molecules. We found that in the simulant molecule where both ester and thioester side chains are present, the thioester side reaction has a lower forward barrier and smaller endothermicity. The thioester side reaction is therefore favoured both kinetically and thermodynamically. In the current paper, we extend this work to study the reaction of VX itself. To the best of our knowledge, this study is the first attempt on calculating reaction rate constants, especially quantum tunneling contribution to rate constants for gas-phase VX thermal decompositions via the pericyclic H transfer route. The major difficulties in performing such calculation is the computational costs due to the size of the molecule, which for VX is 42 atoms. The challenge comes in two folds. On one hand the ab initio quantum chemistry calculations performed for a single geometry point on the PES for these relatively large systems are already expensive, on the other hand, in order to construct a reasonable representation of the PES to analyse the quantum tunneling effect adds another dimension into the cost.

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,²⁴ and the basis set jul-cc-pVTZ²⁵ 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³⁰ and number/density of states³² 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,²¹ 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,⁵² 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₂ (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.

Fig. 1 Schematic of the mechanisms of the pericyclic H-transfer reactions of VX.

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.

2 Theory background

2.1 One-dimensional semiclassical transition state theory

In this study, we use one-dimensional semiclassical transition state theory (1D SCTST)^36,54 to calculate the canonical rate constants, k(T), at high pressure limit. We showed in a previous study²¹ that our method is also applicable to the calculation of the microcanonical rate constant, k(E). Recent developments of SCTST, including the 1D approximation can be found in ref. 55. Here we only give a short outline of the method. In 1D SCTST, k(T) is given bywhere β = k_B⁻¹T⁻¹. The reaction probability, P(E_v) is given byand θ is the penetration integral expressed as a function of the energy, E_v.

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.⁵⁸ 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⁵⁸ 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

where

Ω = Im(ħω_R),

Here, the subscript “R” denotes the reaction vibrational mode, and f₃ and f₄ are the third and fourth order derivatives of the potential energy with respect to this vibrational mode. To obtain the values of f₃ and f₄, one can perform numerical differentiation using either the Hessian matrix elements at 2 displaced geometries or single point energies at 4 displaced geometries. In previous studies,^21,22,36 the Richardson extrapolation^59–61 was used to improve the accuracy of the numerical derivatives. However, as our recent study of the reactions with simulant molecules²² showed that these reactions are strongly asymmetric, it is possible that the Richardson extrapolation fails to produce the correct derivatives. In the current study, we therefore use both the second order Richardson extrapolation and fitting a fourth order function to the energies at displacement geometries to calculate the higher order derivatives.

It is known that this VPT2 expression of θ(E_v) is inaccurate for reactant or product asymptotes known as the deep tunneling (DT) region.⁶² In the DT correction proposed by Wagner et al.⁶²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.

2.2 Approximated hindered rotor treatment

In large molecular systems, hindered rotor (HR) vibrations become a prominent feature. It has been shown that it is incorrect to treat them using harmonic oscillator models.^42–45 Although it is possible to calculate the exact 1D^21,40,43 or 2D^50,51 HR partition functions, the treatment to multi-dimensional (MD) HR remains a challenge. For large molecules, the challenge comes from both the computational cost of constructing the MD HR potential energy surface as well as solving the MD HR Schrödinger equation. The multi-structural all structure (MS-T) method developed by Truhlar and co-workers⁴⁵ provides a relatively inexpensive yet sufficiently accurate solution to the problem. Using the MS-T method, the 1D SCTST rate constant can then be written as

For a system that consists a total of J torsional conformers, the conformational rovibrational partition function, Q_con-rovib(T), is given by

where Q_rot,j and Q_HO,j are the rotational and harmonic oscillator vibrational partition functions of the j-th conformer, respectively. U_j is the relative energy between the j-th conformer and the most stable conformer. c_HR,j is the MD HR correction to the vibrational partition function. Details of how this correction is carried out can be found elsewhere.^22,45

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₃ group rotation and can be approximated as having a C₃ symmetry, we then have 8 hindered rotors each have 2–3 conformer geometries to analyse, which translates to between 128 (= 2⁸) and 6561 (= 3⁸) geometries. In our previous study of decomposition of GB,²¹ 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⁴⁵ 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

where the contribution from the i-th atom, C_i, is given by

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₃ 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.

Fig. 2 Plots of atomic locality for the reaction mode vibration of the TSs of a simulant molecule, O,S-diethyl methylphosphonothiolate, with (a) O-side TS and (b) S-side TS. The atomic coordinates for each TS together with the order of atomic numbers can be found in the supporting information.

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₃ 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.

Fig. 3 Illustration of the active and spectator HR modes of O,S-diethyl methylphosphonothiolate for the O-side reaction. The curved red arrows indicate the rotations across the single bonds in the active HR modes, the curved green arrows are for the inactive HR modes.

Recently, Wu et al.⁶⁷ 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₆–C₈ 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,⁶⁷ 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.

3 Practical details

3.1 Computational details

In this study, we used Gaussian09 software⁶⁸ for the quantum chemistry calculations. The geometry optimizations and vibrational analyses were performed using the DFT method M06-2X.²⁴ The jul-cc-pVTZ basis set, developed by Truhlar and co-workers,²⁵ was used. We used an ultrafine integration grid and the “tight” option in Gaussian09 as the optimization condition. The starting low energy conformer geometry was found using the MS-RS method.⁴⁵ Single point energy calculations were then carried out at each of the optimized geometries at the CBS/QB3^69,70 level. Note that we set the geometry optimization threshold to a high value in the CBS/QB3 input to avoid reoptimization of the input geometry that is usually carried out by this method at B3LYP level; and only the resulting single point energy values were used in our calculation. The anharmonic constant calculations for the reaction vibrational mode in the TS were carried out using second order vibrational perturbation theory (VPT2),²⁹ with the information of the third and fourth order derivatives of the potential energy surface with respect to the reaction mode. These higher order derivatives were calculated by numerically differentiating the single point energies obtained at displaced geometries at M06-2X/jul-cc-pVTZ level of theory. A total of 8 displacement geometries for each TS were used in the anharmonic constant calculation, and a step-size of 0.03a₀Da^0.5 was chosen.

3.2 Experimental details for IR spectra of VX

The IR spectrum of VX liquid (ca. 2 μL; >80% purity) was collected on a commercially available spectrometer in ATR mode at room temperature.

4 Results and discussions

4.1 Structures and energetics

We show in Fig. 5 the optimized geometries of the stationary structures of the VX reactions. Similar to our findings in previous studies on selected simulants of VX,²² the TSs of both the S-side and the O-side reaction consist of 6-membered-ring structures. Key bond lengths and bond angles are reported in Table 1, they are also very similar to what we have found in the simulant molecules.²² One main difference between the VX and the simulants is that the TS for the S-side reaction, in addition to the P atom, has a chiral centre on the secondary C atom on the thioester chain, where the H is transferred from (shown on Fig. 5). Therefore there are two stereoisomers of the S-side TS, similar to the TS for the pericyclic reaction of sarin.²¹ It is worth noting that these two TSs will result in the same products. The more stable stereoisomer is shown in Fig. 5.

Fig. 5 Optimized geometries of the stationary structures of the VX reactions.

Table 1 Table of key bond lengths and angles of the O-side and S-side TSs for the VX reactions. Bond lengths are in Angstroms and bond angles are in degrees. The geometries of the stereoisomer of the S-side TS are in parentheses

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
P=O	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⁻¹, approximately 10 kJ mol⁻¹ higher than the barrier for the sarin reaction.²¹ 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⁻¹ and a smaller endothermicity by ∼20 kJ mol⁻¹. It should be noted that, unlike in the sarin case,²¹ 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₃ 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)₂NCH₂CH₂– group.

Table 2 Table of adiabatic forward and reverse barrier heights for the VX reactions and the relative conformer energies for each reaction paths. All values are in the unit of kJ mol⁻¹. The energy barriers for the stereoisomer of the S-side TS are in parentheses

	S-side	O-side
ΔVf	176.24 (186.19)	190.64
ΔVr	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	—

---

4.2 Vibrational analysis

Our calculated IR spectrum of VX is compared to experimental data in Fig. 6. A scaling factor of 0.971⁷¹ is used to correct the M06-2X/jul-cc-pVTZ frequencies. Overall a good agreement is found between the two spectra. Note that the calculation results include hindered rotor modes at a frequency typically below 600 cm⁻¹ and these modes are not shown in the experimental spectrum. The two spectra have the best agreement for vibrations between ∼700 cm⁻¹ to ∼1500 cm⁻¹, where we can see that the calculation manages to catch almost all the peaks in the experimental results in terms of both peak position and relative intensity. However, for the H-containing bond stretching modes, for which the frequencies are around 300 cm⁻¹, it can be seen that the calculation overestimates by ∼50 cm⁻¹.

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.

Table 3 Table of the third and fourth order derivatives of the PES at the TSs (in a.u.). The values for the stereoisomer of the S-side TS are in parentheses. The RSS of the polynomial fittings are 2.36 × 10⁻¹³, 3.72 × 10⁻¹², and 1.58 × 10⁻¹² for the S-side TS, its stereoisomer and the O-side TS, respectively

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₀) 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₀ on reaction rate constant is negligible. The error in the calculation of G₀ 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.

Table 4 Table of the reaction mode vibrational frequencies (in cm⁻¹), anharmonic constants (in cm⁻¹), and anharmonic correction to the energy barriers (in kJ mol⁻¹). The values for the stereoisomer of the S-side TS are in parentheses

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

---

4.3 Rate constants

In Fig. 7(a), the logarithmic plots of the reaction rate constants for both the S- and O-sides are presented. All rate constant data used to construct the plots are given in the Supporting Information. We can see for either reaction sides, the corresponding 1D SCTST and TST rate constants converge at high temperature as expected. At low temperatures, typically below 300 K, we can see significant quantum tunneling contributions. We take the S-side reaction as an example. At 300 K, the 1D SCTST rate constant is ∼17.4 times larger than the corresponding TST result. At temperatures ranged from 400 K to 1000 K, the ratio between 1D SCTST and TST calculated reaction constant reduces from ∼3.89 to ∼1.24. These results show that quantum tunneling still plays an important role in the reaction mechanism. As seen in our previous work on sarin,²¹ only when the quantum tunneling contributions at relatively high temperatures are accounted for, the resulting %survival rate of the molecule agrees with experimental data. A similar quantum tunneling effect can also be found for the O-side reaction. We compare the rate constants for the two reactions between 700 K and 1000 K in Fig. 7(b). The rate constant of the S-side reaction is ∼10.03 times greater than the O-side one at 1000 K and becomes greater as temperature decreases. At 700 K, it is ∼22.49 times greater. This is expected due to the smaller forward barrier for the S-side reaction.

Fig. 7 Calculated reaction rate constants: (a) comparison between the TST and 1D SCTST results for both reaction sides from 200 K to 2000 K and (b) comparison between the two reaction sides from 700 K to 1000 K.

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:

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,²¹ 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.²¹

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⁻¹, is approximately 10 kJ mol⁻¹ higher than the barrier of the GB reaction, which is 166 kJ mol⁻¹.²¹

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,⁷² 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.⁹ 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⁷³ 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,²² while it is difficult to perform further hydrolysis of EA 2192.³

5 Conclusions

In this work, we investigated theoretically the thermal decomposition of the nerve agent VX. The canonical reaction rate constants over a wide range of temperatures were calculated. Quantum tunneling effects were included in the calculation using the one-dimensional semiclassical transition state theory (1D SCTST). A large number of hindered rotor (HR) vibrational modes were observed in VX. In order to study these reactions efficiently while maintaining sufficient accuracy, we developed the reduced-dimensional (RD)-HR method. In this method, only a subset of the HRs, which directly affects the position of the reacting H atom are treated explicitly. In this study, we tested the RD-HR method on a simulant molecule then applied it to VX.

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.²²

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

X. S. acknowledges the use of the University of Oxford Advanced Research Computing (ARC) facility in carrying out the theoretical studies in this work. X. S. and D. C. C. acknowledge support from DSTL Grant No. R1000121682. M. R. S. thanks Dr Linda Lee and the Chemical Sensing team at Dstl (CBR Division) for the provision of spectral data. Content includes material subject to Crown copyright (2018), Dstl. This material is licensed under the terms of the Open Government Licence except where otherwise stated. To view this licence, visit http://www.nationalarchives.gov.uk/doc/open-government-licence/version/3 or write to the Information Policy Team, The National Archives, Kew, London, TW9 4DU, or email psi@nationalarchives.gsi.gov.uk.

References

1. K. Kim, O. G. Tsay, D. A. Atwood and D. G. Churchill, Chem. Rev., 2011, 111, 5345–5403.
2. H. Rice, C. H. Dalton, M. E. Price, S. J. Graham, A. C. Green, J. Jenner, H. J. Groombridge and C. M. Timperley, Proc. R. Soc. A, 2015, 471, 20140891.
3. Y. J. Jang, K. Kim, O. G. Tsay, D. A. Atwood and D. G. Churchill, Chem. Rev., 2015, 115, PR1–PR76.
4. A. R. Wilmsmeyer, W. O. Gordon, E. D. Davis, D. Troya, B. A. Mantooth, T. A. Lalain and J. R. Morris, J. Phys. Chem. C, 2013, 117, 15685–15697.
5. D. E. Taylor, K. Runge, M. G. Cory, D. S. Burns, J. L. Vasey, J. D. Hearn, K. Griffith and M. V. Henley, J. Phys. Chem. C, 2013, 117, 2699–2708.
6. B. N. Papas, I. D. Petsalakis, G. Theodorakopoulos and J. L. Whitten, J. Phys. Chem. C, 2014, 118, 23042–23048.
7. N. Q. Le, C. E. Ekuma, B. I. Dunlap and D. Gunlycke, J. Phys. Chem. C, 2018, 122, 2832–2839.
8. Y.-C. Yang, Acc. Chem. Res., 1999, 32, 109–115.
9. N. B. Munro, S. S. Talmage, G. D. Griffin, L. C. Waters, A. P. Watson, J. F. King and V. Hauschild, Environ. Health Perspect., 1999, 107, 933–974.
10. J. J. Kiddle and S. P. Mezyk, J. Phys. Chem. B, 2004, 108, 9568–9570.
11. E. Gershonov, I. Columbus and Y. Zafrani, J. Org. Chem., 2009, 74, 329–338.
12. G. O. Bizzigotti, H. Castelly, A. M. Hafez, W. H. B. Smith and M. T. Whitmire, Chem. Rev., 2009, 109, 236–256.
13. G. W. Wagner and Y.-C. Yang, Ind. Eng. Chem. Res., 2002, 41, 1925–1928.
14. R. Zboril, M. Andrle, F. Oplustil, L. Machala, J. Tucek, J. Filp, Z. Marusak and V. K. Sharma, J. Hazard. Mater., 2012, 211–212, 126–130.
15. G. W. Peterson and G. W. Wagner, J. Porous Mater., 2014, 21, 121–126.
16. L. Bromberg, W. R. Creasy, D. J. McGarvey, E. Wilusz and T. A. Hatton, ACS Appl. Mater. Interfaces, 2016, 7, 22001–22011.
17. J. Zhao, D. T. Lee, R. W. Yaga, M. G. Hall, H. F. Barton, I. R. Woodward, C. J. Oldham, H. J. Walls, G. W. Peterson and G. N. Parsons, Angew. Chem., Int. Ed., 2016, 55, 13224–13228.
18. R. Gil-Sam-Millan, E. López-Maya, M. Hall, N. M. Padial, G. W. Peterson, J. B. DeCoste, L. M. Rodríguez-Albelo, J. E. Oltra, E. Barea and J. A. R. Navarro, ACS Appl. Mater. Interfaces, 2017, 9, 23967–23973.
19. I. Bandyopdhyay, M. J. Kim, Y. S. Lee and D. G. Churchill, J. Phys. Chem. A, 2006, 110, 3655–3661.
20. T. Ash, T. Debnath, T. Banu and A. K. Das, Chem. Res. Toxicol., 2016, 29, 1439–1457.
21. X. Shan, J. C. Vincent, S. Kirkpatrick, M. D. Walker, M. R. Sambrook and D. C. Clary, J. Phys. Chem. A, 2017, 121, 6200–6210.
22. X. Shan, M. R. Sambrook and D. C. Clary, J. Phys. Chem. A, 2019, 123, 59–72.
23. L. Yang, R. M. Shroll, J. Zhang, U. Lourderaj and W. L. Hase, J. Phys. Chem. A, 2009, 113, 13762–13771.
24. Y. Zhao and D. G. Truhlar, Theor. Chem. Acc., 2008, 120, 215–241.
25. E. Papajak, J. Zheng, X. Xu, H. R. Leverentz and D. G. Truhlar, J. Chem. Theory Comput., 2011, 7, 3027–3034.
26. W. H. Miller, J. Chem. Phys., 1975, 62, 1899.
27. W. H. Miller, Faraday Discuss. Chem. Soc., 1977, 62, 40–46.
28. W. H. Miller, Annu. Rev. Phys. Chem., 1990, 41, 245–281.
29. W. H. Miller, R. Hernandez, N. C. Handy, D. Jayatilaka and A. Willetts, Chem. Phys. Lett., 1990, 172, 62–68.
30. M. J. Cohen, N. C. Handy, R. Hernandez and W. H. Miller, Chem. Phys. Lett., 1992, 192, 407–416.
31. R. Hernandez and W. H. Miller, Chem. Phys. Lett., 1993, 214, 129–136.
32. T. L. Nguyen, J. F. Stanton and J. R. Barker, J. Phys. Chem. A, 2011, 115, 5118–5126.
33. J. R. Barker, T. L. Nguyen and J. F. Stanton, J. Phys. Chem. A, 2012, 116, 6408–6419.
34. R. E. Weston, Jr., T. L. Nguyen, J. F. Stanton and J. R. Barker, J. Phys. Chem. A, 2013, 117, 821–835.
35. T. L. Nguyen and J. F. Stanton, J. Chem. Phys., 2017, 147, 152704.
36. S. M. Greene, X. Shan and D. C. Clary, J. Chem. Phys., 2016, 144, 244116.
37. C. Aieta, F. Gabas and M. Ceotto, J. Phys. Chem. A, 2016, 120, 4853–4862.
38. C. Aieta, F. Gabas and M. Ceotto, J. Chem. Theory Comput., 2019, 15, 2142–2153.
39. X. Shan and D. C. Clary, Philos. Trans. R. Soc., A, 2017, 376, 20170147.
40. T. A. H. Burd, X. Shan and D. C. Clary, Chem. Phys. Lett., 2018, 693, 88–94.
41. T. A. H. Burd, X. Shan and D. C. Clary, Phys. Chem. Chem. Phys., 2018, 20, 25224.
42. Y.-Y. Chuang and D. G. Truhlar, J. Chem. Phys., 2000, 112, 1221–1228.
43. B. A. Ellingson, V. A. Lynch, S. L. Mielke and D. G. Truhlar, J. Chem. Phys., 2006, 125, 084305.
44. J. Zheng and D. G. Truhlar, Phys. Chem. Chem. Phys., 2010, 12, 7782–7793.
45. J. Zheng, T. Yu, E. Papajak, I. M. Alecu, S. L. Mielke and D. G. Truhlar, Phys. Chem. Chem. Phys., 2011, 13, 10885–10907.
46. J. Zheng and D. G. Truhlar, J. Chem. Theory Comput., 2013, 9, 1356–1367.
47. K. S. Pitzer and W. D. Gwinn, J. Chem. Phys., 1942, 10, 428–440.
48. K. S. Pitzer, J. Chem. Phys., 1946, 14, 239–243.
49. J. E. Kilpatrick and K. S. Pitzer, J. Chem. Phys., 1949, 17, 1064–1075.
50. A. Fernández-Ramos, J. Chem. Phys., 2013, 138, 134112.
51. L. Simón-Carballido and A. Fernández-Ramos, J. Mol. Model., 2014, 20, 2190.
52. L. Simón-Carballido, J. L. Bao, T. V. Alves, R. Meana-Pañeda, D. G. Truhlar and A. Fernández-Ramos, J. Chem. Theory Comput., 2017, 13, 3478–3492.
53. J. Zheng, S. L. Mielke, K. L. Clarkson and D. G. Truhlar, Comput. Phys. Commun., 2012, 183, 1803–1812.
54. S. M. Greene, X. Shan and D. C. Clary, J. Chem. Phys., 2016, 144, 084113.
55. X. Shan, T. A. H. Burd and D. C. Clary, J. Phys. Chem. A, 2019, 123, 4639–4657.
56. T. L. Nguyen and J. R. Barker, J. Phys. Chem. A, 2010, 114, 3718–3730.
57. S. M. Greene, X. Shan and D. C. Clary, J. Phys. Chem. A, 2015, 119, 12015–12027.
58. J. F. Stanton, J. Phys. Chem. Lett., 2016, 7, 2708–2713.
59. L. F. Richardson and J. A. Gaunt, Philos. Trans. R. Soc., A, 1927, 226, 299–361.
60. C. Ridders, Adv. Eng. Software, 1982, 4, 75–76.
61. W. H. Press, S. A. Teukolsky, W. T. Vetterling and B. P. Flannery, Numerical Recipes in C: The Art of Scientific Computing, Cambridge University Press, Cambridge, 2nd edn, 1992.
62. A. F. Wagner, J. Phys. Chem. A, 2013, 117, 13089–13100.
63. J. M. Bowman, Adv. Chem. Phys., 1985, 61, 115–167.
64. J. M. Bowman, J. Phys. Chem., 1991, 95, 4960–4968.
65. J. M. Bowman, Theor. Chem. Acc., 2002, 108, 125–133.
66. P. T. Panek and C. R. Jacob, J. Chem. Phys., 2016, 144, 164111.
67. J. Wu, H. Ning, X. Xu and W. Ren, Phys. Chem. Chem. Phys., 2019, 21, 10003–10010.
68. M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Robb, J. R. Cheeseman, G. Scalmani, V. Barone, G. A. Petersson, H. Nakatsuji, X. Li, M. Caricato, A. Marenich, J. Bloino, B. G. Janesko, R. Gomperts, B. Mennucci, H. P. Hratchian, J. V. Ortiz, A. F. Izmaylov, J. L. Sonnenberg, D. Williams-Young, F. Ding, F. Lipparini, F. Egidi, J. Goings, B. Peng, A. Petrone, T. Henderson, D. Ranasinghe, V. G. Zakrzewski, J. Gao, N. Rega, G. Zheng, W. Liang, M. Hada, M. Ehara, K. Toyota, R. Fukuda, J. Hasegawa, M. Ishida, T. Nakajima, Y. Honda, O. Kitao, H. Nakai, T. Vreven, K. Throssell, J. A. Montgomery, Jr., J. E. Peralta, F. Ogliaro, M. Bearpark, J. J. Heyd, E. Brothers, K. N. Kudin, V. N. Staroverov, T. Keith, R. Kobayashi, J. Normand, K. Raghavachari, A. Rendell, J. C. Burant, S. S. Iyengar, J. Tomasi, M. Cossi, J. M. Millam, M. Klene, C. Adamo, R. Cammi, J. W. Ochterski, R. L. Martin, K. Morokuma, O. Farkas, J. B. Foresman and D. J. Fox, Gaussian 09, Revision E.01, Gaussian, Inc., Wallingford CT, 2016.
69. J. A. Montgomery, Jr., M. J. Frisch, J. W. Ochterski and G. A. Petersson, J. Chem. Phys., 1999, 110, 2822–2827.
70. J. A. Montgomery, Jr., M. J. Frisch, J. W. Ochterski and G. A. Petersson, J. Chem. Phys., 2000, 112, 6532–6542.
71. H. Zhang, X. Zhang, D. G. Truhlar and X. Xu, J. Phys. Chem. A, 2017, 47, 9033–9044.
72. R. M. Black and J. M. Harrison, in The Chemistry of Organophosphorus Chemical Warfare Agents, ed. F. R. Hartley, John Wiley & Sons, Ltd, 2006, vol. 4, ch. 10, pp. 781–840.
73. M. R. Gravett, F. B. Hopkins, A. J. Self, A. J. Webb, C. M. Timperley and J. R. Riches, Anal. Bioanal. Chem., 2014, 406, 5121–5135.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: 10.1039/c9cp05109k

---