Local Energy Decomposition Analysis and Molecular Properties of Encapsulated Methane in Fullerene (CH4@C60)

Methane has been successfully encapsulated within cages of C 60 fullerene, and it is an appropriate model system to study confinement effects. Its chemistry and physics is also relevant for theoretical model descriptions. Here we provided insights into intermolecular interactions and predicted spectroscopic responses of the CH 4 @C 60 complex and compared with results from other methods and with literature data. Local energy decomposition analysis (LED) within the domain-based local pair natural orbital coupled cluster singles, doubles, and perturbative triples (DLPNO-CCSD(T)) framework was used, and an efficient protocol for studies of endohedral complexes of fullerenes is proposed. This approach allowed us to assess energies in relation to electronic and geometric preparation, electrostatics, exchange, and London dispersion for the CH 4 @C 60 endohedral complex. The calculated stabilization energy of CH 4 inside the C 60 fullerene was −13.5 kcal/mol and its magnitude was significantly larger than the latent heat of evaporation of CH 4 . Evaluation of vibrational frequencies and polarizabilities of the CH 4 @C 60 complex revealed that the infrared (IR) and Raman bands of the endohedral CH 4 were essentially “silent” due to dielectric screening effect of the C 60 , which acted as a molecular Faraday cage. Absorption spectra in the UV-Vis domain and ionization potentials of the C 60 and CH 4 @C 60 were predicted. They were almost identical. The calculated 1 H/ 13 C NMR shifts and Methane has been successfully encapsulated within cages of C 60 fullerene, and it is an appropriate model system to study conﬁnement eﬀects. Its chemistry and physics is also relevant for theoretical model descriptions. Here we provided insights into intermolecular interactions and predicted spectroscopic responses of the CH 4 @C 60 complex and compared with results from other methods and with literature data. Local energy decomposition analysis (LED) within the domain-based local pair natural orbital coupled cluster singles, doubles, and perturbative triples (DLPNO-CCSD(T)) framework was used, and an eﬃcient protocol for studies of endohedral complexes of fullerenes is proposed. This approach allowed us to assess energies in relation to electronic and geometric preparation, electrostatics, exchange, and London dispersion for the CH 4 @C 60 endohedral complex. The calculated stabilization energy of CH 4 inside the C 60 fullerene was − 13.5 kcal/mol and its magnitude was signiﬁcantly larger than the latent heat of evaporation of CH 4 . Evaluation of vibrational frequencies and polarizabilities of the CH 4 @C 60 complex revealed that the infrared (IR) and Raman bands of the endohedral CH 4 were essentially “silent” due to dielectric screening eﬀect of the C 60 , which acted as a molecular Faraday cage. Absorption spectra in the UV-Vis domain and ionization potentials of the C 60 and CH 4 @C 60 were predicted. They were almost identical. The calculated 1 H/ 13 C NMR shifts and spin-spin coupling constants were in very good agreement with experimental data. In addition, reference DLPNO-CCSD(T) interaction energies for complexes with noble gases (Ng@C 60 ; Ng = He, Ne, Ar, Kr) were calculated. The values were compared with those derived from supermolecular MP2/SCS-MP2 calculations and estimates with London-type formulas by Pyykkö and coworkers [ Phys. Chem. Chem. Phys. , 2010, 12, 6187 − 6203], and with values derived from DFT-based symmetry-adapted perturbation theory (DFT-SAPT) by Hesselmann & Korona [ Phys. Chem. Chem. Phys. , 2011, 13, 732 − 743]. Selected points at the potential energy surface of the endohedral He 2 @C 60 trimer were considered. In contrast to previous theoretical attempts with the DFT/MP2/SCS-MP2/DFT-SAPT methods, our calculations at the DLPNO-CCSD(T) level of theory predicted the He 2 @C 60 trimer to be thermodynamically stable, which is in agreement with experimental observations.


