Electron-beam-promoted fullerene dimerization in nanotubes: insights from DFT computations

Introduction

Transmission electron microscopy (TEM) is a technique that has been used for a long time to provide images of molecules, but also to monitor the reactions triggered by the energy transfer of the electron beam to the atoms that build the molecules. In particular, the advances in TEM as well as in methods to anchor molecules on surfaces like graphene or carbon nanotubes have allowed the scientific community to visualize at atomic resolution the structural changes of molecules in situ by single-molecule atomic-resolution real-time TEM imaging (SMART-TEM) [1-5]. Since the initial discovery [2], many movies have been published that record the dynamic behavior of a wide range of molecules and chemical reactions. One such process was the dimerization of C₆₀ fullerene in a carbon nanotube peapod, i.e., hybrid structures consisting of fullerene molecules encapsulated in single-walled carbon nanotubes (SWCNT) [6-8]. Different stages of dimerization of C₆₀ molecules inside a peapod, a reaction confined within a one-dimensional nanoscale space, have been detected in the last decade [3,4,7,9]. Nakamura and co-workers termed ‘phase 1’ the stage with reversible bond formation and ‘phase 2’ the stage with irreversible C–C fusions [7]. In phase 1, a [2 + 2] cycloadduct C₁₂₀ dimer is formed, which was initially proposed to be with C_s symmetry, in contrast to the X-ray structure for the C₁₂₀ dimer that shows D_2h symmetry [3,9]. In phase 2, irreversible structural rearrangements occur leading to a nanotubular-shaped fullerene inside the peapod. Kinetic analysis with the variable-temperature (VT) SMART-TEM method for the aforementioned C₆₀ dimerization has also been reported by these authors [9]. They concluded that the SWCNT, which accumulates energy by the interaction with the electron beam, activates the reaction either via singlet excitation or via radical cation formation (Scheme 1). Estimation of the activation barrier for the [2 + 2] cycloaddition when the nanotube acts as a sensitizer is 33.5 ± 6.8 kJ mol⁻¹. This value agrees with computational predictions for the reaction via an excited singlet state [10].

Although only a few analyses of the reaction mechanisms have been studied due to the complexity of the system, several intermediates inside the CNT have been proposed [3,11,12], which may be different from those proposed to take place in the gas phase or in the solid state at high pressures and high temperatures [13-15]. We aim to shed light in these reaction mechanisms and energy profiles by using complementary methodologies as standard density functional theory (DFT) calculations and first-principles Car–Parrinello molecular dynamics (CPMD) simulations. Firstly, we have analyzed the interaction between C₆₀ and the nanotube within the peapod. Next, we have found that some dimeric C₆₀–C₆₀ fullerene structures inside the carbon nanotube are thermodynamically favorable. Experiments indicate that, besides C₆₀ sensitization via a singlet excited state, the [2 + 2] cycloaddition can also be activated through the formation of C₆₀^+• radical cation [3,9]. This mechanistic proposal for phase 1, which to our knowledge has not yet been explored in detail inside a carbon nanotube, is analyzed here and compared to the non-activated C₆₀ dimerization. Finally, some intermediates for the subsequent irreversible C–C fusions occurring in phase 2 are proposed with the help of accelerated Car–Parrinello MD simulations.

Results and Discussion

Nanotube-C₆₀ interaction: stabilization of the peapod

