A Simple Model for Relative Energies of All Fullerenes Reveals the

Jan 31, 2019 - Graduate School of Engineering, Nagasaki University , Bunkyo 1-14, Nagasaki-shi , Nagasaki 852-8521 , Japan. ‡ RIKEN Center for ...
0 downloads 0 Views 2MB Size
Article Cite This: J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

pubs.acs.org/JCTC

A Simple Model for Relative Energies of All Fullerenes Reveals the Interplay between Intrinsic Resonance and Structural Deformation Effects in Medium-Sized Fullerenes Bun Chan,*,† Yukio Kawashima,‡ William Dawson,‡ Michio Katouda,‡ Takahito Nakajima,‡ and Kimihiko Hirao‡ †

Graduate School of Engineering, Nagasaki University, Bunkyo 1-14, Nagasaki-shi, Nagasaki 852-8521, Japan RIKEN Center for Computational Science, 7-1-26 Minatojima-minami-machi, Chuo-ku, Kobe, Hyogo 650-0047, Japan



J. Chem. Theory Comput. Downloaded from pubs.acs.org by TULANE UNIV on 02/09/19. For personal use only.

S Supporting Information *

ABSTRACT: Fullerenes are sheets of sp2 carbon atoms wrapped around to form spheres. With this simple consideration, we have in the present study devised and (with over 3600 DFT data points) successfully validated a simple model, termed R+D, for estimating the relative energies of fullerenes. This model contains a resonance component to account for the intrinsic differences between the πenergies of different fullerenes, and a deformation component for treating the distortions from planarity. Notably, we find that both terms (and they alone) are required to obtain good relative energies, which lends support to the formulation of the R+D model. An interesting finding is that for some medium-sized IPR fullerenes, their isomers show similar variations in the two components. We deduce that these fullerenes may represent a good opportunity for tuning molecular properties for practical applications. We hope that the promising results of the present study will encourage further investigations into fullerenes from a fundamental perspective.



INTRODUCTION Fullerenes1,2 are simple molecules. They contain just one single element, namely carbon, with all atoms formally having the same orbital hybridization, namely sp2. Every single fullerene consists of just pentagons and hexagons formed by interconnected sp2 carbon atoms, leading to a single sheet folded into a hollow sphere. Furthermore, all fullerenes have exactly 12 pentagons. From a simple two-dimensional perspective, they differ from one another by the number of atoms and the locations of pentagons on the sheet of carbon, and this leads to different shapes of sphere in three dimension. Given the simple picture of bonding in fullerenes, what do we know about the chemistry of these molecules? From a practical standpoint, much is known about the production,3 manipulation,4,5 safety,6,7 and applications of fullerenes.8−10 However, less is known from a fundamental scientific perspective. For instance, the debate about the formation mechanism of fullerenes is still ongoing despite strong evidence for a dynamic process involving both shrinkage as well as growth of the cage.11−15 Determination of their basic chemical properties is limited, and there is a disparity in the amount of attention given to C60 versus everything else in a universe of fullerenes (see refs 16−19 for some limited examples). In recent years, we embark on an endeavor to uncover more information from this void by obtaining accurate heat of © XXXX American Chemical Society

formation values for C60 and larger fullerenes using computational quantum chemistry.20−22 Our latest investigation focused on selected isomers of C84, specifically the 24 that satisfy the so-called “isolated pentagon rule” (IPR),23 among a total of 51592 general isomers.22 We have also developed qualitative explanations for the observed patterns in heats of formation that we computed in another study.21 With the simplicity of our explanations, we are hopeful that a simple and intuitive model can be developed for rationalizing and perhaps predicting heats of formation with a high level of confidence. With decades of scientific research into fullerenes, it comes as no surprise that such a desirable model has been attempted before, with IPR itself being one of the earliest. Models based on the motifs24−26 and local curvatures27,28 of fullerenes have been proposed; related isodesmic-type reactions have been explored;29 correlations between topological indicators and relative isomer energies have also been examined.30 These previously proposed models appear to be useful for the systems for which they are designed, but it is less clear how well they will perform for fullerenes outside of their intended range. To the best of our knowledge, no unified empirical model has been developed and demonstrated to be predictive for Received: September 29, 2018

A

DOI: 10.1021/acs.jctc.8b00981 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Theory and Computation

Refined single-point energies were computed at the B3LYP/6-31G(d) and DSD-PBE-PBE41/cc-pVTZ levels using geometries as specified in Table 1. For the sake of simplicity, in the following sections, we abbreviate B3-LYP/6-31G(d) as DFT and DSD-PBE-PBE/cc-pVTZ as DH-DFT. The FULLERENE program was used to obtain topological resonance energies (ER)42,43 derived from Hückel theory with the algorithm of Babić et al.,44 and estimated heats of formation according to the motif models of refs 24−26. We used a local implementation of the local-curvature model of ref 27 (written in the Python scripting language) to obtain energies of deformation (from planarity) of the fullerenes (ED). A number of methods were used to estimate the relative energies of fullerenes, these include the “R+D” model devised in the study. The R+D model comprises the resonance energy of the fullerene, and its deformation energy, summed in a 1:1 ratio, i.e., ER+D = ER + ED. The calculations of these components are essentially cost-free when compared with DFT computations. For example, the total time required for B3-LYP/6-31G(d) single-point calculations of all IPR isomers of C96 was about 6.6 million seconds. In comparison, obtaining the corresponding ER+D values consumed less than 2700 s in total, with the cost associated mostly with our use of a Python script rather than a native (compiled) program for obtaining ED . To compare energies of fullerenes with different number of carbon atoms, we normalize all energies on a per-carbon basis. While we have obtained DFT energies with B3-LYP/631G(d), we also employ literature energies30,31 calculated using different methods. To enable a direct comparison of all DFT energies, we obtain “quasi-B3-LYP” values using the energies of Ih C60 as reference with the formula: E(Cn,quasi-B3LYP)/n = E(Cn,method)/n + [E(Ih-C60,B3-LYP) − E(IhC60,method)]/60. The relative energies computed in the present study are vibrationless electronic energies. The investigation of the structural effects on the zero-point vibrational energies and on thermal corrections for finite temperatures is beyond the scope of the present study. We however note that, in our previous studies, we find that it is feasible to obtain adequate values for these quantities from scaled harmonic vibrational frequencies computed with fairly low-level of theories.21,45,46

estimating the relative energies for the entire family of fullerenes. In this age of powerful computing when the calculation of icosahedral (Ih) C6000 with a (hybrid) density functional theory (DFT) method is feasible,31 one might perceive that we no longer need approximate models for evaluating fullerene relative energies. This is still far from being practical reality. A hybrid-DFT single-point energy for I h C 6000 took approximately 24 h to compute with full exploitation of its 120 symmetry operators.31 A generic isomer with low or no symmetry would take much longer. Compounding on this difficulty is the rapidly growing number of fullerene isomers with the size; even a moderately sized fullerene, C120, has over 10000 IPR isomers and nearly 1.7 million general isomers. To screen a large number of moderately sized species using DFT would be inefficient to say the least, and in many cases it would be impractical without a nearly cost-free model. In the present study, we use computational quantum chemistry to obtain accurate relative energies for a wide range of fullerenes covering IPR and non-IPR structures, with sizes on both sides of C60. These benchmark values are then used to guide our development toward a simple and chemically intuitive universal model for fullerene relative energies. We hope that the present study will also facilitate and encourage further studies into the fundamentals of fullerenes, which spawn a very rich range of chemistry despite being simple molecules.



COMPUTATIONAL DETAILS Quantum chemistry computations were carried out with Gaussian 1632 and NTChem 2013.33 Unless otherwise noted, lists of fullerenes, as defined by their ring spiral pentagon positions, were generated using FULLERENE program.34 This program was also used to obtain initial geometries of selected isomers by surface projection.35 The revised Wu force field36,37 was used for subsequent preliminary optimizations. Refined geometries were then obtained with the UFF,38 AM1,39 or B3-LYP40/6-31G(d) methods (Table 1, Results and Discussion). Vibrational frequencies were computed at these levels to verify the optimized structures as minima.



