Elucidating the Trends in Reactivity of Aza-1,3-Dipolar Cycloadditions

This report describes a density functional theory investigation into the reactivities of a series of aza-1,3-dipoles with ethylene at the BP86/TZ2P level. A benchmark study was carried out using QMflows, a newly developed program for automated workflows of quantum chemical calculations. In total, 24 1,3-dipolar cycloaddition (1,3-DCA) reactions were benchmarked using the highly accurate G3B3 method as a reference. We screened a number of exchange and correlation functionals, including PBE, OLYP, BP86, BLYP, both with and without explicit dispersion corrections, to assess their accuracies and to determine which of these computationally efficient functionals performed the best for calculating the energetics for cycloaddition reactions. The BP86/TZ2P method produced the smallest errors for the activation and reaction enthalpies. Then, to understand the factors controlling the reactivity in these reactions, seven archetypal aza-1,3-dipolar cycloadditions were investigated using the activation strain model and energy decomposition analysis. Our investigations highlight the fact that differences in activation barrier for these 1,3-DCA reactions do not arise from differences in strain energy of the dipole, as previously proposed. Instead, relative reactivities originate from differences in interaction energy. Analysis of the 1,3-dipole–dipolarophile interactions reveals the reactivity trends primarily result from differences in the extent of the primary orbital interactions.


Introduction
In 1938, Smith first suggested the concept of 1,3-dipolar cycloadditions. [1] These elegantly simple reactions, however, stood by the wayside for many years until Huisgen came along and established the generality of this cycloaddition. [2] Traditionally, these reactions required high reaction temperatures and lead to a mixture of regioisomers. [2] The groups of Sharpless [3] and Meldal [4] independently showcased the large acceleration of the reaction of azides with terminal alkynes via copper catalysis. The copper-catalyzed reaction ("click reaction") is regioselective reactions. The BP86/TZ2P method produced the smallest errors for the activation and reaction enthalpies. Then, to understand the factors controlling the reactivity in these reactions, seven archetypal aza-1,3-dipolar cycloadditions were investigated using the activation strain model and energy decomposition analysis. Our investigations highlight the fact that differences in activation barrier for these 1,3-DCA reactions do not arise from differences in strain energy of the dipole, as previously proposed. Instead, relative reactivities originate from differences in interaction energy. Analysis of the 1,3-dipole-dipolarophile interactions reveals the reactivity trends primarily result from differences in the extent of the primary orbital interactions.
for the 1,4-regioisomer and is high-yielding. [5] Over the years, this 1,3-dipolar cycloaddition has undoubtedly become one of the most important reactions and has found utility in a wide range of areas in chemistry, including heterocyclic chemistry, [6] materials chemistry, [7] drug discovery, [8] novel sensor development, [9] radiochemistry, [10] and chemical biology. [11] In addition to azides, various other 1,3-dipoles have found application in synthesis, [12] while the victory march of the copper-catalyzed azide-alkyne cycloaddition has provided additional stimulus to this field. Diazo compounds have been proposed as agents for bioorthogonal cycloadditions, as they show orthogonality to commonly used azide cycloadditions [13] and provide high reactivity. [14] Nitrones have been found to be valuable building blocks for the preparation of potential antiviral agents [15] and agrochemicals, [16] while nitrile oxides have found the application in the preparation of nucleoside and nucleotide analogues [17] and as bioorthogonal agents. [18] From a mechanistic view, the 1,3-dipolar cycloaddition is a pericyclic reaction [19] and is proposed to proceed via a concerted (sometimes asynchronous) mechanism, [6,20] although a competing stepwise mechanism has also been suggested. [21] 1,3-Dipoles are typically described by the zwitterionic resonance structures X = Y + -Z -↔ X --Y + = Z. Recently, 1,3-dipoles have also been suggested to exist, in some cases, as pseudodiradicals, pseudoradicals, or carbenoids. [22] 2π electrons from the dipolarophile and 4π electrons from the dipole flow in a cyclic fashion through a π4 s + π2 s thermal six-electron Hückelaromatic transition state (Scheme 1).
Various qualitative models for understanding the reactivity and selectivity of 1,3-dipolar cycloadditions have been pro-Scheme 1. The concerted mechanism for a generic 1,3-dipolar cycloaddition reaction.
posed. Frontier molecular orbital (FMO) theory has undoubtedly been the most successful of these theories. FMO theory relates molecular reactivity in terms of the energy and shape of the highest occupied molecular orbital (HOMO) and lowest unoccupied molecular orbital (LUMO) of reactants and has been extensively used to provide insight into and predict the reactivity and selectivity of 1,3-dipolar cycloadditions. [20,23] However, FMO analysis has failed to predict the kinetics and selectivity of some cycloaddition reactions, [24] which justifies the need to utilize another, more encompassing predictive model. For that purpose, we turn to the quantitative activation strain model (also known as the distortion/interaction model), [25] which has been used to investigate the reactivity of pericyclic reactions like azide-alkyne cycloadditions and various Diels-Alder-type reactions. [24c,26] Within early versions of this analysis, activation barriers were dissected into strain energy, or deformation energy, the energy needed to distort the reactants from their equilibrium geometry into transition state geometry, and interaction energy between these distorted fragments. The reactivity of 1,3-dipolar cycloadditions has often be described by the required strain energy at the transition state geometry, leading to the concept of "strain controlled" reactions. [23d,24a,27] However, analysis of the strain and interaction energy for a series of reactions at the transition state alone can be misleading if these transition states occur earlier and later on the reaction coordinate. Such single-point analyses may, for example, result in skewed conclusions, since the position of the transition state can influence the magnitude of the interaction and strain terms. Therefore, the analysis should be conducted along the entire reaction coordinate, which allows for an accurate comparison between different systems and does not suffer from the previously mentioned issues. [25,28] An activation strain analysis along the reaction coordinate has provided quantitative insight into the physical factors controlling Diels-Alder reactions, [29] 3 + 2 cycloadditions, [30] as well as 1,3-dipolar cycloadditions. [31] The quantitative analysis is further supplemented by decomposition of the interaction energy using energy decomposition methods provided within the Kohn-Sham density functional theory, [32,33] providing even more insight into the reactivity and selectivity. Such an analysis has already been applied to Diels-Alder cycloadditions, [34] but to the best of our knowledge has yet to be employed to investigate 1,3-dipolar cycloadditions.
the TZ2P basis set in reproducing the G3B3 computed activation and reaction enthalpies [36] for 24 1,3-dipolar cycloaddition reactions between dipoles 1-12 ( Figure 1) and both ethylene (e) and acetylene (a). The Amsterdam Density Functional (ADF) [37] suite was used for all calculations. G3B3 is a composite ab initio method developed by Pople and co-workers based on B3LYP geometries. [38] Due to the high computational cost of the G3B3 method, applications are generally limited to small systems. Grimme's DFT-D3 [39] dispersion corrections were included in our analysis to determine an optimal method that can be used when dispersion interactions are expected to be present. QMflows is a Python package capable of automatically generating input files for various state-of-the-art quantum chemistry platforms, including ADF, ORCA, CP2K, DIRAC, and GAMESS-US. QMflows efficiently distributes these input files in a supercomputer, executes the calculations with the chosen software package and performs an automated post-analysis to retrieve the relevant data. For our purpose, we utilized a Python-based workflow (see SI for details) to benchmark various GGA-DFT functionals in ADF, both with and without explicit dispersion corrections, to assess their accuracies and to determine an affordable method for larger systems. This highly automated workflow explores the potential energy surface (PES) for the reactions between dipoles 1-12 and dipolarophiles, ethylene (e) and acetylene (a). This first involves optimizations of the two reactants (dipole and dipolarophile) and the cycloadduct. Next, the PES is scanned by the symmetric elongation of the newly forming bonds of the cycloadduct, with a series of constrained geometry optimizations. The highest point on this PES is then used as the input for the transition state calculation. Vibrational frequencies are then calculated for all stationary points in order to assess whether they are energy minima or first-order saddle points. Upon completion of these steps, the results are then automatically processed and the requested energy values, in our case the activation (ΔH ‡ ) and reaction (ΔH rxn ) enthalpies are printed. The energetic results are compared to the reference energies obtained from the high-level G3B3 data. [36] The computed G3B3, reproduced from a study by Houk and co-workers, [36] and XC/TZ2P (where XC = BP86, BP86-D3, PBE, PBE-D3, BLYP, BLYP-D3, and OLYP) activation enthalpies (ΔH ‡ ) and reaction enthalpies (ΔH rxn ) for reactions of 1-12 ( Figure 1) with ethylene (e) and acetylene (a) are presented in Table 1 and  Table 2, respectively. All enthalpies are computed with respect to the separate, optimized reactants. The mean deviations (MD), mean absolute deviations (MAD), standard deviations (SD), and maximum errors (both negative and positive) relative to the Table 1. Calculated activation (ΔH ‡ ) and reaction (ΔH rxn ) enthalpies (in kcal mol -1 ) for cycloaddition between 1-12 and ethylene (e). [a] Dipole   [a] All computations based on TZ2P basis set. Table 3. Summary of statistical analysis (in kcal mol -1 ): mean deviations -ΔH MD , mean absolute deviations ΔH MAD , standard deviation -ΔH SD , maximum negative ΔH max(-) and positive ΔH ‡ max(+) error relative to G3B3-computed enthalpies of activation and reaction energies for the 1,3-dipolar cycloaddition reactions of 1-12 with both ethylene (e) and acetylene (a).

