Rationalizing the Regioselectivity of the Diels-Alder Bis- Cycloaddition of Fullerenes

The physical factors governing the regioselectivity of the double functionalization of fullerenes have been explored by means of Density Functional Theory calculations. To this end, the second Diels-Alder cycloaddition reactions involving 1,3-butadiene and the parent C60-fullerene as well as the ion-encapsulated system Li@C60 have been selected. In agreement with previous experimental findings on related processes, it is found that the cycloaddition reaction, involving


Introduction
The chemistry of fullerenes has arguably experienced a remarkable development since their discovery in 1985. 1 This development can be ascribed to the numerous applications of these species in materials science, photovoltaic solar cells or even in biomedicine. 2 Consequently, chemists have designed a good number of synthetic procedures toward novel fullerene derivatives in order to control or tune their associated properties. 2,3 Among them, 1,3-dipolar and Diels-Alder cycloaddition reactions 4 or the so-called Bingel-Hirsch cyclopropanation reaction 5 should be specially highlighted.
In this sense, multiple functionalization of both empty 6 and endohedral fullerenes 7 has been pursued since the very beginning, as the corresponding derivatives are, for instance, very active as acceptors in organic photovoltaic solar cells 6 or can be used, as recently reported by Martín and co-workers, to produce giant glycofullerenes as inhibitors of Ebola virus infection. 8 Despite that, the factors controlling the regioselectivity of multiple additions to fullerenes remain poorly understood. This is of crucial importance because the properties of the different regioadducts may differ significantly. For instance, it was confirmed that some regioisomerically pure fullerene bisadducts show a better performance in photovoltaic cells than the corresponding isomeric mixtures. 9 Recently, we have successfully applied the combination of the so-called Activation Strain Model (ASM) 10 of reactivity and the Energy Decomposition Analysis (EDA) 11 methods to gain quantitative insight into the selectivity of different Diels-Alder cycloaddition reactions. 12 This approach has been particularly useful to understand the factors controlling the [6,6]-over [5,6]regioselectivity of the Diels-Alder reaction involving the parent C60-fullerene 13 and strongly 3 related fullerene fragments. 14 In addition, this method has been really helpful to not only understand but also to predict the regioselectivity of the Diels-Alder reaction involving the related aza[60]-fullerene, C59NH. 15 Given the good performance of this method, we decided herein to investigate in detail the so far poorly understood factors governing the regioselectivity of multiple cycloadditions to fullerenes. To this end, we have selected the second Diels-Alder reaction between 1,3-butadiene and the corresponding [6,6]-cycloadduct (1) formed upon reaction of C60 and 1,3butadiene (Scheme 1). In addition, the influence of the encapsulation of a lithium cation inside the C60-cage on the process shall be considered as well. Scheme 1. Diels-Alder successive cycloaddition reactions between C60 or Li + @C60 and 1,3butadiene.

Computational Details
Geometry optimizations of the molecules were performed without symmetry constraints using the Gaussian03 16 optimizer together with Turbomole 6.6 17 energies and gradients at the BP86 18 /def2-SVP 19 level of theory using the D3 dispersion correction suggested by Grimme et al. 20 and the resolution-of-identity (RI) approximation. 21 This level is denoted RI-BP86-D3/def2-SVP and has been selected because it provided very good results for Diels-Alder reactions involving related fullerenes. [13][14][15]22 Reactants and cycloadducts were characterized by frequency calculations, and have positive definite Hessian matrices. Transition states (TSs) show only one negative eigenvalue in their diagonalized force constant matrices, and their associated eigenvectors were confirmed to correspond to the motion along the reaction coordinate. Connections between 4 reactant complexes and products were confirmed using the Intrinsic Reaction Coordinate (IRC) method. 23 Single-point energy refinements were carried out using the D3-corrected metahybrid M06-2X 24 functional in conjunction with the triple-ζ-quality def2-TZVPP basis sets. 19 This level is therefore denoted M06-2X-D3/def2-TZVPP//RI-BP86-D3/def2-SVP.