First, we estimated the size of the stabilizing interaction that holds the peapod, that is, the interaction between the C₆₀ surface and the walls of the armchair (10,10) CNT. The interaction or encapsulation energy of a fullerene inside the CNT, defined as E_encap = E_Fuller@CNT – E_Fuller – E_CNT, amounts to −3.23 eV for C₆₀ at our present computational settings (PBE/plane waves, see section Computational methods). This significant amount of energy, which could be overestimated, comes from the π–π interactions between the C₆₀ surface and the CNT wall and is modelled in a first approximation using Grimme’s corrections to the dispersion energy [16]. This interaction energy is comparable to those in similar systems with π–π interactions, as for example a “buckyball catcher” complex with C₆₀ or the interlayer interaction between graphene sheets [16]. It is found to be significantly smaller than for C₆₀@C₂₄₀ [17], but larger than for C₆₀@C₅₄₀ and C₆₀@C₉₆₀ nanoonions [18]. The encapsulation energy of two separated C₆₀ molecules essentially doubles that of a single molecule (−6.48 eV, see Figure S1 in Supporting Information File 1). For different C₆₀ dimers (Figure 1), the interaction energies range between −6.45 and −6.50 eV (Figure S1 in Supporting Information File 1), a range which is around 1% of the total encapsulation energy. As a consequence of the effective π−π interactions, the larger the contact between the fullerene dimer surface and the CNT wall the larger the encapsulation energy. This significant interaction is also apparent from an inspection of the electronic structure of the peapod, where we can observe some molecular orbitals with non-spurious contributions from each of the CNT and C₆₀ fragments (Figures S2–S4 in Supporting Information File 1).

Relative stabilities of C₆₀–C₆₀ dimers

We computed the reaction energies for the dimerization of two C₆₀ molecules in the one-dimensional space within the CNT and compared the results with those for the same reaction in the gas phase. We assumed that the cation radical mechanism takes place, that is, the ionized CNT generates a radical cation C₆₀^•+ that reacts with a C₆₀ molecule to yield different C₁₂₀^•+ dimers. The energies for the dimerization of two neutral C₆₀ molecules were also computed for comparison (Table 1). As products, we have considered dimer 1-D_2h and dimer 1-C_s, which are products of reversible [2 + 2] cycloadditions (phase 1) between two [6,6]-bonds in the former case and, between a [6,6]-bond and a [6,5]-bond in the latter (Figure 1). Dimer 1-C_s is at our computational settings (PBE/PW), more than 15 kcal mol⁻¹ higher in energy than dimer 1-D_2h, the one characterized by X-ray crystallography in the solid state, both in the gas phase and inside the CNT. Similar lower stabilities for dimer 1-C_s are also found for the radical cation products, especially in the gas phase (second column in Table 1). The reaction energies are, however, notably more negative for the radical cation products. Besides, radical cations 1-C_s^•+ and 1-D_2h^•+ show comparable energies. Finally, for the nanotubular-shaped C₁₂₀-NT-D_5d isomer, the formation energies with respect to C₆₀ + C₆₀ are up to −300 kcal mol⁻¹, which reflects the high degree of C–C bond reorganization needed to obtain this fullerene isomer that satisfies the so-called isolated pentagon rule (IPR) [19]. We corroborated our results for the gas-phase products using somewhat different computational settings with a non-periodic electronic structure code as ADF, see values in parenthesis in the two first columns of Table 1 (also at PBE level, but using atomic orbitals instead of plane waves as basis sets, see Computational methods). Although reaction energies are predicted to be somewhat more exothermic in most of the cases, relative energies between isomers are very similar, both for the neutral as well as for the radical cation dimers, what confirms the reliability of the computational settings used in the periodic VASP code.

Energy profile for the reversible [2 + 2] cycloaddition in phase 1