BP86
BP86 corresponding G3B3-calculated energies are included in Table 3. The BP86/TZ2P method outperforms all other methods for the calculation of activation enthalpies and has a MD and MAD of only -1.9 and 2.3 kcal mol -1 , respectively. The BP86 reaction enthalpies are, on average, 6.1 kcal mol -1 higher than the G3B3 values. Inclusion of Grimme's D3 dispersion corrections (BP86-D3/TZ2P), yields consistently lower activation enthalpies, with MAD of 5.0 kcal mol -1 . Reaction enthalpies using this method are slightly more accurate, having slightly more favorable MD, MAD, and SD values. PBE/TZ2P and PBE-D3/TZ2P are both less accurate (larger MD and MAD) than BP86/TZ2P in the calculation of activation enthalpies, but are only slightly better for calculating reaction enthalpies. The BLYP/TZ2P and OLYP/TZ2P methods both overestimate barriers and perform very poorly in calculating reaction enthalpies. The BLYP-D3/TZ2P method performs well in calculating activation enthalpies (MAD = 1.9), but reaction enthalpies deviate significantly from the G3B3 values. Weighing the performance and cost, we selected the BP86/ TZ2P level for the following analysis. We suggest using this functional when dispersion energies are expected to be negligible. When non-covalent interactions are suspected to play a role, BLYP-D3/TZ2P would likely yield accurate barriers.

General Methods
All stationary points were calculated using ADF [37] at the BP86 [40] /TZ2P [41] level. A vibrational analysis was used to verify energy minima and transition states. [42] Energy minima had zero imaginary frequencies, while transition states had a single imaginary frequency. The intrinsic reaction coordinate (IRC) method was used to follow the imaginary eigenvector towards both the reactant complex and the cycloadduct. Optimized structures were illustrated using CYLview. [43]

Activation Strain Model of Reactivity
Quantitative insight into the activation barriers associated with the studied cycloadditions was obtained by means of the activation strain model (ASM). [25,44] This involves decomposing the potential energy surface ΔE( ) along the reaction coordinate into the strain ΔE strain ( ) associated with deforming the reactants from their equilibrium geometry and the interaction ΔE int ( ) between the deformed reactants [Equation (1)].
ΔE( ) = ΔE strain ( ) + ΔE int ( ) (1) The ΔE strain ( ) is determined by the rigidity of the reactants and by the extent to which they must distort in order to achieve the transition state geometry. The ΔE int ( ) is usually stabilizing and is related to the electronic structure of the reactants and how they are mutually oriented over the course of the reaction.

