CASPT2 Study of the Unimolecular Reactions of Nitromethane. A Look at the Roaming Reactions in the Decomposition of Nitromethane: An Exergonic Route at High Temperatures

1. Introduction

Nitro compounds are a class of molecules that play important roles in several areas and applications, such as atmospheric chemistry, explosives, propellants or fuels. [1,2,3] For this reason, the thermal and photochemical decompositions of the simplest organic nitro compound [nitromethane (CH₃NO₂)] have been extensively studied both theoretically [4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20] and experimentally, [18,19,20,21,22,23,24] due precisely to its molecular simplicity and the applications above mentioned. There is consensus in the literature that the weakest bond in the molecule is C-N. [6] However, in spite of the structural simplicity of nitromethane, it is not clear what is the initial reaction or the key reaction step that occurs in the chemically simplest decomposition processes, for example, when acts as an explosive without the intervention of an external reactant, or when is decomposed in IRMPD experiments where the thermal chemistry occurs in collision-free conditions. [21]

The detonation of an energetic material is a process that occurs under the stimuli of an external action that leads to a local chemical reaction releasing an enormous amount of energy which is transformed in internal energy of the material with the consequent local increasing of the temperature and pressure. [14] According to the definition of detonation, the initial local reaction is expected to be an elementary unimolecular process. Thus, in this context, detonation of CH₃NO₂ acting as energetic material remains a major puzzle because all of the elementary reactions of this molecule are highly endothermic (vide infra). In this work, we propose an elementary reaction, that leads to an exergonic decomposition process at high temperatures that passes through the nitro-nitrite isomerization channel via a roaming reaction and ends with the dissociation of cis-methyl nitrite into CH₃O and NO. The roaming mechanism of nitromethane has been previously studied by other authors [7,8,9]; therefore, we assume, as starting hypothesis, that this process takes place in the decomposition mechanism. For the sake of completeness, all of the elementary reactions of nitromethane decomposition have been studied with the CASPT2 approximation.

A roaming reaction is a class of unimolecular reaction that was reported relatively recently. [25] Roaming reactions occur at the near-dissociation limit of the molecule where radical products are almost formed, that is, when the fragments are separated to 3-4 Å, then roaming reorientation of the fragments becomes feasible as the kinetic energy is low and, consequently, the angular forces may be comparable to the radial forces, in other words, roaming rotation of the fragments is almost a free energy process. This effect may allow the system to access a distinct reactive domain with respect to dissociation, such as abstraction or isomerization, with the consequent formation of unexpected highly vibrationally excited products. [19,26] Roaming reactions may be considered a special case of those reactions that are dynamically controlled and deviate from the minimum energy path or intrinsic reaction coordinate [27,28,29,30] and is now widely accepted as a nearly universal aspect of chemical reactivity that could make a significant contribution to product branching in unimolecular reactions. [30,31,32,33,34,35,36]

2. Methods of Calculation

The multi-configurational calculations have been carried out with the complete active space self-consistent field (CASSCF) [37,38,39,40,41,42,43] method and the multi-state second-order perturbation (MS-CASPT2) approach, [44,45] as implemented in the MOLCAS 8.4 and OpenMOLCAS 22.06 programs. [46,47,48,49] Extended relativistic ANO-RCC basis sets [50,51] have been used throughout this work by applying the contraction scheme [C,N,O]/[H]: [4s3p2d1f]/[3s2p1d]. The IPEA empirical correction has been fixed at 0.25 in all of the MS-CASPT2 calculations, equally, to avoid the inclusion of intruder states in such calculations an imaginary shift set to 0.1 has been applied.

All geometry optimizations have been performed with both CASSCF and CASPT2 methods. The characterization of all the species as minimum or transition state has been done by means of frequency calculations, analytically for CASSCF calculations and numerically for the CASPT2 ones. In addition, the Møller-Plesset (MP2) [52] and DFT/M06-2X [53] methods have been applied as they are implemented in GAUSSIAN-16, [54] in conjunction with def2-TZVPP basis sets. [55,56]

