Correlation between the energetic and thermal properties of C 40 fullerene isomers: An accurate machine-learning force field study

The isomerisation energies and thermal stabilities of the entire set of C 40 fullerene isomers are investigated using molecular dynamics simulations with the recently developed machine-learning-based Gaussian Approximation Potential (GAP-20) force field. Our simulations predict high statistical correlation between the relative iso-merisation energy distribution of all C 40 isomers and their thermal fragmentation temperatures. The most energetically and thermally stable isomer is the symmetric C 40 - D 2 isomer and its fragmentation temperature is 3875 ± 25 K. In contrast, the least stable isomer is the D 5d nanorod-shaped isomer which is 146.8 kcal mol (cid:0) 1 higher in energy and decomposes at 2975 ± 25 K, i.e., around 900 K below the most stable isomer. In addition, we show that most isomers lie in the range between 20 and 60 kcal mol (cid:0) 1 above the most stable isomer, and nearly 60% of isomers decompose in the temperature range of 3425 – 3675 K. These computational results provide valuable insights into the synthesis of C 40 fullerenes at high-temperatures.


Introduction
Fullerenes or buckyballs are molecular allotropes of carbon with geodesic graphene-like structures that contain pentagon and hexagon carbon rings, which bend the graphene sheet into spheres, ellipses or cylinders.Fullerenes were discovered by Kroto et al. [1] and were named after Richard Buckminster Fuller for their resemblance to his geodesic domes.The discovery of fullerenes motivated a research interest in finding new forms of carbon, which had previously been limited to graphite, diamond, and amorphous carbon.In addition to being discovered in the laboratory by experimental techniques, fullerenes also exist in natural environments and in the interstellar medium.For example, Buseck et al. [2] found fullerenes in a family of minerals known as "shungites".Furthermore, Cami et al. [3] observed fullerenes in a cloud of cosmic dust surrounding a star and, recently, Cordiner et al. [4] by using the Hubble Space Telescope, detected interstellar fullerenes in the space between two stars.
Fullerenes consist of at least twelve pentagons due to their closed cage structure and varying number of hexagons [5,6].All carbon atoms have a formal sp 2 hybridisation and are bonded to three other carbon atoms.The molecular structure of a fullerene is edgeless, without any dangling bonds.These characteristics make the fullerenes unique and different from such crystalline structures as graphite and diamond with edges and dangling bonds at the surface.Perhaps the most well known fullerenes are C 60 and C 70 ; nevertheless, fullerenes with smaller numbers of carbon atoms (starting from C 20 ) are also ubiquitous and well studied [7].Fullerenes have unique physical and chemical properties (such as spectroscopic, thermodynamic, electronic, electrical, and optical properties) and a number of notable potential applications (e.g., in organic solar cells, chemical sensors, molecular containers, gas storage/separation, fuel cells, hydrogen storage, storage for radioactive isotopes, lightemitting devices, fullerene switches, superconductor, lubricants, and direct bandgap semiconductors) [7][8][9][10][11].Accordingly, much experimental and theoretical effort has been devoted to studying fullerenes.
Computer simulations and in particular molecular dynamics (MD) simulations of carbon materials are a powerful tool for complementing carbon related experiments [12][13][14][15][16][17][18].In addition, one of the strengths of MD simulations is the simplicity of performing virtual experiments, which might be challenging or even impossible in the laboratory.Furthermore, fullerenes are generally produced at high-temperature processes such as arc burning or laser vaporisation of graphite, and understanding the thermal stability of fullerenes may play a significant role in controlling the experimental synthesis of fullerene molecules.A few theoretical investigations have studied the process of thermal ; MD, Molecular Dynamics; DFT, Density Functional Theory; GAP, Gaussian Approximation Potential.fragmentation of carbon fullerenes.Zhang et al. [19] used tight-binding MD simulations to investigate the thermal disintegration of fullerenes ranging from C 20 to C 90 .They found that the fragmentation temperature is sensitive to the cluster size for the small fullerenes (≤ C 58 ) but approximately constant for the larger ones.Furthermore, in a separate work, Zhang et al. [20] reported MD simulations of the thermal stability of C 20 , C 26 , C 36 , and C 60 fullerenes.They used the second generation of the Reactive Empirical Bond Order carbon force field [21] in their simulations for the carbon-carbon interatomic interactions and showed that the large fullerenes are more stable than the small fullerenes.In the present work, we use the recently developed machine-learning-based Gaussian Approximation Potential (known as the GAP-20) which is specifically developed for carbon materials to investigate the isomerisation energies and thermal stabilities of the entire set of C 40 isomers.To the best of our knowledge, this is the first investigation using a highly accurate machine-learning force field for a systematic investigation of the thermal stability of an entire set of fullerene isomers.

Computational details
All the MD simulations were performed using the machine learning GAP-20 force field [22] to model the thermal and energetic stabilities of the C 40 fullerene isomers.The GAP-20 force field was developed via a high-dimensional fit to reference data computed using the optB88-vdW density functional theory (DFT) method, including dispersion and longrange interactions where kernel-based machine-learning is employed instead of a classic empirical functional form.It should be pointed out that the GAP-20 database covers all the possible crystalline and amorphous phases of the carbon and exotic carbon allotropes that comprise approximately 17,000 various structures, each containing up to 240 atoms in their unit cell.The GAP-20 force field outperforms widely used conventional empirical force fields (e.g., Tersoff [23], the second generation of Reactive Empirical Bond Order [21], Adaptive Intermolecular Reactive Empirical Bond Order [24], and the first version of Long-range Carbon Bond Order Potential [25]) for the calculation of lattice parameters, bond lengths, formation energies, and phonon dispersions of numerous crystalline phase of carbon and amorphous carbon, such as graphite, diamond, nanotubes, and fullerenes.
In the context of the present work, it is important to point out that the GAP-20 force field was recently found to significantly outperform a wide range of empirical carbon force fields for the 1811 isomerisation energies of the C 60 fullerenes [26], the geometries and the carbon-carbon bond lengths of prototypical C 60 fullerenes [27], and for the interaction energy of small double-layered fullerene shells [28].In addition, Qian et al. [29] investigated the performance of the GAP-20 force field for various properties, such as lattice constants, cohesive energies, van der Waals interactions, thermal stabilities, and mechanical properties for different carbon allotropes.They found that GAP-20 is the most accurate force field for studying the thermal stability of the high symmetry C 60 -I h fullerene.
All MD simulations were performed using the Large-scale Atomic/ Molecular Massively Parallel Simulator (LAMMPS) package [30], which was built with the Quantum mechanics and Interatomic Potentials (QUIP) module including GAP routines [31,32] (see www.libatoms.orgfor further information).The geometry relaxation and the energy minimisation simulations were performed using the conjugate gradient scheme, and the convergence of energy and force are set to 10 − 12 eV and 10 − 10 eV/Å, respectively.To analyse the thermal stability of isomers, they were annealed at a fixed temperature for 10 ps.To ensure the length of the annealing simulation is sufficient, we considered a longer timescale of 100 ps for the C 40 -D 2 isomer at the three different fixed annealing temperatures of 30, 300, and 1000 K.These results are depicted in Fig. S1 of the Supporting Information.As expected, the potential energy fluctuation depends on the annealing temperature which is being simulated.However, the fluctuations remain qualitatively unchanged after 10 ps, indicating that this timescale is sufficient.
Annealing simulations were carried out in the NVT ensemble (meaning that the number of particles (N), volume (V), and temperature (T) are conserved quantities).Following the work of Qian et al. [29], we used the Nose-Hoover thermostat [33] to control the temperature.The equations of motion are integrated using the velocity Verlet algorithm, and a time step was 0.0005 ps in the entire simulations process.
The initial isomers geometries which are used in our study are taken from an online tool [34].The periodic boundary conditions were used in all three directions for all simulations.Note that to avoid the selfinteraction of the fullerene structure, we performed the simulations in a cubic box with a side length of 20 Å.The Visualisation for Electronic Structural Analysis (VESTA) software [35] is used to visualize the optimised isomers which are shown in Fig. 1, while the visualisations of the annealing simulations (Fig. 2) were performed using the Open Visualisation Tool (OVITO) software [36].Finally, coordination analysis was performed by counting the number of nearest neighbours within a cutoff 1.85 Å.For analysis and visualisation purposes, carbon atoms were considered to be sp-hybridized if they have two neighbours and sp 2hybridized if they have three neighbours.

