Reactive Atomistic Simulations of Diels-Alder Reactions: the Importance of Molecular Rotations

TL;DR

This study uses reactive atomistic simulations with neural networks to analyze Diels-Alder reactions, revealing that molecular rotations significantly influence reaction pathways and are enhanced by bromine substituents.

Contribution

The paper introduces multisurface adiabatic reactive molecular dynamics combined with PhysNet neural networks to investigate Diels-Alder reactions, highlighting the role of molecular rotations.

Findings

Rotational energy accounts for about 65% of the energy driving the reaction.

Bromine substituents increase the importance of rotational excitation.

Reactions at high energies are predominantly direct and synchronous.

Abstract

The Diels-Alder reaction between 2,3-dibromo-1,3-butadiene and maleic anhydride has been studied by means of multisurface adiabatic reactive molecular dynamics and the PhysNet neural network architecture. This system is used as a prototype to explore the concertedness, synchronicity and possible ways of promotion of Diels-Alder reactions. Analysis of the minimum dynamic path indicates that rotational energy is crucial ($\sim$ 65 %) to drive the system towards the transition state in addition to collision energy ($\sim$ 20 %). Comparison with the reaction of butadiene and maleic anhydride shows that the presence of bromine substituents in the diene accentuates the importance of rotational excitation to promote the reaction. At the high total energies at which reactive events are recorded, the reaction is found to be direct and mostly synchronous.

Equations

$$\sigma =2\pi b_{\max }\frac{1}{N_{\mathrm{tot}}}\sum _{i=1}^{N_{\mathrm{reac}}}{b}_i$$

$$k=\sigma \cdot v_{\mathrm{rel}}$$

$$V_{\mathrm{MS}-\mathrm{ARMD}}=\sum _{i=1}^n{w}_i(x)V_i(x)+\sum _{i=1}^{n-1}\sum _{j=i+1}^n[w_i(x)+w_j(x)]\sum _{k=1}^{n_{ij}}\Delta V_{\mathrm{GAPO},k}^{ij}(x)$$

$$V_{\mathrm{PhysNet}}=\sum _{i=1}^N{E}_i+k_e\sum _{i=1}^N\sum _{j>i}^N\frac{q_i{q}_j}{r_{ij}}$$

Full text

I Introduction

The regio- and stereoselective Diels-Alder reaction in which a diene reacts with a dienophile to form a cyclic product is widely used in synthetic organic chemistry.Diels and Alder (1928); Ishihara and Sakakura (2014) In this reaction, two $\rm\sigma$ bonds and one $\rm\pi$ bond are formed from three $\rm\pi$ bonds as depicted in Scheme 1.

There have been many experimental and theoretical studies aimed at unveiling the Diels-Alder mechanism at a molecular level and its dependence on the geometric and electronic characteristics of the reactants.Houk, González, and Li (1995); Yepes et al. (2013); de Souza et al. (2016); Domingo (2014); Horn, Herek, and Zewail (1996); Diau, De Feyter, and Zewail (1999); Saettel et al. (2002); Singleton et al. (2001); Goldstein, Beno, and Houk (1996); Sakai (2000); Domingo and Saez (2009); Donoghue and Wiest (2006) Since two bonds are formed, questions about concertedness and synchronicity render this reaction also interesting from a theoretical point of view.

The concertedness of a reaction is determined by the topology of the potential energy surface (PES).Minkin (1999) A concerted mechanism has only one transition state (TS) between reactants and products such that the reaction takes place in a single step. In a stepwise mechanism, there is at least one stable intermediate between reactants and products.

The time elapsed between the formation of the first and the second bond defines the synchronicity.Minkin (1999) It is usually thought that symmetric TSs give rise to synchronous processes in which both bonds are formed at the same time while asymmetric TSs lead to asynchronous processes in which first one bond is formed and then the second one follows. This definition of synchronicity has recently been challenged by some authors who argue that it should not be defined from a geometric but from a dynamic point of view since the connection between spatial quantities and temporal concepts may not always hold.de Souza et al. (2016); Black et al. (2012)

There is a long-standing discussion about synchronicity and concertedness in the context of Diels-Alder reactions.Houk, González, and Li (1995) While it is often believed to be a concerted, synchronous reaction involving an aromatic TS governed by the Woodward-Hoffmann rules,Hoffmann and Woodward (1968) experiments and calculations show that this is not true for all cases.Horn, Herek, and Zewail (1996); Diau, De Feyter, and Zewail (1999); Saettel et al. (2002); Singleton et al. (2001); Goldstein, Beno, and Houk (1996); Wilsey, Houk, and Zewail (1999); Hofmann and Schaefer III (1999); Bouchoux, Salpin, and Yáñez (2004) In principle, one can think of three possible mechanisms (see Scheme 1): (a) concerted, (b) stepwise with a short-lived intermediate, and (c) stepwise with a long-lived intermediate whose lifetime allows for the rotation around a C-C bond. Note that when the system cannot rotate around a C-C bond, as in cases (a) and (b), the reaction is stereo-selective and only the s-cis conformer of the diene will yield the cyclic Diels-Alder product. Mechanism (c), on the other hand, is not stereo-selective and the s-trans conformer of the diene could also participate in the reaction.

Computational studies of Diels-Alder reactions are especially sensitive to the choice of method and basis set. Gayatri (2011); Linder and Brinck (2013) Different electronic structure methods give different TS structures and activation energies as was summarized in the introduction of our previous paper (Ref. Rivero, Meuwly, and Willitsch, 2017). This leads to different conclusions regarding the synchronicity and concertedness of the reaction. In general, neutral reactions occur in a concerted fashion while in cationic systems a stepwise mechanism is favored. Saettel et al. (2002); Singleton et al. (2001); Goldstein, Beno, and Houk (1996) However, the border between these two mechanisms is diffuse. Moreover, some studies show that both concerted and stepwise mechanisms are present at the temperatures at which these reactions are usually performed due to the energetic proximity of both pathways in many systems making dynamical studies crucial for the exploration of these reactions. To the best of our knowledge, all theoretical studies on the dynamics of Diels-Alder reactions have started from TS-like structuresde Souza et al. (2016); Black et al. (2012); Tan, Hirvonen, and Paton (2018); Wang, Hirschi, and Singleton (2009); Liu et al. (2016) or have used steered dynamics to drive the reaction.Soto-Delgado, Tapia, and Torras (2016) These procedures will most likely bias the final result and do not allow the direct calculation of reaction rates. Gas-phase reactive molecular dynamics simulations starting from an equilibrated ensemble of a statistically significant number of initial conditions, on the other hand, have recently been shown to provide molecular-level details into reactions relevant to atmospheric chemistry Meuwly (2019); Brickel and Meuwly (2017) and reactions in the hypersonic regime.Koner, Bemish, and Meuwly (2018); Denis-Alpizar, Bemish, and Meuwly (2017)

From an experimental perspective, the most precise data on reaction mechanisms and dynamics can be gained from gas-phase studies performed under single-collision conditions. As the progress in molecular-beam experiments allows the probing of ever-larger systems under precisely defined conditions,Chang et al. (2013); Willitsch (2017) the open questions pertaining to the mechanistic details of Diels-Alder reactions become an attractive target of study. Since only scarce information is available on these important mechanistic aspects,Eberlin and Cooks (1993) we present here a molecular dynamics study of the Diels-Alder reaction between 2,3-dibromo-1,3-butadiene (DBB) and maleic anhydride (MA) (see Scheme 1) which may serve as a guide for future experiments.

We have chosen DBB because it is a generic diene which fulfills the experimental requirements for conformational separation of its isomers by electrostatic deflection of a molecular beam,Chang et al. (2013); Willitsch (2017) thus enabling the characterization of conformational aspects and specificities of the reaction. MA is a widely used, activated dienophile which due to its symmetry simplifies the possible products of the reaction. The reaction of DBB and MA thus serves as a prototypical system well suited for the exploration of general mechanistic aspects of Diels-Alder processes in the gas phase.

In this paper, we use reactive molecular dynamics simulations starting from the two reactant molecules approaching each other in order to simulate a collision experiment and address questions such as whether the reaction is synchronous, whether the mechanism is complex-mediated and how the reaction could be promoted.

II Methods

II.1 Molecular Dynamics Simulations

Atomistic simulations were carried out either with the CHARMM programBrooks et al. (2009) using an initial parametrization from SwissParamZoete et al. (2011) or with the Atomic Simulation Environment (ASE)Larsen et al. (2017) using the PhysNet neural network architecture designed for predicting energies, forces and dipole moments of chemical systems.Unke and Meuwly (2019) All bonds involving hydrogen atoms were flexible and the time step used in the simulations was $\Delta t$ = 0.1 fs to ensure conservation of total energy. The velocity Verlet algorithm was used for the propagation of the equations of motion.Swope et al. (1982)

The initial SwissParam parametrization was modified in order to construct a multisurface adiabatic reactive molecular dynamics (MS-ARMD) Nagy, Reyes, and Meuwly (2014) force field for the present Diels-Alder reaction. Ensembles of structures for the parametrization of the MS-ARMD model were generated with CHARMM as follows: the optimization of the system with the adopted Newton-Raphson method was followed by 50 ps of heating dynamics, 50 ps of equilibration at 500 K, 60 ps of cooling down to 300 K and free NVE (microcanonical ensemble) dynamics. The temperature was only raised up to 400 K for parts of the parametrization of the reactant van der Waals complex to avoid dissociation.

For the PhysNet parametrization, initial ensembles for the different fragmentsHuang and von Lilienfeld (2017) were generated using the PM7 level of theory implemented in MOPAC2016Stewart (2013, 2016) and subsequently augmented using adaptive sampling.

In order to generate the initial conditions for the collision simulations, ensembles of the individual molecules (MA and DBB) at different vibrational temperatures were generated using CHARMM as described above. Heating and equilibration temperatures were modified depending on the desired final vibrational temperature ($T_{\rm vib}$). The reactants were placed at an initial distance of 20 $\rm\AA$ (taking the center of mass of each molecule as reference) with a random relative orientation. In order to tune the collision energy ($E_{\rm coll}$), the atomic velocities along the collision axis were modified. The impact parameter ($b$) was uniformly sampled by displacing one of the molecules along a perpendicular axis. Rotational temperature ($T_{\rm rot}$) was added following calculation of the moment of inertia tensor and assuming equipartition among the three rotational degrees of freedom. Excitation of specific vibrational normal modes was achieved by projecting the initial velocities onto the space of normal modes and modifying the kinetic energy of the desired normal mode.

The trajectories were considered reactive when the C-C distance between the carbon atoms involved in the two new bonds was smaller than 1.6 $\rm\AA$. The reactive cross section $\sigma$ was calculated as

[TABLE]

where $b_{\rm max}$ is the maximum impact parameter at which no reactive events are observed any more, $N_{\rm tot}$ the total number of trajectories, $N_{\rm reac}$ the number of reactive trajectories and $b_{i}$ the impact parameter of the reactive trajectory $i$. The reaction rate $k$ can be calculated from

[TABLE]

where $v_{\rm rel}$ is the relative center of mass velocity of the colliding molecules.

The total kinetic energy of the minimum-dynamic-path trajectories was decomposed either into the translational energy of the center of mass of the reactant molecules, the rotational energy corresponding to their angular momentum and vibrational energy as stated in equations S1 - S6 of the supplementary information or by projecting it onto the degrees of freedom of the system. The reference degrees of freedom were calculated as the eigenvectors of the Hessian matrix of the isolated, reactant molecules with geometries corresponding to the last point of each trajectory.

II.2 Parametrization of MS-ARMD

The force fields for the reactant and product states were parametrized to reproduce DFT reference energies at the M06-2X/6-31G* level of theory which was found to yield the best accuracy at the DFT level for this type of reaction by Linder and BrinckLinder and Brinck (2013); Zhao and Truhlar (2008) and following our previous work.Rivero, Meuwly, and Willitsch (2017) All single point calculations were performed using Gaussian09.Frisch et al. (2009) Starting from parameters retrieved from SwissParam,Zoete et al. (2011) first ensembles for reactants and products were generated as described in section II.1. Reference energies were calculated and the individual force fields were fitted using a downhill simplex algorithm.Nelder and Mead (1965) The standard CHARMM harmonic potentials for the description of C-C, C=C, C-O and C-Br bonds were replaced by Morse potentials (C-H and C=O bonds were kept as harmonic). Furthermore, in the reactant force field, the Lennard-Jones potentials between the four carbon atoms involved in the Diels-Alder reaction as well as between Br/O atoms were replaced by a generalized Lennard-Jones potential.Nagy, Reyes, and Meuwly (2014) The remaining terms were parametrized as in the standard CHARMM force field.

With the initial parametrization of the individual force fields, another ensemble of structures was generated for reactants and products with the new parameters and added to the training structures. The reference energies were calculated and the parametrization was further refined. This iterative procedure continued until the root-mean-square deviation (RMSD) of the newly generated ensemble was approximately the same as that of the previous iteration of the parametrization. For the reactant force field, 401 and 2064 structures where required for the parametrization of MA and DBB, respectively. The non-bonded terms of the van der Waals complex were parametrized with 2783 additional structures. For the parametrization of the product force field, 2589 structures were required.

The crossing region between the two force fields was smoothed following the internal reaction coordinate (IRC) by combining the final force fields for reactants and products with GAussian times POlynomial functions (GAPOs).Nagy, Reyes, and Meuwly (2014) A genetic algorithm was used for the fitting of the GAPOs. The global PES is given by Eq. 3 where $V_{i}(x)$ are the individual force fields, $w_{i}(x)$ their weights and $\Delta V^{ij}_{{\rm GAPO},k}(x)$ the GAPOs.Nagy, Reyes, and Meuwly (2014) The product of the Diels-Alder reaction has two possible connectivities (see Fig. S3 of the supplementary information). In order to make the parametrization of the product permutation invariant, two force fields that describe these two possible connectivities were used.

[TABLE]

II.3 Parametrization of PhysNet

Reference data for training PhysNet (energies, forces and dipole moments) were calculated at the DFT M06-2X/6-31G* level of theory using Gaussian09.Frisch et al. (2009) All possible “amon” fragmentsHuang and von Lilienfeld (2017) for DBB, MA, and their reaction product were generated (378 in total) and different geometries for all fragments were sampled by running Langevin dynamics at 1000 K at the PM7 level of theory. After training PhysNet models on this initial dataset, additional structures were generated by adaptive sampling:Behler (2014, 2015) an ensemble of 4 PhysNet models was used to run Langevin dynamics at 1000 K and new ab initio data was calculated for geometries for which the energy predictions between the different models differed by more than 0.5 kcal/mol. For further details on the adaptive sampling method, see Ref. Unke and Meuwly, 2019. The dataset was iteratively augmented in this fashion until no significant deviations between the predictions of individual PhysNet models could be observed (the final dataset contains 224483 structures which is approximately 50 times larger than it was required for the MS-ARMD model). PhysNet was then trained on 200000 structures of the final dataset (with 20000 additional structures used for validation) by minimizing the mean absolute error between neural network predictions and the reference data using the procedure given in Ref. Unke and Meuwly, 2019 (all hyperparameters of the neural network architecture and the training procedure were set to the values recommended in Ref. Unke and Meuwly, 2019). The global PES is given by

[TABLE]

where $N$ is the total number of atoms, $k_{e}$ is the Coulomb constant, $r_{ij}$ is the distance between atoms $i$ and $j$, and $E_{i}$ and $q_{i}$ are atomic energy contributions and partial charges (corrected to guarantee charge conservation)Unke and Meuwly (2019) predicted by PhysNet. Here, the Coulomb potential is damped at short distances to avoid numerical problems (see Ref. Unke and Meuwly, 2019). The PhysNet architecture guarantees that Eq. 4 is invariant with respect to translations, rotations and permutation of atoms sharing the same element type.Unke and Meuwly (2019)

III Results and Discussion

III.1 Parametrization of the Reactive Force Fields

Fig. 1 shows stationary points on the PES for the Diels-Alder reaction between DBB and MA at the M06-2X/6-31G* level of theory. Note that DBB is in a gauche conformation (the C=C$-$C=C dihedral angle has a value of 50°) since the s-cis geometry is not a minimum on the PES. Due to the fact that both reactant molecules are symmetric, there are only two possible paths for the Diels-Alder reaction. These are referred to as “endo” and “exo” depending on the relative orientation of the reactants. The endo product (P-endo) is taken as the zero of the energy scale in Fig. 1. Normally, the intermolecular interactions favor the endo path which M06-2X correctly describes.Alder and Stein (1936); Fernández and Bickelhaupt (2014); Rivero, Meuwly, and Willitsch (2017) The dissociation of the van der Waals complexes (R-endo and R-exo) towards the reactants is more favorable than the reaction over the barrier towards the Diels-Alder products (see Fig. 1). Judging from the PES and the geometries of the TSs, the reaction should be concerted and symmetric.