Introduction
Carbon displays a rich chemistry and physics with a variety of molecular allotropes, including common graphite and diamond, but also fullerenes, carbon nanotubes, and graphene. [1][2][3][4] The most investigated fullerene is the C 60 "Buckminsterfullerene" composed of 20 hexagons and 12 pentagons of sp 2 -hybridized carbon atoms fused into a pseudosphere with a ∼7 Å diameter, as displayed in Figure 1. C 60 occurs in trace amounts on Earth in carbon-rich rocks and soot. 5,6 It has also been observed in the micrometeorite impact residue on the NASA Long Duration Exposure vealed encapsulated atoms of noble gases with a 3 He/ 4 He isotope ratio of clearly extraterrestial origin. 9 Moreover, analyses of the 2019 data collected by the NASA/ESA Hubble Space Telescope confirmed spectral features of the ionized C + 60 species in diffuse interstellar bands making it the largest molecule observed in space and indicating that fullerenes might play an important role in interstellar chemistry. 10,11 The properties and chemistry of C 60 have been studied; for example, the wave-particle duality was experimentally observed for C 60 . 12 Methods for preparation and separation have been established, 13,14 and the possibilities to encapsulate atoms and molecules inside the fullerene cage was recognized soon after its discovery. 15 C 60 endohedral complexes with metal ions, noble gases, H 2 , N 2 , and H 2 O have been prepared by high-energy collisions of ionized fullerene species, harsh conditions of high temperature and pressure, electric arc, or by organic synthesis methods called as molecular surgery. [16][17][18][19] The successful synthesis of the endohedral complex with CH 4 was reported in 2019 by Whitby and coworkers. 20 Methane is the largest, and the first organic molecule encapsulated in the C 60 fullerene, and such a complex denoted as CH 4 @C 60 is the main object of this study.
There is no obvious direct route to measure the stabilization energy in fullerene endohedral complexes and obtain insights into the interaction mechanisms. Experimental observations in this respect have been limited to assessing the efficiency/probability of the given complex to be formed, and the main focus has been on spectroscopic and diffraction studies in relation to unusual physical properties of the encapsulated species. 21 Substantial theoretical efforts have directed to studies of C 60 endohedral complexes and associated intermolecular interactions. Pioneering ab initio studies by Jerzy Cioslowski [22][23][24] at the Hartree-Fock (HF) level of theory were expanded by studies of Bühl et al. 25,26 , Patchkovskii et al. 27 , Darzynkiewicz et al., 28 and Autschbach et al. 29 among others, [30][31][32][33][34][35] where density functional theory (DFT) and secondorder Møller-Plesset perturbation theory (MP2) were used. However, within the supermolecular approach with the interaction energy being the arithmetic relation of related energies (E int = E AB − (E A + E B )), DFT is essentially blind to long-range dispersion. This limitation has typically been addressed by using empirical correction schemes for the dispersion contributions. [36][37][38][39][40] MP2 is the lowest ab initio method that accounts for "real" dispersion effects but it is unbalanced and its performance for weak, noncovalent interactions are modest and system dependent. 41 Symmetry-adapted perturbation theory (SAPT) with monomers description at the DFT level (DFT-SAPT) developed by Krzysztof Szalewicz and coworkers is an alternative approach to account for dispersion contributions. [42][43][44] Within SAPT, the interaction energy is obtained as a sum of physical contributions, free from basis set superposition error (BSSE). 45 Hence, the most reliable stabilization energies for C 60 endohedral complexes so far have been obtained with DFT-SAPT by Hesselmann and Korona. 46,47 In parallel developments, approximate London-type formulas have been derived by Pyykkö and coworkers for the estimation of dispersion interaction in endohedral systems. 48,49 Attractive dispersion interactions between nonpolar species such as C 60 and CH 4 (or noble gases) are purely quantum me-chanical and originate from instantaneous effects of dynamical electron correlation. 50 For systems of chemical interest that can be correctly described by a single reference wave function, the most robust (and still tractable) way of introducing electron correlation is the coupled cluster singles, doubles, and perturbative triples CCSD(T) method. 51 It is the "gold standard" of quantum chemistry. However, a canonical implementation of the CCSD(T) model exhibits a seventh-order scaling with the system size, which results in tremendous computational expenses when considering systems larger than 15−25 atoms. Frank Neese and coworkers have recently developed an efficient implementation of the domain-based local pair natural orbital coupled cluster method (DLPNO-CCSD(T)). [52][53][54] Briefly, in the DLPNO-CCSD(T) approach the correlation energy is expressed as a sum of electron pair correlation energies, which enables the distinction between the "weak pairs" with negligible contribution to the total correlation energy and the "strong pairs" that constitute the dominant, desired part. In this way the "weak pairs' can be treated with a computationally more efficient second-order perturbation theory, whereas only the essential "strong pairs" are subjected to an accurate coupled cluster treatment, which greatly reduces the computational complexity. With appropriately selected pairselection thresholds, this model is capable to recover 99.9 % of the correlation energy of its canonical counterpart. It reproduces the CCSD(T) results within a chemical accuracy at substantially reduced computational efforts. 55,56 This approach extends the possibility of obtaining accurate ab initio energies to systems for which only DFT has been applicable so far. 57,58 Moreover, using a local energy decomposition (LED) protocol allows for a physical meaningful decomposition of the interaction energy within the DLPNO-CCSD(T) framework. 50,59,60 In this study, the goal was to providing an accurate interaction energy decomposition for the CH 4 @C 60 complex and the encapsulation energy barrier using the DLPNO-CCSD(T) method. This approach is not biased by the parametrization inherent to the DFT models, including type of exchange-correlation approximation and dispersion correction scheme. The reference interaction energies for endohedral complexes with noble gases were provided and compared with results by Pyykkö et al. and Hesselmann and Korona. [47][48][49] Methods The counterpoise-corrected interaction energy of molecular fragments X and Y can be expressed as: 59 where the E B A (C) notation denotes energy of fragment A calculated at the energy optimized coordinates of B and using a basis set of system C. The ∆E int term is the "electronic interaction", whereas ∆E geo−prep is the geometric preparation contribution that accounts for the differences between equilibrium molecular ge-ometries of isolated fragments and those in a complex ("deformation energy").
The electronic interaction energy ∆E int is decomposed within the following DLPNO-CCSD(T)-LED decomposition scheme: The interaction energy ∆E int is decomposed into that from the Hartree-Fock level of theory ∆E HF int and the corrections due to inclusion of electron correlation ∆E C int . The latter is decomposed further into the interaction energy contribution at the CCSD level of theory ∆E C−CCSD int and that resulting from the perturbative triple excitations ∆E C−(T) int . The Hartree-Fock interaction energy ∆E HF int is decomposed into the electronic preparation contribution ∆E HF el−prep , which correspond to the energy needed to bring the electronic structures of the isolated fragments into the one optimal for the interaction ("energy investment") and into attractive electrostatic E elstat and exchange E exch contributions. The CCSD correlation interaction energy is partitioned further into the genuine London dispersion interaction energy E C disp and the nondispersive correlation contribution ∆E C non−disp . The latter provides (dynamical) corrections to the Hartree-Fock polarization effects, "dynamic charge polarization". We refer the reader to the original articles for a detailed description of the method and implementation. 50,59,60 Computational details All calculations were performed with the ORCA code 61,62 using a very tight convergence tolerance of 1 × 10 −9 E h . The evaluation of Coulomb and exchange integrals was accelerated with the RIJCOSX approximation 63 with def2/J Coulomb-fitting basis set 64 and tightened grid (GridX5; further increase was verified to have negligible effect). Geometry optimizations were converged to very tight thresholds (VeryTightOpt setting) using the revised PBE "revPBE" exchange-correlation DFT approximation 65,66 together with atom-pairwise dispersion correction based on tight binding partial charges (D4). 40 The polarization-consistent segmented pcseg-1 basis set, 67 and the increased DFT integration grid (Grid5 NoFinalGrid) were used. The choice of the revPBE model was based on its performance in a recent and thorough benchmark study (best among the computationally efficient gradient-corrected GGA functionals). 68 To confirm the global energy minima at the potential energy surfaces, and to evaluate vibrational IR and Raman spectra, Hessians and dipole polarizabilities were calculated. The transition state search involved many consecutive computations of the Hessians towards the firstorder saddle point. Thereby computational efforts were reduced by using the smaller pcseg-0 basis set for atoms of aryl groups for the open-cage models. Cartesian coordinates of the models are provided in the Electronic Supplementary Information (ESI) † . DLPNO-CCSD(T) calculations were performed with a Foster-Boys localization scheme, a full MP2 guess, and employed correlationconsistent cc-pVXZ (X = D, T, Q, 5) orbital basis sets [69][70][71] together with the corresponding cc-pVXZ/C auxiliary basis sets. 72 The chosen PNO truncation thresholds are discussed in the results and discussion section. Computations were performed on a cluster node equipped with two Intel R Xeon R Gold 6126 CPUs (2.6 GHz; 12-core) and 256 GB of RAM.