Activation Strain Analyses of Reaction Profiles
The so-called activation strain model, 10 also known as distortion/interaction model, 25 has allowed us to gain quantitative insight into the physical factors controlling different fundamental processes. 26 The activation strain model is a systematic development of the Energy Decomposition Analysis (EDA) method (see below) to understanding chemical reactions, in which the height of reaction barriers is described and understood in terms of the original reactants. In this approach, the potential energy surface ΔE() is decomposed, along the reaction coordinate , into two contributions, namely the strain energy ΔEstrain(), associated with the structural deformation suffered by the reactants along the reaction, and the actual interaction ΔEint() between these increasingly deformed reactants: It is the interplay between these two contributions which determines if and at which point along  a barrier arises, namely, at the point where dΔEstrain()/d = -dΔEint()/d. For further details of the theoretical background and different applications of the ASM method, we refer readers to the review articles that were published recently. 10 In the cycloaddition reactions considered herein, the reaction coordinate is defined as the projection of the IRC onto the shortest forming C···C distance between the carbon atom of the fullerene and the carbon atom of the diene. This reaction coordinate  undergoes a well-defined change in the course of the reaction from the initially formed reactant complexes to the equilibrium C···C distance in the corresponding TSs.

5
The Energy Decomposition Analysis (EDA) 11 method can be used to further decompose the interaction energy into the following chemically meaningful terms: The term ΔEelstat corresponds to the classical electrostatic interaction between the unperturbed charge distributions of the deformed reactants and is usually attractive. The Pauli repulsion ΔEPauli comprises the destabilizing interactions between occupied orbitals and is responsible for any steric repulsion. The orbital interaction ΔEorb accounts for electron pair bond formation (interaction between singly occupied MOs in each fragment), charge transfer (interaction between occupied orbitals on one moiety with unoccupied orbitals on the other, including HOMO-LUMO interactions) and polarization (empty-occupied orbital mixing on one fragment due to the presence of another fragment). Finally, the ΔEdisp term takes into account the interactions that are due to dispersion forces.
Moreover, the NOCV (Natural Orbital for Chemical Valence) 27,28 extension of the EDA method has been also used to further partitioning the ∆Eorb term. The EDA-NOCV approach provides pairwise energy contributions for each pair of interacting orbitals to the total bond energy. The EDA calculations reported herein were carried out at the dispersion corrected BP86-D3 23 /TZ2P 29 level using the optimized RI-BP86-D3/def2-SVP geometries with the ADF 2017 program package. 30 This level is therefore denoted BP86-D3/TZ2P//RI-BP86-D3/def2-SVP.