Energy Decomposition Analysis
The interaction ΔE int ( ) between the deformed reactants is further decomposed into three physically meaningful terms, in the conceptual framework provided by the Kohn-Sham molecular orbital (KS-MO) model [Equation (2)]. [45] ΔE int ( ) = ΔV elstat ( ) + ΔE Pauli ( ) + ΔE oi ( ) (2) The ΔV elstat ( ) term corresponds to the classical electrostatic interaction between unperturbed charge distributions ρ A (r) + ρ B (r) of the deformed fragments A and B and is usually attractive. The Pauli repulsion ΔE Pauli ( ) comprises the destabilizing interactions between occupied orbitals and is responsible for any steric repulsion. The orbital interaction ΔE oi accounts for charge transfer (interaction between occupied orbitals on one fragment with unoccupied orbitals of the other fragment, including HOMO-LUMO interactions) and polarization (emptyoccupied orbitals mixing on one fragment due to the presence of another fragment).
In the activation strain and energy decomposition plots, the IRC is projected onto the average distance of the newly forming C···X bonds (X = C or N depending on the 1,3-dipole). This resulting reaction coordinate undergoes a well-defined change in the course of the reaction from the reactant complex to the average C···X distance in the transition state and cycloadducts. The activation strain analysis (ASA) was performed with the aid of the PyFrag program [46] along the reaction coordinate.

Results and Discussion
The transition state structures, corresponding activation energies (ΔE ‡ ), and reaction energies (ΔE rxn ) for the 1,3-dipolar cycloaddition reactions between dipoles 1-4, 7-9 and ethylene (e) are shown in Figure 2. The trends in reactivity for these cycloadditions are identical when considering total electronic activation energy (ΔE ‡ ), enthalpy of activation (ΔH ‡ ), or Gibbs activation free energies (ΔG ‡ ) (see Table S1 in the Supporting Information). These series of aza-1,3-dipoles were selected due to their relevance in chemical synthesis, whereas ethylene was chosen for the sake of simplicity. In the dipolar cycloaddition reactions of dipoles 1-4, nitrile ylide (1) is the most reactive compound and proceeds with a low reaction barrier of only 5.4 kcal mol -1 . The activation energy for nitrile imine (2) is slightly higher than for 1. The barrier energies of dipolar cycloaddition reactions for diazomethane (3) and hydrazoic acid (4) are 13.1 and 16.6 kcal mol -1 , respectively, which is 2-3 times larger than that of 1. For the dipolar cycloaddition reactions of dipoles 7-9, azomethine ylide (7) has the lowest reaction barrier (2.1 kcal mol -1 ), which is even lower than that of 1. Azomethine imine (8) and azonium imine (9) are less reactive than 7 and proceed with barriers of 7.0 and 14.8 kcal mol -1 , respectively. The lengths of the newly forming bonds and dipole angle in the transition states are shown in Figure 2. From 1TS-e to 4TSe, the length of the shorter forming bond, as well as the longer forming bond decreases. This shift towards a later transition state is concomitant with a higher reaction barrier and a less exothermic reaction energy. Similar trends are observed for transition states 7TS-e to 9TS-e. Specifically, the reactions for 1 and 7 proceed via the earliest transition states, while 4 and 9 have a much later transition state. Earlier transition states lead to less distorted reactants in the transition state, thus the angles of dipoles 1 and 7 at the transition states are the largest ones among dipoles 1-4 and 7-9, respectively.