As discussed in the introduction, in order to answer questions such as whether the reaction is synchronous or complex mediated, the dynamics of the system must be studied. Therefore, we constructed a computationally efficient PES that allows us to run a statistically significant number of trajectories, such that diverse initial conditions can be sampled. The quality of the parametrized MS-ARMD model is shown in Fig. 2(a). Points generated during the parametrization of the product force field (2589 structures) and the parametrization of the non-bonded interactions of the reactant force field (2783) are shown, as well as the IRCs for the endo (81) and the exo (59) paths. The exo IRC was used for the parametrization of the GAPOs and thus it is better described by the model than the endo path. The energy of the endo IRC is slightly overestimated by the MS-ARMD model. The total RMSD is 1.47 kcal/mol over a range of 80 kcal/mol which we asses to be sufficient for an adequate characterization of the dynamics of the system. It is important to note that the s-trans conformer has not been included and is not stable in this model. Thus, the reactants will not sample the entire conformational space leading to a possible overestimation of the reaction rates since the s-trans conformer does not contribute to the reaction i.e. the number or unreactive trajectories increases when including the s-trans conformer which leads to a reduction of the reaction flux and hence the rate. For the validation of the MS-ARMD model the reactant, transition state and product structures were minimized and compared to the DFT geometries (see Fig. S1 of the supplementary information). The harmonic frequencies computed by DFT and the model were also compared (Fig. S2 of the supplementary information). The parametrized model is given in Tables SI-SVII of the supplementary information.

For validation and direct comparison, a PhysNet model was parametrized for the same system. Fig. 2(b) reports on the quality of this PES by showing the correlation between the reference energies and the PhysNet predictions on the training structures of the MS-ARMD model with a total RMSD of 0.25 kcal/mol. The PhysNet model is significantly more accurate than MS-ARMD, but it has a computational cost 200 times higher.

III.2 Minimum Dynamic Path

A trajectory starting at the TS geometry without kinetic energy (the total energy of the trajectory is equal to the potential energy of the TS) follows the minimum dynamic path (MDP).Unke, Brickel, and Meuwly (2019) The projection of the total kinetic energy along the MDP towards the reactants onto the degrees of freedom of DBB and MA is shown in Fig. 3(a) as sums of translational, rotational and vibrational energies. At $t$ = 0 fs, the system is at the TS and at $t$ = 100 fs it has arrived at the reactant state. At the beginning of the trajectory, the system contains no kinetic energy and moves only slowly, until at around $t$ = 50 fs more potential energy is converted into kinetic energy. By projecting the total kinetic energy onto the different degrees of freedom of DBB and MA, the active degrees of

freedom can be identified. Those degrees of freedom to which most kinetic energy is imparted will be important for driving the system towards the transition state. Fig. 3(a) shows that rotations contain the largest amount of kinetic energy for both DBB and MA. The rotational energy of DBB and MA together accounts for 63 % of the total kinetic energy while translational energy accounts for 18 %. Certain vibrations are also activated (see Fig. S4 of the supplementary information for the individual contributions). An identical result has been obtained from the direct decomposition of the total kinetic energy (Fig. S5 of the supplementary information). Using the sudden vector projection method Guo and Jiang (2014) we also arrive at the conclusion that rotations are the most important degrees of freedom. This is surprising since one may naively assume that the reaction coordinate is mainly a translation and not a rotation.

In order to further validate this result, the MDP was also calculated using the PhysNet PES. The results are shown in Fig. 3(b) and Fig. S6 of the supporting information and qualitatively agree with those of MS-ARMD (68 % of the total kinetic energy is imparted into rotations and 20 % into translations), although less kinetic energy is acquired by the vibrations with PhysNet: vibrational energy accounts for 12 % of the total kinetic energy while for MS-ARMD it is 19 %.

The fact that the MDP between TS and reactant for MS-ARMD extends over 100 fs while for the PhysNet PES it spans for 175 fs is due to differences in the shape of the two PESs near the TS (see Fig. S7 of the supplementary information) which define the initial accelerations of the particles.

In order to investigate whether rotations are important only due to the presence of the heavy bromine atoms in the system (the large mass of bromine atoms potentially enhances the torque on DBB along the MDP), the MDP of the reaction of MA with butadiene (BD) was calculated with the PhysNet model. This is possible since the dataset on which the model was trained contains the necessary information to describe the system with hydrogen atoms instead of bromine atoms. The results are shown in Fig. 3(b) (individual contributions are shown in Fig. S8 of the supplementary information). It was found that the important degrees of freedom for MA remain approximately the same for the reaction with DBB and the reaction with BD. When comparing the distribution of kinetic energy of BD with that of DBB, it can be seen that, even though rotations are still active, translations seem to be imparted the most kinetic energy. The rotational energy of BD and MA together accounts in this case for 48 % of the total kinetic energy while translational energy accounts for 40 %. This indicates that the heavy bromine atoms accentuate the importance of rotational excitation for driving the reaction.

III.3 Cross section for the formation of van der Waals complexes in the entrance channel

The formation of van der Waals complexes in the entrance channel was studied in order to determine whether the reaction is direct (without the formation of complexes) or complex-mediated. Initial ensembles were generated as described in section II and the impact parameter $b$ was uniformly sampled in intervals of 1 $\rm\AA$ until the maximum impact parameter was reached. For each set of initial conditions ($E_{\rm coll}$, $T_{\rm vib}$, $T_{\rm rot}$, $b$) 500 trajectories were run for 10 ps. If at the end of a trajectory, the center of mass distance between the two molecules was $<$ 15 $\rm\AA$, it was considered that a van der Waals complex had been formed. Fig. 4 shows the cross section for the formation of complexes as a function of the collision energy. It can be seen that the cross section diminishes as the collision energy increases. In fact, the cross section is close to zero for $E_{\rm coll}$ $>$ 20 kcal/mol.

The influence of vibrational and rotational temperature is also shown in Fig. 4. Increasing the energy of the system in either vibrational or rotational degrees of freedom reduces the cross section for complex formation. These findings can be explained by the stabilization of the well of the van der Waals complex of around 12 kcal/mol (Fig. 1). Thus when the system has collision energies above 15 kcal/mol or sufficiently high rotational and vibrational excitation, it can dissociate or not get trapped in the potential well at all.

III.4 Reactive collisions

For studying the Diels-Alder reaction itself, approximately $10^{7}$ reactive molecular dynamics simulations were carried out. The sets of simulations that yield reactive events are summarized in Table SVIII of the supplementary information and Fig. 5(a) displays trajectories of two reactive events. The vibrational temperature was 100 K for most of the trajectories although some tests at higher vibrational temperatures have also been performed. Collision energies between 15 and 100 kcal/mol were sampled. The impact parameter was varied from 0 to 6 $\rm\AA$. The rotational temperature was 0, 1000, 2000 and 4000 K.

In order to calculate a reactive cross section, we ran $5\cdot 10^{6}$ trajectories for initial conditions $E_{\rm coll}$ = 100 kcal/mol, $T_{\rm rot}$ = 4000 K and $T_{\rm vib}$ = 100 K at which the largest number of reactive events were observed. The maximum impact parameter $b_{\rm max}$ was 5 $\rm\AA$ which yields a cross section $\sigma=2.13\cdot 10^{-3}~{}\rm\AA^{2}$, corresponding to a rate $k=7.53\cdot 10^{-14}$ cm3 molecule*-1* s*-1*. A collision energy of 100 kcal/mol corresponds to a relative velocity of 3533 m/s which is at the very upper limit of what may be experimentally achievable.Rivero, Meuwly, and Willitsch (2017)

As the MDP shows and Fig. 5(b) confirms, the rotational energy promotes the reaction, although reactive events are still rare (1 in 104). In order to further test this result, a rotational temperature of $T_{\rm rot}$ = 8000 K was used for collisions at $E_{\rm coll}$ = 15, 20 kcal/mol but few reactive events (1 - 2) were observed, indicating that collision energy is also needed for the reaction to take place. It should be noted that only 5 reactive events out of 482 are recorded without rotational excitation (see Table SVIII of the supplementary information).

With a vibrational temperature of 100 K, no reactive events could be observed with $E_{\rm coll}$ $<$ 75 kcal/mol even though, in principle, with a collision energy of 10 kcal/mol the system would have enough energy to overcome the reaction barrier (Fig. 1). This might indicate that the reaction rate with such low $E_{\rm coll}$ is too small combined with the fact that at these energies the reaction could be partly complex mediated and take much longer times than simulated. For a reactive event to take place, the molecules need to collide in a suitable relative orientation in order to overcome steric constraints and with the right distribution of energy, such that the TS can be reached. From Fig. 4, we know that there is no complex formation for $E_{\rm coll}$ $>$ 20 kcal/mol and thus the reactions at the energies in Fig. 5 must be direct.

The influence of vibrational energy was tested by raising the vibrational temperature to $T_{\rm vib}$ = 1000 K and by individualy exciting some of the vibrational normal modes that appear to be active in the MDP. However, reactive events were so rare (0 - 2) that it can be concluded that vibrational activation of the reaction is weak, at least in this energy range.

The time $\delta t$ elapsed between formation of the first and second bond was calculated for all the reactive events in order to determine the synchronicity of the reaction (Fig. S9 of the supplementary information). Out of 482 reactive events (see Table SVIII of the supplementary information), only two are slightly asynchronous with time lapses of 40.4 fs and 43.9 fs which are larger than the ca. 30 fs corresponding to a C-C bond vibrational period in cyclohexene. de Souza et al. (2016); Diau, De Feyter, and Zewail (1999) The difference between a synchronous and a slightly asynchronous process can be seen in Fig. 5(a). We have not found a correlation between the (a)symmetry of the TS structure in the dynamics and the (a)synchronicity of the process (Fig. S9 of the supplementary information) meaning that the elapsed time does not seem to depend on the (a)symmetry of the TS structure. Here, the TS geometry along the reactive trajectory is defined as the structure at the maximum in potential energy before the systems crosses from the reactant surface to the product surface. When recrossings occur, the last crossing point is taken for the definition of the TS structure.

IV Conclusions

We have studied the Diels-Alder reaction between MA and DBB using reactive molecular dynamics. The trajectories start with the two reactant molecules approaching each other as in a collision experiment. The minimum dynamic path indicates that rotations are important to drive the system towards the transition state. Furthermore, this finding has been confirmed by the fact that the majority of reactive collisions occur with rotational excitation. The presence of bromine substituents in the system accentuates the importance of rotations, but they were also found to be important for reactions of non substituted dienes. At the energies at which reactive events were observed, the cross section for the formation of van der Waals complexes in the entrance channel is almost zero and thus the reaction cannot be complex-mediated. Most of the observed reactive events are synchronous.

One of the fundamental aspects in reaction dynamics concerns the question which form of energy (translation, vibration or rotation) is most efficacious for the system to reach and surmount the transition state.Polanyi (1972) Back in the 1970s, when studying model atom plus diatom reactions, Polanyi formulated rules which relate the nature of the transition state (early or late) with the type of energy that promotes the reaction (translational or vibrational energy). ApplicationZhang et al. (2012) and generalizationGuo and Jiang (2014); Jiang, Li, and Guo (2014) of these rules to polyatomic molecules remains a challenging undertaking, see for example the sudden vector projection (SVP) model.Guo and Jiang (2014); Jiang, Li, and Guo (2014) Analysis of a number of atom plus diatom (H+H2, F+H2, F+HCl) or atom plus triatom (H+H2O, F+H2O) reactions highlighted the cases under which rotations may play an important role in promoting or inhibiting a reaction.Jiang, Li, and Guo (2014) The strength of the SVP model is that it requires only information about the normal modes and their directions in the reactant and transition state structures. On the other hand, the “sudden approximation” will not be applicable to situations in which internal vibrational relaxation (IVR) occurs or for collision energies much higher than the reaction threshold, as is the case in the present work. The present work suggests that rotations can play an important role for reactions involving large excess of translational energy and the implications for reaction dynamics involving polyatomic molecules are exciting.Li and Guo (2014a); Song et al. (2016); Kilaj et al. (2018); Jiang, Li, and Guo (2014); Li and Guo (2014b)

As a next step, we intend to perform reactive molecular dynamics simulations for the cationic reaction between DBB and MA that is expected to be faster and in which concerted and stepwise mechanisms are anticipated to coexist.Rivero, Meuwly, and Willitsch (2017)

To the best of our knowledge, no simulation study had been performed in Diels-Alder reactions starting from the beginning of the reaction without steered dynamics. The present results indicate that the need of high collision energies together with rotational excitation and the low reaction rate of reaction will make the study of this reaction in single-collision experiments challenging.

Acknowledgements

We acknowledge funding from the Swiss National Science Foundation, grant nr. BSCGI0_157874, and the University of Basel.