We analyzed in detail the energy profile for the first step in the dimerization process, that is, the reversible [2 + 2] cycloaddition to obtain dimers 1-C_s and 1-D_2h. We initially considered dimerization in the gas phase to check the reliability of our methodology (periodic calculations using VASP code) by comparing with other more standard procedures (non-periodic calculations using ADF code; see Computational methods for details) when computing the energy profile for the radical cation C₆₀ + C₆₀⁺. Besides, we computed the profile for the neutral dimerization C₆₀ + C₆₀ to assess the effect ionization produced by the electron beam has in the process. Before the reaction takes place, a stabilizing van der Waals complex between the two C₆₀ molecules was formed (Figure 2), with an interdimer distance around 3 Å (slightly larger for the neutral profile). For dimer 1-C_s, the stabilization of this complex is significantly more important for the radical cation (around 20 kcal mol⁻¹) than for the neutral complex (less than 10 kcal mol⁻¹). To reach intermediate I-1, the singly-bonded dimer, a transition state TS-1 has to be overcome. The barrier for the radical cation is much smaller (6–7 kcal mol⁻¹) than for the neutral dimer (28.8 kcal mol⁻¹), so we confirm that the process is activated for the cation. The interdimer C···C distance is smaller than 2 Å for the radical cation and for the neutral species. Once intermediate I-1 is formed, the interdimer C–C distance is around 1.60–1.70 Å. Formation of the second C–C bond to yield dimer 1-C_s requires to overcome a second transition state TS-2 with energy barriers that range between 10–13 kcal mol⁻¹ from the immediate intermediate depending on the profile. The interdimer C···C distance of the forming bond is slightly larger than 2 Å in all three cases. Finally, formation of dimer 1-C_s is exothermic for the radical cation profile, but moderately endothermic for the neutral dimer, as previously observed. As general trends, we find that the radical cation profiles are qualitatively parallel with some minimal differences, which corroborates the validity of our periodic approximation to study the dimerization of the radical cation inside the CNT (see below). The rate-determining transition state for the radical cation and for the neutral dimer is the same (TS2) as well as the lowest-energy intermediate (van der Waals complex), so the two profiles are not that different, according to the energetic spam model by Kozuch and Shaik [20]. In any case, the mechanism via radical cation is the most favorable, both kinetically and thermodynamically. For dimer 1-D_2h^•+, the energy profile up to I-1 is very similar to that of dimer 1-C_s^•+ (see Supporting Information File 1, Figure S9). However, TS-2 is significantly lower in energy. Therefore, in gas phase, dimer 1-D_2h^•+ is predicted to be the thermodynamic and the kinetic product.

Once our methodology was validated, the energy profile for the formation of dimer 1-C_s inside the CNT was analyzed. In contrast to the gas phase reaction, we now only found the initial van der Waals complex and a single transition state TS for both the neutral and the radical cation profiles (Figure 3). The van der Waals complex shows an interdimer C···C distance of 2.86 Å, slightly shorter than the one observed for the gas phase reaction. The transition state TS corresponds now to the formation of the second interdimer C···C distance, around 2.15 Å (Figure 3), once the first C–C bond is already formed. All the attempts to obtain an intermediate with a single C–C interdimer bond failed. The energy barrier to overcome the TS is predicted to be 29.1 kcal mol⁻¹ for the neutral dimer, but only 2.0 kcal mol⁻¹ for the radical cation. Therefore, the first step of the C₆₀ dimerization is significantly faster for the radical cation than for the neutral system. Formation of dimer 1-C_s@CNT is appreciably endothermic for the neutral profile (>10 kcal mol⁻¹), but fairly exothermic (−25 kcal mol⁻¹, Figure 3) for the radical cation dimer (C₆₀-C₆₀)^•+@CNT. Albeit some differences are present, the general features of the energy profiles inside the CNT are not that different from those in the gas phase; the C₆₀ + C₆₀ dimerization is slightly endothermic with a barrier around 30 kcal mol⁻¹, whereas dimerization for the radical cation is exothermic with a tiny (or almost nil) barrier, which makes this initial step of the reaction drastically faster. Therefore, the main function of the CNT is to constrain the translation of C₆₀ molecules to a one-dimensional space to maximize the rate of collisions between them. Finally, we would like to point out that in the reversible phase 1, radical cation (1-D_2h)^•+@CNT, which shows slightly lower energy than (1-C_s)^•+@CNT (2 kcal mol⁻¹ with our settings) is also predicted to be formed. Free energy difference between the two isomers, estimated in the gas phase, is somewhat reduced when increasing temperature (see Supporting Information File 1, Table S1). Therefore, the molar fraction of (1-C_s)^•+@CNT would increase slightly at higher temperatures.

