Energy transfer, pre-reactive complex formation and recombination reactions during the collision of peroxy radicals

Abstract

In this paper we study collisions between polyatomic radicals – an important process in fields ranging from biology to combustion. Energy transfer, formation of intermediate complexes and recombination reactions are treated, with applications to peroxy radicals in atmospheric chemistry. Multi-reference perturbation theory, supplemented by coupled-cluster calculations, describes the potential energy surfaces with high accuracy, including the interaction of singlet and triplet spin states during radical recombination. Our multi-reference molecular dynamics (MD) trajectories on methyl peroxy radicals confirm the reaction mechanism postulated in earlier studies. Specifically, they show that if suitable pre-reactive complexes are formed, they will rapidly lead to the formation and subsequent decomposition of tetroxide intermediates. However, generating multi-reference MD trajectories is exceedingly computationally demanding, and we cannot adequately sample the whole conformational space. To answer this challenge, we promote the use of a novel simplified semi-empirical MD methodology. It assumes the collision is governed by two states, a singlet (S₀) and a triplet (T₁) state. The method predicts differences between collisions on S₀ and T₁ surfaces, and qualitatively includes not only pre-reactive complex formation, but also recombination processes such as tetroxide formation. Finally, classical MD simulations using force-fields for non-reactive collisions are employed to generate thousands of collision trajectories, to verify that the semi-empirical method is sampling collisions adequately, and to carry out preliminary investigations of larger systems. For systems with low activation energies, the experimental rate coefficient is surprisingly well reproduced by simply multiplying the gas-kinetic collision rate by the simulated probability for long-lived complex formation.

---

1 Introduction

Peroxy radicals (also referred to as peroxyradicals, as well as peroxide or peroxyl radicals), with the general structure RO₂, are ubiquitous intermediates formed in the oxidation of almost any organic compound in the presence of molecular oxygen. Their reactions are therefore important in diverse contexts, including atmospheric chemistry, combustion, and biochemistry.

Due to their lack of reactivity with O₂, as well as the absence of rapid (i.e. ≫1 s⁻¹ at 298 K) unimolecular reaction channels,¹ RO₂ in the gas phase are relatively long-lived compared to most other free radicals, such as alkyl or alkoxy radicals. Their chemistry is also more versatile, with a given RO₂ typically having several different competing radical propagation and termination channels. While unimolecular isomerization reactions have recently^2–4 been shown to be important for highly functionalized RO₂ especially in clean atmospheric conditions, bimolecular reactions form the main sink for most RO₂ species in the atmosphere. The main bimolecular reaction partners for atmospheric RO₂ are nitrogen monoxide (NO) in polluted conditions, and the hydroperoxyl radical (HO₂) in clean conditions.¹ Self- and cross-reactions are rarely the dominant sink for gas phase RO₂ outside of laboratory studies. Nevertheless, this reaction class has recently received increasing attention, as it is one of the only known gas-phase routes for forming low-volatility accretion products relevant to atmospheric aerosol formation.^3,5

Russell⁶ already suggested in 1957 that the reaction of two peroxy radicals proceeds via the formation of a tetroxide (RO₄R′) intermediate. A schematic of the tetroxide formation is shown for two methyl peroxy radicals in Fig. 1. In the original Russell hypothesis, the reaction channels available to a given RO₂ + R′O₂ pair correspond to different structural rearrangements of the RO₄R′ tetroxide. Later computational studies by us and others^7–10 have nuanced this picture somewhat, as they indicate that rather than rearranging directly to the products, all tetroxides inevitably decompose, yielding two alkoxy (RO) radicals and O₂.

Fig. 1 Schematic of the formation of CH₃O₄CH₃ from two colliding methyl peroxy radicals, proceeding via the formation of a prereactive complex.

Individual peroxy radicals are open-shell systems, and as such radical pairs can couple as either a ground-state singlet (denoted S₀) or in an excited triplet state (T₁). However, the tetroxide can only be formed on the S₀ potential energy surface, and the O₂ must be formed in its ground triplet state for the reaction to be energetically feasible. Therefore, the two alkoxy radicals must accordingly be coupled as a triplet. Different RO₂ + R′O₂ reaction channels then correspond to different fates of this coupled alkoxy radical pair, including dissociation to yield the free RO + R′O radicals, hydrogen shifts forming alcohol and carbonyl products (R = O + R′OH), and intersystem crossing (ISC) to the singlet surface, followed by recombination to form the peroxide accretion product ROOR′ (the O₂ co-product is assumed to be ejected, as it is only very weakly bound to the alkoxy radical pair).

Previous computational studies^7–10 have focused only on the energetics of the possible reaction channels, and computed e.g. activation energies (barrier heights) corresponding to energy differences between various minima and transition states (saddle points) between them. However, to our knowledge the actual dynamics of peroxy radical self- and cross-reactions have not been studied before. Using molecular dynamics (MD) to study this system is very challenging. The reactions themselves require the use of computationally expensive multi-reference methods, possibly combined with coupled-cluster methods, in order to be described correctly. This makes the simulation of long trajectories of a sufficient number to be statistically representative prohibitively time-consuming. On the other hand, empirical force field models can be used to generate a large number of long trajectories, but they cannot model bond making and breaking processes. Semi-empirical methods can in principle describe the reactions, even including their multi-reference character correctly, and they are much more computationally efficient than the purely ab initio wavefunction-based methods. However they may not be accurate enough or sufficiently well parameterized for practical use. In this study, we use multiple different methods to investigate the dynamics of the interaction of two peroxy radicals. We focus primarily on the initial steps of the reaction, including both collisional energy transfer, and the possible formation, lifetime and fate of weakly bound pre-reactive complexes. These all affect the pre-exponential part of the overall reaction rate coefficient, as expressed e.g. by the Arrhenius equation, and thus complement previous studies focusing on activation energies.

We employ three complementary approaches to compute molecular dynamics trajectories for the interaction and reaction of two peroxy radicals. First, our multi-reference MD trajectories on methyl peroxy self-reactions confirm the previously postulated reaction mechanism involving the rapid formation and decomposition of tetroxide intermediates, and further illustrate the key role played by the pre-reactive complexes. Then, to sample many more trajectories while taking into account the spin state of the system, we promote the use of a simplified semi-empirical MD methodology, which relies on singlet (S₀) and triplet (T₁) single-reference potential energy surfaces. While the focus of the present paper is on peroxy radical self-reactions, we note that this methodology is applicable more generally, to any collisions of polyatomic radicals.

Finally, force-field simulations verify that the semi-empirical method is sampling collisions adequately. We also use force-field based MD to lay the basis for investigations of larger peroxy systems, t-butyl peroxy ((CH₃)₃COO), and acetonyl peroxy (CH₃COCH₂OO). The choice of the two larger systems is motivated by the fact that the overall self-reaction rate of acetonyl peroxy radicals is high, while that of t-butyl peroxy radicals is exceptionally low: over 20 times higher, and 160 000 times lower, respectively, than that of the methyl peroxy system.¹ The acetonyl peroxy radical self-reaction is also the smallest known system able to produce ROOR accretion products in the gas phase.⁵

The outline of the paper is as follows: in the Computational methods section we discuss in detail and separately the four methods that we used for computing the potential energy surface of the two radicals, as well as the details of the MD simulations. The Results and discussion section presents the MD simulations of the radicals obtained using three different methods, while highlighting the similarities and differences between the different simulations. In the Concluding remarks section we compare our findings resulting from the different methodologies and present some new insights into the dynamics of the peroxy radical reactions.

2 Computational methods

2.1 Multi-reference perturbation theory benchmarked by coupled-cluster calculations

The gas phase dynamics of the colliding CH₃OO radicals are highly dependent on the spin state of the system, as the outcomes of triplet state (T₁) collisions may differ significantly from collisions in the singlet state (S₀). Furthermore, the electronic structure of the system, whether in the S₀ or T₁ state, is of a complex nature due to the interactions of different possible electronic configurations. Such interactions are non-negligible in open shell systems such as the system of interest. Thus, the method of choice for computing the potential energy surface of the two radicals must be an electronic structure-based method that accounts for the different possible spin states, and that also incorporates configuration interaction to account for the complex electron correlations.

Our method of choice which can correctly describe all of the relevant configurational space is the extended quasi-degenerate 2nd-order multi-reference perturbation theory (XMC-QDPT2).¹¹ This theory takes into account both static and dynamic electronic correlation. The minimum-energy structure of the pre-reactive complex of CH₃OO⋯CH₃OO on the singlet surface was obtained using the XMC-QDPT2 method, an active space with 10 electrons in 8 molecular orbitals (MOs), and the 6-311++G(d,p) basis set. This active space includes the bonding σ-MOs of O–O bonds, two lone pair MOs, two MOs for the π-bonding of O–O bonds and two anti-bonding σ*-MOs located on each radical. This active space is quite stable during the geometry optimisation and also can be applied to determine the transition state for tetroxide formation, which according to our previous study lies around 1.5 kcal mol⁻¹ above the pre-reactive complex.⁹ Subsequent MD calculations were carried out using the same approach, as described in Section 2.4.

Note that calculations using smaller active spaces (for example, 4 electrons in 4 MOs) were not stable, leading to swapping orbitals between active and inactive spaces and finally to incorrect gradient values and molecular dynamics. Also, we should note that only including static electronic correlation in the gradient calculation of the molecular dynamics scheme did not lead to RO₄R formation, presumably since it does not correctly include the dispersion interaction between two radicals.

To benchmark the XMC-QDPT2 results, especially with respect to the energy of the tetroxide intermediate relative to the isolated free peroxy radicals (a quantity difficult to compute accurately with multi-reference approaches), coupled-cluster calculations were conducted with the explicitly correlated R(O)CCSD(T)-F12 level of theory. We used the diagonal 3C ansatz with fixed amplitudes, and geminal exponent β = 0.9 with the cc-pVDZ-F12 basis set.^12,13 The underlying Hartree–Fock/cc-pVDZ-F12 wavefunction was perturbed by generating 100 input files with 15 random orbital rotations between occupied and unoccupied orbitals, and subsequently solving the SCF solution. Orbital rotations were done in order to drive the Hartree–Fock (HF) solution from a possible local minimum to the global minimum.¹⁴ The lowest energy HF wavefunction was then used as the starting point for CCSD(T)-F12 calculations. The F12a solution was used instead of F12b, as suggested for basis sets smaller than triple-ζ.¹⁵ All orbital rotation and coupled-cluster calculations were done with MOLPRO software.¹⁶

Molecular geometries for coupled-cluster calculations were generated by thoroughly investigating the conformational space of different molecular species and selecting the global minimum conformer. The general approach is one we have developed in our previous studies of larger peroxy radical complexes, where considerable conformational sampling is necessary.^9,10 First, a list of plausible conformers was generated by scanning over relevant torsional degrees of freedom as a function of the MMFF94 energy,¹⁷ using Spartan 16 software.¹⁸ Then, single-point energies were calculated using B3LYP density functional theory (DFT) with the 6-31+G(d) basis set. The number of relevant conformers was narrowed down by uniqueness filtering and a 5 kcal mol⁻¹ energetic cutoff from the apparent global minimum. The remaining conformers were then further optimized at the B3LYP/6-311+G(d) level, and the conformers were further narrowed down by filtering and a 2 kcal mol⁻¹ cutoff. Finally, the remaining conformers were optimized and vibrational analysis was carried out with the ωB97XD level of theory and the aug-cc-pVTZ basis set. All DFT calculations were done with Gaussian 16 software.¹⁹

2.2 Semi-empirical calculations of the collision dynamics and recombination of the peroxy radicals

2.2.1 Potential energy calculations. The XMC-QDPT2 method, as described in Section 2.1, has been demonstrated to give quantitatively and qualitatively reliable results for peroxy radicals and tetroxides, and the transition state between them.⁹ However, this first-principles method is very computationally demanding, and therefore more computationally efficient methods which can still describe at least some of the relevant states are desirable so that enough trajectories can be generated for statistical analysis.

Hence, to complement the XMC-QDPT2 results and to extend their scope, we investigated the use of the orthogonalization corrected methods including Grimme's dispersion corrections, including the multi-reference configuration interaction (ODMx/MRCI). We calculated the potential energy surfaces using different active spaces and different sets of reference states. The ODMx methods are based on orthogonalization-corrected MNDO theory with additional dispersion correction including three-body terms,^20–22 but ODM2 and ODM3 each use different parameterization procedures to determine the best fit semi-empirical parameters.^23–26 Multi-reference effects can also be included via a GUGA-CI approach.²⁴ Potential energy surfaces generated with these methods should, in principle, yield qualitatively correct results assuming a large enough active space and the relevant reference states are used. It may also allow for generation of a sufficient number of trajectories for statistical analysis of the results.

As seen in our results for potential energies of different states of two peroxy radicals (available in the ESI†), including MRCI in the semi-empirical calculations does have advantages over the single-reference ODMx methods. However, longer trajectories generated using MRCI suffered from instabilities, such as discontinuities in the energy or failure at the self-consistent field (SCF) part of the calculation.

We also compared the results from the ODM2 and ODM3 methods. Compared with high-level results, ODM3 predicts the tetroxide minimum in reasonable agreement, but with a higher barrier for formation. ODM2 predicts a deeper tetroxide minimum, but a similar small barrier for recombination ∼1 kcal mol⁻¹ in good agreement with the XMC-QDPT2 value (∼1.5 kcal mol⁻¹). The barrier height is a crucial quantity for observing the tetroxide formation. Therefore, in this study we focused our molecular dynamics simulations on the ODM2 single-reference method, since it has been shown to describe peroxides and related Criegee intermediates well,^27,28 and because it is computationally cheap.

Another consideration is the choice between restricted and unrestricted Hartree–Fock wavefunctions (RHF, UHF) for the singlet surface. As shown in the ESI,† RHF wavefunctions are not a good description for the separated radicals, or the pre-reactive complex. Therefore, we have used UHF wavefunctions for all of our ODM2 simulations on both S₀ and T₁ surfaces.

2.2.2 A semi-empirical model using single reference wavefunctions for T₁ → S₀ transitions. A single-reference potential energy surface cannot describe all of the active processes in collisions between two radicals. One obvious difficulty is with the possible occurrence of non-adiabatic transitions. On the other hand, more general and rigorous methods such as the XMC-QDPT2 algorithm discussed in Section 2.1 are extremely difficult to apply for all but the simplest chemical systems. In view of this we propose here a computationally simple two state model which is meant to provide a practical approach to the collision problem which can describe all of the various reactive and non-reactive processes addressed here.

The ingredients of this new model are as follows: (1) we describe the system by two potential energy surfaces (electronic states) and these states will be constructed on the basis of the semi-empirical quantum chemical method ODM2. The lower of the two potential energy surfaces will have essentially the character of the ground singlet (S₀) state between the colliders at least in the ODM2 framework. The second “excited state” has a triplet character (T₁) in the framework of ODM2. The calibration of these states will be discussed later. (2) Non adiabatic transitions between the two states are allowed in this limited framework by assuming it happens at a critical distance between the respective radical groups of the collision partners. This certainly significantly oversimplifies the treatment and must be viewed as no more than an empirical replacement of proper surface hopping. We will refer here to the two potential energy surfaces as S₀ and T₁ retaining only these two states.

This geometric condition (transitions happen with a probability of 1.0 once a critical distance between radical groups is reached) was chosen because the spin–orbit coupling between the S₀ and T₁ states is expected to be most sensitive to r_OO, the distance between the terminal oxygen atoms of the colliders. A strong spin–orbit coupling is necessary for the transition to occur, and the coupling increases for short distances between the radicals. Minaev and Yashchuk²⁹ studied the spin–orbit coupling between the singlet and triplet states of O₂, and found sufficiently strong spin–orbit coupling for distances of r_OO ≤ 2.5 Å between the oxygen atoms. We chose this as the critical distance between terminal oxygen atoms in each radical for triplet to singlet transition to occur in our model. To test our choice we carried out simulations for different assumed values of the critical distance. The results are shown in Fig. S3 of the ESI.† If the critical distance is set at r_OO ≤ 2.3 Å the transition probability is negligibly small. If the critical distance is taken to be r_OO ≥ 2.8 Å the transition probability obtained seems unphysically large. We conclude that the value of r_OO = 2.5 Å is compatible with our model.

It may be argued that the T₁ → S₀ transitions that are simulated by this model are somewhat artificial because they do not correspond to a detailed microscopic treatment of these non-adiabatic transitions, e.g. as in the surface hopping approach. However, they do mimic the real physics of the transitions because the model is grounded on physical data (spin–orbit coupling) that reflects the occurrence and rate of electronic transitions directly. Since we strive for a model which will have both significant predictive power as well as computational efficiency, we found this model satisfactory for our purposes.

In principle, the non-adiabatic transitions S₀ → T₁ can be treated by our model in a similar way to T₁ → S₀ transitions. We chose to ignore the “excitation” S₀ → T₁ transitions altogether. These were neglected because our main interest in this study was focused on the possibility of the unreactive triplet surface transitioning over to the reactive singlet surface. In addition, it turns out that for our system the number of recombination events is very small, therefore the inclusion of all possible reactive events is crucial. On the other hand, transitions from S₀ to the unreactive T₁ surface are relatively insignificant. We carried out an equal number of MD trajectories using the model for S₀ and T₁ initial conditions, as in both cases the free energy of the colliders (at the beginning of the interaction, i.e. at large distances) is the same.

2.3 Empirical fixed-charge force fields

Although simple fixed-charge force fields cannot describe bond making and breaking processes, they may still be adequate to describe the formation of pre-reactive complexes. Their simplicity makes the generation of a large number of trajectories convenient, and they may also be easily extended to modelling larger peroxy systems which are not tractable using explicitly quantum mechanical or even semi-empirical methods.

We base the empirical force fields we use for peroxy radicals on the OPLS-AA parameter set.^30–32 We used the LigParGen webserver^33,34 to automatically generate the LAMMPS data file. Since the OPLS-AA parameter set does not explicitly include peroxy radicals, we have also added partial charges, Lennard-Jones parameters, bonds, and angular terms from other work on peroxy radicals.³⁵

2.4 Molecular dynamics

To the greatest extent possible we have used the same methods for propagating the MD trajectories, regardless of the very different methods used to calculate energies.

XMC-QDPT2 trajectories were initiated, starting from the geometry of the minimum-energy pre-reactive complex, by generating random atomic velocities chosen from a Boltzmann distribution at T = 300 K using the Newton-X code.^36,37 The gradient of XMC-QDPT2 was calculated numerically using Firefly software³⁸ and then was used in the integration of Newton's equations using the Newton-X software for each point on the potential energy surface (PES). The integration time step was set to 2 fs, due to the high computational cost of the gradient calculation. Although a smaller time step would be preferable to ensure correct treatment of the XH vibrations, we found no significant difference in the time required for tetroxide formation when using a time step of 0.5 fs.

The adiabatic molecular dynamics simulation was done for 1 ps for each trajectory on the S₀ surface. Note that the T₁ surface is repulsive at short distances according to XMC-QDPT2 calculation and thus reactive collisions occur only on the S₀ surface in the absence of surface-hopping. Therefore we limited our XMC-QDPT2 calculations to the state specific CASSCF (SS-CASSCF).

We also carried out XMC-QDPT2 MD trajectory simulations for the decomposition of the CH₃–O₄–CH₃ tetroxide, using the same size of active space: 10 electrons in 8 MOs. However, the details of the active space are subtly different from the previous set of trajectories, as it contains the bonding σ-MOs and two anti-bonding σ*-MOs of the O–O bonds, and the two π-MOs of nascent O₂. As the main focus of the present article is on the initial collision of two peroxy radicals, we calculated only one trajectory for the tetroxide decomposition, using the same MD parameters as described above. The minimum-energy CH₃–O₄–CH₃ structure, optimized at the XMC-QDPT2(10,8)/6-311++G(d,p) level of theory, was used as the starting point for this trajectory.

For semi-empirical ODM2-based MD simulations we made use of the MNDO99 software package,³⁹ and for OPLS-AA we used a recent version of the LAMMPS simulation package.⁴⁰ In both the semi-empirical ODM2 and OPLS-AA force field simulations, two radicals were initially placed some distance apart (10 Å for ODM2, 40 Å for OPLS-AA), with randomised orientations, and with each atom given an initial velocity consistent with a total molecular temperature of 300 K. We used a timestep of 0.1 fs for the ODM2 simulations and 0.2 fs for the OPLS-AA simulations.

After some equilibration using a Nosé–Hoover thermostat (OPLS-AA) or simple velocity scaling (ODM2), one of the radicals was given a center-of-mass velocity in the direction of the center-of-mass of the other radical and the remainder of the trajectory was simulated within the NVE ensemble. A collision velocity of v_z = 520 m s⁻¹ was used to correspond with the average relative velocity of methyl peroxy radicals at 300 K. This velocity was lowered for simulations of larger peroxy radicals to keep the kinetic energy the same. We also did some simulations with lower collision velocities to measure the impact of temperature on the formation and lifetime of the pre-reactive complex.

The simulations were run for a total time of 10 ps in the case of ODM2, and for 200 to 400 ps in the case of OPLS-AA. This was sufficient time to observe the initial collision of the two radicals, and except in extremely rare cases the eventual dissociation and separation of the pre-reactive complex as well.

The three different MD methodologies have vastly different computational demands, both in the number of processors required and in the computational time. A comparison of the methods is shown in Table S3 of the ESI.† The OPLS-AA simulations require only 6.3 μs per timestep, while the XMC-QDPT trajectories require ∼3 minutes per timestep and multiple processors. ODM2 simulations are ∼200× slower than OPLS-AA but also ∼15 000× faster than XMC-QDPT.

3 Results and discussion

We benchmarked the potential energy surfaces of the different methods by computing potential energies for various relevant states of our systems, from two separate peroxy radicals through the formation of the pre-reactive complex, the RO₄R complex, and where possible, the 2RO + O₂ system. These data are available in the ESI.†

A comparison of the different energetics for recombination shows that ODM2 with UHF wavefunctions is in reasonable accord with the high level XMC-QDPT2 potential. Most importantly, both methods predict that the reaction path for recombination has a chemically insignificant barrier, and thus we expect ODM2 to yield a valid prediction of the recombination rate.

In the rest of this Section, we focus on presenting our molecular dynamics results. Table 1 collects information about the overall number of trajectories computed using all of our different methods.

Table 1 Summary of the number of collision trajectories n_traj generated for different peroxy radical systems, models and initial conditions

System	v z /m s^−1	n traj, XMC-QDPT2	Initial state
Methyl	—	100	Pre-reactive
	—	1	RO4R

---

System	v z /m s^−1	n traj, ODM2	n traj, OPLS-AA
Methyl	520	1000(S0), 1000(T1)	30 000
	400		15 000
	200		15 000
Acetonyl	378		20 000
t-Butyl	378		10 000

---

We begin in Section 3.1 by discussing the multi-reference ab initio results using the XMC-QDPT2 method, since this approach has been previously established as a benchmark method able to accurately simulate the formation of RO₄R from the pre-reactive complex.

The main limitations of our XMC-QDPT2 findings are the relatively small number of trajectories and the difficulty of starting these trajectories from configurations with two radicals far apart. In Section 3.2 we turn our focus to the use of semi-empirical methods and empirical force fields to generate large numbers of collision trajectories from random orientations and study the statistics of these collisions.

In Section 3.3 we compare collisions on the singlet (S₀) surface versus the triplet (T₁) surface and the results of our model for transitions from T₁ → S₀, followed by comparing how the formation of the RO₄R intermediate is described by our semi-empirical approach and by the XMC-QDPT2 method.

We end this section by considering our results using empirical force fields for the collisional energy transfer and binding times of larger peroxy systems in Section 3.4.

3.1 XMC-QDPT2 trajectories

When starting from the pre-reactive S₀ complex CH₃OO⋯CH₃OO, collision of the two radicals, and subsequent formation of the CH₃–O₄–CH₃ tetroxide, occurs within 1 ps in all XMC-QDPT2 trajectories. By choosing random atomic velocities chosen from a Boltzmann distribution at 300 K, we have initiated multiple trajectories from the same initial geometry. Over all 100 trajectories with a timestep of 2.0 fs, the average time for tetroxide formation (defined as the first time when the distance between the terminal oxygens of each peroxy radical r_OO < 1.6 Å) was 863 ± 118 fs. One test trajectory with a shorter timestep of 0.5 fs formed the tetroxide in 875 fs. Snapshots showing examples of the XMD-QDPT2 MD trajectories leading to tetroxide formation, as well as the subsequent decomposition into 2CH₃O + O₂, are shown in Fig. 2.

Fig. 2 Snapshots taken from two MD trajectories showing the recombination reaction to form RO₄R, and the decomposition reaction to form 2CH₃O + O₂, both calculated using the XMC-QDPT2 method.

The XMC-QDPT2 trajectory results validate the mechanism postulated in previous computational studies: peroxy radical self-reactions indeed proceed via the formation of a tetroxide intermediate, which is – at least for the methyl peroxy system – extremely efficient provided that the initial collision forms a pre-reactive complex with the correct geometry. The overall reaction rate is thus likely to be primarily governed by the probability that collisions lead to pre-reactive complexes in the first place, as well as by the lifetime of these complexes, which in turn determines the probability that they reach the specific geometry allowing for tetroxide formation. As a simple order-of-magnitude estimate, the reaction rate coefficient in this framework should then correspond to the product of the gas kinetic collision rate and the simulated probability that a collision forms a sufficiently long-lived complex.

Our single trajectory starting from the optimized structure of the CH₃–O₄–CH₃ tetroxide shows that the breaking of two O–O bonds occurs within 200 fs, followed by the release of O₂, and leaving a complex of CH₃O radicals in the triplet state. Though we caution that further studies with multiple trajectories are likely needed for quantitative confirmation, this result agrees with earlier computational studies,^7–10 and suggests that the modification of the original Russell hypothesis⁶ (i.e. branching between the different reaction channels occurs from CH₃O⋯CH₃O rather than from CH₃–O₄–CH₃) is likely valid.

3.2 Collisional energy transfer and association lifetimes

Using the ODM2 potential energy surface we calculated 1000 trajectories on the excited triplet state (T₁) surface, and 1000 trajectories on the ground state (S₀). As was expected, recombination occurs only on the S₀ surface, while the T₁ state is repulsive and produces only inelastic collisions. This calls for a comparison between the excited T₁ state collisions and the unreactive portion of the S₀ collisions, as well as between the non-reactive collisions generated with the OPLS-AA force field. The comparison also validates that the MD trajectories performed with ODM2 adequately sample the phase space.

Fig. 3 shows our analysis of unreactive collisions, comparing the ODM2 results with the OPLS-AA results. We note that the ODM2 results here do not include the fraction of T₁ trajectories which underwent transitions to the S₀ ground state, according to our model described in Section 2.2.2. The results shown also do not include those trajectories which led to tetroxide formation. These transitions will be described later in Section 3.3.

Fig. 3 A comparison between unreactive collisions of CH₃OO initially in the triplet (T₁) state or in the singlet state (S₀) using ODM2 (black), and collisions of CH₃OO modelled with the OPLS-AA force-field (red). The collision velocity is 520 m s⁻¹ in both cases. (left) Histograms of the energy transfer during the collision, indicated by the ratio between the initial translational energy of the colliding molecules versus the final translational energy. (right) Histograms of the association time t_assoc of the colliding molecules. Note that the ODM2 data does not include those trajectories that underwent T₁ → S₀ transitions, or that formed the tetroxide.

First, Fig. 3 shows the ratio of the translational kinetic energy of the system, after and before the end and the start of the collision event, i.e. E_trans(f)/E_trans(i). Momentum conservation places a lower limit of 0.5 on E_trans(f)/E_trans(i) for a perfectly inelastic collision between two equal masses, while an elastic collision would conserve translational energy so that E_trans(f)/E_trans(i) = 1.0. These outcomes represent idealized cases; in reality, some transfer of energy from or to internal modes of each molecule can be expected, so that E_trans(f)/E_trans(i) may take on values ranging from 0.5 to well above one.

The profile of the translational energy transfer during collisions between radicals modelled using ODM2 agrees quite closely with the spinless collisions of the OPLS-AA simulations. However, as can be inferred from Fig. 3, in around 80% of ODM2 collisions energy was transferred from the translational degrees of freedom to internal degrees of freedom of the colliders, whereas collisions simulated by the OPLS-AA force field are somewhat less likely (ca. 60%) to lose translational energy.

In Fig. 3 we also show the histograms of association times t_assoc, i.e. the time between the initial collision and the re-separation, for methyl peroxy collisions. Defining r_12,min as the smallest inter-molecular atom–atom distance between the two radicals (hydrogen atoms are not included), the association duration spans from the first time r_12,min < 4.0 Å up to the last time r_12,min < 4.0 Å. This allows for transient oscillations of r_12,min above the distance criterion without ending the association event. In Table 2, we also summarise the data for t_assoc for all of the systems and methods we studied by showing the probability for t_assoc to be <1 ps, >8 ps and >100 ps.

Table 2 Probability of different durations of association times t_assoc between two radicals during unreactive collisions with different collision velocities v_z. 520 m s⁻¹ for the methyl system and 378 for the larger systems correspond roughly to the most likely translational velocity from a Boltzmann distribution at 298 K

System (method)	v z /m s^−1	% chance for tassoc to be:
		<1 ps	>8 ps	>100 ps
Methyl(ODM2)	520	68	0.7	N/A
Methyl(OPLS-AA)	200	41	20	1.7
	400	49	15	0.7
	520	56	11	0.4
Acetonyl(OPLS-AA)	378	9.4	62	8.2
t-Butyl(OPLS-AA)	378	43	24	0.9

---

Overall, the OPLS-AA results for the methyl peroxy system are rather similar to the non-reactive ODM2 collisions, which gives us confidence that the fixed-charge force field methods are giving the correct qualitative picture of the dynamics of the pre-reactive complex. As can be seen in Fig. 3, a large proportion of collisions on the ODM2 as well as on OPLS-AA surfaces quickly re-separate, with t_assoc < 1.0 ps. Hence, the most common scenario after these non-reactive collisions is a very short-lived complex or no complex formation at all.

However, a noticeable difference between the ODM2 and OPLS-AA trajectories is the possibility of some very long association times (up to 200 ps) in a small fraction of OPLS-AA trajectories, and these were not observed on the ODM2 T₁ and S₀ surfaces. For illustration, the data represented by Fig. 3, comprises only ∼90% of the OPLS-AA trajectories, with the rest of the trajectories resulting in complex formation where t_assoc > 8 ps. This difference certainly stems from the slightly more attractive character of the OPLS-AA force field compared to ODM2. However, complex formation on the ODM2 potential is also possible, although rare, as a small fraction of collisions (<1%) resulted in a binding of the radicals for the entire simulation after the initial collision.

We also studied the effect of changing the collision velocity using the OPLS-AA model. Although we only have data for different collision velocities of methyl peroxy, it is likely that these findings are quite general and can be relevant to the larger peroxy systems we discuss in Section 3.4.

The association time histograms in Fig. 4 show that the chance of fast re-separation rises as the collision velocity is increased; at v_z = 200, 400, and 520 m s⁻¹, respectively ∼41, 49 and 56% of the collisions re-separate within 1.0 ps. In the rest of the collisions the radicals remain associated for a longer time. For all collision velocities the distribution of t_assoc shows a broad peak in probability for t_assoc between 1.5 and 2 ps, followed by a monotonic decrease as t_assoc increases. The distributions of association times for the same molecular pairs are not affected very much by the collision velocity when t_assoc > 1.0 ps.

Fig. 4 Histograms of t_assoc for methyl peroxy dimer collisions modelled with the OPLS-AA force field with three different initial velocities.

The probability of overall translational kinetic energy loss increases as the collision velocity v_z increases (results not shown). For the highest incident velocity we studied, the majority of the collisions undergo some amount of kinetic energy loss, as some of the energy is transferred into internal degrees of freedom and/or inter-molecular potential energy.

3.3 Recombination and tetroxide formation

The 1000 trajectories that were computed with ODM2 and started on the excited T₁ surface, were allowed to undergo a transition to S₀ following our crude model for surface hopping described in Section 2.2.2. T₁ → S₀ transitions were observed in 63 trajectories (6.3% of trajectories) which occurred, according to the model, when the distance r_OO between the two terminal oxygen atoms became smaller than 2.5 Å.

Only one out of the 63 collisions that started on the T₁ surface and transitioned to S₀ resulted in recombination. This is also the only collision out of the 1000 simulations that started on T₁ that resulted in recombination, as triplet collisions that did not undergo transition to S₀ were completely unreactive. Snapshots taken from this single reactive trajectory are shown in Fig. 5. It can be seen that following transition to S₀, recombination occurred in less than 0.1 ps.

Fig. 5 Snapshots taken from two collision trajectories in which recombination occurred, representing the two possible routes for recombination in the framework of our crude semi-empirical model.

A larger number of recombination events was recorded for collisions that started on the S₀ surface of ODM2. These trajectories remained exclusively on S₀, as transitions to T₁ were not treated by our model. Recombination occurred in 13 out of 1000 collisions (1.3%). In the majority of cases (11 collisions) the peroxy radicals recombined over the course of only 0.5 ps following the initial association of the radicals. In two collision simulations only, 5.0 and 7.0 ps passed after the initial association and before recombination. This result implies that recombination events that are preceded by a formation of a long-lived pre-reaction complex are less common than recombination events during the first pico-second following the association of the radicals. Fig. 5 shows snapshots from one of the reactive trajectories in which recombination was preceded by the formation of a pre-reactive complex.

Strong support for our crude semi-empirical model comes from its agreement with experiments that measured the self-reaction rate coefficient of methyl peroxy radicals. Assuming that the likelihood of a triplet collision is three times higher than the likelihood of a singlet collision, due to the existence of three configurations for the triplet state while only one for the singlet, we conclude the total recombination probability is

0.25 × (P_r(S₀)) + 0.75 × (P_r(T₁)) = 0.4%, (1)

where P_r is the probability for recombination according to our ODM2 simulations. This is indeed quite close to the ratio between the experimental rate coefficient¹ (on the order of 10⁻¹³ cm³ molecule⁻¹ s⁻¹) and the gas-kinetic collision rate for small molecules (on the order of 10⁻¹⁰ cm³ molecule⁻¹ s⁻¹).

3.4 Larger systems

In Fig. 6 we compare the distributions of t_assoc and the pre- and post-collision translational kinetic energy ratio between the three different species we study, for the largest value of incident velocity (v_z = 520 m s⁻¹ for methyl peroxy, v_z = 378 m s⁻¹ for acetonyl and t-butyl). We plot t_assoc over two different time ranges. First, the left side of Fig. 6b shows the distributions of t_assoc over short times (<8 ps). The right side of Fig. 6b shows a broader range of times, to highlight the rare cases where the binding time is longer. Again, we can see that all species have a significant probability for collisions to immediately result in re-separation of the colliding molecules with t_assoc effectively less than 1 ps. However, this probability is much lower for acetonyl peroxy collisions (ca. 9%) versus the other species (ca. 56% and 43% for methyl and t-butyl, respectively).

Fig. 6 Comparison of gas phase collisions between two colliding molecules of the three systems we have studied (methyl, acetonyl and t-butyl peroxy) modelled with the OPLS-AA force field using our highest collision velocity. (a) Histograms of the energy transfer during the collision, indicated by the ratio between the initial translational energy of the colliding molecules versus the final translational energy. (b) Histograms of the complex association time t_assoc of the colliding molecules.

Looking at larger values of t_assoc, we see that the acetonyl peroxy has a significantly higher proportion of collisions leading to longer association times, with ∼8% chance of t_assoc being larger than 100 ps. In contrast, methyl peroxy and t-butyl peroxy collisions tend to be very short, with only ∼20% of collisions remaining associated for longer than 8 ps and less than 1% for longer than 100 ps.

Approximately 60% of the collisions between methyl peroxy and t-butyl peroxy radicals with v = 520 m s⁻¹ result in the loss of translational energy. This proportion is even higher (∼70%) for the acetonyl peroxy radicals with the same translational kinetic energy, which goes some way towards explaining the higher chance of longer association between these, despite the rather similar potential energy landscape for the pre-reactive complexes at this level of description.

Our finding that the acetonyl peroxy radical pair has the longest average t_assoc qualitatively agrees with the fact that the reaction between acetonyl peroxy radicals has a higher overall rate than the other two reactions.¹ The longer association between the acetonyl peroxy radicals may also help explain why this system is able to produce ROOR accretion products in the gas phase.⁵ In contrast, the similarity of our MD results for methyl and t-butyl peroxy collisions suggests that the over 5 orders of magnitude slower rate coefficient of the t-butyl self-reaction¹ is likely caused by higher reaction barriers (e.g. for the formation of the tetroxide, as suggested previously⁹), rather than by substantial differences in the initial collision dynamics.

It would seem that the main factor driving the increased association times and higher chance of translational energy loss in acetonyl peroxy is the presence of somewhat more attractive interactions between the carbonyl groups and the methyl groups (as reflected by the interaction energies in Table S1, ESI†). In methyl and t-butyl peroxy systems, on the other hand, there is a higher chance that the initial collision results in repulsive interactions between methyl groups, which quickly lead to re-separation.

4 Concluding remarks

Our study focuses on extending our previous work elucidating gas phase peroxy radical interactions and reactions by including molecular dynamics trajectories. We present results from three different methods of computing energies and forces/gradients in MD trajectories: (1) multi-reference perturbation theory (XMC-QDPT2), (2) semi-empirical ODM2 method, and (3) empirical OPLS-AA force-field method.

The XMC-QDPT2 results support our earlier findings from purely energetic considerations,⁹i.e. that the initial reaction between RO₂ radicals rapidly forms a RO₄R tetroxide intermediate on the ground state (S₀) surface, provided that the initial collision forms a pre-reactive complex with the correct geometry.

We have presented a novel, very simplified model to describe surface hopping transitions from the excited triplet (T₁) state to the S₀ state when the distance between terminal oxygen atoms in each radical becomes small. We find that the collisions on the T₁ surface are unreactive, while a small fraction (1.3%) of collisions on the S₀ surface lead to recombination, forming the tetroxide intermediate. Only 0.1% of the total number of trajectories that start on T₁ result in recombination after a transition to the S₀ surface. Overall, our ODM2 results predict that tetroxide formation should result from approximately 0.4% of all methyl peroxy collisions, in good agreement with experimental results indicating that the rate coefficient is roughly three orders of magnitude smaller than the gas-kinetic collision rate.¹

Finally, we present the results of a large number of collision trajectories generated with the empirical OPLS-AA force field. We find that overall, the results for these inherently unreactive collisions agree well with unreactive collisions on the semi-empirical ODM2 surfaces. Furthermore, we note that analogous to the ODM2 case, multiplying the gas-kinetic rate coefficient for the collision of two methyl peroxy radicals (on the order of 10⁻¹⁰ cm⁻³ s⁻¹ at 298 K) by the probability that they form a complex with a lifetime on the order of tens of picoseconds (between 0.1 and 1 percent according to OPLS-AA simulations) leads to a result fairly close to the experimentally measured rate coefficient for methyl peroxy self-reactions (on the order of 10⁻¹³ cm⁻³ s⁻¹ at 298 K).¹ This demonstrates the utility of combining a small number of high-level ab initio trajectories with a much larger number of lower-level simulations – even unreactive simulations – sampling the initial collision dynamics.

We have also been able to use the empirical models to study some larger peroxy systems. We find that, compared to the methyl peroxy system, the acetonyl peroxy radical collisions have a roughly 20 times higher probability of forming pre-reactive complexes with lifetimes longer than 100 picoseconds – in perfect agreement with experimental data indicating a 20 times higher self-reaction rate coefficient for this system.¹ In contrast, the initial collision dynamics of t-butyl peroxy radicals is not very different from that of methyl peroxy radicals, indicating that the large differences in experimentally measured reaction rates must be due to higher energy barriers in the t-butyl system, as already suggested in our earlier work.⁹

Our future work will proceed in two main directions. First, we intend to use high-level ab initio methods to thoroughly explore the final steps of the reaction, where the RO₄R tetroxide breaks up via (at least) three different pathways. One improvement we will consider is the use of relativistic perturbation theory, where energies and gradients are computed using the Dirac equation which includes spin–orbit coupling explicitly.

Second, we will further examine the general utility of our novel simplified semi-empirical method for describing two-state surface hopping during radical collisions. In order to demonstrate the full scope of its utility, we will compare a wider variety of relevant radical systems and their collisions, using both semi-empirical methodologies and empirical force-field models. We are optimistic that this model, despite its crudeness, can be a pathway towards improving our understanding of radical reactivity in the atmosphere.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

T. K., R. V. and C. D. D. are supported by the Jane and Aatos Erkko foundation. T. K. and R. V. are also supported by the Academy of Finland (especially the Centre of Excellence VILMA; grant number 346369). R. B. G. and I. Z. are supported by a grant from the Israel Science Foundation, number 593/20. Computing resources were supplied by Finland's Center for Scientific Computing (CSC).

Notes and references

1. J. J. Orlando and G. S. Tyndall, Chem. Soc. Rev., 2012, 41, 6294–6317.
2. M. Ehn, J. A. Thornton, E. Kleist, M. Sipilä, H. Junninen, I. Pullinen, M. Springer, F. Rubach, R. Tillmann, B. Lee, F. Lopez-Hilfiker, S. Andres, I.-H. Acir, M. Rissanen, T. Jokinen, S. Schobesberger, J. Kangasluoma, J. Kontkanen, T. Nieminen, T. Kurtèn, L. B. Nielsen, S. Jørgensen, H. G. Kjaergaard, M. Canagaratna, M. D. Maso, T. Berndt, T. Petäjä, A. Wahner, V.-M. Kerminen, M. Kulmala, D. R. Worsnop, J. Wildt and T. F. Mentel, Nature, 2014, 506, 476–479.
3. F. Bianchi, T. Kurten, M. Riva, C. Mohr, M. P. Rissanen, P. Roldin, T. Berndt, J. D. Crounse, P. O. Wennberg, T. F. Mentel, J. Wildt, H. Junninen, T. Jokinen, M. Kulmala, D. R. Worsnop, J. A. Thornton, N. Donahue, H. G. Kjaergaard and M. Ehn, Chem. Rev., 2019, 119, 3472–3509.
4. L. Vereecken and B. Nozière, Atmos. Chem. Phys., 2020, 20, 7429–7458.
5. T. Berndt, W. Scholz, B. Mentler, L. Fischer, H. Herrmann, M. Kulmala and A. Hansel, Angew. Chem., Int. Ed., 2018, 57, 3820–3824.
6. G. A. Russell, J. Am. Chem. Soc., 1957, 79, 3871–3877.
7. G. Ghigo, A. Maranzana and G. Tonachini, J. Chem. Phys., 2003, 118, 10575–10583.
8. R. Lee, G. Gryn'ova, K. U. Ingold and M. L. Coote, Phys. Chem. Chem. Phys., 2016, 18, 23673–23679.
9. R. R. Valiev, G. Hasan, V.-T. Salo, J. Kubečka and T. Kurtén, J. Phys. Chem. A, 2019, 123, 6596–6604.
10. G. Hasan, V.-T. Salo, R. R. Valiev, J. Kubečka and T. Kurtén, J. Phys. Chem. A, 2020, 124, 8305–8320.
11. A. A. Granovsky, J. Chem. Phys., 2011, 134, 214113.
12. T. B. Adler, G. Knizia and H.-J. Werner, J. Chem. Phys., 2007, 127, 221106.
13. K. A. Peterson, T. B. Adler and H.-J. Werner, J. Chem. Phys., 2008, 128, 084102.
14. A. C. Vaucher and M. Reiher, J. Chem. Theory Comput., 2017, 13, 1219–1228.
15. D. Feller, K. A. Peterson and J. G. Hill, J. Chem. Phys., 2010, 133, 184102.
16. H.-J. Werner, P. J. Knowles, G. Knizia, F. R. Manby and M. Schütz, Wiley Interdiscip. Rev.: Comput. Mol. Sci., 2012, 2, 242–253.
17. T. A. Halgren, J. Comput. Chem., 1996, 17, 520–552.
18. Wavefunction, Inc., Irvine, CA, Spartan'16.
19. 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. V. 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. J. Bearpark, J. J. Heyd, E. N. Brothers, K. N. Kudin, V. N. Staroverov, T. A. Keith, R. Kobayashi, J. Normand, K. Raghavachari, A. P. 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 16 Revision C.01, Gaussian, Inc., Wallingford, CT, 2016.
20. S. Grimme, J. Comput. Chem., 2004, 25, 1463–1473.
21. S. Grimme, J. Comput. Chem., 2006, 27, 1787–1799.
22. P. O. Dral, X. Wu and W. Thiel, J. Chem. Theory Comput., 2019, 15, 1743–1760.
23. W. Weber and W. Thiel, Theor. Chem. Acc., 2000, 103, 495–506.
24. A. Koslowski, M. E. Beck and W. Thiel, J. Comput. Chem., 2003, 24, 714–726.
25. P. O. Dral, X. Wu, L. Spörkel, A. Koslowski, W. Weber, R. Steiger, M. Scholten and W. Thiel, J. Chem. Theory Comput., 2016, 12, 1082–1096.
26. D. Tuna, Y. Lu, A. Koslowski and W. Thiel, J. Chem. Theory Comput., 2016, 12, 4400–4422.
27. D. Shemesh and R. B. Gerber, Mol. Phys., 2012, 110, 605–617.
28. R. B. Gerber, D. Shemesh, M. E. Varner, J. Kalinowski and B. Hirshberg, Phys. Chem. Chem. Phys., 2014, 16, 9760–9775.
29. B. F. Minaev and L. B. Yashchuk, Opt. Spectrosc., 2003, 95, 553–559.
30. W. L. Jorgensen, D. S. Maxwell and J. Tirado-Rives, J. Am. Chem. Soc., 1996, 118, 11225–11236.
31. G. A. Kaminski, R. A. Friesner, J. Tirado-Rives and W. L. Jorgensen, J. Phys. Chem. B, 2001, 105, 6474–6487.
32. W. L. Jorgensen and J. Tirado-Rives, Proc. Natl. Acad. Sci. U. S. A., 2005, 102, 6665–6670.
33. L. S. Dodda, I. C. de Vaca, J. Tirado-Rives and W. L. Jorgensen, Nucleic Acids Res., 2017, 45, W331–W336.
34. L. S. Dodda, J. Z. Vilseck, J. Tirado-Rives and W. L. Jorgensen, J. Phys. Chem. B, 2017, 121, 3864–3870.
35. J. Garrec, A. Monari, X. Assfeld, L. M. Mir and M. Tarek, J. Phys. Chem. Lett., 2014, 5, 1653–1658.
36. M. Barbatti, Newton-X, 2.2, http://newtonx.org.
37. M. Barbatti, M. Ruckenbauer, F. Plasser, J. Pittner, G. Granucci, M. Persico and H. Lischka, Wiley Interdiscip. Rev.: Comput. Mol. Sci., 2014, 4, 26–33.
38. A. A. Granovsky, Firefly, 8.2.0, http://classic.chem.msu.su/gran/firefly/index.html.
39. MNDO99, version mndo99_20170427_intel64_ifort-13.1.3.192_mkl-11.1.4.211. Available by request: https://www.kofo.mpg.de/en/research/molecular-theory-and-spectroscopy.
40. S. J. Plimpton, J. Comput. Phys., 1995, 117, 1–19.

---

Footnotes

† Electronic supplementary information (ESI) available: Calculations of benchmark potential energy surfaces using ab initio, semi-empirical ODMx methods, and empirical force fields. Movies of example collision trajectories. Computational time for different MD methods. See DOI: 10.1039/d1cp04720e

‡ These authors contributed equally to this work

---