Activation Strain Analysis
Application of the activation strain model (ASM) has yielded quantitative insights into the origins of the reactivity differences for reactions 1e (black), 2e (green), 3e (blue), 4e (red) (Figure 3).  The activation strain diagram (ASD) shows our results for the 1,3-dipolar cycloaddition of 1-4 with e (reactions 1e-4e). The strain, ΔE strain , and interaction, ΔE int , curves vary greatly for these reactions. The lowest barrier associated with reaction 1e, proceeds with the largest strain energy along the entire reaction coordinate. This can be understood by the fact that 1 begins bending earlier on the reaction coordinate, thus resulting in an earlier, more destabilizing strain energy in contrast to 4, which bends at a later point along the coordinate. Interestingly, the strain energy associated with the in-plane bending away from the dipole's equilibrium geometry is nearly the same for dipoles 1-4 (see Figure S1 in the Supporting Information). The origin of these trends can be traced back to the orbital overlap between the FMOs of 1-4 and e. Specifically, the overlap for 1e occurs much earlier than for 4e due to larger porbitals on C compared to N, thus resulting in both a more destabilizing strain energy, but also an increasingly stabilizing interaction energy. [47] The reactivity differences of dipoles 1-4 originate primarily from the differences in the strengths of the interactions between deformed reactants along the reaction coordinate. Figure 3b shows the energy decomposition analysis (EDA) of the interaction energy for reactions 1e-4e. Analysis of the EDA terms reveals that the trends in the destabilizing Pauli repulsion, ΔE Pauli , are more or less offset by the trends in the stabilizing electrostatic interaction, ΔV elstat . The orbital interaction, ΔE oi , on the other hand, plays an important role in determining the trends in interaction energy, and thus, the height of the reaction barriers.
The factors controlling the reactivity of dipoles 7-9 are similar to what has been discussed above for the related dipoles 1-4 ( Figure 4). Inspection of the ASD (Figure 4a) reveals that the reaction with the lowest barrier, 7e, proceeds with the most unfavorable ΔE strain along the reaction coordinate. Interestingly, the ΔE strain is overcome by an even more stabilizing ΔE int , which pulls the barrier down significantly. Reactivity differences between dipoles 7-9 result due to differences in ΔE int . The energy decomposition of the ΔE int for reactions 7e-9e is shown in Figure 4b. Analysis of the EDA terms shows that the trends in ΔE Pauli are offset to a large degree by stabilization through electrostatic interaction. Again, ΔE oi strongly influences the trends in ΔE int and are mainly responsible for determining the reactivity in reactions 7e-9e.

Frontier Molecular Orbital (FMO) Analysis
The key factor in determining the reactivity of these 1,3-dipolar cycloaddition reactions are the orbital interactions. In order to provide a quantitative evaluation of ΔE oi , the molecular orbital (MO) diagrams and overlaps were assessed at the BP86/TZ2P level. Our Kohn-Sham MO analyses were performed on consistent geometries with average C···X bond forming distances of 2.30 Å for dipoles 1-4 and 7-9 to ensure that our results are not skewed by the fact that the transition states shift significantly along the reaction coordinate. Frontier molecular orbital (FMO) analysis reveals the importance of both the normal electron demand (HOMO 1-4 -LUMOe) and inverse electron demand (HOMOe-LUMO 1-4 ) orbital interactions for reactions 1e-4e ( Figure 5). For 1e to 4e, the calculated HOMO 1-4 -LUMOe energy gaps and orbital overlaps range from 3.3 to 5.4 eV and 0.29 to 0.10, respectively. These trends are consistent with the ΔE oi trends, which, as we know from the EDA results, diminish from 1 to 4. A second normal demand FMO interaction, present in reactions 2e and 4e, involving HOMO-1 2,4 and LUMOe occurs with a slightly larger FMO gap. The HOMO-1 1-4 are perpendicular to the HOMO 1-4 and have essentially no overlap and therefore no interaction with LUMOe. However, due to the tilted approach of the two reactants in reactions 2e and 4e, this interaction leads to favorable and productive overlap. This interaction, present in 2e and absent in 3e, compensates for the weaker primary orbital interactions for the former reaction and leads to an overall more stabilizing Figure 6. (a) MO diagrams with calculated energy gaps and orbital overlaps for the normal demand HOMO 7-9 -LUMO e interaction and (b) the inverse demand HOMOe-LUMO 7-9 interaction from the cycloaddition reactions between dipoles 7 (black), 8 (blue), 9 (red), and ethylene (e). All data computed at the BP86/ TZ2P level at a consistent geometry with the average C···X bond forming distances of 2.30 Å. orbital interaction. The inverse electron demand energy gaps (HOMOe-FMO 1-4 ) range from 4.0 to 3.3 eV. The orientation of the LUMO of 3, being perpendicular to the HOMOe, prevents it from overlapping in a favorable manner and thus the HOMOe-LUMO+13 interaction predominates. The computed overlaps are greatest for 1 (S = 0.24) and diminish as the number of nitrogens increases, for example in 4 (S = 0.17), which can be explained by the contracted nature of the p-orbital of nitrogen compared to carbon. [47] The important FMO interactions for reactions 7e to 9e involve both the HOMO 7-9 -LUMOe and HOMOe-LUMO 7-9 (Figure 6). The energy gaps for the normal demand (HOMO 7-9 -LUMOe) interaction become larger as the number of nitrogen atoms increases in the 1,3-dipole, from 7 to 9. The smallest HOMO7-LUMOe gap for 7, only 2.3 eV, also has the most favorable orbital overlap of S = 0.28. The overlap becomes less favorable as the number of nitrogen atoms increase from 7 to 9, again due to the more contracted nature of nitrogen p-orbitals compared to those of carbon. The HOMOe-LUMO 7-9 energy gaps for 7e to 9e range from 4.6 to 3.9 eV. The increased orbital interactions for 7 are reinforced by significant orbital overlap (S = 0.24). Overlap is less efficient for 9 and the value decreases to only S = 0.15. Interestingly, there is a strong correlation (R 2 = 0.98) between the activation barrier ΔE ‡ and S 2 /Δε for the HOMO 7-9 -LUMOe and HOMOe-LUMO 7-9 interactions at the transition state for reactions 7e-9e (See Figure S2 in the Supporting Information), further highlighting the significant role of FMO interactions in determining the computed cycloaddition reactivities for these 1,3-DCA reactions.