Mechanistic insights on the irreversible phase 2: first steps

Once phase 1 was analyzed in detail, we studied the irreversible formation of C–C bonds in the interdimer region, known as phase 2 in the dimerization process. Due to the large number of possible pathways to be investigated, we made use of a simulation technique that allowed us to explore the free energy surface in a fast and efficient way. In particular, we performed Car–Parrinello metadynamics simulations by choosing collective variables (CVs) that describe the formation and breaking of C–C bonds in the interdimer region (see Computational methods for more details). It is not the goal of the present work to explain all the steps up to the formation of fullertube C₁₂₀, but the steps that follow the formation of dimers 1-C_s^•+ and 1-D_2h^•+. We first run a metadynamics simulation for dimer 1-C_s^•+ using as CV the coordination number of nine carbon atoms of one C₆₀ molecule (those of contiguous hexagon and pentagon) with respect to nine carbon atoms in the other C₆₀, all of them in the interdimer region (see Supporting Information File 1). Choosing this single CV, we aim to rapidly and efficiently explore the region of the free energy surface that describes the formation of irreversible C–C bonds. Molecular dynamics simulations were done in the radical C₁₂₀ dimer alone; we did not consider the interaction with the CNT. After a 4 ps metadynamics, we observed the sequential formation of C–C bonds between the two C₆₀ cages up to a number of six in structures that we have called dimer 3B-C_s^•+ (three bonds), dimer 4B-C_s^•+ (four bonds), dimer 5B-C_s^•+ (five bonds) and finally dimer HPR-C_s^•+ (six bonds, see Figure 4). In the latter, one hexagon from each C₆₀ face each other with six interdimer bonds formed between the C atoms at the vertexes of each hexagon forming a hexagonal prism (HPR). Several authors have already proposed this HPR structure for the dimerization of two neutral C₆₀ cages [13-15]. We characterized this structure as a minimum of the potential energy surface at 49.2 kcal mol⁻¹ (PBE/TZP) higher than dimer 1-C_s^•+ (Table 2). The C–C distances of the six new bonds are 1.597 Å and the six C–C distances within each hexagon have been elongated from 1.41–1.49 in C₆₀ to 1.55–1.57 Å (Table S2). Dimers 3B-C_s^•+ and 4B-C_s^•+ are also characterized as minima at 18.7 and 34.3 kcal mol⁻¹ with respect to dimer 1-C_s^•+ (Table 2). Structure of dimer 5B-C_s^•+ found in the metadynamics, however, is not a minimum and leads to dimer 3B-C_s^•+ upon geometry optimization. C–C distances of the new bonds are within the range 1.53–1.71 Å. For dimer 1-C_s^•+, the two distances are equivalent (1.587 Å) as they are the six distances for HPR-C_s^•+ (1.597 Å); for dimers 3B-C_s^•+ and 4B-C_s^•+, however, more asymmetry in the distances is found. Similar C–C distances within the interacting hexagons exist, all of them within the range 1.55–1.60 Å (Supporting Information File 1, Table S2). A second metadynamics run for dimer 1-D_2h^•+ was also done using as CV the coordination number of ten carbon atoms of one C₆₀ molecule (two contiguous hexagons) with respect to ten carbon atoms in the other C₆₀, analogous to the first run. After 7 ps, we also observed the sequential formation of C–C bonds between the two cages in structures now called dimer 3B-D_2h^•+ (three bonds), dimer 4B-D_2h^•+ (four bonds), dimer 5B-D_2h^•+ (five bonds) and finally dimer HPR-D_2h^•+ (six bonds, see Figure S10). Dimers 3B-D_2h^•+, 4B-D_2h^•+ and HPR-D_2h^•+ are characterized as minima at 32.7, 39.3 and 63.3 kcal mol⁻¹ with respect to dimer 1-D_2h^•+ (Table 2). The C–C bond distances are very similar to those of the corresponding C_s dimers (Supporting Information File 1, Table S2).

