Reaction Mechanism and Regioselectivity of the Bingel- Hirsch Addition of Dimethyl Bromomalonate to La@C2v-C82

We have quantum chemically explored the thermodynamics and kinetics of all 65 possible mechanistic pathways of the Bingel-Hirsch addition of dimethyl bromomalonate to the endohedral metallofullerene La@C2v-C82 using dispersion-corrected DFT that result from the combination of 24 nonequivalent carbon atoms and 35 different bonds present in La@C2v-C82. Experimentally, this reaction leads to four singly-bonded derivatives and one fulleroid adduct. Out of these five products, only the singly-bonded derivative on C23 could be unambiguously identified. Our calculations show that La@C2v-C82 is not particularly regioselective under Bingel-Hirsch conditions. From the obtained results, however, it is possible to make a tentative assignment of the products experimentally observed. We propose that the observed fulleroid adduct results from the attack to bond 19 and that the singly-bonded derivatives correspond to the C2, C19, C21, and C23 initial attacks. However, other possibilities cannot be ruled out completely.

Scheme 1. a) The 24 nonequivalent carbon atoms and the 35 different bonds of La@C2v-C82. For bonds, numbers denote [6,6] bonds and lower-case letters denote [5,6] bonds. Bond types are distinguished by color code (type A: blue, type B: black, type C: green, type D: red). Labels as assigned in a previous study. [18] b) General mechanism of the nucleophilic [2+1] Bingel-Hirsch reaction in La@C2v-C82. Nomenclature used in the current work is also provided. The attack on a bond γ formed between adjacent carbon atoms Cα and Cβ leads to the formation of an intermediate, either ICα,γ or ICβ,γ; thus proceeding to the respective transition state, either TSCα,γ or TSCβ,γ, which leads to the formation of a methanofullerene *PCαCβ,γ or fulleroid PCαCβ,γ.
Chemical functionalization of EMFs [5,6] is commonly achieved through cycloaddition reactions, being the most important: Diels-Alder, [18,29] 1,3 dipolar (or Prato), [30] and nucleophilic [2+1] Bingel-Hirsch (BH) additions. [31][32][33] The BH addition is a cyclopropanation reaction where a fullerene and diethyl bromomalonate in the presence of a strong base such as 1,8diazabicycloundec-7-ene (DBU) or sodium hydride react to produce a methanofullerene or a fulleroid. [34] In the first step of the mechanism, the base abstracts the acidic proton of the malonate derivative to generate a carbanion or enolate. Then this carbanion nucleophilically attacks the fullerene, thus generating a new carbanion with charge localized at the cage. In the final step, bromide is displaced in a nucleophilic substitution SN2 reaction causing an intramolecular three-membered ring closure. Scheme 1b illustrates the mechanism of the BH reaction in La@C2v-C82, in which we have introduced the nomenclature used through the whole manuscript. Accordingly, the initial attack on bond γ formed between adjacent carbon atoms Cα and Cβ leads to the formation of a singly-bonded derivative which is labeled as the intermediate structure, either ICα,γ or ICβ,γ, followed by the corresponding transition state structure, either TSCα,γ or TSCβ,γ, leading to the formation of a closed-cage Bingel product methanofullerene *PCαCβ,γ or an open-cage fulleroid PCαCβ,γ. Moreover, since there are 19 [6,6] and 16 [5,6] different bonds in La@C2v-C82, γ can take 35  The first successful functionalization of EMFs through the BH reaction was achieved by Alford and co-workers, wherein Gd@C60[C(COOC2H5]n (n = 1-10) was the main product. [35] The same procedure applied to 212 Pb@C60 led to the formation of 212 Pb@C60[C(CO2H)2]x. [36] In 2005, Echegoyen et al. [37] reported the selective BH formation of Y3N@Ih-C80 and Er3N@Ih-C80 [6,6] monoadducts under mild conditions. [38] In some cases [6,6]-products could be further converted to [5,6]-adducts by heating. The same authors also noticed that the BH reaction did not take place in Sc3N@Ih-C80 and Lu3N@Ih-C80 under the same reaction conditions. [37] X-ray structure of Y3N@Ih-C80[C(CO2CH2Ph)2] confirmed that the BH monoadduct was the result of the attack to a [6,6] bond that leads to the formation of a fulleroid structure with one of the yttrium atoms directly pointing to the attacked open bond. [39] DFT calculations indicated that the open [6,6] adduct is more stable than the closed one and that the rotation of the metallic cluster inside the cage is partially hindered once the adduct is formed. [39] Later, in 2010, Echegoyen and coworkers reported that Sc3N@Ih-C80 and Lu3N@Ih-C80 can undergo BH additions to yield open [6,6] adducts by changing the common BH reaction conditions. NaH was used as the base and a 4:1 mixture of o-DCB and N,N-dimethylformamide as the solvent. [40] Interestingly, the BH reaction of TiSc2N@Ih-C80 afforded two unconventional singly bonded monoadducts, revealing a change 5 in the addition pattern and an improved reactivity of TiSc2N@Ih-C80 when compared to Sc3N@Ih-C80. [41] A study of the BH reaction of Gd3N@C80, Gd3N@C84, and Gd3N@C88 performed by Echegoyen et al. led to the conclusion that EMFs reactivity is reduced when the size of the fullerene cage increases. [42] BH additions to non-IPR cages have been experimentally achieved for Sc3N@D3(6140)-C68, [43] Gd3N@Cs(39663)-C82, and Gd3N@Cs(51365)-C84 EMFs. [42,44] For the former, according to theoretical studies [45] and experiments, [43] the addition takes place at a [6,6] bond close to a [5,5] bond and to the Sc atom. For the other systems, [42,44] the BH addition occurred on a [5,6] bond adjacent to the unique [5,5] bond. Unexpectedly, the addition never took place on the a priori more reactive [5,5] bond. This result was attributed to the larger aromaticity of the adducts of the [6,6]-and [5,6]-additions to the Sc3N@D3(6140)-C68 and Gd3N@Cs(51365)-C84 EMFs, respectively. [46] It is worth mentioning that theoretical studies of the reaction mechanism of BH additions to Sc3N@Ih-C80 and Sc3N@D3(6140)-C68 [45] and to Y3N@Cs(39663)-C82 and Y3N@Cs(51365)-C84 [47] were performed by Poblet and co-workers.
Their results show that the BH thermodynamic and kinetic products do not coincide, being the kinetically controlled product the one observed under experimental conditions. Very recently, some of us in collaboration with Echegoyen's group have studied experimentally and theoretically the BH to Sc3N@D5h-C80. [48] Our results show that the addition takes place under kinetic control in a [6,6] bond and confirm that the most stable thermodynamic and kinetic products differ.
The BH addition to La@C2v-C82 was experimentally explored by Nagase, Akasaka, and coworkers [20,49,50] using diethyl bromomalonate as a reactant. The BH reaction was monitored by multistage high-performance liquid chromatography (HPCL), which allowed the identification of different monoadducts. They could isolate four ESR-inactive adducts, which were assigned to singly-bonded products, and one ESR-active species that was considered to be a BH adduct. Xray crystallographic analysis of one ESR-inactive product confirmed a singly-bonded derivative of La@C2v-C82 that corresponds to an oxidized BH intermediate formed at C23 position (see Scheme 1a). Moreover, when heated to 80 ºC in anhydrous o-dichlorobenzene (o-DCB), all the synthesized singly-bonded products decompose to give the parent La@C2v-C82 as the major product. On the contrary, the BH product showed higher thermal stability. Based on the fact that a more positively charged carbon atom could be a more favored reaction site towards a nucleophilic attack, the authors proposed C18, C14, and C21 as the possible addition sites for the three unknown ESR-inactive products. The authors made their estimations based on Mulliken charge distribution analysis as determined by density functional theory (DFT) calculations. [49] We consider that the theoretical study of the BH addition to La@C2v-C82 is timely and relevant for three main reasons: i) because the relatively high abundance of La@C82 among endohedral 6 monometallofullerenes makes the reactivity of this EMF one of the most studied (although not fully understood yet); ii) because of the radical character of the fullerenic cage in La@C2v-C82, a non-conventional reaction mechanism for the BH addition to La@C2v-C82 can be expected (this is confirmed by the presence of singly-bonded derivatives in the observed products); and iii) because experiments were unable to definitely identify all adducts generated in the BH addition to La@C2v-C82. In this contribution, we have carried out an extensive quantum chemical (DFT) exploration of the various mechanistic pathways of the BH reaction between the anion dimethyl bromomalonate (dmbm -) and La@C2v-C82 through different reaction sites with two main objectives: i) to determine the reaction mechanism of the different reaction pathways; and ii) to evaluate the regioselectivity of the BH reaction of dmbmto La@C2v-C82 to provide a comprehensive description for the preference of different positions to react. For the particular case of the Diels-Alder reaction of (1,2,3,4,5-pentamethyl)cyclopentadiene to La@C2v-C82, [51] some of us [18] demonstrated by means of DFT calculations that this cycloaddition occurs regioselectively at bond o (shown in Scheme 1a) which corresponds to both, the thermodynamically and kinetically most favored attack, independently of the diene used for the reaction. We anticipate here that our results show that the BH reaction of dmbmto La@C2v-C82 is not as regioselective as the Diels-Alder reaction.

Computational Details
All DFT calculations were performed by using the Amsterdam Density Functional (ADF) program. [52] An uncontracted set of Slater-type orbitals (STOs) of double-ζ (DZP) and triple-ζ (TZP) quality containing diffuse functions and one set of polarization functions were used to expand the molecular orbitals (MOs). The frozen-core approximation (FCA) [52] was employed during the SCF procedure for the core orbitals of C, O, and La (1s for C and O and 1s2s2p3s3p4s3d4p for La). It was shown that the FCA has a negligible effect on optimized geometries. [53,54] Scalar relativistic corrections were included self-consistently by using the zeroth order regular approximation (ZORA). [55] The local density approximation (Slater exchange) with non-local corrections for exchange (Becke88) [56] and correlation (Perdew86) [57] (i.e. the BP86 functional) were used to selfconsistently calculate energies and gradients. Although it is well documented that standard DFT functionals like BP86 underestimate energy barriers, [58] this underestimation should be similar for all the BH transition states we encounter here and should not affect the main conclusions.
Calculations were executed under the unrestricted formalism. Moreover, dispersion energy corrections as developed by Grimme et al. [59,60] including the so-called Becke and Johnson damping (D3-BJ) were added in the calculation of DFT energies and gradients. It has been shown 7 that dispersion corrections are essential for a correct description of the thermodynamics and kinetics of reactions with fullerenes, nanotubes and other systems and reactions involving nonrepulsive steric contact. [61] Geometry optimizations were performed without symmetry constraints in the gas phase. All stationary points were characterized by analytical frequency calculations. Electronic energies were obtained in solution with the TZP basis set by single-point energy calculations at the geometries optimized with the DZP basis (i.e. BP86-D3(BJ)/TZP//BP86-D3(BJ)/DZP). Solvent effects were introduced using the conductor-like screening model (COSMO) [62] as implemented in ADF, performing single point energy calculations on the gas-phase optimized structures and using toluene as a solvent. In the Supporting Information (see Tables S1 and S2) we show that the gas-phase geometry is practically identical to the geometry optimized under the implicit presence of toluene.
In the search of minima and first-order saddle points the QUILD code (quantum regions interconnected by local descriptions) [63] was used. QUILD works as a wrapper around the ADF program; it creates input files for ADF, then executes this program and collects energies and gradients. QUILD uses adapted delocalized coordinates and constructs model Hessians with the appropriate number of eigenvalues. [64] This latter feature is particularly useful for the search of transition state structures. corrections by zero-point energies, thermal contributions to the internal energy, and the entropy term were determined in the gas phase at the BP86-(BJ-D3)/DZP level of theory using the statistical thermodynamics expressions of an ideal gas in standard conditions at 298 K. Lastly, to account for the condition change from 1 atm to 1 M concentration related to the phase change from gas to solution, we added to the Gtol values a concentration correction of 1.89 kcal/mol. [65][66][67]

Results and Discussion
This section is divided as follows. Firstly, the 65 possible reaction pathways are defined to introduce the reader to the many different stationary points. Then, based on a complete thermodynamic study of the system, a classification of available energy profiles is performed so as to identify the thermodynamic products and, in turn, the potentially accessible kinetic products. Such a classification is designed with the purpose of predicting the experimentally observed products, either a singly-bonded or a cyclopropanated derivative. The kinetic study 8 therefore is conducted through all those reaction pathways that, in fact, lead to a cyclopropanated derivative prone to be experimentally observed.

General aspects of the 65 possible reaction pathways
In the initial step of the BH addition, a singly-bonded derivative is formed when dmbmis attached to a specific carbon atom of La@C2v-C82. In the resulting structure, bromine is oriented toward an adjacent 6-MR or five-membered ring (5-MR) so that the reactive carbon atom of dmbmfaces the carbon atom of the cage that will be attacked (see Fig. 1). The rotation of the dihedral angle φ by ±120º in the resulting intermediate structure (φ is a structural parameter indicating the alignment between the bromine atom and a bond of the EMF cage; φ is equal to 0º or ±180º when bromine is exactly aligned with a [6,6] bond), generates three orientational isomers for each intermediate structure leading to different products. This is general, except when dmbmis linked to carbon atoms C2, C3, C5, C8, C17, C20, and C24, in which only two orientational isomers can be distinguished because of the symmetry of the cage. Every single nonequivalent carbon atom represents one reaction site, and the formation of distinct orientational isomers of a singly-bonded derivative gives rise to 65 different reaction pathways, some of them leading to the same product (e.g., cyclopropanation to bond 1 can be achieved from either C1 or C2, see Scheme 1a).

Figure 1.
Structural parameters of a singly-bonded derivative. The distances dn are geometrical descriptors of the cyclopropanation reaction happening at the bond highlighted in black (to the left). Angle θ and dihedral angle φ (this latter formed by the black-marked carbon atoms to the right) describe the position and orientation of the bromine atom (see text for details).
As mentioned in the introduction, it is well known that under typical conditions and in the case of EMFs, the BH reaction occurs under kinetic control. [47,68] Therefore, a complete DFT-based study must be performed through the 65 different reaction pathways available in La@C2v-C82.
However, it is important to keep in mind that there are two possible reaction pathways leading to PCαCβ,γ (or *PCαCβ,γ if the case). One of them follows the consecutive formations of ICα,γ and TSCα,γ; and analogously the other goes through ICβ,γ and TSCβ,γ. In this regard, based on the energetic stability of ICα,γ, ICβ,γ, and PCαCβ,γ (or *PCαCβ,γ), in Scheme 2 we draw a picture with the intention of representing the three situations susceptible to occur. Observe that profile type I is defined as the energy profile in which both reaction pathways lead to a Bingel product that is more stable than the two possible intermediate precursors (ICα,γ and ICβ,γ). If one reaction pathway could be reverted, then profile type II is defined, where the BH product is more stable than only one of the intermediate precursors. Finally, profile type III corresponds to the energy profile involving a destabilized Bingel product with respect to ICα,γ and ICβ,γ. For profiles II and III we expect that the corresponding BH retro-reaction will happen and the BH product will not be accumulated. Accordingly, the following subsections are conducted regarding the profile classification of the 65 reaction pathways with the aim of describing all the most thermodynamically and kinetically favored BH reactions.

Scheme 2.
Representation of the two reaction pathways leading to the formation of a Bingel-Hirsch adduct as a fulleroid or a methanofullerene through the attack on bond γ with adjacent reaction centers Cα and Cβ. The energy profile in which the retro-Bingel-Hirsch addition is hampered is classified as profile type I. On the contrary, if one or two reaction pathways can be regressed to a singly-bonded derivative then the resulting profiles are respectively classified as II or III. For every diagram, the most likely structure to be experimentally accumulated is enclosed in a rectangular box.

Classification of the 65 possible reaction pathways
Gibbs energies in toluene (ΔGtol) relative to separated reactants for the formation of the 65 different ICα,γ structures are reported in Table 1, as well as ΔGtol for the 35 most stable BH adducts (mostly PCαCβ,γ adducts). Based on Scheme 2 and the values of ΔGtol for the structures related to a single pathway, it is possible to define the profile type in which that reaction pathway can be classified (e.g. for the attack on bond 1, PC1C2,1 is more stabilized than IC1,1 and IC2,1; therefore the situation is described by profile type I). After analyzing the profile type for all the 65 available reaction pathways, as reported in Table 1, we find that only 15 attacks out of 35 are of profile type I. These 15 attacks can occur through 26 different pathways (C = C for bonds 9, 16, 18, and a) and generate relatively stable BH products. To determine the kinetics of the entire BH addition to La@C2v-C82, we focus on reaction pathways with profile type I because only in this situation the BH adduct is going to be experimentally observed. Furthermore, in principle the 65 different intermediates ICα,γ must be completely characterized. However, in a few cases exhibiting type II and III profiles, we did not optimize all possible orientational isomers of intermediate ICα,γ as they exhibit similar energetics (see superscript a in Table 1). To sum up, by strategically defining three different profiles of the BH addition to La@C2v-C82 based on ΔGtol of all the possible intermediates, methanofullerenes, and fulleroids, we can avoid the study of a number of reaction pathways and still provide a full description of this reaction; therefore in the following subsection the kinetic study is described for the 26 reaction pathways related to profile type I. .87 II a ΔGtol for these intermediates was estimated from the energy of the most stable orientational isomer for that specific reaction site. For instance, the energy of IC9,6 is not calculated, but it is taken as the energy of IC9,f since orientational isomers are energetically very similar. Moreover, from the energy of IC12,6 it is possible to anticipate that bond 6 is not indeed involved in a profile type I. b Only in these four cases, *PC7C10,4, *PC9C9,g, *PC16C19,12, and *PC22C22,16, the most stable Bingel-Hirsch adduct is a methanofullerene instead of a fulleroid structure.

Analysis of structural parameters through the reaction coordinate
Relevant structural parameters (d1, d2, d3, θ, and φ as defined in Figure 1) are summarized in Table 2 for the 26 reaction pathways under consideration. The distance d1 describes the bond length of the attacked C-C bond γ of the cage, d2 corresponds to the CC82-Cdmbm-that is formed during the ring closure, and d3 refers to the Cdmbm--CBr bond length. The angles θ and φ are useful parameters for the examination of the exact localization of bromine through the reaction coordinate. In addition to that, from the reaction sites Cα involved in those 26 reaction pathways, there are seven 666 positions (carbon atoms surrounded by three hexagons) and eleven 566 (carbon atoms surrounded by one pentagon and two hexagons) that lead to eleven [6,6] and four [5,6] BH adducts (one bond is type A, six type B, four type C, and four type D). Therefore, there are fifteen potential PCαCβ,γ (or *PCαCβ,γ if more stable) structures to be observed. The other parameters, d3, θ, and φ, account for the position of bromine through the reaction coordinate. The bond length d3 between bromine and the reactive carbon atom of dmbmis ca.
2.0 Å for the initial structure ICα,γ; and it is elongated by 0.7-0.9 Å while d2 is simultaneously shortened to reach the transition state structure in a SN2-like process. In the final product, bromine is completely dissociated as bromide, d3→∞. Moreover, the reactive carbon atom of dmbmhas a tetrahedral arrangement in every ICα,γ structure since θ is nearly 109.5º. When TSCα,γ is reached, θ is decreased to 90º as the reaction evolves toward the ring closure. On the other hand, we observe that ICα,γ, ICα,δ, and ICα,ε (the latter if available) formed at the center Cα are structurally very similar and they differ only in their dihedral angle φ. If Cα is a 566 center then bromine can be oriented ±180º with respect to a [6,6] bond γ, and the other [5,6] bonds δ or ε result in φ ≈ ±180º ± 120º ≈ ±60º (see for instance φ for intermediates at C1 in Table 2). In the case of a 666 center, we selected the bond having the label with the smallest number as the reference point (e.g. in Table 2, for IC15,9 φ ≈ ±180º and the others, IC15,10 and IC15,11, have φ ≈ ±60º).  Figure 1. ΔG ‡ tol (kcal/mol) is the Gibbs energy barrier in toluene for the ring-closure step; this is estimated as the energy difference between TSCα,γ and ICα,γ if the intermediate is more stable than separated reactants; otherwise it is the energy difference between TSCα,γ and separated reactants. All reaction pathways are of profile type I as explained in Scheme 2. b For the dihedral angle φ, if Cα is a 566 center then bromine can be oriented ±180º with respect to a [6,6] bond γ, and the other [5,6] bonds result in φ ≈ ±180º ± 120º ≈ ±60º. On the other hand, if Cα is a 666 center, we select the bond having the label with the smallest number as the reference point (e.g., for IC15,9 φ ≈ ±180º and the others, IC15,10 and IC15,11, have φ ≈ ±60º). c See Figure S1 in the Supporting Information for the discussion of the transition states for these reaction pathways. ΔGtol ‡ values are estimated from linear transits. d See Figures S2 and S3 in the Supporting Information for a discussion of the movement of the La atom in these structures. e Only in the case of *PC22C22,16 the most stable Bingel-Hirsch adduct is a methanofullerene instead of a fulleroid structure.

Fulleroids vs methanofullerenes
We have studied the relative stability of open-cage/closed-cage (fulleroid/methanofullerene) monoadducts. In all cases in which it was possible to optimize both the closed-and open-cage adducts, the fulleroid structure was found to be more stable than the methanofullerene adduct; except for the case of *PC9C9,g. In general, one may optimize both structures but the fulleroid is usually the most stable (see Table S3 in the Supporting Information for a complete analysis of all optimized closed-and open-cage BH adducts). [68] In the cases of *PC7C10,4, *PC16C19,12, and *PC22C22,16, we were not able to optimize an open-cage adduct and only the methanofullerene BH adducts were located. The reaction barrier for the conversion of the cyclopropanated adduct to the fulleroid structure is found to be around 1 kcal/mol. To illustrate this point, Figure 2 depicts the gas-phase linear transit from methanofullerenes *PC1C2,1 and *PC21C23,o to fulleroids PC1C2,1 and PC21C23,o, respectively. For the former, the fulleroid is the most stable while, in the latter, the cyclopropanated structure is marginally more stable than the analogous fulleroid in the gas phase but not in solution. Solvent effects stabilize open-cage adducts more than closed-cage ones (see Table S3). To discuss the origin of the higher stability of the open-cage adducts we performed an activation strain analysis (ASM) [69][70][71] of the reactivity (also known as the distortion/interaction model) [72][73][74] along the reaction pathways leading to the open-cage PC1C2,19 (see Figure 3). In the ASM analysis, the relative energy with respect to separated reactants, ΔE, along the reaction coordinate 15 is decomposed into the strain ΔEstrain associated with deforming the individual reactants and the actual interaction ΔEint between the deformed reactants:

ΔE = ΔEstrain + ΔEint
The open-cage structure PC1C2,19 is found to be more distorted by 33.2 kcal/mol in terms of ΔEstrain than the *PC1C2,19 closed one. This large structural deformation in the open-cage structure is compensated by the interaction term, ΔEint, which is 43.6 kcal/mol more stabilizing than the interaction in its corresponding closed-cage adduct. The resulting binding energy given as ΔEdef First, the breaking of the attacked C-C bond in the fullerene cage stabilizes the LUMO and destabilizes the HOMO of the carbon cage, making frontier orbitals interactions more favorable. [75] Moreover, in open-cage adducts all carbon atoms keep their sp 2 hybridization forming homoaromatic rings, maintaining their π-delocalization. This π-homoconjugation has been shown to be crucial for the final stabilization of open-cage EMFs BH adducts. [19,46] In addition, we have applied the ASM analysis to other reaction pathways and, in all cases, we arrive at the same conclusion that, regardless the higher structural distortion, the open-cage adduct tends to be more stabilized than the closed-cage one because of the higher interaction between the reactants (see Table S4 in the Supporting Information).

The most favored reaction pathways
We were able to optimize 24 TSCα,γ out of the 26 possible transition states of reaction pathways with profile type I. In Figure S1 in the Supporting Information, linear transit calculations of the two missing TSCα,γ structures show that they are very high in energy. From the experimental evidence provided by Nagase et al., [20,49] five different structures were isolated revealing the following product distribution (in percentage): 55.4, 22.1, 5.5, 5.1, and 11.9. From the reaction energies ΔGtol reported in Table 1, we have determined that the most stable intermediate is generated at the C2 position with the orientational isomer IC2,19 very close in energy to IC23,o by 1.2 kcal/mol. The orientational isomers of IC2 and IC23 are connected to the corresponding addition products via relatively high energy transition states. Indeed, TSC2,1, TSC2,19, TSC23,17,   TSC23,o, and TSC23,p are associated with energy barriers of 13.1, 15.8, 13.5, 21.5 and 15.2 kcal/mol, respectively (the three latter barriers are not shown in Table 2 because they were not classified into profile type I; however, they are reported in Table S5  to the non-IPR Sc2@C66; wherein reaction pathways with the lowest-energy barriers but destabilized intermediates were discarded to produce the experimentally observed BH adduct. [68] Based on these results we conclude that the most favored reaction pathway leads to the formation of a fulleroid derivative on bond 19 rather than on bond 11 (although the formation of this latter cannot be totally discarded). These results are in agreement with the experimental evidence since only one of the five characterized structures corresponds to a BH product (P). In fact, from the 13 C NMR spectrum of anion P, and on the basis of the similar NMR spectra of previously reported cycloadditions leading to open-cage adducts, [19,49] the authors concluded that the observed signals correspond to sp 2 carbon atoms of an open-cage adduct attacked at a [6,6] position. For instance, in the case of the formation of the Y3N@C80 fulleroid, X-ray analysis confirmed the formation of an open-cage monoadduct on a [6,6] bond. [39] Moreover, based on NMR and UV spectra and DFT calculations, a [6,6] bond was suggested to be the attacked site for the formation of Sc3N@C80 BH derivative. [76] We conclude, therefore, that the isolated La@C2v-C82 monoadduct structure corresponds to PC2C3,19, a fulleroid formed on a [6,6] bond.
There is also another competitive reaction pathway leading to *PC22C22,16, a methanofullerene, for which the respective energy barrier is only 0.5 kcal/mol higher than the one for the formation of PC2C3,19 (see Figure 4). Nevertheless, the formation of a BH product on bond 16 is discarded because the obtained closed-cage monoadduct is not in agreement with the experimentally observed fulleroid structure. It is interesting to notice that additions with the lowest energy barriers produce well-stabilized fulleroids. Indeed the most stable fulleroid is PC2C3,19. The great stability of PC2C3,19 is in line with the experimental higher stability of the fulleroid as compared to the singly-bonded derivatives when heated at 80 ºC. [49]  Three [6,6] bonds are determined to be the most reactive for the formation of the final BH adduct: 11 (type B), 16 (type A), and 19 (type B). For 16 and 19, the initial intermediate formation occurs on a 566 carbon atom (C3 and C22), and for 11 at a 666 carbon atom (C15).
The closure and final BH product formation takes place always on a [6,6] bond. The attack to bonds of type C (or pyrene), which only involve 666 carbons, results in relatively high energy barriers. Additionally, all the reaction pathways on pyrene bonds are classified with profile type I. Type D (or corannulene) bonds can be considered relatively unreactive since many reaction pathways happening on corannulene positions exhibit type II and III energy profiles, and those additions with profile I involving type D bonds have energy barriers higher than 10 kcal/mol. In general, we observe that the ring-closure between dmbmand a [6,6] bond during the second step of the BH reaction is favored when a 5-MR is adjacent to the initially functionalized carbon atom, as shown for bond types A and B. Interestingly, bonds 11, 16, and 19 are among the ten most reactive bonds for the Diels-Alder cycloaddition to La@C2v-C82, but they are less reactive than bond o, which is the most reactive in this type of reactions. [18] 18 Based on the ESR, NMR, and X-ray analysis carried out by Nagase et al., [49] four adducts were assigned to be singly-bonded derivatives. Those structures should correspond to the thermodynamically most stable intermediates: 666 positions C2 and C7 and 566 carbons C11, C13, C18, C19, C21, and C23 (see Table 1). The reaction pathways involving the intermediates formed at C23 are classified into profiles II or III: intermediates IC23,17, IC23,o, and IC23,p are at least 3 kcal/mol more stable than the corresponding fulleroids PC23C24,17, PC21C23,o, and PC22C23,p.
In addition, their formation is associated with activation energies twice or three times larger than the lowest ones (see Table S5). Consequently, we can conclude that the formation of the three products involving the formation of a C23 intermediate are kinetically hampered, resulting in an accumulation of the intermediate IC23 as the major product instead of the cyclopropanated derivative. This is in good agreement with the experimental observations since IC23 was the only structure confirmed by X-ray and had the highest yield (55.4%). The other structure with the second highest yield (22.1%) is also a singly-bonded derivative. We assign this structure to be IC2 because it is actually the thermodynamically most stable singly-bonded structure (1.2 kcal/mol below the energy of IC23), and the reaction pathways through the bonds around C2 are associated with activation energies higher than 13 kcal/mol. In view of that, one might expect the same yield for IC2 and IC23. Nonetheless, unlike C23, reaction pathways around C2 do actually produce well-stabilized fulleroids and are classified as type I. In fact the most stable structure is PC2C3,19. Consequently, a portion of IC2 may be consumed to produce PC2C3,19 during the BH reaction, thus being the second most abundant singly-bonded adduct. under Bingel-Hirsch conditions. As a final remark, let us mention that we did not find any correlation between the predicted adducts and the position of the La atom in the La@C2v-C82 cage in terms of Voronoi-deformation-density (VDD) charge distribution [77] and spin density, as reported in Table 3. These results are actually highlighting that these quantities cannot be used to make predictions about possible BH reaction sites. 20

Conclusions
We have quantum chemically identified the five experimentally observed but hitherto uncharacterized products in the Bingel-Hirsch (BH) addition of dimethyl bromomalonate to the La@C2v-C82. Only one of the experimentally observed products is a fulleroid. Our computations indicate that this is the fulleroid at bond 19 which emerges as the kinetically and thermodynamically most favorable Bingel-Hirsch adduct; however, our analysis made through all possible cyclopropanated derivatives suggests that the fulleroid on bond 11 cannot be totally discarded.
Based on experimental results by Nagase et al., [49] four adducts can be assigned to be singlybonded derivatives. We assign two of the observed products with the highest yields to the thermodynamically most stable intermediates, IC2 and IC23. According to our computations, these intermediates are traps on the Bingel-Hirsch reaction pathways. Formation of the final Bingel-Hirsch adduct is kinetically hampered leading to accumulation of these intermediates, especially in the case of IC23. Likewise, the remaining two experimentally reported singly-bonded derivatives with a reduced product distribution are most likely IC19 and IC21. Note, however, that the other well-stabilized intermediates formed at C7, C11, C13, and C18 may be in direct competition.
Finally, our computational exploration also shows that the Bingel-Hirsch addition to La@C2v-C82 is not particularly regioselective. It differs in this respect from the Diels-Alder cycloaddition that occurs exclusively at bond o. Thus, in the case of the Bingel-Hirsch reaction one may anticipate the formation of at least 10 products even under relatively mild reaction conditions. 21 The authors are also grateful for the computer resources, technical expertise, and assistance provided by the HPC centers Barcelona Supercomputing Center -Centro Nacional de Supercomputación, SURFsara in Amsterdam, and the Centre de Serveis Científics i Acadèmics de Catalunya (CESCA). Support for the research of M.S. was received through the ICREA Academia 2014 prize for excellence in research funded by the Generalitat de Catalunya.

Supporting Information
A comparison between transition state geometries optimized in the gas phase and under solvent environment is shown in Tables S1 and S2. Figure S1 contains linear transit calculations for the two missing transition state structures; that is TSC10,i and TSC17,14. Figures S2 and S3 Table S3 the energies of methanofullerenes and fulleroids are compared. Table S4 reports the activation strain model for some selected reaction pathways. In Table S5 the structural parameters and Gibbs energy barriers in toluene are gathered for some selected reaction pathways classified into profile type II or III. Table S6 contains the absolute energies for all the structures under study (from reactants to products). Table S7 contains the Cartesian coordinates for all the 141 DFT-optimized structures located in this work.