Isomerisation energy
In the first part of our study, we investigate the prediction of the isomerisation energy for forty C 40 fullerene isomers.It should be pointed out that for the geometry relaxation and the energy minimisation, we follow the same methodology as used before in our previous evaluation of the performance of the GAP-20 and empirical force fields for the isomerisation energies of C 60 isomers relative to hybrid-DFT at the PW6B95-D3/Def2-QZVP level of theory (see Ref. [26] for further details).most and least energetically and thermally stable isomers) are also illustrated at the top of Fig. 1.The values of the relative isomerisation energies and the point-group symmetry for all isomers are provided in Table 1.
The most stable isomer is a spherical fullerene with D 2 symmetry (isomer#1).Relative to the energy of this isomer, the isomerisation energies of the higher-energy isomers range from 8.25 kcal mol − 1 for the second most stable C s symmetry isomer up to 146.79 kcal mol − 1 for the least stable nanorod-shaped isomers of C 40 -D 5d (isomer#40).As shown in Fig. 1, the distribution energy of all isomers increases in a nonlinear fashion.Similar to the energy distribution of the C 60 and C 80 fullerenes [16], the relative energy curve of the C 40 isomers has a logit-type distribution.It should be pointed out that the average value for the relative energy of the entire set of forty isomers is 56.15 kcal mol − 1 , and the results show that most isomers (~53%) lie in the range between 20 and 60 kcal mol − 1 above the most stable isomer.Finally, it is of interest to note that the two lowest nanorod-shaped isomers, with the D 2 (#39) and D 5d (#40) symmetry, are separated by 12.16 kcal mol − 1 .

Thermal stability
We proceed to use MD simulations with the machine-learning-based GAP-20 force field to examine the thermal stability and thermal disintegration of the C 40 isomers.For this purpose, we annealed the optimized isomers for 10 ps at 35 different temperatures (ranging from 30 to 4000 K).As mentioned in the methodology section, the length of the annealing simulations is sufficient for the temperature and potential energy to be equilibrated.Observing the annealing process indicates that, at sufficiently high temperatures, the carbon-carbon bonds break due to the large thermal fluctuations of the atoms.In particular, the atoms with dangling bonds that are located outside the cage-like surface fluctuate much faster than the fully bonded carbon atoms.This results in more cleavages of carbon-carbon bonds until the whole structure is completely decomposed.Fig. 1 illustrates the fragmentation temperatures (T frag , light blue circles) for all the C 40 isomers.The fragmentation temperatures and corresponding error bars are calculated as follows; if an isomer decomposes at a temperature T 2 but not at a lower temperature T 1 , then the fragmentation temperature is taken as T frag = (T 1 + T 2 )/2 and the error bar is simply taken as T 2 -T 1 .The fragmentation temperature (T frag ) for each isomer is provided in Table 1.As a general trend, our results predict that the fragmentation temperature of C 40 isomers decrease non-linearly (with a distribution curve that follows the shape of the relative energy distribution).In addition, as can be seen from the figure, the results reveal that the most spherical-shaped isomers decompose at higher temperatures, while the nanorod-shaped isomers decompose at lower annealing temperatures.
The fragmentation temperatures of the C 40 isomers range from 3875 ± 25 K for the thermally and energetically most stable C 40 -D 2 isomer (#1) to 2975 ± 25 K for the least stable D 5d symmetric nanorod isomer (#40).In particular, the least stable nanorod-shaped isomer decomposes at a temperature ~900 K lower than the most stable D 2 isomer.As can be seen from Fig. 1, the average value of the fragmentation temperatures of all isomers is 3520 K.In addition, it is worth noting that the fragmentation temperature of most isomers (i.e., 60% of isomers) lies in a range of 3425 to 3675 K, which is nearly 200 to 400 K below the most stable isomer.
Atomistic details of the time evaluation of the fullerene decomposition process illustrates that the carbon atoms experience more pronounced fluctuations at the fragmentation temperature than at lower annealing temperatures.Therefore, the carbon-carbon bonds are broken due to drastic local fluctuations.Once the bond breaking occurs, the sphybridized atoms with dangling bonds fluctuate much faster than the sp 2 -hybridized carbons, and finally, the entire fullerene structure is decomposed.Snapshots of the atomistic structures of the most and least energetically and thermally stable isomers of C 40 and seven other selected isomers between them as a function of annealed temperature are illustrated in Fig. 2. For example, our results predict that C 40 -D 2 isomer (which is the most energetically and thermally stable isomer) somewhat deformed at 3800 K but remains intact; however, the structure is decomposed and transformed to an amorphous-like carbon at 3900 K, and at the higher annealing temperature the structure converts to a long carbon chain.Finally, it should be emphasised that the machine-learning GAP-20 isomerisation energies and the fragmentation temperatures are consistent with each other for the entire set of C 40 isomers.That is, they follow the same order of energetic stability and thermal instability.This qualitative correlation between the energy of the isomers and their corresponding thermal fragmentation temperatures is an important observation which increases our confidence in the GAP-20 results.To the best of our knowledge, this correlation has not previously been reported for a systematic set of fullerene isomers.Our developed methodology can be easily extended to investigate the energy and thermal stabilises of other fullerenes.Albeit, one has to be conscious of the fact that the number of isomers increases exponentially with the number of carbon atoms; for example, C 60 has 1812 isomers [6].

Conclusions
We perform extensive molecular dynamics simulations using the recently developed machine-learning GAP-20 force field to predict the isomerisation energies and thermal stabilities of the entire set of C 40 fullerene isomers.Our simulations predict a systematic correlation between the relative isomerisation energy distribution and their thermal fragmentation temperatures of the forty C 40 isomers.The results show that the energetically most stable isomer is the spherical D 2 isomer and the least stable one is the D 5d symmetric nanorod isomer, where the energy difference between them is 146.8 kcal mol − 1 .Most isomers (~53%) lie in the energy range of 20 to 60 kcal mol − 1 above the most stable isomer.Notably, our simulations predict the thermal stability of C 40 isomers and provide reliable values for the fragmentation temperatures distribution for the forty isomers, where the maximum fragmentation temperature is 3875 ± 25 K for the energetically most stable C 40 -D 2 isomer, the average fragmentation temperature is 3520 K, and nearly 60% of isomers fall in a fragmentation temperature window of 3425 to 3675 ± 25 K.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.A. Aghajamali and A. Karton

Fig. 1 Fig. 1 .
Fig. 1.Correlation between relative isomerisation energy (pink circles) and fragmentation temperature (T frag , light blue circles) for the forty C 40 fullerene isomers.The error bars around the fragmentation temperatures indicate the uncertainty in the fragmentation temperature (see main text).Optimized structures of the energetically and thermally most and least stable C 40 fullerenes are shown along with an intermediary stable structure.(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 2 .
Fig. 2. Atomic structure snapshots of nine different C 40 isomers after annealing for 10 ps at five different temperatures.Red, green, and blue circles denote sp, sp 2 , and sp 3 hybridisation, respectively; yellow carbons have just one neighbour.Note that a cutoff bond distance of 1.85 Å is used to determine the coordination number of the carbon atoms.(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Table 1
Point-group symmetries, isomerisation energies relative to the most stable isomer (ΔE isomer )*, and the thermal fragmentation temperatures (T frag ) for the C 40 isomers.