Changes in Structure and Reactivity of Ng2 Encapsulated in Fullerenes: A Density Functional Theory Study

Noble gas can be no noble in certain situations from the perspective of structure, bonding, and reactivity. These situations could be extreme experimental conditions or others. In this contribution, we systematically investigate the impact of fullerene encapsulation on molecular structure and chemical reactivity of noble gas dimers (Ng2) in a few fullerene molecules. To that end, we consider He2, Ne2, and Ar2 dimers encapsulated in C50, C60, and C70 fullerenes. We unveil that bond distances of Ng2 inside fullerene become substantially smaller and noble gas atoms become more electrophilic. In return, these noble gas dimers make fullerene molecules more nucleophilic. Using analytical tools from density functional theory, conceptual density functional theory, and information-theoretic approach, we appreciate the nature and origin of these structure and reactivity changes. The results and conclusions from this work should provide more new insights from the viewpoint of changing the perspectives of noble gas reactivity.


INTRODUCTION
Regarded as the most unreactive elements in the periodic table, noble gas atoms could indeed be reactive and no noble at all under certain circumstances, [e.g., high temperature, high pressure, special conditions (e.g., confinement, etc.)]. The first successful experimental realization was the synthesis of xenon hexafluoro platinate XePtF 6 in 1962 (Bartlett, 1962;Graham et al., 2000). Since then many more compounds with noble gas elements included have been reported (Turner and Pimentel, 1963;Nelson and Pimentel, 1967;Bondybey, 1971;Stein, 1973;Holloway and Hope, 1999) and a brand new discipline of noble gas chemistry has emerged thereafter.
In this contribution, we explore the possibility of fullerene encapsulation as another feasible way to make noble gas elements no longer noble. In specific, we consider the impact of the encapsulation by fullerene molecules on geometrical structure and chemical reactivity for noble gas species. Previously, by heating fullerenes at 650 • C under 3,000 atmosphere, helium (Saunders et al., 1993), neon (Saunders et al., 1993), argon (DiCamillo et al., 1996), krypton (Yamamoto et al., 1999), and xenon (Syamala et al., 2002) noble gas atoms were experimentally introduced into the fullerene cage. This experimental realization of fullerene encapsulation demonstrates the feasibility, and there were reports of experimental realization of noble gas dimers encapsulation in fulerenes Laskin et al., 1998;Peres et al., 2001;Popov et al., 2013;Saha et al., 2019). Also, earlier, Krapp and Frenking performed a computational study on Ng 2 @C 60 (Ng = He, Ne, Ar, Kr, and Xe) systems from the bonding perspective and concluded that He 2 @C 60 and Ne 2 @C 60 were weakly bonded van der Waals complexes (Krapp and Frenking, 2007). On the other hand, Solà et al. explored the same systems from the reactivity perspective and they attributed the changes in reactivity to the stabilized LUMO, increased fullerene strain energy, and compressed Ng 2 unit (Osuna et al., 2009). There are other theoretical and computational studies on the encapsulation of either noble gas or other species in the literature (Osuna et al., 2011;Cheng and Sheng, 2013;Dolgonos and Peslherbe, 2014;Khatua et al., 2014;Bickelhaupt et al., 2015;Kryachko et al., 2015;Nikolaienko and Kryachko, 2015;Jalife et al., 2016;Srivastava et al., 2016;Nikolaienko et al., 2017;Foroutan-Nejad et al., 2018;Jaroš et al., 2018;Gómez and Restrepo, 2019;Das et al., 2020).
In this work, we reconsider the likelihood of these species from the computational perspective using a few newly developed theoretical tools. To that end, we examine noble gas dimers, He 2 , Ne 2 , and Ar 2, inside a few fullerene systems, C 50 , C 60 , and C 70 , as illustrative examples in this work. In addition, to pinpoint the nature and origin of these structure and reactivity changes, we make use of a few analytical tools available in the literature, including density functional theory, conceptual density functional theory, informationtheoretic approach, energy decomposition analysis, bonding energy analysis, natural population analysis, and noncovalent interaction analysis. In information theory, Shannon introduced a measure of information content or a lack of it as an entropy through an analogy with the statistical mechanical variant of thermodynamic entropy as proposed by Boltzmann. Hence the stability of a system and spontaneous evolution toward an equilibrium state can be understood by following the variation of the entropy during the thermodynamic process. These analyses should provide better understanding from the viewpoint of changing the perspectives of noble gas reactivity.