The choice of PNO truncation thresholds and basis sets
To facilitate accuracy control in a user friendly manner, the authors of the DLPNO-CCSD(T) method have implemented three levels of predefined PNO truncation thresholds. These levels converge towards the method limit at increasing computational cost: LoosePNO, NormalPNO, and TightPNO. 55 The first two offer sufficient accuracy for most applications (<0.5 kcal/mol deviation for the evaluations of total energy with respect to the CCSD(T) reference), 55 but for analysis of weak intermolecular interactions, the TightPNO setting should be used. It ensures that the electron pairs that dominate the interaction are being treated at the coupled cluster level. 59 However, fullerenes are challenging for local coupled cluster methods. The large number of (long-range) π − π interactions in the highly delocalized π-system of fullerenes render calculations with the TightPNO setting and accurate basis sets very demanding. 56,73,74 The DLPNO-CCSD(T) in its current implementation is formally a linear-scaling method when considering the iterative part, but the RI-PNO integral transformations on large systems add substantial prefactors to the total computation times, which in turn limit the feasibility, also due to the memory and disk space requirements. Therefore, based on test calculations for the CH 4 · · · C 6 H 6 dimer ( Figure 2), which is expected to exhibit similar nature of noncovalent interactions to those in the CH 4 @C 60 complex, we used a multilevel approach as proposed by Sparta et al. 57 Within the multilevel DLPNO approach the test CH 4 · · · C 6 H 6 system was divided into CH 4 and C 6 H 6 fragments. By this division, the intrafragment electron pairs with their orbitals entirely localized on one molecular fragment could be separated from those that gave rise to intermolecular interaction. 57 In Table 1, dispersion interaction and time used for the DLPNO-CCSD(T)-LED calculations for the CH 4 · · · C 6 H 6 dimer are presented. We monitored the convergence of the energy component for the dispersion interaction (E C disp ) since it depends solely on the treatment of the electron correlation, and therefore is critically sensitive to the truncation thresholds of the PNO and the (in)completeness of the basis set. By using a TightPNO setting for both intrafragment and interfragment pairs, smooth convergence towards the limit of a complete basis set (CBS) was observed (see Fig. 2). The dispersion interaction energy essentially converged at the TightPNO/cc-pVQZ level. Unfortunately, this setup would involve prohibitive computational costs when applied to endohedral fullerene complexes. Therefore, for routine applications we propose a more tractable multilevel DLPNO scheme. In this scheme, the truncation thresholds for the intrafragment PNO Table 1 Dispersion contribution (E C disp ; kcal/mol) to the interaction energy in the CH 4 · · · C 6 H 6 dimer calculated within the DLPNO-CCSD(T)-LED scheme using different PNO truncation settings a and basis sets, as well as extrapolated to the complete basis set limit (CBS) b Inter CH 4 · · · C 6 H 6 Intra CH LoosePNO: T CutPairs = 10 −3 , T CutDO = 2 × 10 −2 , T CutPNO = 1.00 × 10 −7 , T CutMKN = 10 −3 ; see ref. 55 b Estimated from cc-pVXZ (X = T, Q, 5) results using extrapolation scheme by Helgaker et al. 75 for the guest (in this case CH 4 ) and for a troublesome delocalized π-system (C 6 H 6 , C 60 ) are reduced to the NormalPNO and LoosePNO, respectively, whereas the critical interfragment pairs are subjected to an accurate TightPNO evaluation. Together with the combination of the cc-pVQZ/cc-pVTZ basis sets, the multilevel scheme proposed in the last row in Table 1 offers massive computational savings without compromising accuracy to any significant extent. For a test conducted for the CH 4 · · · C 6 H 6 dimer, this approach recovered >95 % of the dispersion interaction energy when compared to the TightPNO/CBS reference, while being only twice as expensive as the TightPNO/cc-pVDZ calculation. The contribution from weak pairs to the E C disp was less than 4 % throughout the calculations presented in Table 1. Moreover, the related interaction energies compared well with previously reported accurate ab initio calculated energies for the CH 4 · · · C 6 H 6 dimer. When using the proposed multilevel DLPNO-CCSD(T) setup, a total interaction energy ∆E int =−1.32 kcal/mol was calculated, which agreed very well with the CCSD(T)/aug-cc-pVTZ result of −1.39 kcal/mol by Ringer et al. 77 The estimated dispersion contribution of −2.03 kcal/mol at the SAPT2/aug-cc-pVDZ level of theory by Ringer et al. 77 was close to our value of −2.22 kcal/mol. Above indicated that the multilevel DLPNO-CCSD(T) setup tailored for substantially more demanding calculations on endohedral complexes of fullerenes is robust.