Conclusions
We show quantum chemically that the reactivity of two series of aza-1,3-dipoles in cycloadditions with ethylene decreases as the number of nitrogen atoms in the dipole increases. The two series of aza-1,3-dipoles comprise: (a) nitrile ylide (1), nitrile imine (2), diazomethane (3) and hydrazoic acid (4); and (b) azomethine ylide (7), azomethine imine (8) and azonium imine (9). Our computations were performed using the BP86/TZ2P method, which was, for the first time, benchmarked using the free and highly-automated QMflows software. Compared to G3B3 benchmark data, the computationally efficient BP86 functional provides relatively accurate activation and reaction enthalpies.
Our activation strain analyses reveal that the vast differences in reactivity of the 1,3-dipoles originates from significant differences in interaction energies. These reactions do not operate under strain control. In fact, the fastest reactions have the most unfavorable strain curves! Trends in orbital interactions are the decisive factor controlling the reactivity of aza-1,3-dipolar cycloadditions. They become less stabilizing as the number of hetero atoms in the dipole increases along 1-4 and 7-9. One reason is the stabilization of the dipole's HOMO which leads to an increase in the normal-electron demand HOMOdipole-LUMOdipolarophile energy gap. The other reason is the more compact nature of the nitrogen 2p orbitals, which overlap less efficiently with the carbon 2p orbitals of ethylene.
These results have important implications for the rational design of novel 1,3-dipolar cycloaddition reactions. One must have a firm grasp of the underlying physical mechanism by which reactivity is controlled, prior to designing elaborate chemical syntheses. Performing the activation strain analysis along the reaction coordinate is highly advised when the position of the transition states shift from early to late. Analysis at the transition state alone can be skewed and leads to incorrect interpretations of the factors governing reactivity. Strain energy is undeniably the main contributor to the activation barrier, but the differences in interaction energy for these 1,3-dipolar cycloadditions determine the trends in reactivity. Our study highlights the ability for the activation strain model, coupled with quantitative Kohn-Sham MO theory, to pinpoint the factors governing reactivity. This robust methodology can be extended to any number of novel chemistries, a feat we hope to see in the near future.