Results and Discussion
We first focused on the regioselectivity of the cycloaddition reaction involving the [6,6]cycloadduct 1 and 1,3-butadiene (Scheme 1). Eight different [6,6]-pyracylennic bonds in 1 may serve as a dienophile in the reaction with 1,3-butadiene (see Figure 1). In addition, two possible conformers (two possible conformations for the formed six-membered ring) per [6,6]-bond can be produced in the Diels-Alder reaction. Our calculations indicate that the barrier and reaction energy 6 differences between both approaches is in most cases negligible (< 0.5 kcal mol -1 ) and therefore, below we only refer to the most favored conformer for each [6,6]-bond. Figure 1. [6,6]-bonds in cycloadduct 1 considered in this study (highlighted in green). The usual nomenclature used for these bonds is given in parentheses.
As readily seen in Figure 2, which shows the computed reaction profile for the formation of regioisomer CA4, there are no significant differences between the first and second cycloaddition reactions, which agrees with previous computational studies. 31 Thus, they both proceed concertedly via a highly synchronous transition state from an initial van der Waals reactant complex, which lies ca. -3.0 kcal mol -1 below the separate reactants, and lead to the corresponding cycloadduct in a highly exothermic process (see Figure 2). The other [4+2] cycloadditions forming the regioisomers CA1-8 proceed similarly (for a representation of the associated transition states TS1-8, see Figure S1 in the Supporting Information).  Table 1 gathers the computed activation and reaction energies (BP86-D3/TZVP//RI-BP86-D3/def2-SVP level) of the Diels-Alder reactions involving the initially formed [6,6]-cycloadduct 1 and 1,3-butadiene which lead to cycloadducts CA1-8. Our calculations indicate that, for all bonds 1-8, formation of bis-adducts is thermodynamically and kinetically less favored (by at least 0.5 kcal mol -1 ) than the production of the monoadduct. Moreover, we find that the formation of cycloadducts CA1 and CA2 is kinetically much more difficult than that of the rest of cycloadducts CA4-8 (∆∆E  ca. 8 kcal mol -1 ). A similar result is found for cycloadduct CA3, although in this case the activation barrier difference is lower (∆∆E  ca. 4 kcal mol -1 ). Moreover, these cycloadducts (CA1-3) are formed in comparatively less exothermic processes (∆∆ER ca. 6-8 kcal mol -1 ), which suggests that these regioisomers are selectively not produced in the reaction. This computational prediction is fully consistent with the previous experimental results reported by 8 Kräutler and co-workers on the bis-adducts formed in the strongly related reaction between C60 and anthracene which leads to a mixture of the corresponding cycloadducts CA4-8 (in similar reaction yields) with no trace of regioisomers CA1-3. 32 The lack of formation of cycloadducts CA1-3 observed experimentally can be initially ascribed to the bulkiness of the anthracene moiety.
However, our calculations involving the much smaller 1,3-butadiene show the same reactivity trend, thus indicating that the steric hindrance is not the exclusive factor responsible for the low reactivity of [6,6]-bonds 1-3. Finally, although shorter C-C bonds in fullerenes are typically considered as the most reactive bonds, our calculations indicate that there is no clear relationship between reactivity and the computed C-C bond lengths in cycloadduct 1 (see Table 1).  We applied next the Activation Strain Model (ASM) of reactivity to gain more quantitative insight into the physical factors governing the different reactivity of the [6,6]-bonds 1-8 which ultimately defines the regioselectivity of the selected bis-Diels-Alder reaction. To this end, we 9 considered one of the least reactive bonds (bond 1), one of the most reactive bonds (bond 6) and also bond 3, which has an activation barrier intermediate between those computed for the former bonds. The corresponding activation strain diagrams (ASDs) for the cycloadditions involving these bonds from the corresponding reactant complexes up to the respective transition states are depicted in Figure 3. As readily seen in this figure, all processes exhibit rather similar ASDs in the sense that the interaction energy between the deformed reactants, measured by the ∆Eint term, becomes slightly destabilizing at the beginning of the cycloaddition and then inverts at a certain point along the reaction coordinate (i.e., at forming C···C distances of ca. 2.5 Å) thus becoming more and more stabilizing as approaching the corresponding transition state. Differently, the strain energy (∆Estrain term) becomes monotonically more and more destabilizing along the reaction coordinate and is able to compensate the stabilizing effect of the interaction energy at the transition state region. A similar behavior has been found not only for related Diels-Alder reactions 12-15 but also for other pericyclic reactions. 33 Despite that, significant differences are found when comparing the different ASDs. Thus, the much higher activation barrier computed for bond 1 (leading to cycloadduct CA1) can be clearly ascribed to the much higher deformation energy required for this particular process along the entire reaction coordinate and in particular, at the transition state region. This is consistent with the computed rather low pyramidalization angle of 9.82º associated with the carbon atoms of this bond (for comparison, the pyramidalization angle of C60 is much higher: 11.64º). 34 At variance, bond 6 (leading to cycloadduct CA6) benefits from the lowest strain energy (pyramidalization angle of 11.68º) as well as the highest interaction energy between the deformed reactants practically along the entire reaction coordinate. As a result, bond 6 becomes the most reactive bond in the selected series. In the case of bond 3 (leading to cycloadduct CA3), although the corresponding cycloaddition presents a nearly identical (albeit slightly less destabilizing) strain energy than the process involving bond 6, the computed interaction energy between the deformed reactants is significantly weaker. For instance, at the same C···C forming distance of 2.35 Å, a value of ∆Eint = -4.5 kcal mol -1 was computed for the reaction involving bond 6 whereas a lower (i.e. less stabilizing) value of ∆Eint = -1.8 kcal mol -1 was computed for the cycloaddition involving bond 3 ( Figure 3). As a consequence, the activation barrier computed for the cycloaddition involving the latter bond is intermediate between those involving bonds 1 and 6. According to the EDA method, which further decomposes the interaction energy into chemically meaningful contributions (see above), the stronger interaction energy computed for the process involving bond 6 (as compared to bond 3) is exclusively due to stronger orbital interactions between the reactants along the reaction coordinate as the rest of the contributions (namely, ∆EPauli, ∆Velstat and ∆Edisp) are identical for both cycloaddition reactions (Figure 4). The origins of the stronger orbital interactions in the process involving 1 can be found by using the NOCV (Natural Orbital for Chemical Valence) 27,28 extension of the EDA method. The main orbital interaction identified by the EDA-NOCV method corresponds to the π(diene)π*(fullerene) interaction, which confirms the normal electronic demand nature of these [4+2]-cycloaddition reactions (see Figure 5, charge flow takes place in the direction redblue).
Interestingly, this molecular orbital interaction is clearly stronger for the reaction involving bond 6 (see the data computed for the processes at the same C···C bond forming distance of ca. 2.4 Å), which is translated into the computed higher orbital attractions (measured by the ∆Eorb term) discussed above. Once the regioselectivity of the second cycloaddition reaction involving the parent C60-fullerene has been analyzed, we next considered the analogous transformation involving the ionencapsulated counterpart Li + @C60 in order to investigate the influence of the charge on the regioselectivity. This endohedral system has proven to be much more reactive (ca 2400 times faster) than C60 in the otherwise same Diels-Alder reactions with cyclopentadiene or 1,3hexadiene. 35 The reasons for this enhanced reactivity have been analyzed in detail by us recently using the ASM-EDA(NOCV) method. 22 In addition, Li + @C60 has been also found to have interesting properties in materials science. 36 Despite that, nothing is known about the multiple functionalization of this species.
Although the scenario in this endohedral system is somewhat more complicated than that found in the parent C60-fullerene due to the mobility of the lithium cation inside the fullerenic cage, we focused, once again, on the same eight chemically different [6,6]-bonds for comparison reasons with the lithium cation close to center of the cage. Similar to the processes involving C60, the successive cycloaddition reactions involving Li + @C60 also proceed concertedly from an initial reactant complex lying below the separate reactants (see Figure S2 in the Supporting Information for a representation of the associated transition states TS1-8-Li). Table 2 gathers the computed activation and reaction energies of the Diels-Alder reactions involving the initially formed [6,6]-cycloadduct 1-Li and 1,3-butadiene and leading to cycloadducts CA1-8-Li together with the values for the initial reaction involving Li + @C60. As expected, the computed activation barriers for the processes encapsulating the lithium cation are systematically lower (∆∆E ‡ ca. 4 kcal mol -1 ) than the analogous processes involving the hollow C60, which is consistent with the observed acceleration of the reaction induced by the cation. 35,37 At variance with the cycloaddition of 1 and 1,3-butadiene, now, for some bonds, bis-addition is favored thermodynamically (bonds 1, 2, 4, 6, and 7) and kinetically (bonds 6 and 7) as compared to monoaddition. Regarding the regioselectivity of the second cycloaddition, the same trend to that found for cycloadduct 1 is observed for 1-Li. Thus, the formation of cycloadducts CA1-Li, CA2-Li and CA3-Li (the latter in a lesser extent) is kinetically hampered in comparison to the rest of cycloadducts CA4-8-Li, which exhibit much lower activation barriers (∆E ‡ ca. 10-11 kcal mol -1 ).
Therefore, it can be concluded that the encapsulation of a lithium cation inside the C60-cage, although induces a remarkable enhancement of the exohedral Diels-Alder reactivity of the system, has a negligible influence on the regioselectivity of the successive cycloaddition reaction. Li + @C60 -4.
The ASM method was also applied to confirm the factors controlling the regioselectivity of the second cycloaddition reaction between cycloadduct 1-Li and 1,3-butadiene. Once again, we focused on the most representative bonds 1, 3, and 6, whose ASDs from the initial reactant complexes up to their corresponding transition states are depicted in Figure 6. Similar conclusions to those discussed above for cycloadduct 1 can be drawn for this cationic fullerene. Thus, the high value of the deformation energy (∆Estrain) computed for the reaction involving bond 1 is again responsible for its high activation barrier which results in its lack of formation. This can be qualitatively attributed to the computed rather low pyramidalization angle of 9.77º of the carbon atoms of this particular bond. Similarly, the most reactive bond 6 exhibits a much lower strain energy (pyramidalization angle of 11.68º) and a stronger interaction energy between the reactants, particularly at the transition state region. Although bond 3 requires a nearly identical deformation energy to the process involving bond 6, the associated interaction energy is comparatively much weaker, and as a result, the activation barrier computed for this process is intermediate between those computed for the cycloaddition involving bonds 1 and 6. According to the EDA method, the stronger interaction computed for the most reactive bond 6 derives almost exclusively from the much higher (i.e. more stabilizing) orbital interactions between the reactants occurring at this bond along the entire reaction coordinate (see Figure 7).  We finally used the EDA-NOCV method to confirm that, once again, bond 6 benefits from a higher π(diene)π*(fullerene) interaction which is responsible for the computed stronger orbital interactions between the reactants. Indeed, as readily seen in Figure 8, the associated orbital energy (∆E()) is clearly higher (i.e. more stabilizing) for the process involving bond 6 (data computed for the processes at the same C···C bond forming distance of ca. 2.4 Å). Moreover, the ∆E() values computed for 1-Li are about 4 kcal mol -1 more stabilizing than those computed for 1, which mainly derives from the stabilization of the LUMO (i.e. the π* molecular orbital) of the fullerene induced by the Li + cation. Interestingly, the computed energy difference ∆∆E() = 2.0 kcal mol -1 is rather similar to that computed for the cycloaddition involving the parent cycloadduct 1 (∆∆E() = 2.2 kcal mol -1 , see Figure 5), which further confirms the negligible role of the encapsulated cation on controlling the regioselectivity of the transformation.

Conclusions
The second Diels-Alder reaction involving 1,3-butadiene and the cycloadduct initially formed in the reaction between C60 and 1,3-butadiene also proceeds concertedly through highly synchronous transition states. Among the eight different [6,6]-bonds which can serve as dienophile, the cycloaddition occurs selectively at bonds 4-8 in view of their computed lower activation barriers. According to the ASM of reactivity, it is the interplay between the strain and interaction energies which controls the regioselectivity of the transformation. For instance, the cycloadditions involving bonds 1-3 are kinetically disfavored due to the much higher deformation energy required by the reactants to adopt the geometries of the corresponding transition states.
Differently, one of the most reactive bonds, bond 6, benefits not only from a lower strain energy but also from a comparatively stronger interaction between the deformed reactants. According to the EDA method, this is exclusively due to much higher orbital attractions between the reactants, mainly as a result of a stronger π(diene)π*(fullerene) interaction. The same results are found for the analogous process involving the ion-encapsulated system Li + @C60, therefore indicating that the endohedral cation has a negligible influence on the regioselectivity of the considered Diels-Alder cycloaddition.

Supporting Information Available.
Figures S1 and S2, Cartesian coordinates (in Å) and total energies of all the stationary points discussed in the text. This material is available free of charge via the Internet at http://pubs.acs.org.