RESULTS AND DISCUSSION General Considerations. The list of fullerenes that we have examined is shown in Table 1. The smallest isomeric fullerenes that we have examined are of C38 and the largest of C200. For each of these fullerenes, we divide its isomers into IPR and non-IPR sets. Within each set, we include all isomers when the total number is smaller than 200, and otherwise we randomly select 200 different isomers. An exception to this practice is C60, for which relative energies have been previously investigated for all 1812 isomers with DFT, and 10 among these with higher-level protocols.30 In this case, we adopt the number of isomers and energies that correspond to an equivalent methodology used in the present study. For icosahedral fullerenes from C60 to C6000, the DFT relative energies per carbon have previously been computed,31 and these are also included in the present study. In general, we obtain the relative energies with two different protocols. The higher-level approach (DFT//AM1) involves DFT [specifically, B3-LYP40/6-31G(d)] energy calculations on geometries optimized with the AM139 semiempirical proce-

Table 1. Fullerenes and Numbers of Isomers Examined in the Present Study, and the Methods Useda DH-DFTb//DFTc fullerene C38 C48 C60 C86 C96 C120 C150 C200 largere

IPR

non-IPR

DFTc//AM1 IPR

17 199 10d 19 187

19 187 200 200 200 9f

non-IPR 17 199 1812d 200 200 200 200

R+D//UFF IPR

19 187 200 200 200 9

non-IPR 17 199 1812 200 200 200 200

a The notation A//B signifies energies obtained using method A on geometries optimized with B. bDSD-PBE-PBE/cc-pVTZ. cB3-LYP/631G(d). dInclude icosahedral C60; energy values obtained from ref 30. e Icosahedral C240, C540, C960, C1500, C2160, C2940, C3840, C4860 and C6000. f Energy values obtained from ref 31.

B

DOI: 10.1021/acs.jctc.8b00981 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Theory and Computation

quantities of both microscropic (Hü ckel β value) and macroscopic (graphene elastic constants) sp2 carbon materials. The potential success of the model would lend further support to the notion of a continuous variation between fullerenes and graphene (and by extension, similar materials such as nanotubes) in terms of relative energies. We note that, the resonance and strain of fullerenes have been considered previously (for example, refs 51−54). Conceptually similar approaches to the R+D model have been proposed for charged fullerenes,55,56 with a distinction being that, in those protocols, “strain” is estimated implicitly without direct treatment of structural distortions. In this regard, a conceivable extension of R+D to charged systems may provide an account of their fundamentals in greater details. Validating the R+D Model for Relative Energies of Fullerenes. Before we inspect the performance of the R+D model, let us first be reassured that the relative energies that we will use as the benchmark are of good accuracy. Figure 1 shows

dure. The lower-level method constitutes the model (termed R +D) that we devise in the present study, applied to structures obtained with the UFF38 molecular mechanics method. For selected isomers of the smaller systems (up to C96), we have also obtained their relative energies computed with a doublehybrid (DH) DFT method (specifically, DSD-PBE-PBE41/ccpVTZ) using DFT-optimized structures. We have in previous studies20,21 found this protocol to be highly accurate for the computations of heats of formation for fullerenes. The R+D Model Based on Resonance and Deformation Energies of Fullerenes. Let us now consider a thought process for building fullerenes. We connect the carbon atoms according to the graph of a fullerene,35,47 on a hypothetical plane where all atoms can attain their ideal bonding without structural distortion. We may reasonably consider that energy variations for different fullerenes are due mainly to distinct resonance energies of different π-bond patterns. These can be estimated from Hückel48 energies at the simplest level. We use a value of −0.111 atomic units (derived from experimental properties of C60) for the β parameter in the Hückel method.34 While we are not aware of a β value determined specifically for graphene, we note that the fullerene value is similar to the estimated ones based on polyenes (−0.096 to −0.124)49 and the free-electron model (−0.112).50 Together with appropriate reference systems, resonance energies (ER) can be obtained from the Hückel energies by ER = EHückel − Eref. In the present study, we use the topological resonance energy in which Eref is the energy for the matching polynomial of the molecular graph. Full details of the underlying theory are given in refs 42−44 and references therein. We now turn the hypothetical system into fullerenes by folding the plane into a sphere. One could argue that the effects of specific structural deformations at the molecular level, be it the strain in σ-system or the weakening of π-bonds due to nonideal orbital overlap, would all be displayed as a whole in the experimentally determined elastic properties of graphene, which are similar to our hypothetical planar fullerenes. Indeed, this is the rationale behind the local curvature model proposed in ref 27, which applies the continuum elasticity theory to graphene by using its flexural rigidity and Poisson ratio. Specifically and briefly, the distortion energy (ED) is obtained by the formula ED = DA∑[2k2 − (1 − α)G] in which D is the flexural rigidity, A is the area per carbon atom, and α is the Poisson ratio. In addition to these three characteristic constants of graphene, the distortion at any particular atomic site is defined by k and G, which are the local mean curvature and local Gaussian curvature, respectively. A thorough discussion of the theoretical background of this approach is provided in ref 27. Herein we simply note that k is given by k = 1/R in which R is the radius of curvature, as defined by the sphere that contains a particular atom and its three nearest neighbors. This detail will be relevant to our later discussion. From a practical point of view, we also note that this model for ED is demonstrated to be applicable to finite systems such as fullerenes and nanotubes, yielding reasonable deformation energies for these systems.27,28 We combine the resonance and deformation energies, and term this protocol “R+D”. We consider a simple 1:1 ratio to be appropriate because the parameters that are involved are all derived experimentally from fullerene and graphene. In other words, it has no fitted parameter from our data set. It uses

Figure 1. DH-DFT [DSD-PBE-PBE/cc-pVTZ] relative energies per carbon atom [kJ mol−1, w.r.t. C96(181)] for C38 to C96 isomers (over 400 data points, Table 1) versus the ones obtained with DFT [B3LYP/6-31G(d)].

a plot of DH-DFT relative energies (on a per-carbon basis) against the corresponding ones obtained with the lower-level DFT method for the (over 400) data points (Table 1) for which DH-DFT energies are available. In this figure, we take the energies for the C96 isomer 181 (numbered according to ref 47) as zero on our relative energy scale because this species is found to have the lowest energy per carbon at the DH-DFT level (as well as with DFT energies). We can see that there is a good correspondence between the two levels of relative energies, not only qualitatively but almost quantitatively. The coefficient of determination (R2) is 0.999, and the mean absolute deviation (MAD) between the two sets of data is 0.9 kJ mol−1 over a range of ∼50 kJ mol−1 of relative energies. Somewhat more noticeable deviations can be seen for some large relative energies but the actual deviations are not excessive (a few kJ mol−1). The results lend us confidence in using the DFT relative energies to assess the performance of the R+D model. A plot of DFT relative energies versus the corresponding ones estimated by the R+D model for all (over 3600) data points is given in Figure 2, with the energies for icosahedral (Ih) C6000 taken as zero. The R2 and MAD between the DFT and R+D values are 0.995 and 1.1 kJ mol−1, respectively, over a range of ∼80 kJ mol−1 of relative energies. If we consider just C

DOI: 10.1021/acs.jctc.8b00981 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Theory and Computation

has requirements for a fullerene with a particular energy and a certain size, there is a reasonable possibility that one will find a suitable candidate. Our R+D model provides a rapid and reliable means for the screening of f ullerenes with certain sizes and thermochemical stabilities. Let us now compare the R+D model with some of the previously proposed models for estimating fullerene energies. Table 2 shows the R2 values between the DFT relative energies Table 2. Correlations (R2) between Relative Energies of Fullerene Isomers Obtained with DFT [B3-LYP/6-31G(d)] and Those Computed with Various Models Figure 2. DFT [B3-LYP/6-31G(d)] relative energies per carbon atom [kJ mol−1, w.r.t. Ih C6000] for all (∼3600) fullerenes examined versus the ones obtained with the R+D model.

the IPR isomers, the corresponding values are 0.995 and 0.7 kJ mol−1 for an ∼40 kJ mol−1 range. The MAD values are associated with both systematic and nonsystematic deviations. Thus, it is also useful to look at the standard deviation (SD) between the deviations, which would give us some insights into the nonsystematic component of the deviations. For the general and IPR sets of data, the SD values are 1.2 and 0.3 kJ mol−1, respectively. The R+D model relies on adequate geometries for estimating the deformation component, and these geometries were obtained with the UFF method in the present study. It is thus instructive to examine how the R+D relative energies fare against those calculated with UFF. The MAD from benchmark values for the full set of per-carbon relative energies for UFF is 9.1 kJ mol −1 , which is substantially larger than the corresponding R+D value of 1.1 kJ mol−1. Beyond examining statistical measures R2, MAD, and SD, it may also be relevant to ask whether the model is capable of, for example, identifying the lowest-energy isomer for a given sized fullerene. We have examined eight sets of fullerene isomers in the present study. In five cases, the R+D model identifies the lowest-energy isomer. In each of the other three cases, it selects an isomer that is within 0.2 kJ mol−1 atom−1 of the lowestenergy one. Is such a difference in an acceptable range for practical applications? Consider, as an example, that for C84, experimentally, isomers up to ∼0.6 kJ mol−1 atom−1 above the lowest-energy one have been observed (see ref 22 and references therein). Thus, we deem a difference of 0.2 kJ mol−1 atom−1 to be within the margin for practical use. While an accuracy of ∼0.2 kJ mol−1 atom−1 for the nearly cost-free R+D model enables efficient preselection of a small proportion of low-energy fullerene isomers from a much larger number, we caution that even a small deviation on a percarbon basis could amount to a sizable deviation for the total energy. Thus, a precise distinction of isomers that are very close in total energies remains a challenge. In this regard, highlevel quantum chemistry computations, such as those employed in refs 20−22, represent a powerful tool for bridging the last mile toward highly accurate theoretical thermochemistry of fullerenes. As an aside, it is noteworthy that a given value of energy percarbon often corresponds to multiple fullerenes with different numbers of carbon atoms. This is by no means a surprising finding given the diversity of fullerenes. However, its implication may be of practical relevance, such that if one

fullerenea

R+D

gen-07b

IPR-00c

IPR-17d

C38 C48 C60 C86 C96 C120 C150 all C86 (IPR) C96 (IPR) C120 (IPR) C150 (IPR) C200 (IPR) Ih (IPR)e all (IPR)

0.933 0.961 0.960 0.987 0.994 0.989 0.988 0.995 0.816 0.890 0.890 0.939 0.958 0.999 0.995

0.844 0.751 0.735 0.788 0.826 0.733 0.667 0.960 0.180 0.089 0.163 0.074 0.140 0.992 0.949

n/a n/a n/a n/a n/a n/a n/a n/a 0.790 0.798 0.774 0.602 0.511 0.439 0.002

n/a n/a n/a n/a n/a n/a n/a n/a 0.836 0.686 0.682 0.570 0.533 0.999 0.986

a

Sets contain all of the examined isomers (IPR and/or non-IPR), unless otherwise noted; see Table 1. bReference 25. cReference 24. d Equation 15 of ref 26. eIh C60 to Ih C6000.

and those obtained with various other models, which include R +D and the structural motif models of refs 24−26 (according to eq 15 in that paper). In the present study, we will refer to the three motif models as gen-07 (ref 25), IPR-00 (ref 24) and IPR-17 (ref 26) as the first model can be applied to all fullerenes whereas the latter two are only applicable to the IPR ones. In addition to evaluating the accuracy for all fullerenes to which a particular model is applicable, we have also examined the performance for subsets of fullerenes. For each fullerene of a given number of carbon atoms, we have assessed the correlations for two sets of isomers. The general set comprises both IPR and non-IPR isomers, and we have also separately assessed the IPR subset, as these isomers are typically lower in energies and may be more relevant to practical applications. For the R+D model, we can see that the correlations with the benchmark DFT energies are good (R2 ∼ 0.85 or better) for all subsets, which is not surprising given the smooth linear correlation that we see in Figure 2. Notably, it compares favorably with the other three models. As an example, Figure 3 shows a plot of the correlations for IPR C96 that has a considerable number of 187 isomers, that for which the R2 value is just 0.855 and to which all models are applicable. Curiously, the IPR-17 model has a good R2 value of 0.986 when it is applied to all IPR fullerenes examined, but its R2 values for isomers of any given sized fullerene are significantly worse. A plot of its relative energies in comparison with the DFT benchmark values reveals the cause of this observation (Figure 4). The IPR-17 model gives a good correlation for fullerenes with different number of carbon atoms, but the D

DOI: 10.1021/acs.jctc.8b00981 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Theory and Computation

Figure 5. DFT [B3-LYP/6-31G(d)] relative energies per carbon atom [kJ mol−1, w.r.t. Ih C6000) for all (∼3600) fullerenes examined versus the R+D and its component values.

Figure 3. DFT [B3-LYP/6-31G(d)] relative energies per carbon atom [kJ mol−1, w.r.t. C96(181)] for all C96 IPR isomers versus the ones obtained with the various models.

two components. This again shows the potential drawback of masking finer details by using a single overall metric. Table 3 shows the R2 values for the correlation with the DFT relative energies for R+D and its resonance and Table 3. Correlations (R2) between Relative Energies of Fullerene Isomers Obtained with the DFT Protocol [B3LYP/6-31G(d)] and Those Computed with the R+D Model and Its Components

Figure 4. DFT [B3-LYP/6-31G(d)] relative energies per carbon atom [kJ mol−1, w.r.t. Ih C6000] for all IPR isomers versus the ones obtained with the R+D and IPR-17 models.

relative energies for isomers of a particularly sized fullerene aggregate together and show poorer correlations. The good size-dependency appears to be sufficient to yield a good R2 when assessing with all IPR isomers, but this single number masks the finer details provided by each subset of fullerene isomers. In comparison, we can see that the R+D model gives good correlations in general and for all subsets of fullerene isomers. Figure 4 also shows that the R+D model yields across-theboard good quantitative agreements in addition to good qualitative trends. We provide in the Supporting Information additional validations and comparisons with the specific systems given in ref 26. It is noteworthy that, in the gen-07, IPR-00, and IPR-17 protocols, the motif connectivity and structural distortion effects are not explicitly treated by these models. It is quite remarkable that they already provide reasonable relative energies in a limited fashion. Importance of Resonance and Deformation Contributions. The R+D estimates consist of two components, and it is instructive to examine the contributions of the two factors to R+D. The relative energies obtained from the resonance and deformation energy terms are shown in Figure 5 for all (∼3600) fullerenes considered in the present study. We can see that both of these components show reasonable correlations with the DFT relative energies, though quantitative agreement can only be achieved by combining the two. As we shall see, the correlations for certain subsets of fullerenes can be substantially worse when we consider just one of the

fullerenea

R+D

resonance

deformation

C38 C48 C60 C86 C96 C120 C150 all C86 (IPR) C96 (IPR) C120 (IPR) C150 (IPR) C200 (IPR) Ih (IPR)b all (IPR)

0.933 0.961 0.960 0.987 0.994 0.989 0.988 0.995 0.816 0.890 0.890 0.939 0.958 0.999 0.995

0.444 0.480 0.282 0.275 0.463 0.248 0.250 0.958 0.421 0.001 0.019 0.014 0.001 0.990 0.885

0.917 0.961 0.955 0.971 0.976 0.972 0.975 0.994 0.006 0.519 0.660 0.814 0.898 0.999 0.985

a

Sets contain all of the examined isomers (IPR and/or non-IPR), unless otherwise noted; see Table 1. bIh C60 to Ih C6000.

deformation components for all and subsets of fullerenes. While the two components give reasonable R2 values of 0.958 (resonance) and 0.994 (deformation) when all fullerenes are used in the correlation, (significantly) smaller values can be seen for the subsets. In particular, the R2 values obtained using resonance energies are typically smaller than 0.5 for fullerene isomers with the same number of carbons. The corresponding ones for the deformation estimates are usually notably better, often with values that are larger than 0.95. For the IPR fullerene isomers, the correlations are generally poorer when compared with those for the corresponding general isomer sets. In particular, for some medium-sized fullerenes, neither the resonance nor the deformation component correlates well with DFT. For example, for C86, the R2 values computed using the resonance and deformation energies independently are, respectively, 0.366 and 0.020. Likewise, the corresponding R2 values for C96 are 0.001 E

DOI: 10.1021/acs.jctc.8b00981 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Theory and Computation

convergence toward the limit. In this case, the SDs for the deformation energies are notably larger than those for the resonance energies, and they dominate the total R+D values. In comparison, for the IPR isomers, while the SD values for the resonance energies become smaller from 0.6 kJ mol−1 for C86 to 0.1 kJ mol−1 for C200, there is no clear trend for the corresponding SDs for the deformation energies. The absence of a clear trend in the SDs for deformation energies is reflected in the absence of one in the SDs for the overall R+D relative energies. This is reasonable because the deformation component, while not dominating in all cases, contributes at least to similar extents to those from the resonance component for the IPR fullerenes. Why is there not a clear trend for the deformation-energy SDs for these IPR isomers? We have previously shown that the convergence of fullerene heats of formation to the graphene limit is rather slow,21 and we have put forward a hypothesis that this is related to a presumed slow convergence of deformation energies in the system. On the basis of these previous results, we speculate that the lack of observable trend in the deformation-energy SDs for IPR isomers is due to two factors: (1) the range of IPR fullerenes (C86 to C200) is simply not sufficiently large for such a convergence to be observed, and (2) the IPR isomers lack a difference in the number of pentagon−pentagon junctions (none of them have any), which is a major differentiating factor in generic fullerene isomers. We reiterate that, for medium-sized fullerenes such as C86 and C96, both the resonance and deformation energies are important contributors to the variations in the combined R+D values. Fullerenes undergo a wide range of reactions, among which some are controlled by the intrinsic π-energies (for example, pericyclic reactions),57 whereas in other cases bond strains and factors such as steric hindrance also contribute to reaction outcomes (for example, single and multiple atom additions).58−60 The magnetism of fullerenes are related to their π-connectivities,61 whereas their mechanical properties are associated with strain.62 A delicate balance between the resonance and deformation may represent an opportunity for tuning molecular properties for practical applications. A Case Study with C86 and C96. To further examine the interplay between the resonance and distortion effects in medium-sized fullerenes, let us use IPR C86 and IPR C96 for a succinct case study (Table 5). For C86, the isomers with the

(resonance) and 0.482 (deformation). In such cases, both components are essential to adequately describe the differences in energies between fullerene isomers (see for example Figure 6). Nonetheless, we see a tendency for the R2 value for

Figure 6. DFT [B3-LYP/6-31G(d)] relative energies per carbon atom (kJ mol−1) for all C96 IPR isomers versus the ones obtained with the R +D model and its components.

deformation energies to become larger with increasing size of the fullerene, which suggests that this factor dominates the difference in energies for isomers of larger fullerenes. We now turn to another indicator in order to better understand the contributing factors to the variations in the relative energies for fullerene isomers. For a particular fullerene size, we determine the variations in their relative energies using the standard deviation (SD), and these values are shown in Table 4. For the general set, the variations in R+D relative Table 4. Standard Deviations (kJ mol−1) between the perCarbon Relative Energies Obtained Using the R+D Model and Its Components for Various Fullerene Isomers fullerenea

R+D

resonance

deformation

C38 C48 C60 C86 C96 C120 C150 C86 (IPR) C96 (IPR) C120 (IPR) C150 (IPR) C200 (IPR)

3.7 3.4 3.1 2.9 3.5 2.3 1.9 0.5 0.6 0.6 0.7 0.6

0.8 0.7 0.6 0.4 0.6 0.3 0.2 0.6 0.5 0.3 0.2 0.1

3.1 2.9 2.8 2.7 3.2 2.2 1.8 0.6 0.8 0.7 0.7 0.6

Table 5. C86 and C96 IPR Isomersa with the Lowest and Highest R+D Relative Energies, and Those with the Lowest and Highest Relative Energies for Only the Resonance and Only the Deformation Components fullerene isomer

a

Sets contain all of the examined isomers (IPR and/or non-IPR), unless otherwise noted; see Table 1.

C86 C86 C96 C96

energies between isomers tend to become smaller as the size of the fullerene grows. At the two ends of the scale, the SD values for C38 and C150 are 3.8 and 1.9 kJ mol−1, respectively, though the progression is by no means completely monotonic. As fullerenes become larger, the number of hexagons increases while the number of pentagons remains a constant, thus all isomers should then converge to the same graphene limit. At that stage we would expect an SD value of zero. Let us now look at the component values for the general isomer set. For each of the two terms, there is a trending

a

lowest highest lowest highest

R+D

resonance

deformation

17 9 181 33

2 9 33 187

19 2 163 33

Numbered according to ref 47.

lowest and highest DFT relative energies are isomers 17 and 9, respectively, numbered according to ref 47. The R+D model is consistent with DFT in identifying these two at the extremes of the scale. Likewise, for C96, both DFT and R+D find isomer 181 to be the lowest in energy and that 33 to be the highest. What about the isomers that give the lowest and highest relative energies for their individual resonance and deformation F

DOI: 10.1021/acs.jctc.8b00981 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Theory and Computation components? Interestingly, we find that a different isomer of C86, isomer 2, has the lowest energy for resonance but it corresponds to the highest relative energy for the deformation component. We also find this to be the case for C96, with isomer 33 having the lowest relative resonance energy but it is highest in energy for the deformation component. We now examine these isomers [C86(2) and C96(33)] in greater detail to deduce whether such an observation is more likely to be coincidental or there is a straightforward rationalization. The structures of these isomers are shown in Figure 7. We can see that they are characterized by having a

which is a direct measure of deviation from planarity, with R being the radius of a local sphere at a particular carbon atom. Accordingly, a structurally more inhomogeneous system would have an overall larger deformation energy, with benefits of maintaining flatness over a large area overshadowed by the quadratically increased energy penalties in other parts. We note that such a rationalization for the deformation component in the R+D model is in spirit similar to the maximum pentagon separation rule proposed for anionic fullerenes.63 In that study, both anionic fullerenes and the corresponding neutral ones were examined; it has been found that the energies of neutral fullerenes estimated by applying the maximum pentagon separation rule are poorer correlated with their benchmark energies than those for the anions. Overall, for these two medium-sized fullerenes, we again see a delicate balance between the two factors. The lowest energy isomer for C86(17) does not correspond to the one with the most favorable resonance energy (2) nor is it the one with least overall deformation (19). This is also the case for C96, with isomer 33 being the overall most favorable, 163 for resonance and 181 for deformation. At the other end of the scale, the C86 isomer with the highest R+D relative energy (9) is the one with the least favorable resonance energy. However, for C96, the highest-energy isomer is associated with the highest deformation energy. There may be more of such interesting cases for the many medium-sized fullerenes, but an exhaustive investigation is beyond the scope of the present study. Convergence from C60 to C6000. Finally, we examine the convergence of the R+D energies and its components by zooming in on the icosahedral fullerenes from C60 to C6000, by looking at the actual per-atom resonance, deformation, and R +D energies, as opposed to the ones relative to C6000 (Figure 8). We can see that the resonance energies are much larger in

Figure 7. Isomers for IPR C86 and IPR C96 with the lowest relative resonance energies and highest relative deformation energies, viewed from x-, y-, and z-axis.

dish shape, with relatively large areas that are quite flat, linked by sharper edges. This is especially apparent for C96(33). We postulate that the inverse relationship between the two factors in these cases is a consequence of a tendency to minimize the penalties due to structural deformations on their intrinsically favorable resonance energies. We recall that the resonance energy component is for a hypothetical planar fullerene. We can qualitatively deduce that the favorable π interactions in these medium-sized planar species are associated with relatively large areas of hexagons without interruption by pentagons. This is easiest to see when we look at C96(33) from the z-axis. A tendency to sustain these interactions in the actual fullerenes would then lead to the dish-shape structures with less distortion to large areas of uninterrupted hexagons, at the expense of larger distortions to atoms along the edges of the dish. Let us now consider the equation for calculating the deformation energy ED = DA∑[2k2 − (1 − α)G]. It is proportional to the square of the local mean curvature k = 1/R,

Figure 8. Convergence of per-carbon R+D energies (kJ mol−1) and their resonance and deformation components for icosahedral fullerenes from C60 to C6000, both the left and right y-axes span a range of 45 kJ mol−1, with the R+D and resonance energies shown on the left, and the deformation energies on the right.

magnitude than those for deformation energies. As a result, the sign and general magnitudes of the R+D energies resemble those for the resonance energies. However, the overall convergence behavior is determined by the deformation energies, such that the shapes of the R+D and deformation curves are similar. Specifically, we can see a more rapid convergence of the resonance contribution, with the limit of −167 kJ mol−1 being essentially reached at C960. This essentially converged value corresponds to 0.573 β, which is G

DOI: 10.1021/acs.jctc.8b00981 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Theory and Computation

An interesting finding from our analysis is that, for some medium-sized IPR fullerenes (e.g., C86 and C96), their isomers show similar variations in the two constituents of the R+D relative energies. By contrast, for larger fullerenes structural strain alone dominates the relative energies. The medium-sized fullerenes may therefore be more versatile candidates than the larger ones for tuning molecular properties for practical applications. The promising results of the R+D model show that, for fullerenes with simple chemical structures, we can apply simple concepts to estimate their relative thermochemical stability with good confidence. While we shall not speculate that other properties of fullerenes can also be described in such straightforward manners, we believe that this target is worth attempting.

practically the value obtained for graphite (0.575 β).64 In contrast, the deformation energy is still changing at C6000, though only slowly at this point. While the deformation energy has not reached zero at C6000, it has decreased to a fairly small value of 1.1 kJ mol−1 at this point. If we fit a curve of the general form E = aNb + c to each set of data points for the resonance and deformation energies per carbon, where N is the number of carbon atoms, we obtain exponent (b) values of −0.97 and −0.71 for the resonance and deformation curves, respectively. Let us compare the convergence behavior of the resonance energies for these Ih fullerenes to those for other conjugated systems.65 For polyenes and polyacenes, the b values are −0.98 and −0.99, respectively, which are similar to the value for the hypothetical planar fullerenes. For the variation in per-carbon deformation energies, if we assume Ih fullerenes to be perfect spheres with all carbon atoms being evenly distributed, we can use eq 1 of ref 27 to derive a simple formula E = 4πD(1 + α)/N. In other words, it corresponds to a power decay with b = −1. The fact that the actual exponent of −0.71 deviates substantially from this value can be attributed to the nonspherical shape of Ih fullerenes, where deformations occur locally at pentagon vertexes and the connecting hexagon-chain edges. We emphasize that, in order to quantitatively reach an accurate graphene limit (−37 kJ mol−1 from C60),21 both the resonance (−6.3) and the deformation (−30.5) components are essential. While deformation is not the major contributor to the total R+D energy in absolute terms, it accounts for a large portion of the energy differences. Even though the variation in resonance energies is small compare to that for deformation energies, there is still a considerable intrinsic difference between the π-system for graphene, which contains only hexagons, and that for C60, in which every single atom is part of a pentagon. The R+D model incorporates the two factors, intrinsic resonance energies and deformation f rom planarity, that together quantif y most of the thermochemical stability of f ullerenes.



ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jctc.8b00981. Additional comparisons between benchmark relative energies and those for more approximate methods (PDF) Relevant raw energies (XLSX)



AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. ORCID

Bun Chan: 0000-0002-0082-5497 Takahito Nakajima: 0000-0002-0229-3666 Kimihiko Hirao: 0000-0002-9545-7569 Funding

We gratefully acknowledge financial support from Japan Society for the Promotion of Science (Grants 15K17816, 16H07074, and 17H05169), Ministry of Education, Culture, Sports, Science and Technology, Japan, and computer time from RIKEN Advanced Center for Computing and Communication, Japan (Project Q18266), High Performance Computing Infrastructure (Project HP170163), Japan, and Research Center for Computational Science, the Institute for Molecular Science, Japan.



CONCLUDING REMARKS Fullerenes consist of sheets of sp2 carbon atoms wrapped around to form spheres. The simplicity in their molecular description enables us to use an equally simple model to account for the relative energies of a universe of different fullerenes. Our R+D model comprises just two nonempirical constituents that together describe the energies of fullerenes. The resonance component (derived from Hückel energies) accounts for the intrinsic differences between the π-energies of different fullerenes associated with their distinct atomic connectivity patterns. The deformation component (obtained with local curvatures derived from continuum elasticity theory) takes the distortions from planarity into account. We have assessed the R+D model using validated DFT energies for a variety of fullerenes of different sizes and isomeric structures, covering both the IPR and non-IPR types. The results show that the model reproduces relative energies for all fullerenes examined with good accuracy. Thus, the R+D model enables rapid and reliable screening for pinpointing fullerenes with certain sizes and thermochemical stabilities. Importantly, both the resonance and the deformation components are required in order to obtain good correlations and numerical agreements between R+D and DFT relative energies. This supports the fundamentals that underpins the formulation of the model.

Notes

The authors declare no competing financial interest.



REFERENCES

(1) Kroto, H. W.; Fischer, J. E.; Cox, D. E. The Fullerenes; Pergamon Press: Oxford, 2012. (2) Tománek, D. Guide through the Nanocarbon Jungle: Buckyballs, Nanotubes, Graphene, and Beyond; Morgan & Claypool Publishers: San Rafael, 2014. (3) Osawa, E. Perspectives of Fullerene Nanotechnology; Kluwer Academic Publishers: Dordrecht, 2002. (4) Taylor, R. Lecture Notes on Fullerene Chemistry: A Handbook for Chemists; Imperial College Press: London, 1999. (5) Yamada, M.; Nakahodo, T.; Wakahara, T.; Tsuchiya, T.; Maeda, Y.; Akasaka, T.; Kako, M.; Yoza, K.; Horn, E.; Mizorogi, N.; Kobayashi, K.; Nagase, S. Positional Control of Encapsulated Atoms Inside a Fullerene Cage by Exohedral Addition. J. Am. Chem. Soc. 2002, 127, 14570−14571. (6) Jia, G.; Wang, H.; Yan, L.; Wang, X.; Pei, R.; Yan, T.; Zhao, Y.; Guo, X. Cytotoxicity of Carbon Nanomaterials: Single-Wall Nanotube, Multi-Wall Nanotube, and Fullerene. Environ. Sci. Technol. 2005, 39, 1378−1383. H

DOI: 10.1021/acs.jctc.8b00981 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Theory and Computation

(27) Guan, J.; Jin, Z.; Zhu, Z.; Chuang, C.; Jin, B.-Y.; Tománek, D. Local Curvature and Stability of Two-Dimensional Systems. Phys. Rev. B: Condens. Matter Mater. Phys. 2014, 90, 245403−1−6. (28) Chuang, C.; Guan, J.; Witalka, D.; Zhu, Z.; Jin, B.-Y.; Tománek, D. Relative Stability and Local Curvature Analysis in Carbon Nanotori. Phys. Rev. B: Condens. Matter Mater. Phys. 2015, 91, 165433−1−11. (29) Sun, C. H.; Lu, G. Q.; Cheng, H. M. Comparative Studies of Hyperhomodesmotic Reactions for the Calculation of Standard Heats of Formation of Fullerenes. J. Phys. Chem. A 2006, 110, 299−302. (30) Sure, R.; Hansen, A.; Schwerdtfeger, P.; Grimme, S. Comprehensive Theoretical Study of All 1812 C60 Isomers. Phys. Chem. Chem. Phys. 2017, 19, 14296−14305. (31) Noël, Y.; De La Pierre, M.; Zicovich-Wilson, C. M.; Orlando, R.; Dovesi, R. Structural, Electronic and Energetic Properties of Giant Icosahedral Fullerenes Up to C6000: Insights from an Ab Initio Hybrid DFT Study. Phys. Chem. Chem. Phys. 2014, 16, 13390−13401. (32) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Scalmani, G.; Barone, V.; Petersson, G. A.; Nakatsuji, H.; Li, X.; Caricato, M.; Marenich, A. V.; Bloino, J.; Janesko, B. G.; Gomperts, R.; Mennucci, B.; Hratchian, H. P.; Ortiz, J. V.; Izmaylov, A. F.; Sonnenberg, J. L.; Williams-Young, D.; Ding, F.; Lipparini, F.; Egidi, F.; Goings, J.; Peng, B.; Petrone, A.; Henderson, T.; Ranasinghe, D.; Zakrzewski, V. G.; Gao, J.; Rega, N.; Zheng, G.; Liang, W.; Hada, M.; Ehara, M.; Toyota, K.; Fukuda, R.; Hasegawa, J.; Ishida, M.; Nakajima, T.; Honda, Y.; Kitao, O.; Nakai, H.; Vreven, T.; Throssell, K.; Montgomery, J. A., Jr.; Peralta, J. E.; Ogliaro, F.; Bearpark, M. J.; Heyd, J. J.; Brothers, E. N.; Kudin, K. N.; Staroverov, V. N.; Keith, T. A.; Kobayashi, R.; Normand, J.; Raghavachari, K.; Rendell, A. P.; Burant, J. C.; Iyengar, S. S.; Tomasi, J.; Cossi, M.; Millam, J. M.; Klene, M.; Adamo, C.; Cammi, R.; Ochterski, J. W.; Martin, R. L.; Morokuma, K.; Farkas, O.; Foresman, J. B.; Fox, D. J. Gaussian 16, Revision A.03; Gaussian, Inc.: Wallingford CT, 2016. (33) Nakajima, T.; Katouda, M.; Kamiya, M.; Nakatsuka, Y. NTChem: A High-Performance Software Package for Quantum Molecular Simulation. Int. J. Quantum Chem. 2015, 115, 349−359. (34) Schwerdtfeger, P.; Wirz, L. N.; Avery, J. Fullerene − a Software Package for Constructing and Analyzing Structures of Regular Fullerenes. J. Comput. Chem. 2013, 34, 1508−1526. (35) Schwerdtfeger, P.; Wirz, L. N.; Avery, J. The Topology of Fullerenes. WIREs Comput. Mol. Sci. 2015, 5, 96−145. (36) Wu, Z.; Jelski, D. A.; George, T. F. Vibrational Motions of Buckminsterfullerene. Chem. Phys. Lett. 1987, 137, 291−294. (37) Wirz, L. N.; Tonner, R.; Hermann, A.; Sure, R.; Schwerdtfeger, P. From Small Fullerenes to the Graphene Limit: A Harmonic ForceField Method for Fullerenes and a Comparison to Density Functional Calculations for Goldberg−Coxeter Fullerenes Up to C980. J. Comput. Chem. 2016, 37, 10−17. (38) Rappé, A. K.; Casewit, C. J.; Colwell, K. S.; Goddard, W. A., III; Skiff, W. M. UFF, a Full Periodic-Table Force-Field for Molecular Mechanics and Molecular-Dynamics Simulations. J. Am. Chem. Soc. 1992, 114, 10024−10035. (39) Dewar, M. J. S.; Zoebisch, E. G.; Healy, E. F.; Stewart, J. J. P. AM1: A New General Purpose Quantum Mechanical Molecular Model. J. Am. Chem. Soc. 1985, 107, 3902−3909. (40) Stephens, P. J.; Devlin, F. J.; Chabalowski, C. F.; Frisch, M. J. Ab Initio Calculation of Vibrational Absorption and Circular Dichroism Spectra Using Density Functional Force Fields. J. Phys. Chem. 1994, 98, 11623−11627. (41) Kozuch, S.; Martin, J. M. L. Spin-Component-Scaled Double Hybrids: An Extensive Search for the Best Fifth-Rung Functionals Blending DFT and Perturbation Theory. J. Comput. Chem. 2013, 34, 2327−2344. (42) Aihara, J. A New Definition of Dewar-Type Resonance Energies. J. Am. Chem. Soc. 1976, 98, 2750−2758. (43) Gutman, I.; Milun, M.; Trinajstic, N. Graph Theory and Molecular Orbitals. 19. Nonparametric Resonance Energies of Arbitrary Conjugated Systems. J. Am. Chem. Soc. 1977, 99, 1692− 1704.

(7) Lewinski, N.; Colvin, V.; Drezek, R. Cytotoxicity of Nanoparticles. Small 2008, 4, 26−49. (8) Langa, F.; Nierengarten, J. F. Fullerenes: Principles and Applications; Royal Society of Chemistry: Cambridge, 2007. (9) Cataldo, F.; da Ros, T. Medicinal Chemistry and Pharmacological Potential of Fullerenes and Carbon Nanotubes; Springer Science +Business Media: Berlin, 2008. (10) Morinaka, Y.; Zhang, R.; Sato, S.; Nikawa, H.; Kato, T.; Furukawa, K.; Yamada, M.; Maeda, Y.; Murata, M.; Wakamiya, A.; Nagase, S.; Akasaka, T.; Murata, Y. Fullerene C70 as a “Nano-Flask” that Reveals the Chemical Reactivity of Atomic Nitrogen. Angew. Chem., Int. Ed. 2017, 56, 6488−6491. (11) Irle, S.; Zheng, G.; Wang, Z.; Morokuma, K. The C60 Formation Puzzle “Solved”: QM/MD Simulations Reveal the Shrinking Hot Giant Road of the Dynamic Fullerene Self-Assembly Mechanism. J. Phys. Chem. B 2006, 110, 14531−14545. (12) Saha, B.; Irle, S.; Morokuma, K. Hot Giant Fullerenes can Eject and Capture C2 Molecules: QM/MD Simulations with Constant Density. J. Phys. Chem. C 2011, 115, 22707−22716. (13) Zhang, J.; Bowles, F. L.; Bearden, D. W.; Ray, W. K.; Fuhrer, T.; Ye, Y.; Dixon, C.; Harich, K.; Helm, R. F.; Olmstead, M. M.; Balch, A. L.; Dorn, H. C. A Missing Link in the Transformation from Asymmetric to Symmetric Metallofullerene Cages Implies a TopDown Fullerene Formation Mechanism. Nat. Chem. 2013, 5, 880− 885. (14) Dunk, P. W. Bottom-Up Formation of Endohedral MonoMetallofullerenes is Directed by Charge Transfer. Nat. Commun. 2014, 5, 5844−1−8. (15) Qian, H.-J.; Wang, Y.; Morokuma, K. Quantum Mechanical Simulation Reveals the Role of Cold Helium Atoms and the Coexistence of Bottom-Up and Top-Down Formation Mechanisms of Buckminsterfullerene from Carbon Vapor. Carbon 2017, 114, 635− 641. (16) Boltalina, O. V.; Dashkova, E. V.; Sidorov, L. N. Gibbs Energies of Gas-Phase Electron Transfer Reactions Involving the Larger Fullerene Anions. Chem. Phys. Lett. 1996, 256, 253−260. (17) Boltalina, O. V.; Ioffe, I. N.; Sidorov, L. N.; Seifert, G.; Vietze, K. Ionization Energy of Fullerenes. J. Am. Chem. Soc. 2000, 122, 9745−9749. (18) NIST Chemistry WebBook, NIST Standard Reference Database Number 69; Linstrom, P. J.; Mallard, W. G., Eds; National Institute of Standards and Technology: Gaithersburg, MD, (retrieved 26 June 2018). (19) Katz, E. A. In Nanostructured Materials for Solar Energy Conversion; Soga, T., Ed.; Elsevier: Amsterdam, 2006; pp 361−443. (20) Karton, A.; Chan, B.; Raghavachari, K.; Radom, L. Evaluation of the Heats of Formation of Corannulene and C60 By Means of High-Level Theoretical Procedures. J. Phys. Chem. A 2013, 117, 1834−1842. (21) Chan, B.; Kawashima, Y.; Katouda, M.; Nakajima, T.; Hirao, K. From C60 to Infinity: Large-Scale Quantum Chemistry Calculations of the Heats of Formation of Higher Fullerenes. J. Am. Chem. Soc. 2016, 138, 1420−1429. (22) Waite, S. L.; Chan, B.; Karton, A.; Page, A. J. Accurate Thermochemical and Kinetic Stabilities of C84 Isomers. J. Phys. Chem. A 2018, 122, 4768−4777. (23) Kroto, H. W. The stability of the Fullerenes Cn, with n = 24, 28, 32, 36, 50, 60 and 70. Nature 1987, 329, 529−531. (24) Cioslowski, J.; Rao, N.; Moncrieff, D. Standard Enthalpies of Formation of Fullerenes and Their Dependence on Structural Motifs. J. Am. Chem. Soc. 2000, 122, 8265−8270. (25) Alcamí, M.; Sánchez, G.; Díaz-Tendero, S.; Wang, Y.; Martín, F. Structural Patterns in Fullerenes Showing Adjacent Pentagons: C20 to C72. J. Nanosci. Nanotechnol. 2007, 7, 1329−1338. (26) Wang, Y.; Díaz-Tendero, S.; Alcamí, M.; Martín, F. Generalized Structural Motif Model for Studying the Thermodynamic Stability of Fullerenes: From C60 to Graphene Passing Through Giant Fullerenes. Phys. Chem. Chem. Phys. 2017, 19, 19646−19655. I

DOI: 10.1021/acs.jctc.8b00981 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Theory and Computation (44) Babić, D.; Brinkmann, G.; Dress, A. Topological Resonance Energy of Fullerenes. J. Chem. Inf. Comput. Sci. 1997, 37, 920−923. (45) Chan, B.; Radom, L. Frequency Scale Factors for Some Double-Hybrid Density Functional Theory Procedures: Accurate Thermochemical Components for High-Level Composite Protocols. J. Chem. Theory Comput. 2016, 12, 3774−3780. (46) Chan, B. Use of Low-Cost Quantum Chemistry Procedures for Geometry Optimization and Vibrational Frequency Calculations: Determination of Frequency Scale Factors and Application to Reactions of Large Systems. J. Chem. Theory Comput. 2017, 13, 6052−6060. (47) Fowler, P. W.; Manolopoulos, D. E. An Atlas of Fullerenes; Dover Publications: New York, 2006. (48) Coulson, C. A.; O’Leary, B.; Mallion, R. B. Hückel Theory for Organic Chemists; Academic Press: London, 1978. (49) Bahnick, D. A. Use of Hückel Molecular Orbital Theory in Interpreting the Visible Spectra of Polymethine Dyes: An Undergraduate Physical Chemistry Experiment. J. Chem. Educ. 1994, 71, 171−173. (50) Taubmann, G. Calculation of the Hückel Parameter β From the Free-Electron Model. J. Chem. Educ. 1992, 69, 96−97. (51) Bühl, M.; Hirsch, A. Spherical Aromaticity of Fullerenes. Chem. Rev. 2001, 101, 1153−1183. (52) Cyranski, M. K.; Howard, S. T.; Chodkiewicz, M. L. Bond Energy, Aromatic Stabilization Energy and Strain in IPR Fullerenes. Chem. Commun. 2004, 7, 2458−2459. (53) Salcedo, R.; Fomina, L. Homodesmotic Reaction for Fullerenes. Tetrahedron Lett. 2007, 48, 3949−3951. (54) Suresh, C. H.; Lincy, T. L.; Mohan, N.; Rakhi, R. Aromatization Energy and Strain Energy of Buckminsterfullerene from Homodesmotic Reactions. J. Phys. Chem. A 2015, 119, 6683−6688. (55) Wang, Y.; Díaz-Tendero, S.; Alcamí, M.; Martín, F. Cage Connectivity and Frontier π Orbitals Govern the Relative Stability of Charged Fullerene Isomers. Nat. Chem. 2015, 7, 927−934. (56) Wang, Y.; Díaz-Tendero, S.; Alcamí, M.; Martín, F. TopologyBased Approach to Predict Relative Stabilities of Charged and Functionalized Fullerenes. J. Chem. Theory Comput. 2018, 14, 1791− 1810. (57) Śliwa, W. Cycloaddition Reactions of Fullerenes. Fullerene Sci. Technol. 1995, 3, 243−281. (58) Haddon, R. C. Chemistry of the Fullerenes: The Manifestation of Strain in a Class of Continuous Aromatic Molecules. Science 1993, 261, 1545−1550. (59) Addicoat, M. A.; Page, A. J.; Brain, Z. E.; Flack, L.; Morokuma, K.; Irle, S. Optimization of a Genetic Algorithm for the Functionalization of Fullerenes. J. Chem. Theory Comput. 2012, 8, 1841−1851. (60) Wang, Y.; Díaz-Tendero, S.; Alcamí, M.; Martín, F. Relative Stability of Empty Exohedral Fullerenes: π Delocalization versus Strain and Steric Hindrance. J. Am. Chem. Soc. 2017, 139, 1609−1617. (61) Haddon, R. C.; Pasquarello, A. Magnetism of Carbon Clusters. Phys. Rev. B: Condens. Matter Mater. Phys. 1994, 50, 16459−16463. (62) Giannopoulos, G. I.; Georgantzinos, S. K.; Kakavas, P. A.; Anifantis, N. K. Radial Stiffness and Natural Frequencies of Fullerenes via a Structural Mechanics Spring-based Method. Fullerenes, Nanotubes, Carbon Nanostruct. 2013, 21, 248−257. (63) Rodríguez-Fortea, A.; Alegret, N.; Balch, A. L.; Poblet, J. M. The Maximum Pentagon Separation Rule Provides a Guideline for the Structures of Endohedral Metallofullerenes. Nat. Chem. 2010, 2, 955− 961. (64) Barriol, J.; Metzger, J. Application de la Méthode des Orbitales Moléculaires au Réseau du Graphite. J. Chim. Phys. Phys.-Chim. Biol. 1950, 47, 432−436. (65) Kutzelnigg, W. What I Like About Hückel Theory. J. Comput. Chem. 2007, 28, 25−34.

J

DOI: 10.1021/acs.jctc.8b00981 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX