A novel partitioning scheme for the application of the distortion/interaction - activation strain model to intramolecular reactions

The distortion/interaction or activation strain model, developed by Houk and Bickelhaupt, relates chemical reactivity to the reagents deformations and reciprocal electronic influences. However, in its original formulation, it struggles to elucidate the mechanistic insights of intramolecular reactions, those unimolecular processes in which two parts of a molecule, the reaction centers, linked by a connector, are brought together to yield a different chemical species. Here we present a modification of the distortion/interaction procedure for its application on intramolecular reactions. This new procedure allows the calculation of the influence exerted by the connector over the reaction pathway in an indirect way, from the distortions of the two reaction centers and their interaction energy. This procedure does not include additional, undesired interactions and offers the possibility of calculating very large connectors in a computationally inexpensive way. We applied this methodology in the normal electron-demand Diels–Alder reaction of 1,3,8-nonatriene derivatives, with different functionalizations and connector lengths. In-depth analysis of the IRC showed that the reaction pathway can be subdivided in three main regions, what we called the oncoming, conversion and relaxation phases, each of them characterized by different evolutions of the distortion and interaction energies, and with clear geometry changes. We suggest that this new formulation can provide additional information for intramolecular reactions, especially to those processes for which the connector is said to play a crucial role in the observed reaction rates.


Introduction
A thorough comprehension of any given chemical reaction requires detailed knowledge about the energy changes associated with the sequential transformations that take place. While we are nowadays proficient in the study of kinetics and thermodynamics, the origins of these properties often lie in the dark. In this regard, different theories and models have been developed throughout the last decades to provide straightforward approaches to the reaction energetics based on simple chemistry-based concepts. In 1952, Fukui and coworkers presented Frontier Molecular Orbital (FMO) theory, that relates reactivity to the interaction of the highest occupied and lowest unoccupied molecular orbitals of the reagents [1]. Four years later, in 1956, Marcus theory (MT) emerged as an explanation for the reaction pathways of electron transfer [2]. In 1971, Morokuma introduced energy decomposition analysis (EDA), that divides the total interaction energy between two molecular entities into different, chemically meaningful terms [3]. About 30 years later, the groups of Bickelhaupt and Houk proposed in parallel an innovative approach to reactivity, called activation strain model (ASM) by the former and distortion/interaction model (DIM) by the latter. This model relies on the assumption that the reaction energy arises from two contributions. Using the DIM terminology (see Fig. 1), the distortion energy or ΔE d , the energy (usually loss) associated with the geometry deformations experienced by the reagents along the reaction coordinate; and the interaction energy or ΔE i , the energy (usually gain) that comes from the changes in the electronic structure to yield the final product [4]. Another set of definitions is that ΔE d represents the energy needed to give the reagents the geometry they show at the TS while keeping them at infinite distance (Fig. 1, blue and green arrows), while ΔE i is the energy released when the distorted reagents form the TS (Fig. 1, yellow arrow). This model also assumes that ΔE d is additive (the total distortion is equal to the sum of the distortions of each reagent), since the individual distortions are intrinsic to each corresponding reagent and independent on the other(s).
Since its proposal, the ASM/DIM has been applied to many different reactions, even to those that involve transition metal catalysts [5][6][7]. Nevertheless, if there is one unique type of organic reactions for which this model is particularly successful, it is for pericyclic reactions. In those cases, the ΔE d and ΔE i terms can be readily related to the bond reorganizations and orbital interactions, respectively. Among all pericyclic reactions, the Diels-Alder [8][9][10][11][12][13][14][15] and the 1,3-dipolar cycloaddition [16][17][18][19][20] have been thoroughly studied, with an endless variety of alkenes and dienes/ dipoles. Other pericyclic reactions such as the Alder-ene [21] and the [2 + 2] cycloaddition [22] have been also investigated. Another strength of this model is its ability to be easily combined with previous reactivity models, like EDA [23] or FMO [5, 6, 9, 10, 14-17, 20, 22, 24], so a more complete analysis of the reaction profile can be provided.
However, the original formulation of the ASM/DIM has one clear limitation. As it needs all reagents to be separable to calculate ΔE d , it is difficult to apply to intramolecular reactions, in which one only has one reagent that contains different reaction centers, which will evolve together toward the final adduct. Given that intramolecular pericyclic reactions are common [25][26][27], especially those oriented to the synthesis of natural products [28], it becomes of high interest to formulate an alternative to ASM/DIM to describe such processes. In 2011, Houk and coworkers presented a possible solution to this problem [29]. Suppose that, for an intramolecular reaction, the molecule can be separated in three components (see Fig. 2), the two reaction centers, R1 and R2, and the tether or connector, C, that joins them. In their approach (Fig. 2a), they stated that, if one assumes that the distortion energy is still additive, then it is possible to detach each component of the molecule, replace the loose ends with a chemically similar functional group (X) and compute all the distortion energies. While this should be a safe procedure for R1 and R2 as long as X and X′ are chosen carefully, it may induce errors in the calculation of ΔE d (C) as it gives rise to a new X-X′ interaction not included in the original molecule. Alternatively, Bickelhaupt and coworkers described another fragmentation procedure to separate the reaction centers [30][31][32][33]. This second procedure carries out a homolytic bond cleavage in the connector (Fig. 2b), assuming that the reaction centers can be considered as "strongly bound reagents." This way, ΔE d (R1) and ΔE d (R2) can be calculated just carrying out single-point energy calculations considering the independent fragments as openshell doublets (if the original molecule was a closed-shell singlet). However, this second approach cannot disentangle the contribution of ΔE d (C) to the activation energy, as it is included inside ΔE d (R1) and ΔE d (R2) and thus no value can be provided.
In this work, we present a novel approach for the application of the ASM/DIM on intramolecular reactions, which can quantitatively calculate the value of ΔE d (C) as a linear combination of all the remaining energy terms. To do so, we propose an alternative pathway for the calculation of ΔE i exclusively from the two reaction centers (Fig. 2c). In   [29]. b Bickelhaupt and coworkers' approach [30]. c Our approach presented in this work. In cases a and c, the fourth remaining term can be calculated from the difference with respect to the total energy; case b does not calculate ΔE d (C). R1: First reaction center. R2: Second reaction center. C: Connector this case, we took the normal electron-demand Diels-Alder reaction as model, by placing a methoxy group at the diene and different electron-withdrawing groups (EWGs) at the alkene (see Fig. 3). As connectors, alkyl chains of different length were tested.

Computational details
Minimum and transition state (TS) optimizations, frequency analysis, intrinsic reaction coordinate (IRC) calculations and Born-Oppenheimer molecular dynamics (BOMD) simulations have been carried out within the framework of density functional theory (DFT). In particular, the exchange-correlation M06-2X hybrid functional was employed [34], in combination with the cc-pVDZ basis set [35]. All minima were characterized as regions of the potential energy surface with zero gradient and zero imaginary frequencies, whereas TSs showed one single imaginary frequency, that corresponds to the vibration toward the acyclic structure and the Diels-Alder adduct. For some of the procedures applied in this communication, a step-by-step description is provided in the Supporting Information (SI). These calculations were performed using the Gaussian 09 package [36].

Results and discussion
In our first step, provided that we have established the difficulties of directly finding the value of ΔE d (C), we decided to explore the feasibility of calculating the value of ΔE i by other methods and, with ΔE d (R1) and ΔE d (R2) obtained by the same methodology previously described [29], then it should be possible to calculate indirectly ΔE d (C) from the sum of those three terms and the difference with respect to the total energy, ΔE. The most critical point of this approach is to find a way to accurately calculate ΔE i . In this regard, let us assume that this interaction energy can be decomposed into four different terms, ΔE i = ΔE i (R1-R2) + ΔE i (R1-C) + ΔE i (R2-C) + ΔE disp , where the first three terms are associated to the pair-wise interactions and the last one to non-local, dispersion reorganizations. In the event of a chemical reaction between R1 and R2, where C is an inert spectator (in the sense that it is not chemically modified in any way), the first term should be clearly predominant and To corroborate this hypothesis, we took a closer look to our particular systems. In this case, the condition of C being an inert connector weakly interacting with the two reaction centers is clearly fulfilled, as it is as simple as an alkyl chain. In fact, the influence of C over the total interaction energy should be remarkably close to that of a methane molecule (see Fig. 4). Therefore, we decided to explore the intermolecular Diels-Alder reaction between 1-methoxy-1,3-pentadiene and 1-nitro-1-propylene in the presence of methane, given that, in this situation, all distortion and pair-wise interaction energies (see SI for the procedure for the calculation of these latter terms) can be calculated very easily as all reagents are separate entities. We found that the diene-dienophile interaction energy is not only the clear predominant term but recovers up to 99.2% of the total interaction energy (see SI for further information about energy contributions), thus allowing us to infer that ΔE i ≈ ΔE i (R1-R2) is an accurate approximation. Performing this study using the same DFT functional but with the augmented aug-cc-pVDZ basis set (using M06-2X/cc-pVDZ geometries), which should describe better long-range interactions, yielded a similar result (100.0%). Substitution of NO 2 by CF 3 led to a mild overestimation of the total interaction energy (105.0%), due to a larger, positive value of ΔE i (R1-C) (see SI). Before moving on, it is worth mentioning that this approximation has one noticeable flaw. If there exists a strong interaction between C and R1/R2 (or both), then the corresponding pair-wise terms will be no longer negligible, and the approximation fails. Indeed, substituting the methane molecule by 1,3-propanediol, and placing it in a position where it can form hydrogen bonds with both OMe and NO 2 groups (see SI), led to the recovery of 83.3% of the total interaction energy by just considering ΔE i (R1-R2), which slightly increases to 84.7% with the aug-cc-pVDZ basis set.
With a method to obtain all ΔE d contributions and ΔE i , we calculated the reaction profiles for the Diels-Alder reaction that yields the structures shown in Fig. 3. Entries 1-8 vary in the strength of the EWG, entry 9 extends the π system of the dienophile to study the effect of the EWG at long distances, entries 10-11 increase the size of the connector to yield larger cycles, and entry 12 is a transannular reaction, a special type of intramolecular reactions in which the reaction centers are present at opposite places of a macrocycle, and the process leads to the formation of polycyclic scaffolds in a single step [37-39]. Entries 9-12 contain an aldehyde group as EWG as the best compromise between relative strength and size (in terms of number of atoms and electrons), offering a consistent electron-withdrawing effect with little computational effort. Thus, their results should be compared with entry 5 only. For the calculation of all energies, we took as reference the acyclic conformer which is directly connected to the TS of the Diels-Alder reaction, without regard of it being the most stable conformer or not. For all entries, we found that the reaction follows a concerted but asynchronous mechanism: the δ1-β2 bond (see Fig. 3) is formed in first place, quickly followed by the formation of the α1-α2 bond and the cyclohexene motif. This is in perfect agreement with the Lewis model, in which the OMe group can increase the electron density at the δ1 position of the diene while the EWG decreases it at the β2 position of the dienophile. Table 1 shows the activation (ΔE ≠ ) and reaction energies (ΔE reac ) for entries 1-12, and the values of ΔE d for each component of the molecule, the total ΔE d , and ΔE i , calculated at the TS points taking as reference the acyclic structures (reagents). For entries 1-8, the value of σ − of the EWG is also provided. This value is a variation of the Hammett constant, more accurate to describe the effect of the functional group over the reactivity when it is conjugated with a negative charge/partial charge during the reaction pathway [40]. For the calculation of ΔE d (R1), ΔE d (R2) and ΔE i , the connector was replaced by methyl groups (see Fig. 5 and SI), which are nearly identical to the connector (a -CH 2 -group has been replaced by a -CH 3 ); thus, any possible X-X′ interaction that might arise for the calculation of ΔE i (see Fig. 2C) was present in the original molecule. This holds for the main difference with respect to Houk and coworkers' approach [29]. In their procedure to calculate ΔE d (C), they substituted a nitrone, one of the reaction centers, by a primary amine. While uncapable of evolving via 1,3-dipolar cycloaddition, it gives rise to a spurious aminealkene interaction that contaminates the results of both ΔE d (C) and ΔE i . In our case, replacing the alkyl chain that forms the connector by methyl groups does not add new interactions and therefore can lead to more accurate results.
In all cases, the largest contributor to the total ΔE d is R1, and ΔE d (C) takes much smaller values. This is not surprising, given that R1 is larger than R2 and contains more bonds to be reorganized during the reaction, and that the connector only experiences a small compression and/ or light conformational torsions. If we take a closer look at the values of ΔE d (C), its negative value for entry 9 can be explained with the IRC analysis, which is presented later on. Compared to entry 5 (EWG = CHO, formation of a 5-membered ring), entry 11 (formation of a 7-membered ring) shows a larger value of ΔE d (C), which can be related to a larger connector; on contrary, entry 10 (formation of a 6-membered ring) shows a lower value of ΔE d (C) despite having a larger connector than entry 5. This can be explained by the greater stability of the arising cyclohexane with respect to both cyclopentanes and cycloheptanes. Finally, ΔE d (C) for entry 12 (transannular) is 2.4 times greater than that of entry 5, in good agreement with the fact that it contains the same connector twice and therefore a value of ΔE d (C) twice of that with only one connector could be expected. However, given that the "total" ΔE d (C) is calculated, there is no systematic way with this procedure to obtain the individual distortion energies of each connector. For entry 5 (EWG = CHO), we compared the results for ΔE d (C) and ΔE i obtained with our approach and confronted them to the results using Houk's methodology [29]. To do so, we calculated ΔE d (C) directly by substituting R1 and R2 for other functional groups. As inferred from these results, shown in Table 2, this procedure is highly reliant on the chosen replacements for the reaction centers, with ΔE d (C) values that oscillate from 1.05 to 31.40 kcal/mol. Entry 5b represents the closest example to the Houk and coworkers' approach, R2 has been simplified by removal of the formyl moiety and R1 was replaced by a vinyl group, electronically similar to the original structure but unable to undergo Diels-Alder cycloaddition. The strong π-π interaction overestimates ΔE d (C), thus yielding a much larger value for ΔE i . Moreover, such large values of ΔE d (C) are not consistent with the mild deformations experienced by the connector from the initial structure to the TS. While entry 1b (X = X′ = H) produces a similar value of ΔE d (C) that the one obtained with our methodology, we cannot ensure that this is not a fortuitous coincidence and it will be observed for other reactions and/or connectors.
Entries 6b and 7b show the results obtained using Bickelhaupt's approach [30]. Given that this procedure involves a homolytic bond breaking in the connector to separate the reaction centers, the value of ΔE d (C) is not calculated explicitly but it is included in ΔE d (R1) and ΔE d (R2). Therefore, one may expect that these values, obtained with Bickelhaupt's approach, should be slightly higher than the same values obtained with our approach, for which ΔE d (C) is considered separately. However, homolytic cleavage of C a bond (see Fig. 3) leads to lower values of ΔE d (R1) and ΔE i , while ΔE d (R2) gives a nearly identical result. We propose that this arises from the conjugation of the generated CH 2 radical with the diene moiety (revealed by analysis of the  spin density, see SI), that perturbs the electronic structure of the π system compared to the original, non-disconnected molecule. On the other hand, C b bond cleavage, for which the conjugation of the CH 2 radical with the dienophile is less efficient and does not provoke a substantial perturbation of the electron density, leads to the expected trend for the values of the distortion energies. Moreover, all energy terms of entries 5 and 7b are in particularly good agreement. This suggests that, as long as an adequate bond is selected for its homolytic cleavage, both Bickelhaupt's approach and ours yield identical results. While Bickelhaupt's approach is less computationally demanding (requires a lower amount of calculations), our approach allows to quantitatively assign the effect of the connector over the activation barrier. It also provides additional evidence to confirm that the approximation we used to calculate ΔE i is reliable. By taking entries 1-8, we can analyze the effect of the strength of the EWG on the activation, distortion and interaction energies. As can be seen from Fig. 6, there is a linear correlation (R 2 = 0.91) between σ − and ΔE ≠ . Stronger EWGs modify the electron density of the dienophile in a more pronounced way, making it more susceptible to undergo normal electron-demand Diels-Alder reactions, decreasing ΔE ≠ . This relation between the strength of a functional group and its influence in the reaction rates is precisely the usefulness of Hammett constants, and it has been observed in previous studies of this kind [9,22]. There also exists a linear correlation (R 2 = 0.87) between σ − and ΔE i . Following the same reasoning, a stronger EWG should provide the dienophile a larger positive partial charge that, combined with the negative partial charge of the diene (provided by the OMe group), leads to greater values of ΔE i due to simple electrostatic interactions. On the other hand, ΔE d is less affected by the strength of the EWG (R 2 = 0.56). As the main geometry deformations (pyramidalization of the α1, δ1, α2 and β2 positions, elongation of the α1-β1, γ1-δ1 and α2-β2 bonds, and shortening of the β1-γ1 bond) are a constant in all variations of the Diels-Alder reaction, then one should expect a weak dependence of these distortions with respect to the substitution, as we observed.
We conducted an in-depth study by analyzing the evolution of all these energy parameters along the reaction pathway. For this purpose, we calculated the IRC for all reactions and projected them onto the variation of the central CCC angle of the connector (see Fig. 3), which provides a description of how the connector must be distorted to bring the two reaction centers closer in space. Only for entries 10 and 11 was it necessary to project the IRC onto a different parameter, the C-C distance between the two methylene groups directly attached to R1 and R2 (see also Fig. 3). While there are many possible CCC angles in those structures, this parameter also provides a description of the connector compression during the course of the reaction. Figure 7 shows the results for entry 5 (EWG = CHO), all the remaining profiles can be found on the SI (IRC analysis using Houk's and Bickelhaupt's approaches for entry 5 can also be found on the SI). The early stage of the reaction is characterized by a positive ΔE i , even larger than the total ΔE d (Fig. 7, inset on left). This can be related with the increasing steric repulsions between the two reaction centers. On the other hand, the required ΔE d (C) is nearly energetically free (or even favorable), as it takes very low values (Fig. 7, right). We call this first stage (0° to − 5° for entry 5) the oncoming phase of the reaction, no relevant distortions take place in the reaction centers and steric repulsions predominate. After this stage, ΔE i quickly starts to decrease to take negative values and becomes the driving force. While the distortion energies increase at a higher rate than before, it is not enough to overcome ΔE i . This is also true for ΔE d (C), which also experiences a steady increase. We refer to this second stage as the conversion phase of the reaction (− 5° to − 9.5° for entry 5), where all the chemical transformations (bond breakings and formations) take place, and therefore the TS is located here. Finally, the reaction enters a stage in which the growths of ΔE i and all ΔE d are slowed down, and the reaction energy gets just a mild stabilization. This is the relaxation phase of the reaction, in which conformational rearrangements are carried out to reduce steric repulsions (− 9.5° to − 11° for entry 5). These rearrangements are less energetically demanding that any process of bond formation/ breaking, which explains the energy growth reduction.
Interestingly, ΔE d (C) does not grow slowly in the relaxation phase, but in fact it starts to decrease. While this was quite subtle for some systems (like the one presented here, Fig. 7 right), it was more evident in others like entry 2 (EWG = CF 3 , see SI). To provide further insights into the behavior of ΔE d (C) at the conversion and relaxation phases, we performed a 100 fs BOMD calculation, so the  Due to computational cost, this analysis was carried out exclusively for entry 5 (EWG = CHO). Free evolution (no additional excitation energy) from the TS structure toward the Diels-Alder adduct led to the formation of the first C-C bond at 18 fs [41], while the second C-C bond is formed at 37 fs, demonstrating the asynchronous mechanism. At this stage, the HCCH dihedral angle in the alkyl chain that forms the connector has been reduced with respect to the TS (see Fig. 8). In fact, it reaches its minimum value at 72 fs, which explains why ΔE d (C) keeps increasing after surpassing the TS, as shown in Fig. 7. Above 72 fs, the dihedral angle starts to increase and reaches its maximum value, 48.6°, at 93 fs, higher than that of the TS (39.9°). Adopting a more staggered conformation softens steric repulsions, leading to the reduction of ΔE d (C) previously characterized. Therefore, the BOMD simulation is in good agreement with the reaction profile extracted from the IRC. While we did not perform BOMD calculations for the rest of the systems, we corroborated that their relaxation stage, where ΔE d (C) starts to decrease, is associated with an increase of the HCCH dihedral angle to get a nearly staggered conformation.

Conclusions
As a final remark, we have developed in this work a novel way to tackle intramolecular reactions in the framework of the distortion/interaction -activation strain model. The analysis of the variation of the energy components along the IRC allowed us to differentiate three phases in the evolution of the reaction pathway. A similar division of the reaction profile has been proposed before according to variations in the ELF [42], here we reported a complementary analysis based on the changing behaviors of ΔE d and ΔE i throughout the cyclization. Compared to previous applications of ASM/ DIM for intramolecular reactions [29,30], our procedure shows several advantages. First, no interactions absent in the original structure are introduced, and the approximation used for the calculation of ΔE i was tested to be reliable. Second, given that the connector is never calculated on its own, our approach allows the quantitative calculation of ΔE d (C) with little computational cost, with absolutely no regard of how large the connector could be. However, we have found two limitations in our procedure. First, it fails to describe situations where there exists a strong interaction between the connector and the reaction centers. Second, it cannot disentangle the individual participation of the two connectors in a transannular reaction. Nevertheless, we believe that this procedure can be a step forward in the comprehension of the chemical evolution of intramolecular reactions. One of the  Fig. 8 Snapshots of the BOMD calculation of entry 5, from the TS structure toward the Diels-Alder adduct, at different simulation times. The analyzed HCCH dihedral angle is marked in blue. Only the connector atoms are shown in ball-and-wire model for the sake of clarity possible targets we would like to highlight are those intramolecular reactions in which the connector plays a crucial role in the reaction rates. For instance, in those affected by the Thorpe-Ingold or gem-dimethyl effect, the disubstitution at a methylene position of the connector to increase the reaction rate in several orders of magnitude [43,44]. Understanding the role of disubstitution into ΔE d (C) can help to elucidate the origins of this rate enhancement.