Conclusion

Two different computational methodologies, molecular approach versus periodic solid-state approach provide analogous and consistent results, which gives strong support and reliability to the predictions. In general, the reaction of two C₆₀ molecules to yield different dimers is exothermic and the exothermicity is enhanced when radical cation dimers are formed. The thermodynamics of dimerization inside the nanotube is not that different as found in the gas phase; the nanotube mainly provides a 1D confinement for the reaction to proceed. The barriers for the reversible phase 1 are found to be easily surmountable at ambient temperature for both the forward and the inverse processes, especially for the dimerization of radical cation fullerenes. For the initial steps of irreversible phase 2, up to six C–C bonds in the interdimer region are formed leading to a hexagonal prism shared by the two cages, with significantly larger barriers than in phase 1, and the processes, consequently, slower. Next steps of radical cation dimerization inside the CNT (phase 2), which are under study, combined with new SMART-TEM images and movies will provide us a deeper understanding of fullerene coalescence processes in peapods.

Computational Methods

The Amsterdam Density Functional (ADF) code [21,22] was used for the electronic structure calculations and to optimize reactants, products, intermediates and transition states. The Perdew, Burke and Ernzerhoff (PBE) functional provided the electronic density [23]. Electrons were described with Slater-type basis functions of triple-ζ + polarization quality. We have included scalar relativistic corrections using the zeroth-order regular approximation (ZORA) formalism. Grimme3 BJ-DAMP dispersion corrections have also been included in all calculations [24]. Stationary points were fully characterized by computing the Hessian matrix.

Car–Parrinello molecular dynamics (CPMD) simulations were performed by means of the CPMD program [25,26]. The description of the electronic structure was based on the expansion of the valence electronic wave functions in a plane wave (PW) basis set, which was limited by an energy cutoff of 40 Ry. The interaction between the valence electrons and the ionic cores was treated through the pseudopotential (PP) approximation (Martins–Troullier type) [27]. The PBE functional was selected as the density functional. Dispersion corrections were also considered in the calculations. We used a fictitious electron mass of 800 a.u. The simulations were carried out using periodic boundary conditions in a tetragonal cell (side length of 11 Å and height of 20 Å) and a time step of 0.144 fs. We used the metadynamics technique to analyze the dimerization reaction mechanism [28-30]. The collective variable (CV) considered for the exploration of the free-energy surface was the coordination number of nine C atoms of one C₆₀ (those that are involved in a hexagon fused with a pentagon) with respect to the nine C atoms (hexagon-pentagon) on the other C₆₀ molecule (see Supporting Information File 1 for details). The armchair-type single-wall carbon nanotube (SWCNT) (10,10) with a length of 22.14 Å and a diameter of 13.65 Å was employed in the calculations.

The fullerene dimerization mechanism inside a SWCNT was also studied by density functional theory applied to periodic systems through VASP code [31]. The exchange-correlation functional used was PBE with the zero damping DFT-D3 method of Grimme et al. [16]. Inner electrons were replaced by PAW pseudopotentials [16] while valence electron density was expanded in plane waves with a maximum kinetic energy of 400 eV. The model of SWCNT was 29.61 Å long and with a diameter of 13.65 Å, embedded in a box of 30 × 30 Å of vacuum in the plane perpendicular to the nanotube axis. The gas-phase structures were in a cubic box of 30 Å side. The k-point sampling was performed through the Gamma–Pack scheme with one point. Transition states were obtained through the climbing image version of the nudged elastic band algorithm [32] and dimer method [33]. These structures showed a single imaginary frequency or with some negligible imaginary frequencies (under 10 cm⁻¹). The projected density of states was calculated and own code estimated the bands that belong to fullerene or SWCNT. The excited states of the fullerene dimer inside the SWCNT were studied with a single-point energy calculation of a given band structure where one electron of the highest-occupied band of the fullerene is excited to the lowest-unoccupied band of SWCNT.