DLPNO-CCSD(T)-LED analysis of the CH 4 @C 60
In Table 2 the corresponding bond lengths of energy optimized geometries of the CH 4 , C 60 , and CH 4 @C 60 endohedral complex are shown. As a consequence of I h symmetry the C 60 fullerene molecular structure is defined by the two distinct C−C distances r 1 and r 2 that originate from the bonds between fused pentagons and hexagons (r 1 ) and the shorter ones between two hexagons (r 2 ). The energy optimized model of C 60 exhibited excellent agreement with experimental bond lengths estimates. For r 1 the deviation was <0.005 Å, and r 2 coincided with the empirical C−C distance. This indicated that the revPBE-D4/pcseg-1 level of theory was capable to deliver robust models of fullerene systems. For comparison, reported geometries of C 60 optimized at the (ab initio) HF and MP2 levels of theory have revealed considerable deviations of C−C bond lengths, 33,34 whereas previous tests of different DFT approximations have not included dispersion corrections. 33,34 Noteworthy, at the revPBE-D4/pcseg-1 level of theory the geometries of both CH 4 and C 60 remained essentially unchanged upon formation of the CH 4 @C 60 endohedral complex. We note however that in our DFT optimized model of CH 4 the C−H bond length was overestimated by ∼0.01 Å, and did not cor- respond to an energy minimum at coupled cluster level of theory, for which the equilibrium C−H distance was shorter, and closer to experimental estimate. The encapsulation of CH 4 in C 60 is associated with a tiny shortening (<0.001 Å) of the C−H bond lengths. Therefore, the total DLPNO-CCSD(T) energy of the CH 4 molecule calculated at molecular geometry corresponding to that in the CH 4 @C 60 complex was lower compared to that for the isolated CH 4 . Although this effect was very small and had no implications on the evaluation of the electronic interaction energy ∆E int in the CH 4 @C 60 complex, it would lead to an unphysical lowering of the "deformation energy", ∆E geo−prep term in Equation 1. Therefore, to provide the most realistic values the evaluation of ∆E geo−prep included only the contribution from the deformation of the C 60 cage.
In Figure 3, interaction energy contributions are presented from the DLPNO-CCSD(T)-LED analysis of the CH 4 @C 60 endohedral complex. The total interaction energy for the complex is represented as a sum of seven physical contributions: , and ∆E geo−prep (according to Equations 1 and 2). The large and positive electronic preparation term ∆E HF el−prep = +81.22 kcal/mol is counteracted by attractive contributions due to electrostatics and exchange (E elstat = −39.73 kcal/mol, E exch = −21.00 kcal/mol). However, the summed components of the interaction energy at the Hartree-Fock level (∆E HF el−prep + E elstat + E exch ) resulted in substantially repulsive interaction of ∆E HF int = +20.49 kcal/mol. This value was basically identical to the value calculated by Pyykkö and coworkers at the HF/def2-QZVPP level of theory (the same as the +20.50 kcal/mol). 49 This summation can be regarded as an estimate of the extent of steric repulsion. 48 Noteworthy is that the CH 4 @C 60 is predicted to be unstable also by DFT if empirical dispersion corrections are not used. 30,31 The non-dispersive correction due to electron correlation was small and repulsive (∆E C non−disp = +1.00 kcal/mol). As expected, London dispersion is the dominant intermolecular interaction mechanism. The magnitude of the E C disp term of −29.96 kcal/mol was larger than the substantially repulsive Hartree-Fock interaction and the ∆E C non−disp correction, resulting in an endohedral complex stabilization by −8.47 kcal/mol. Yet, further attractive correction came from the contribution from perturbative triple excitations ∆E C−(T) int that stabilized the complex by an additional estimated contribution of −5.03 kcal/mol. The correction from perturbative triples was important, given that it increased the net binding energy in the complex by nearly 60 % (from −8.47 to −13.50 kcal/mol; see inset in Fig. 3). Therefore, the final electronic interaction energy at the DLPNO-CCSD(T)/cc-pVQZ(CH 4 )/cc-pVTZ(C 60 ) level of theory for the CH 4 @C 60 endohedral complex was −13.50 kcal/mol. The geometry preparation ("deformation energy") term was very small (∆E geo−prep = +0.03 kcal/mol). C 60 is a rigid molecule and encapsulation of the CH 4 guest had a negligible effect on its geometry (see Table 2). Our reference stabilization energy of ∆E = −13.47 kcal/mol was compared with best reported estimates. For the CH 4 @C 60 complex, the most robust results have been reported by Pyykkö and coworkers. 49 In that study interaction energies were obtained with supermolecular MP2 and its spin component scaled counterpart (SCS-MP2). Calculations were performed with the def2-TZVPP and def2-QZVPP basis sets, interaction energies were corrected for basis set superposition error and extrapolated to the complete basis set limit. The obtained MP2 interaction energy of −21.37 kcal/mol was clearly overestimated. The value obtained with SCS-MP2 (−11.97 kcal/mol) was closer to the coupled cluster reference, but underestimated. These calculated interaction energies followed the pattern observed in previously reported benchmark calculations. For the CH 4 · · · C 6 H 6 dimer, MP2 overesti-mates the CCSD(T) interaction energy, as was shown by Ringer et al., 77 and for the endohedral CH 4 @C 60 complex, this overestimation seems to be even more pronounced. The same trend of deviation was observed by Pyykkö and coworkers for dimers composed of atoms of noble gases and benzene (Ng· · · C 6 H 6 ), where MP2 overestimated the CCSD(T) reference interaction energies significantly, whereas SCS-MP2 was generally much closer to coupled cluster results, but consistently underestimated the interaction. 49 Pyykkö and coworkers have also developed London-type formulas to estimate dispersion interaction energies in endohedral systems. 48,49 The input parameters to these formulas such as ionization potentials and polarizabilities can be readily computed at the DFT level. Using the data from the study of Pyykkö and coworkers 49 (Table 17, equations 69+86) the dispersion energy estimate of −17.33 kcal/mol was obtained for the CH 4 @C 60 complex. This energy was much smaller than for the DLPNO-CCSD(T)-LED (−29.96 kcal/mol), and would not overcome the steric repulsion estimate of +20.49 kcal/mol, and the complex would not be estimated to be stabilized in that description. In Figure 4, energies are shown in relation to the energy barrier of CH 4 insertion through the orifice of the open-cage C 60 .
The model of open-fullerene is designed to match the molecule used by Whitby and coworkers 20 in their successful synthesis of the CH 4 @C 60 endohedral complex. The insertion energy barrier was calculated at the DLPNO-CCSD(T)/cc-pVQZ(CH 4 )/cc-pVDZ(open-fullerene) level of theory and was significant (∆E = +19.81 kcal/mol), and in agreement with experimental obser-vations that high pressure and elevated temperature conditions (1645 atm, 190 • C for 22 h) are necessary to achieve a high degree of CH 4 insertion. 20 Noteworthy is that the electronic repulsive interaction at the orifice ∆E int = + 9.19 kcal/mol amounted to only less than half of the insertion energy barrier. The remaining geometry preparation term corresponded to the energy needed to deform the open-cage fullerene from its equilibrium geometry to the one optimal for CH 4 insertion (∆E geo−prep = +10.62 kcal/mol). After insertion, the CH 4 molecule is predicted to be stabilized inside the open-fullerene cage by ∆E = −10.16 kcal/mol.

Spectroscopic properties of the CH 4 @C 60 in IR and UV-Vis
Harmonic vibrational frequencies of the endohedral CH 4 @C 60 complex have been computed at the level of GGA and hybrid DFT approximations without using dispersion corrections, 30,31 and at the Hartree-Fock level of theory. 32 Hence, those frequencies were evaluated on structures corresponding to energy minima in situation where London dispersion interactions had not been accounted for. In those studies the intensities in the resulting calculated IR and Raman spectra were not discussed as well. Therefore, we calculated the harmonic vibrational frequencies together with the respective IR absorption coefficients and Raman scattering factors for the CH 4 , C 60 , and the CH 4 @C 60 complex at the revPBE-D4/pcseg-1 level of theory. Related frequencies, absorption coefficients and scattering factors are presented in Table 3.
The calculated vibrational frequencies of the C 60 were in very good agreement with experimental data and virtually not changed upon CH 4 encapsulation. This situation was in agreement with the experimental IR spectrum of the H 2 O@C 60 , where the vibrational frequencies of the fullerene cage were the same as those of the free C 60 . 19 The calculated IR absorption coefficients and Raman scattering factors (for the fullerene cage) were predicted to be slightly affected by CH 4 encapsulation and resulted on average in a <5 % loss in spectral intensity. Frequencies of the encapsulated CH 4 on the other hand were blue shifted with respect to the free CH 4 molecule, and were in agreement with previous theoretical predictions. 31,32 Our results suggested that the triple degeneracy of the asymmetric bending and stretching IR modes (1287 and 3107 cm −1 ) of CH 4 was partially removed due to the interaction with the cage. However, what was the most important, for both IR and Raman a substantial loss in spectral intensities for the encapsulated CH 4 was revealed. This intensity loss was in line with experimental IR spectra of the CH 4 @openfullerene complex, where vibrations of the CH 4 could not be observed. 35 In addition, vibrational features of the H 2 O were very weak in the experimental IR spectrum of H 2 O@C 60 , and the potential screening effect of the fullerene cage was indicated. 19,83 Dielectric measurements conducted at low temperature and IR spectra of H 2 O@C 60 collected at liquid helium temperature have revealed that the dipole moment of the encapsulated H 2 O was 0.5±0.1 D, 84,85 which is about 25 % of free the H 2 O molecule (in agreement with theoretical predictions 86 ). A very similar extent of dipole moment reduction has been observed for HF in HF@C 60 as well. 87 Table 3 Harmonic vibrational frequencies (ν; cm −1 ), IR absorption coefficients (A; 10 5 cm/mol) a and Raman scattering factors (S; Å 4 /amu) a calculated at the revPBE-D4/pcseg-1 level of theory b ; experimental values for the free CH 4  The fullerene cage protects encapsulated species from influence of the external electric field, and acts as a molecular Faraday cage. 88,89 Such a screening effect can be assessed by calculating a difference between the dipole polarizability (α) of an endohedral complex, and the sum of polarizabilities of an isolated guest and the empty fullerene: 90 Negative ∆α corresponds to the polarizability depression that results from the decrease of polarizability of the encapsulated guest. Therefore, the dielectric screening coefficient can be evaluated: 90 To inspect these effects for the CH 4 @C 60 endohedral complex, the dipole polarizabilities of CH 4 , C 60 , and CH 4 @C 60 were calcu-lated at the CAM-B3LYP/Sadlej-pVTZ level of theory. The rangeseparated and Coulomb-attenuating method called as the CAM-B3LYP DFT approximation 91 was shown to deliver accurate polarizabilities, 92 and a balanced description of electronic excited states. 93 The Sadlej-pVTZ basis set was specifically developed for calculations of polarizabilities and other electric molecular properties. [94][95][96] The calculated polarizabilities and the respective values obtained from Equations 3 and 4 are shown in Table 4. The values for the polarizabilities for CH 4 and C 60 were in good agreement with experimental data. The polarizability of the CH 4 @C 60 complex was predicted to be almost the same as that of the empty C 60 , which in turn was reflected in a substantial polarizability depression of ∆α=−2.48. The value for the dielectric screening coefficient (c = 0.97) indicated a particularly strong effect for the CH 4 @C 60 endohedral complex. The polarizability of the encapsulated CH 4 molecule was essentially quenched. The electronic excited states energies and absorption in the UV-Vis domain were calculated at the CAM-B3LYP/cc-pVTZ level of theory with a time-dependent DFT (TD-DFT) within the efficient sTD-DFT implementation of Bannwarth and Grimme. 99 Transition energies, wavelengths, and oscillator strengths for C 60 and CH 4 @C 60 are shown in Table 5. Table 5 Transition energies (E; eV), corresponding wavelengths (λ ; nm), and oscillator strengths ( f osc ) a calculated at the CAM-B3LYP/cc-pVTZ level of theory The predicted spectroscopic characteristics of C 60 and CH 4 @C 60 were almost identical in the UV-Vis range. Both entities exhibited the same transition energies. The oscillator strengths were marginally lower for the complex compared to the empty fullerene. These results were in agreement with experimental observations for C 60 and H 2 O@C 60 . They display close to identical UV-Vis spectra, despite revealing a slightly lower absorption for the complex. 19 The data in Table 5 compare favorably with the experimental UV-Vis spectrum of C 60 in the gas phase. 100 The set of three most intense bands at 187/231/300 nm that reveal oscillator strengths of 1.907/1.922/0.384 correspond to the experimentally observed transitions at 205/257/330 nm that exhibit extinction coefficients of 4.8/3.5/0.9 (10 5 L mole −1 cm −1 ). 100 Hence, the pattern of UV-Vis bands correlated well between theoretical predictions and the experimentally observed spectra, however, transition energies were overestimated (too short wavelengths) at the TD-DFT level of theory.
Because of the rich chemistry of the ionized fullerene species (C n+ 60 ; n = 1, 2, 3), it is interesting to compare the ionization potentials of the C 60 and its endohedral complex. Ionized fullerenes exhibit a diversity of ionization mechanisms and a variety of reactions with potential implications to chemistry in the interstellar medium. 101 The vertical ionization energies (V IE) for C 60 and CH 4 @C 60 were calculated at the CAM-B3LYP/Sadlej-pVTZ level of theory according to: where E n+ cation and E 0 neutral denote the total energies of the ionized and neutral species, respectively. They were calculated at the equilibrium geometry of the ground state. The calculations for ionized species involved an unrestricted (UDFT) formalism due to higher than singlet multiplicities. The obtained results indicated that the ionization potentials of the CH 4 @C 60 complex were almost identical to those of empty C 60 . The calculated V IE for C 60 agreed very well with experimental data, as can be seen from values in Table 6.

NMR properties of the CH 4 @C 60
Fullerenes constitute the only known allotrope of carbon that can be dissolved in organic solvents at room temperature. 103 Therefore, high resolution liquid-state NMR spectra of C 60 and its endohedral complexes can be measured in common NMR solvents. 104,105 Whitby and coworkers collected and analyzed 1 H and 13 C NMR spectra of the CH 4 @C 60 complex dissolved in 1,2dichlorobenzene. 20 For prediction of such spectra with quantum chemistry methods, a robust model including subtle interactions of the fullerene cage with solvent molecules is important. Hence, we constructed systems composed of the C 60 and CH 4 @C 60 explicitly solvated by 25 molecules of 1,2-dichlorobenzene, whereas solvent effects at the outer sphere were accounted implicitly by a polarizable continuum model (PCM) assuming a dielectric con-stant ε = 9.93. The initial configuration was obtained with the Packmol software, 106 and the coordinates were energy optimized at the revPBE-D4/pcseg-1(CH 4 , C 60 )/pcseg-0(C 6 H 4 Cl 2 ) level of theory up to the energy change of < 5 × 10 −6 E h ; see Figure 5. The 1 H and 13 C nuclear magnetic shielding tensors 107,108 (σ σ σ ; assessed with the GIAO approach 109 ) and indirect nuclear spinspin coupling constants (J) were calculated at the level of PBE0 DFT approximation. 110 The calculations were performed with segmented pcS-1 and pcJ-1 basis sets, which have been specifically developed to provide fast convergence towards the Kohn-Sham limit for NMR shieldings and spin-spin couplings. 111,112 The chosen setup represents a reasonable compromise for the calculations of 1 H/ 13 C NMR parameters when compared to more accurate methods given the size of the system. 113,114 For solvent molecules, the pcseg-0 basis set was used. The calculations of NMR observables involved all electrons (no frozen core) and very tight grids (GridX8 Grid7). The calculated isotropic 1 H/ 13 C NMR shielding: was converted into isotropic NMR chemical shift δ according to: where δ j and σ j correspond to the chemical shift and shielding of an atom of interest j, whereas σ re f , calc CH 4 (gas) and δ re f , exp CH 4 (gas) represent the calculated shielding and experimental shift of the reference, respectively. The CH 4 molecule was used as a reference, since its proton and carbon chemical shifts measured in the gas phase are available. 115 Shifts of chemically equivalent atoms were averaged. Spin-spin coupling constants were represented as a sum of four physical contributions: the Fermi contact (FC), spin-dipole (SD), paramagnetic spin-orbit (PSO), and diamagnetic spin-orbit (DSO) terms.
a From ref. 20,115,116 The calculated 1 H and 13 C NMR chemical shifts as well as the 1 H− 13 C spin-spin couplings presented in Table 7 revealed a very good agreement with experimental data. For protons in the CH 4 molecule, the change in chemical shift (∆δ ) upon encapsulation in C 60 was predicted to be −7.54 ppm, which compared very well to −7.88 ppm observed in an experiment. 20 This change on encapsulation is associated with the NMR shielding inside the fullerene cage, where the locally induced magnetic fields counteract the applied external field of the NMR instrument. The corresponding effect for the 13 C in CH 4 was slightly smaller, as revealed by the calculated and experimental ∆δ values of −5.87 and −4.98 ppm, respectively. For carbon atoms of the fullerene cage, the presence of endohedral CH 4 results in a deshielding of the 13 C NMR signal. Therefore, this effect is opposite to that observed for the 13 CH 4 inside the cage. The calculated 13 C deshielding of the cage of +0.51 ppm was in an excellent agreement with the experimental value of +0.52 ppm. 20 The spin-spin coupling constant 1 J HC in the CH 4 was dominated by the Fermi-contact (FC) mechanism. The calculated coupling strength of 1 J HC = 126.4 Hz for the free CH 4 was close to the experimental value of 125.3 Hz. 116 Encapsulation of CH 4 in C 60 had little effect on the 1 J HC , and the calculated value of 124.1 Hz was very close to the experimentally determined 124.3 Hz. 20 Therefore, our theoretical model predicts correctly the sign and the small magnitude of the ∆J. The small change of the coupling was consistent with negligible deformation of the CH 4 geometry upon encapsulation.
Complexes with noble gases Ng@C 60 (Ng = He, Ne, Ar, Kr) The efficient DLPNO-CCSD(T) setup designed for probing intermolecular interactions in the CH 4 @C 60 complex was also used to obtain reference interaction energies for endohedral complexes with atoms of noble gases (Ng@C 60 ; He, Ne, Ar, Kr). Hence, in this section we compared our coupled cluster results to previously reported estimates at lower levels of theory. Table 8 Interaction energies (∆E int ) calculated at the DLPNO-CCSD(T)/cc-pVQZ(Ng)/cc-pVTZ(C 60 ) level of theory for Ng@C 60  In Table 8, interaction energies (∆E int ) calculated at the DLPNO-CCSD(T)/cc-pVQZ(Ng)/cc-pVTZ(C 60 ) level of theory are presented together with values from the MP2/SCS-MP2 calculations by Pyykkö and coworkers, 49 and from the DFT-SAPT calculations by Hesselmann and Korona. 47 The MP2 method exhibited the most unbalanced performance. The associated interaction energies for He was too low with this method, for Ne they were quite reasonable (accidentally), but those for Ar and Kr were severely overestimated. These trends resembled those observed for values derived with the MP2 method for the CH 4 @C 60 complex. The spin component scaled variant (SCS-MP2) displayed an improved description. Although interaction energies for He and Ne were too low, the result for Ar was close to the coupled cluster reference, and the overestimation for Kr was not as severe as with the MP2. The situation with the DFT-SAPT results was more complex and somewhat difficult to judge. On the one hand, all interaction energies for these systems obtained with the DFT-SAPT were consistently underestimated when compared to the DLPNO-CCSD(T) values. On the other hand, the performance of the DFT-SAPT was consistent and well balanced across the He, Ne, Ar, and Kr complexes. This complexity became apparent when the results of DLPNO-CCSD(T) were plotted against the DFT-SAPT counterpart in Figure 6a. The correlation among the two approaches was very good, despite the consistently underestimated interaction energies by the latter. For comparison, the correlation with the SCS-MP2 data shown in panel b was much worse and significantly less convincing.
In Table 9, estimations of dispersion contributions to the interaction energy for the Ng@C 60 complexes are presented for different methods of calculations. The dispersion contributions (E C disp ) from the DLPNO-CCSD(T)-LED were in good agreements with the dispersion components from the DFT-SAPT for all four considered complexes. Hence, it was concluded that the DFT-SAPT compared favorably to the DLPNO-CCSD(T)-LED for prediction of the dispersion interaction in fullerene complexes. This favorable comparison for DFT-SAPT is illustrated in Figure 7a where results of DLPNO-CCSD(T) and DFT-SAPT are plotted against each other. Therefore, the too low interaction energies obtained with the DFT-SAPT for the complexes with noble gases did not originate from inappropriate description of the dispersive part of the in- , and yet, for small guests (He, Ne) for which repulsive interaction at the Hartree-Fock level (∆E HF el−prep + E elstat + E exch ) is small, the non-dispersive corrections due to electron correlation (∆E C non−disp ) were attractive. The dispersion interaction as predicted with London-type formulas by Pyykkö and coworkers were substantially underestimated when compared to the DLPNO-CCSD(T)-LED results, although a linear trend with the latter was revealed; see Table 9 and Figure 7b. This underestimation indicates that the sum of the dipole-dipole and quadrupole-quadrupole terms in the "Pyykkö model 2010" (equations 69+72 from ref. 49 ) was not sufficient, and that the higher order multipole-multipole contributions are necessary to include to obtain better agreement with high-level quantum chemistry methods. 47,49   The existence of the He 2 @C 60 trimer, where two helium atoms are encapsulated inside the C 60 was discovered with 3 He NMR by Rabinovitz and coworkers. 117 The observed He 2 @C 60 :He@C 60 ratio of 1:200 was 10 times smaller than that for the He 2 @C 70 :He@C 70 (1:20). This reduction suggested that the smaller cavity of C 60 was significantly less suited for the encapsulation of two He atoms as compared to the C 70 fullerene.
The stability of the He 2 @C 60 trimer was studied with quantum chemistry methods by Darzynkiewicz & Scuseria, 28 Krapp & Franking, 118 and Hesselmann & Korona. 47 However, all methods applied, including DFT, MP2, SCS-MP2, and DFT-SAPT predicted repulsive interaction in the range from +1.13 to +10.23 kcal/mol, depending on the method and the basis sets used. Hence, theoretical investigations reported so far suggest that the He 2 @C 60 is thermodynamically unstable towards loss of the noble gas atom, in stark contrast to the experimental observation.
To gain insight into this challenging system and confront discrepancy between theory and experiment with the DLPNO-CCSD(T) approach, a relaxed potential energy surface scan at the revPBE-D4/pcseg-1 level of theory was performed for the He 2 @C 60 . The He−He distance was sampled with a 0.02 Å increment, and for each step all other coordinates were subjected to unconstrained optimization. Subsequently, single-point DLPNO-CCSD(T) calculations were performed on the resulting geometries to locate the "true" energy minimum; see Figure 8.
At the DLPNO-CCSD(T)/cc-pVQZ(He)/cc-pVTZ(C 60 ) level of theory, the equilibrium He−He distance inside the C 60 was 1.94 Å. This distance was not only substantially shorter than that of 3.00 Å calculated for the free He 2 dimer, but it corresponded to a clearly repulsive regime for the latter. Our result was close to the value of 1.95 Å obtained by Krapp & Franking 118 at the DFT BP86/TZVPP level, and slightly shorter than 1.98 Å estimated by Kryachko et al. at the DFT M062X/6-31G(d) level. 119 In agreement with the analysis in the latter work, we observed that the He 2 dimer inside the C 60 was fractionally ionized by +0.0141 |e| (Mulliken charge). The corresponding effect for the case of a single He atom in the He@C 60 complex was +0.0036 |e|; hence, comparably smaller. Note that the interaction well for the isolated He 2 dimer was relatively shallow, whereas the relation of potential energy change upon He−He distance variation was very steep for the He 2 @C 60 trimer. For the equilibrium distance r He−He = 1.94 Å inside C 60 the stabilization energy ∆E int of the He 2 @C 60 trimer: evaluated at the DLPNO-CCSD(T)/cc-pVQZ(He 2 )/cc-pVTZ(C 60 ) level of theory was −1.43 kcal/mol. Noteworthy was that the stabilization energy for the He 2 @C 60 trimer at the DLPNO-CCSD(T) level of theory was almost as high as for the complex with a single He atom with the DFT-SAPT. The ab initio calculations predicted the He 2 @C 60 trimer to be stable, which is in agreement with experimental observations.