The selection of the active space to study the reactions of nitro-nitrite derivatives (X-NO_2; X-=NO) consists of 14 electrons distributed in 11 orbitals, [57,58,59,60,61,62] or 16 electrons and 13 orbitals when proton migration is under study. Molecular dynamics calculations at the CASSCF level were performed with the velocity algorithm of Verlet [63,64] implemented in MOLCAS.

The analysis of molecular orbitals and molecular geometries have been done with the programs Molden [65] and Gabedit. [66] The analysis of vibrational normal modes have been performed with the MacMolplt program. [67]

Construction of the 1D-potential energy curves (PECs) have been performed by applying a linear interpolation method that uses the full space of non-redundant internal coordinates. [68,69,70,71,72]

Mapping of 2D-Potential Energy Surfaces

Mapping of 2D-potential energy surface of a reaction was sketched out in a previous work. [73] Here, we describe the procedure in detail (Figure 1). First of all, to start the mapping, three points (geometrical configurations) are required for defining the reaction domain, such points usually content critical points of the system under study, e.g., minima, transition states, surface crossings or a combination of them. Secondly, we define a common set of non-redundant internal coordinates (${\overrightarrow{\mathit{R}}}_{\mathit{n}}=[{\mathit{r}}_{\mathit{i}},\dots ,{\mathit{\alpha }}_{\mathit{j}},\dots ,{\mathsf{\Gamma }}_{\mathit{k}},\dots ])$ for the three reference points, where ${\mathit{r}}_{\mathit{i}}\mathit{\text{'}}\mathit{s}$ correspond to internuclear distances and ${\mathit{\alpha }}_{\mathit{j}}\mathit{\text{'}}\mathit{s}$ and ${\mathsf{\Gamma }}_{\mathit{k}}\mathit{\text{'}}\mathit{s}$ are valence bond angles and dihedral angles, respectively. Third, one of the selected points is chosen as reference (0 in Figure 1A). Then, we built the linear interpolation vectors ${\overrightarrow{\mathit{R}}}_{\mathbf{01}}$ and ${\overrightarrow{\mathit{R}}}_{\mathbf{02}}$ by subtraction of the i, j non-redundant internal coordinates:

$${\overrightarrow{\mathit{R}}}_{\mathbf{01}}={\overrightarrow{\mathit{R}}}_{\mathbf{1}}-{\overrightarrow{\mathit{R}}}_{\mathbf{0}}$$

(1)

$${\overrightarrow{\mathit{R}}}_{\mathbf{02}}={\overrightarrow{\mathit{R}}}_{\mathbf{2}}-{\overrightarrow{\mathit{R}}}_{\mathbf{0}}$$

(2)

The coordinates of ${\overrightarrow{\mathit{R}}}_{\mathbf{0}\mathit{i}}$ correspond to the following elements: non-dimensional internuclear distances ${\mathit{r}}_{\mathbf{0}\mathit{i}}=({\mathit{r}}_{\mathbf{0}}-{\mathit{r}}_{\mathit{i}})/{\mathit{r}}_{\mathbf{o}}$ (${\mathit{r}}_{\mathit{k}}$ internuclear distance in Å and ${\mathit{r}}_{\mathbf{o}}=1\mathit{\AA }$) and the differences between the ${\mathit{\alpha }}_{\mathbf{0}\mathit{j}}\mathit{\text{'}}\mathit{s}$ and ${\mathsf{\Gamma }}_{\mathbf{0}\mathit{k}}\mathit{\text{'}}\mathit{s}$, respectively, which correspond to non-dimensional valence bond angle coordinates or dihedral angle coordinates of ${\overrightarrow{\mathit{R}}}_{\mathbf{0}}$ and ${\overrightarrow{\mathit{R}}}_{\mathit{i}}$ in radians.

Given that we are working in internal coordinates, the relative orientation of the ${\overrightarrow{\mathit{R}}}_{\mathbf{01}}$ and ${\overrightarrow{\mathit{R}}}_{\mathbf{02}}$ vectors can be arbitrarily chosen. Thus, in this case, the relative orientation is selected in such a way that the bisector line of the angle $\widehat{\mathbf{102}}$ is aligned with the Q_s-axis and is orthogonal to the Q_R-axis (Figure 1A.1).

The next step is the normalization of the vectors ${\overrightarrow{R}}_{01}$ and ${\overrightarrow{R}}_{02}$:

$${\overrightarrow{r}}_{01}={\overrightarrow{R}}_{01}/\left|{\overrightarrow{R}}_{01}\right|$$

(3)

$${\overrightarrow{r}}_{02}={\overrightarrow{R}}_{02}/\left|{\overrightarrow{R}}_{02}\right|$$

(4)

Thus, the addition and subtraction of the unit vectors ${\overrightarrow{r}}_{01}$ and ${\overrightarrow{r}}_{02}$ -with equal lengths- yields vectors ${\overrightarrow{v}}_S$ and ${\overrightarrow{v}}_R$ that are mutually orthogonal and are in line with the Q_s-axis and Q_R -axis, respectively (Figure 1A.2). Normalization of ${\overrightarrow{v}}_S$ and ${\overrightarrow{v}}_R$ yields ${\overrightarrow{q}}_S$ and ${\overrightarrow{q}}_R$ that form an orthonormal basis

$${\overrightarrow{v}}_S={\overrightarrow{r}}_{01}+{\overrightarrow{r}}_{02}$$

(5)

$${\overrightarrow{q}}_S={\overrightarrow{v}}_S/\left|{\overrightarrow{v}}_S\right|$$

$${\overrightarrow{v}}_R={\overrightarrow{r}}_{01}-{\overrightarrow{r}}_{02}$$

(6)

$${\overrightarrow{q}}_R={\overrightarrow{v}}_R/\left|{\overrightarrow{v}}_{SR}\right|$$

In order to span the configurational domains of species 1 and 2, the orthonormal vectors ${\overrightarrow{q}}_S$ and ${\overrightarrow{q}}_R$ are scaled by the norm $\left(f_m\right)$ of the largest vector, ${\overrightarrow{R}}_{01}$ or ${\overrightarrow{R}}_{02}$. Thus, we obtain two vectors ${\overrightarrow{Q}}_s$ and ${\overrightarrow{Q}}_R$, which are mutually orthogonal and with the same norm (Figure 1A.3) as follows:

$$f_m\ge \left|{\overrightarrow{R}}_{0X}\right|X=\left(1or2\right)$$

(7)

$${\overrightarrow{Q}}_S=f_m·{\overrightarrow{q}}_s$$

(8)

$${\overrightarrow{Q}}_R=f_m·{\overrightarrow{q}}_R$$

(9)

$$f_m=\left|{\overrightarrow{Q}}_s\right|=\left|{\overrightarrow{Q}}_R\right|$$

(10)

Once the vectors ${\overrightarrow{Q}}_s$ and ${\overrightarrow{Q}}_R$ are built, in order to map the potential energy surface, it is necessary to select a step size s,

$$s=f_m/n_{step}$$

(11)

where $n_{step}$is an arbitrary integer number that will determine the grid size of the 2D-surface. Thus, the position of any point $(k,l)$, which represents a distortion of the reference geometry 0 ($k·s$ times along the vector ${\overrightarrow{Q}}_S$ and $l·s$ times along the vector ${\overrightarrow{Q}}_R$), is given by the position vector ${\overrightarrow{P}}_{kl}$ (Figure 1A.4).

$${\overrightarrow{P}}_{kl}=k·s·{\overrightarrow{Q}}_S/\left|{\overrightarrow{Q}}_S\right|+l·s·{\overrightarrow{Q}}_R/\left|{\overrightarrow{Q}}_R\right|$$

$$(k,l=\pm 0,\pm 1,\dots ,\pm n_{step})$$

(12)

In practice, we run the mapping in horizontal mode, that is, following the interpolation vector ${\overrightarrow{Q}}_S^l$ (Figure 1A.4). As an example of the application of the mapping procedure described above, Figure 1B,C depict the 2D and 3D-representations of the potential energy surfaces of the ternary system [nitromethane:trans-methylnitrite:cis-methylnitrite]. The mapping of such a system identifies two saddle points, labelled as pTS1 and pTS4, that would connect (i) nitromethane with trans-methyl nitrite and (ii) trans-methyl nitrite with cis-methyl nitrite, respectively. Thus, each one of these saddle points has been successfully used as good starting geometries to optimize the true transition states (TS1 and TS4; vide infra) for the reactions CH₃NO₂ $\to$ trans-CH₃ONO and trans-CH₃ONO $\to$ cis-CH₃ONO, respectively.

To finish this section, it must be remarked that the points used to mapping the PESs satisfy the condition of being contained in a plane because they are constructed from the basis vectors ${\overrightarrow{q}}_S$ and ${\overrightarrow{q}}_R$ that form a subspace of $R^n$. In contrast, the points that would be obtained by the so-called constrained optimization method would not satisfy the above condition because they would not be in the same vector subspace.

3. Results and Discussion

3.1. Unimolecular Reactions of Nitromethane and Methyl Nitrite

The elementary unimolecular reactions of nitromethane and methyl nitrite have been widely studied by many authors with a variety of theoretical approximations,[4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20] mostly with single-determinantal representation for the electrons. In this section, with the objective of showing the accuracy and reliability of the CASPT2 theory, whose reference wave function is multi-determinantal, we have studied some of the main unimolecular reactions of nitromethane and methyl nitrite with such a theory (energetics, frequencies and geometries of stable species, intermediates and transition states), as well as where possible, CASPT2 results are compared with experimental data [74,75,76] and with MP2 and DFT calculations obtained in this work. To be specific, the reactions reported are: (i) dissociation of nitromethane into methyl and nitrogen dioxide radicals; (ii) trans-methyl nitrite dissociation into methoxy and nitric oxide radicals; (iii) nitromethane dissociation into nitrosomethane and atomic oxygen; (iv) nitro-nitrite isomerization of nitromethane leading to formation of trans-methyl nitrite; (v) proton migration in nitromethane to give aci-nitromethane [CH₂N(O)OH]; (vi) proton migration in methyl nitrite to give formaldehyde and nitrosyl hydride (HNO); and (vii) trans-cis isomerization of methyl nitrite. The thermodynamics properties (enthalpies and Gibbs free energies) of such reactions computed with the expressions of Statistical Thermodynamics and obtained from three different theoretical approaches (CASPT2, MP2/HF and DFT) are compared with experimental values in Tables 1,2. In general, there exists an excellent agreement between experimental and CASPT2 calculated values for the reaction enthalpies; the major source of deviation from experimental data arises from the harmonic CASPT2 frequencies. Comparison of the dissociation and activation Gibbs energies for the tabulated reactions shows that the C-N bond cleavage is the most energetically favourable process of nitromethane. As another probe/prove of the accuracy of the CASPT2 method to dealt with this kind of compounds, the comparison of experimental [77,78,79] and calculated geometrical parameters of the species studied in this work for nitromethane, intermediates and products, is given in ESI (Tables S1-S8). Equally, in order to show that the MS-CASPT2 approach is an appropriate method to deal with dissociation reactions, [80,81,82,83,84] the potential energy surfaces of the dissociation reactions collected in Table 1 are given in Figures S1,S2, where is clearly observed as every potential energy curve reaches the asymptotic limit at the dissociation region.

3.2. A Look at the so-Called Loose Transition State

Two different reaction mechanisms for nitro-nitrite isomerization of nitromethane (CH₃NO₂) to methyl nitrite (CH₃ONO) have been proposed in the literature: one involves a so-called loose transition state (Figure 2a), in which the methyl and NO₂ radicals are almost formed and well separated from each other [4] and the other mechanism passes through a so-called tight transition state in which the internuclear distances among the atoms of the C-NO₂ group are shorter than in the loose one (Figure 2b). The geometries of these two types of transition states obtained in this work are given in Figure 2 and the CASSCF and CASPT2 structural parameters are collected in Table 3. It should be noted that due to the very flexible nature of the transition state given in Figure. 2a, we could not find a precise stationary point; although the convergence criteria for maximum gradient and root mean square of gradient were fulfilled; the stopping criterium for displacements was not reached.

The loose transition state has acquired a special interest in the context of the roaming reaction for the nitro-nitrite isomerization of nitromethane, [7] where is proposed a nitro-nitrite isomerization mechanism which evolves through a loose transition state or roaming saddle point (RSP). In this respect, inclusion of such an RSP seems to deviate from the original statement of the definition of roaming reaction, [25] that is, a reaction that does not pass through a conventional transition state. [25,27] For this reason, we have re-investigated these two types of transition states (loose and tight). To obtain such transition states, two different active spaces must be selected: (i) 10 electrons distributed in 7 orbitals for the loose transition state [CAS(10,7)], and (ii) 14 electrons in 11 orbitals [CAS(14,11)] for the tight one. The CASSCF molecular orbitals of the C_s symmetry equilibrium geometry of nitromethane, which are included in each one of the active spaces, are represented in Figure 3. The reduced active space (Figure 3a) comprises the following orbitals that coincide with the description of Saxon and Yoshimine:⁵ one O σ lone pair (n_σ), one O π lone pair (n_π), two N-O σ bonds (σ_NO), one N-O σ∗ bond (σ∗_NO), one N-O π bond (π_NO) and one N-O π* (π∗_NO). The larger active space (Figure 3b) comprises one O σ lone pair (n_σ), one O π λone pair (n_π), two N-O σ bonds (σ_NO) and the two correlating σ∗ orbitals (σ∗_NO), one N-O π (π_NO) and the correlating N-O π* orbital (π∗_NO), one C-N σ bond (σ_CN) and the correlating σ∗ orbital (σ∗_CN), and the 2s-orbital of nitrogen (2s_N) that has a strong C-N bonding character. Although both types of calculations have been done without symmetry restrictions, the CAS(14,11) wave function of the larger active space keeps the C_s symmetry properties. In contrast, as a result of the unbalanced selection of the CAS(10,7) reduced active space; it is clearly observed in the orbitals that we have obtained a symmetry breaking solution of the wave function. Symmetry breaking has an important impact on the results, generally, such results are unappropriated. For example, Figure 4 represents the 1D-potential energy surfaces that connects the equilibrium geometry of nitromethane with the geometry of the loose transition state (Figure 2a), obtained with the two active spaces, CAS(10,7) and CAS(14,11), and with the linear interpolation method. [68,69,70,71,72] The curve obtained with the CAS(10,7) space shows a discontinuity between points 21 and 22, the reason for this behaviour is orbital rotations (Figure 3) among the inactive, active and secondary spaces, what, in practice, is equivalent to a change in the wave function in passing from one point to the other on the reactive domain, which is, from our point of view, chemically unacceptable to describe any reaction. Therefore, the loose transition state represented in Figure 2a must be considered an artefact of the unappropriated CAS(10,7) active space.

Figure 5a,b depict the 2D and 3D-representations of the potential energy surfaces [CASPT2/CASSCF(14e, 11o)] that comprises the configurational domain of the ternary system [CH₃NO₂:Loose:cis-CH₃ONO], the reference geometries correspond to the CASPT2 minimum geometries of nitromethane and cis-methyl nitrite plus the CASSCF loose transition state. It is clearly observed in such a figure that there is not any signal of a transition state (saddle point) which relates nitromethane with cis-methyl nitrite. In contrast, seems that they are related through a roaming mechanism.

To gain more insights about the roaming structures of the nitro-nitrite isomerization of nitromethane, as in Figure 5, we have built four CASPT2 2D-potential energy surfaces for the ternary system [CH₃NO₂-Roaming-Structure-cis-CH₃ONO] by taking different geometries for the quasi-dissociated species (Figure 6) and using a CASSCF(14, 11) reference wave function. The first surface (Figure 6a), in fact, would not correspond to a roaming intermediate because the analysis of the wave function indicates that the two fragments are well formed radicals, that is, the molecule is completely dissociated. The geometries and relative energies of the reference points (0_n) are included in the graphic. It is shown that a significant change in the geometry of the reference point is accompanied by a small variation of its relative energy, that is, roaming rotation must not be a very hindered process at these points of the potential energy surface, especially if rotation is assumed to be free at the higher energy geometry (0₁). Figure 6b–d represent the topology of the potential surfaces for true roaming structures. What it is observed in these latter figures is that as the roaming complex is more compact with decreasing internuclear distances, the product (cis-CH₃ONO) is closer to the roaming structure, in addition, the slope of the surface becomes more negative, which enhances the formation of the intermediate [cis-CH3ONO]* from the roaming structure. To finish this paragraph, the brown flat regions on the surfaces depicted in Figure 6 correspond to very high energy points arising from geometries with very close separation between atoms

To demonstrate that the roaming process leads to decomposition of nitromethane into CH₃O and NO passing through the methyl nitrite intermediate, we have performed molecular dynamics calculations at the CASSCF level with the algorithm of Verlet. [63,64] The obtained results are shown in Figure 7. These calculations are started at the geometry represented in Figure 6c (0₃). The initial condition imposed on the dynamical calculations consists of addition of kinetic energy on the NO₂ fragment as is shown in Figure 7. The pattern of formation/dissociation of cis-methyl nitrite is similar in the four calculations: (i) the nitrite molecule is formed from the roaming structure, more rapid as higher is the initial kinetic energy; (ii) the molecule lives during approximately 80 fs in the well of cis-methyl nitrite; (iii) the molecule dissociates into NO and CH₃O. (iv) Excess of kinetics energy (E_kin = 80 kcal/mol) leads to a delay in the dissociation process of [cis-CH3ONO]*.

The global process starting at nitromethane is represented in equation 13

CH₃NO₂→ [CH₃^…NO₂]* → [CH₃ONO]* → CH₃O + NO

(13)

Figure 8 collects diagrammatically the Gibbs free energies at four different temperatures of the process represented in equation 13 and Figure 6b–c. In accordance with such diagrams, the activation energy of the first step decreases as temperature increases and is the higher energy step in every case. Therefore, when the system reaches the region of the roaming rearrangement, it has accumulated enough energy to complete reaction 13. Curiously, the activation energies of the process that passes through the roaming molecular arrangement with a C-N internuclear distance of 4.0 Å (Figure 6c) are the lower ones at any temperature (Figure 8b). Gibbs free energies at different temperatures have been calculated with the standard expressions of Statistical Thermodynamics, the thermal corrections to the energies of the roaming structures have computed after projecting out the vibrational Hessian matrix the rotation-translation eigenvectors and gradient vectors. [85,86]

4. Conclusions

It is described a mapping method (orthogonalized 3D-representation) for the PESs by defining an orthonormal basis consisting of two $R^n$ orthonormal vectors that allows to obtain a set of ordered points in the vector subspace defined by the plane that contains the basis vectors.

The main elementary unimolecular reactions of nitromethane -including energetics, geometrical optimizations and calculated vibrational harmonic frequencies- have been studied with the CASPT2 method by using a well-balanced CASSCF reference wave function of 14(16) electrons distributed in 11(13) orbitals. All the kinetically controlled reactions are endergonic, being the C-N bond breaking the process with lower ΔG energy among this class of reactions. In contrast, it is found that the roaming process [CH₃NO₂→ [CH₃^…NO₂]* → [CH₃ONO]* → CH₃O + NO] is exergonic at high temperatures, where the star symbols indicate vibrationally excited species. At this point, it is important to note that this mechanism has been previously proposed by other authors. [7,8,9] Furthermore, all the roaming structures analyzed in this work, the higher energy points in the mechanism given in equation 13, are well below the dissociation limit of nitromethane that leads to CH₃ and NO₂.

The so-called loose transition state has been re-investigated with two CASSCF active spaces: (i) 10 electrons distributed in 7 orbitals and (ii) 14 electrons in 10 orbitals. It is shown that the smaller active space yields artifactual results arising from an unbalanced description of the electronic structure. In fact, the selection of the active space for studying nitromethane would be straightforward after the work of Blahous et al. on NO₂ radical, [87] where they pointed out that the CASSCF wave function to describe the whole configurational domain of the potential energy surface must contain 13 electrons distributed in 10 orbitals.

To finish, it must be remarked that an exhaustive study of a roaming process requires intensive and extensive molecular dynamics calculations, such as those performed in relevant works given in references cited above. [7,8,9,10] However, this aim is beyond the scope of this manuscript.