Models
Nine model systems were built to examine the impact of fullerene encapsulation on structure and reactivity properties for noble gas dimers (Ng 2 ) using three fullerenes C 50 , C 60 , and C 70 and three Ng 2 species, He 2 , Ne 2 , and Ar 2 . They are represented by following notations: He 2 @C 50 , He 2 @C 60 , He 2 @C 70 , Ne 2 @C 50 , Ne 2 @C 60 , Ne 2 @C 70 , Ar 2 @C 50 , Ar 2 @C 60 , and Ar 2 @C 70 . Also, as the reference, the structural and reactivity properties of these three fullerene molecules and three Ng 2 dimers in vacuum were also simulated and compared.

Analyses
Structure, bonding and reactivity properties for these species were analyzed through a number of well-established analytical tools available in the literature. They include analyses of the total energy decomposition (Parr and Yang, 1989;Liu, 2007), interaction energy decomposition (Bickelhaupt and Baerends, 2000), natural population (Glendening et al., 2012), non-covalent interactions (Johnson et al., 2010), conceptual density functional theory (Geerlings et al., 2003(Geerlings et al., , 2020Liu, 2009Liu, , 2013, and informational-theoretic approach (Liu, 2016;Rong et al., 2019). These analyses enable in-depth understanding about the origin and nature of the changes in structure, bonding, and reactivity. Formulations as well as the details of how to perform these analyses are available elsewhere.

Methodologies
Structure optimizations and natural population analyses were performed with the Gaussian 16 (Frisch et al., 2016) package (versionA03), with ultrafine integration grids and tight selfconsistent field (SCF) convergence. To obtain the most stable orientation of Ng 2 dimer in a fullerene cage, quantum molecular dynamics simulations were performed for each species at the DFT B3LYP/3-21G (Lee et al., 1988;Becke, 1993) level of theory at 3,000 K for 300 ps with a step size of 1 fs. A total of 20 structures with the lowest total energy were selected as the possible structure candidates from the trajectory and then optimized at the DFT M06-2X/6-311G(d) (Zhao and Truhlar, 2008) level of theory. The structure with the lowest total energy from these 20 candidates was selected as the global minimum for the species. The full structure optimization was followed by a single-point frequency calculation to ensure that the final optimized structure has no imaginary frequency. As a measure of quality control, we conducted a benchmark test on the impact of different choices of basis sets and approximate DFT functionals (Ma et al., 2014). The Multiwfn 3.6.1 program (Lu and Chen, 2012) was utilized to calculate the information-theoretic quantities by using the checkpoint file from the Gaussian calculations as the input file. The total energy components were obtained from Gaussian calculations with the keyword of iop(5/33 = 1). The ADF (Amsterdam Density Functional) (Lee et al., 1988) package was employed to perform the bond energy decomposition (BED) analysis. The BED analysis was conducted using the previously optimized structure and the DFT BP86 approximate functional with the double zeta basis set, zero order regular approximation for relativistic correction, and Grimme4 dispersion correction (Te Velde et al., 2001). Table 1 is the benchmark results using Ar 2 @C 60 as an example with 14 basis sets and 14 DFT functionals as well as HF and MP2 methods to examine their impact on the Ar-Ar bond distance. It is compared with Ar-Ar length in vacuum. As can be seen from the table, the bond length in Ar2@C 60 is not significantly dependent on the choice of basis sets and functionals, which is unsurprising because the Ng 2 molecule is confined. As a reference, the experimental Ar-Ar length in vacuum is 3.758 Å (Computational Chemistry Comparison Benchmark DataBase, 2019). There are a handful of computational results for the Ar-Ar length (Colburn and Douglas, 1976). In all analyses below, we chose M06-2X/6-311G(d). Figure 1 shows the final structure obtained for Ng 2 @C 50 , Ng 2 @C 60 , and Ng 2 @C 70 with both top and side views. With a given fullerene, the binding mode is the same for different Ng 2 dimers, but for different fullerenes, as shown in Figure 1, the same Ng 2 dimer prefers to bind differently. This is not surprising either. This is because, for instance, C 70 is elongated so the dimer aligns with the long axis, as one would expect from a simple steric argument. These similarity and difference illustrate the fact that the size and nature of fullerene molecules play more important role during the encapsulation process of Ng 2 dimers. Table 2 exhibits a few selected structural and electronic properties of these species. As can be seen from the table, the When basis sets were tested, M06-2X functional was employed. When functionals were tested, 6-31G(d) basis set was employed. Units in Å.

RESULTS AND DISCUSSION
bond length of Ng 2 dimers becomes smaller when encapsulated in fullerenes, and the smaller the fullerene cage the shorter the Ng-Ng distance. This latter trend was resulted from the fact that and E int are the interaction energy between Ng 2 dimer and the fullerene molecule with and without the BSSE (basis set superposition error) correction considered, respectively; Hirshfeld is the Hirshfeld charge based on the stockholder partition of atoms in molecules; and NPA charge is the electron charged based on the natural population analysis.
FIGURE 1 | Optimized structures of noble gas dimers encapsulated in fullerene cages with Ng = Ar as an illustrative example. The point group for these three encapsulated systems are C 2v , D 3d , and D 5h , respectively. a smaller fullerene has a smaller cage inside so Ng 2 dimers are forced to have shorter bond distance. Energetically, both BSSE (basis set superposition error; Boys and Bernardi, 1970;Simon et al., 1996)-corrected and uncorrected interaction energy results show that for He 2 , all three fullerene molecules yield negative (favorable) interaction energies, whereas for larger Ng 2 dimers, (e.g., Ne 2 and Ar 2 , it appears that only larger-sized fullerenes such as C 70 are able to yield attractive interactions between Ng 2 and fullerene, albeit with the limitations of the method of computation used). This is because larger Ng 2 dimers tend to occupy more space within the fullerene molecule, making it harder for smaller-sized fullerenes to compensate the energetic cost from electrostatic and other repulsions. We are, however, interested in analyzing the qualitative trends only.
Shown in Table 2 are two charge results for Ng atoms encapsulated in fullerene molecules. The first one is the NPA FIGURE 2 | HOMO and LUMO contour surfaces of Ar 2 @C 70 molecule. (natural population analysis) charge from NBO (natural bond orbital) analysis (Glendening et al., 2012), and the other is the Hirshfeld charge based on the shareholder partition of atoms in molecules (Hirshfeld, 1977). NPA charges for Ng atoms could be either positive or negative, but the Hirshfeld charge value is always positive for Ng atoms. According to the reference literature on the Hirshfeld charge (Liu et al., 2014), which is a good descriptor of electrophilicity and nucleophilicity, a positive Hirshfeld charge is an indication that the atom is electrophilic. The larger value of a positive Hirshfeld charge indicates that it is more electrophilic, able to accept more electrons from a nucleophile. Hirshfeld charge results in Table 2 displays that Ng 2 in smaller fullerenes tend to be more electrophilic. Shown in Figure 2 are frontier molecular orbitals from both top and side views for Ar 2 @C 70 as an illustrative example, FIGURE 3 | The mapped contour surface of the molecular electrostatic potential onto the van der Waals surface for Ne 2 and Ar 2 dimers encapsulated in C 60 fullerene. Frontiers in Chemistry | www.frontiersin.org from which one can see that they are both delocalized over the entire molecule. Notice that in both orbitals, we do not see discernible contributions from the Ng 2 dimer overlapping with any of the atomic orbitals from the fullerene, suggesting that atomic orbitals from Ng atoms are not involved in the frontier molecular orbital delocalization. In Table 3, using the information from these frontier orbitals, we calculated a few global reactivity descriptors from conceptual DFT, including chemical potential µ (Parr et al., 1978), chemical hardness η (Parr and Pearson, 1983), and electrophilicity index ω (Parr et al., 1999) for the nine systems studied in this work. A larger value in hardness is usually an indication of being more stable, while a larger value in electrophilicity index suggests that the species is more reactive. From the results in Table 2, we can see that numerical values of these reactivity indices are slightly changed after Ng 2 encapsulation. This should be mainly due to the charge transfer between Ng 2 and fullerene. Ng 2 in smaller fullerenes are more reactive with a larger ω value, whereas in larger fullerenes it becomes often more stable with a larger η value. These results agree well with the conventional chemical wisdom. Figure 3 shows the contour surface of the molecular electrostatic potential mapped onto the van der Waals surface for Ne 2 and Ar 2 dimers encapsulated in C 60 fullerene. It shows that the carbon atoms perpendicular to the Ng-Ng bond is more electronegative. This result is consistent with our Hirshfeld charge result, showing that Ng atoms become more positive and the neighboring carbon atoms in fullerene become more negative. Since Hirshfeld charges are reliable descriptors of electrophilicity and nucleophilicity, Figure 3 also shows that carbon atoms of the fullerene molecule become more nucleophilic, capable of donating more electrons to an electrophile.
In conceptual DFT, Fukui functions (Yang et al., 1984) can also be employed as local descriptors to predict nucleophilic and electrophilic attacks. Shown in Figure 4 are these local functions, where both electrophilic (Figures 4A,B) and nucleophilic (Figures 4C,D) Fukui functions for Ar 2 dimer encapsulated in C 60 fullerene are exhibited. From the figure, we found that (i) carbon atoms in fullerene can be both electrophilic (dark blue areas in Figures 4A,B) and nucleophilic (dark blue regions in  Figure 3 (red areas), suggesting that using Hirshfeld charge and Fukui charge, we obtain different results. According our recent findings, this discrepancy is not unusual. As our recent study demonstrated , results from the Hirshfeld charge should be more robust and reliable. Notice that Krapp and Frenking (2007) also observed that Ng 2 molecules were oxidized inside fullerenes, implying that Ng 2 became more electrophilic in return as electrons were transferred to C 60 , and this also made fullerenes more nucleophilic. This result was confirmed by Solà et al. (Osuna et al., 2009(Osuna et al., , 2011, also who unveiled that the combination of the LUMO stabilization, increased fullerene strain energy, and greater compression of the encapsulated Ng 2 contributed to the reactivity change. To qualitatively identify and distinguish weak interactions between Ng 2 and fullerene molecules, noncovalent interaction (NCI) analysis (Johnson et al., 2010) was employed, whose results are shown in Figure 5. The reduced density gradient is large and positive in the regions far from the molecule, but it becomes small in the vicinity of both covalent and non-covalent interactions.
To identify different types of interactions, the sign of Laplacian of the density is utilized and compared. Plotting low-gradient isosurfaces subject to a further low-density constraint enables real-space visualization of non-covalent interactions. Results in Figure 5 show that when He 2 encapsulated from C 50 to C 70 , a "spike" area becomes more apparent from −0.02 to 0 a.u., indicating that there existed more attractive interactions within the complexes as the fullerene cage becomes larger. The same trend is also observed for Ne 2 and Ar 2 . These results are also in qualitative agreement with the total energy difference result as shown in Table 2. Table 4 is the results from the two schemes (Parr and Yang, 1989;Liu, 2007) of the total energy partition for the nine Ng 2 -fullerene systems studied in this work. Different from the traditional binding energy or interaction energy partition, the total energy partition decomposes the total energy E of a molecular system in terms of a few components that are physiochemically meaningful. In DFT (Parr and Yang, 1989), we have E = T S + E xc + E e , where T S , E xc , and E e are the non-interacting kinetic energy, exchange-correlation energy, and electrostatic interaction energy, respectively. An alternative scheme was proposed by one of us (Liu, 2007), where E = E S + E e + E q , where E S , E e , and E q stand for the energetic contribution from three physiochemical effects, steric, electrostatic, and quantum (due to the exchange-correlation effects). There is a common component in these two partition schemes, the electrostatic energy E e . Considering the energy difference E between two chemical processes sharing the same molecular fragments, we have E = T S + E xc + E e , and E = E S + E e + E q . Table 4 shows the values of these five components from the above two partition schemes using Ng 2 and fullerene in vacuum as the reference. From the table, we found that in the first scheme (Parr and Yang, 1989;Liu, 2007), both T s and E xc contribute positively to E < 0, and it is the electrostatic contribution E e that contributes negatively to their stability, in order to make E < 0. Since E e itself is negative in value (Liu, 2007), this result suggests that to put Ng 2 dimer in the fullerene cage, one has to overcome the large electrostatic repulsion between Ng 2 and fullerene. This result is further verified by the second total energy partition scheme (Liu, 2007), where we found that in this case both electrostatic and quantum (due to exchange-correlation effects) contribute negatively to E < 0 but the steric contribution E S is positive. This result shows that after Ng 2 dimer is put into the cage, smaller space is occupied, but to do so, large electrostatic and Fermionic quantum repulsions have to be overcome. Notice that this total energy difference is different from the BSSE corrected interaction energy in Table 2 and our present results are consistent with what we observed elsewhere for other systems. Figure 6 illustrates the total energy partition analysis to elucidate which energy component is dominant over others in contributing to the total energy difference. A strong linear correlation of the total energy difference E with the electrostatic Frontiers in Chemistry | www.frontiersin.org energy difference E e is shown in Figure 6A, suggesting that it is the electrostatic repulsion that dominates the energetic contribution, which agrees well with the numerical results in Table 4. If the strategy of two-variable fitting is employed, shown in Figures 6B,C, we found that E e combining with either E xc and E S yields even stronger correlations, indicating that contributions from the exchange-correlation and steric effects are minor yet indispensable.
To confirm the results in Tables 4, 5 shows the total bonding energy E bond analysis (Bickelhaupt and Baerends, 2000) results using the ADF package. In this analysis, E bond consists of contributions from the Pauli repulsion E Pauli , steric interaction E steric , electrostatic attraction E elstat , orbital interactions E orb , and dispersion Pauli repulsion E disp , E bond = E Pauli + E steric + E orb + E disp + E elstat . From the results in Table 5, it can be seen that the steric effect (even though its definition is substantially different) contributes favorably and the major opposite contribution is from the electrostatic interaction, agreeing well with our results in Table 4.
Lastly, we examine the results of eight quantities from the information-theoretic approach, which are shown in Table 6, including Shannon entropy (Shannon, 1948), Fisher information (Fisher, 1925), Ghosh-Berkowitz-Parr entropy (Ghosh et al., 1984), information gain (Kullback and Leibler, 1951), 2nd and 3rd orders of absolute Rényi and relative Rényi entropy (Rényi, 1970;Nagy and Romera, 2015). These quantities each have their own physio-chemical meaning and thus representing different   yet intrinsic properties of molecular systems. For example, Shannon entropy is a measure of the electron density distribution uniformity, Fisher information gauges the heterogeneity of the same electron density distribution, and information gain is a robust measurement of regioselectivity, electrophilicity, and nucleophilicity, with its first-order approximation yielding the Hirshfeld charge (Liu et al., 2014). To make sense of these ITA results, Table 7 lists the correlation coefficient of these quantities with respect to either Ng 2 or fullerene types, and Figure 7 displays illustrative example of strong linear relationships of the total energy difference E with Shannon entropy, Fisher information, and information gain for He 2 , Ne 2 , and Ar 2 encapsulated in three fullerenes. Strong correlations are seen in most cases, as shown in the figure, but as highlighted in Table 7, these strong correlations are true only for Ng 2 dimers cross different fullerenes. If one fits a fullerene with three Ng 2 dimers, no such strong correlation is often observed. Also, we tried to put all data points in one plot, and no significant correlation was observed either. Strong correlations in Figure 7 show that ITA quantities could be good descriptors of change trends for Ng 2 dimers encapsulated in fullerene cages. For instance, for He 2 encapsulated in from C 50 to C 70 , binding energy increases, so does the Shannon entropy, due to the more delocalization of the electron density distribution. Another example is the information gain, which is a direct measure of the Hirshfeld charge. Again, for He 2 dimer, from C 50 to C 70 , the difference in information gain decreases, leading to the decrease in Hirshfeld charge, which is in good agreement with the result in Table 2. We noticed the abnormality of Shannon entropy for Ar 2 dimer in Figure 7, whose trend is opposite to that of He 2 and Ne 2 . The reason is unknown. More studies are in need.

CONCLUSIONS
To summarize, in this work, we investigated the possibility of making noble gas reactive through the means of fullerene encapsulation. To that end, we built a total of nine model systems, with three Ng 2 dimers (He 2 , Ne 2 , and Ar 2 ) encapsulated in three fullerene (C 50 , C 60 , and C 70 ) cages. Quantum molecular dynamics simulations were employed to explore the potential energy surface of Ng 2 inside the fullerene molecule. To obtain the lowest energy structure it was further examined through a few well-established analytical tools such as conceptual density functional theory, information-theoretic approach, total energy decomposition, bonding energy decomposition, and others.
Our results show that bond distances of Ng 2 inside fullerene become substantially smaller and noble gas atoms become more electrophilic. In return, these noble gas dimers make fullerene molecules more nucleophilic. Using these analytical tools, we appreciate the nature and origin of these structure and reactivity changes. Our results and conclusions drawn from the present study should provide more understanding and new insights from the viewpoint of changing the perspectives of noble gas reactivity.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation, to any qualified researcher.

AUTHOR CONTRIBUTIONS
ML performed the calculations and analyses and prepared the draft. XH and BW helped in calculations and analyses. DZ and CR managed progression of the project. PC initiated the project and revised the draft. SL worked on the plan, assisted the calculations and analyses, and revised the manuscript. All authors contributed to the article and approved the submitted version.