Conclusions
The reference interaction energies for endohedral complexes of the C 60 fullerene with He, Ne, Ar, Kr, and CH 4 were calculated at the DLPNO-CCSD(T) level of theory and decomposed into physical contributions with the LED scheme. An accurate and efficient multilevel DLPNO-CCSD(T) setup was proposed, which was applicable to routine studies of endohedral complexes of C 60 and larger fullerenes. Calculated molecular properties of the CH 4 @C 60 complex revealed that the IR and Raman bands of the endohedral CH 4 were essentially "silent" due to the dielectric screening effect of the C 60 , which acted as a molecular Faraday cage. Absorption spectra in the UV-Vis and ionization potentials of C 60 and CH 4 @C 60 were predicted to be almost the same. Calculated 1 H/ 13 C NMR shifts and spin-spin coupling constants were in very good agreement with experimental data. Lastly, selected points at the potential energy surface of the endohedral He 2 @C 60 trimer were calculated at the DLPNO-CCSD(T) level of theory. In contrast to previous theoretical studies with DFT, MP2, SCS-MP2 and DFT-SAPT, where all these methods predicted the He 2 @C 60 to be thermodynamically unstable towards the loss of the noble gas atom, our calculations predicted the He 2 @C 60 to be stable, which is in agreement with experimental observations. Therefore, the case of the He 2 @C 60 trimer clearly indicated that the DLPNO-CCSD(T) level of theory is indispensable in studies of weakly interacting systems, and should be used whenever applicable.

Contents
Cartesian coordinates of the models: