Comparative Study of Carbon Force Fields for the Simulation of Carbon Onions

Keywords: carbon nano-onions, carbon nanomaterials, fullerenes, machine-learning force fields, empirical force fields, interaction energy.

Carbon onions are nanoforms of carbon materials and members of the fullerene family, with quasi-spherical concentric fullerene-like shells that range from double to multi-layered structures. They are known as onion-like carbon or carbon nano-onions and were first synthesized by Iijima^[1] in arc discharge experiments. The intershell distance in carbon onions is about 3.4 Å which is slightly larger than the distance between graphene layers in graphite.^[2] Carbon onions are an attractive research target due to their unique thermodynamic, electronic, electrical and optical properties and the possibility of conversion into other forms of carbon.^[2,3] For example, they have high electrical conductivities and improved dispersibility relative to graphene, making them promising materials for electrochemical applications. Due to their excellent capacitance retention at high current densities and improved solubility, carbon onions are ideal additives to lithium-ion batteries^[4,5] and supercapacitor electrodes for high power applications.^[6,7] In addition, they have potential applications in cellular imaging,^[8,9] tribology, catalysis, composite materials, and electromagnetic devices.^[3,10,11]

Over the past four decades, in addition to arc discharge, which is the most common approach to generating carbon onions, other synthetic methods have been employed. These include electron beam irradiation,^[12–14] laser irradiation,^[15] high-temperature nanometre-sized diamond annealing,^[16–21] implantation of carbon ions onto metal particles^[22,23] and chemical vapour deposition (CVD). It should be noted that with the exclusion of the CVD method, all the above-mentioned methods are costly due to their high-temperature processes, and they produce a limited number of carbon onions in highly condensed forms. In contrast, CVD is an effective method for generating pure carbon onions in large amounts.^[24]

Atomistic simulations play an important role in modelling carbon onions. Generally, calculating the energy of a system of atoms is a prerequisite requirement for carrying out atomistic simulations. A number of density functional theory (DFT) studies have explored the formation and interaction energies in carbon onions.^[25–33] However, due to the computational cost associated with DFT methods, these studies have been limited to systems with hundreds of carbons. Molecular mechanics, on the other hand, is computationally much more economical and is routinely used for calculating structural properties of nanomaterials. Molecular mechanics also provides a practical avenue for running molecular dynamics (MD) simulations which are of key importance in studying the dynamic behavior of carbon nanomaterials. For examples of force field-based MD simulations of carbon onions see:^[34–41] The stability of the multi-layered carbon onions alongside the interaction energies between their layers are important in a wide range of chemical applications, such as lithium-ion batteries and supercapacitors.

In the present work we use a wide range of commonly used carbon force fields to investigate the interaction energies in a range of carbon onion structures. Namely, the double-layered C₂₀@C₆₀, C₂₀@C₈₀, C₆₀@C₁₈₀, C₆₀@C₂₄₀, C₆₀@C₂₄₀ and C₂₄₀@C₅₄₀ carbon onions, and C₆₀@C₂₄₀@C₅₄₀ and C₈₀@C₂₄₀@C₅₄₀ triple-layered carbon onions. We evaluate the performance of ten commonly used force fields developed for carbon relative to the DFT interaction energies obtained by Hashmi and Lein.^[33] We then proceed to use to most reliable force field to study the stability and the interaction energies in a giant carbon onion with a C₆₀ core (C₆₀@C₂₄₀@…@C₄₈₆₀@C₆₀₀₀).

The first force field for carbon (Tersoff^[42]) was developed in 1988 and modified the bond strengths according to the number of neighbours. It was fitted to experimental data such as lattice constant, binding energies, bulk modulus of diamond, and binding energy of graphite. The energies of carbon–carbon bond energies and simple cubic, body-centred cubic and face-centred cubic phases are fitted to DFT calculations. In 2013, to improve the description of making and breaking bonds, an environment-dependent screening function was added to the Tersoff force field, which is known as the Tersoff-S.^[43]

The second generation of reactive empirical bond-ordered (REBO-II) force field,^[44] an extended version of the Tersoff force field, was developed in 2002. REBO-II improves the description of short-range bonding by modifying the functional form. Its parameters are fitted to experimental data such as binding energies and lattice constants of graphite and diamond, and DFT data such as binding energies and lattice constants of simple cubic and face-centred cubic phases of carbon. In 2008, the description of the bond-making and breaking was improved by adding an environment-dependent screening function to the REBO-II force field (known as the REBO-II-S).^[45]

In 2000, the adaptive intermolecular reactive empirical bond-order (AIREBO)^[46] force field based on the REBO-II force field was developed, in which a Lennard-Jones term with a switching function is used to inactive the long-range terms at short distances. Another version of the AIREBO force field (named AIREBO-M) is derived by replacing the Lennard-Jones form with a Morse expression.^[47] The long-range carbon bond-order potential (LCBOP-I),^[48] which is generally similar to the AIREBO force field, was developed in 2003 to improve the accuracy and transferability of the model by complete reparameterization of the bond order force field.^[49]

The reactive force field potential ReaxFF-15 was developed in 2015 and originally developed for hydrocarbon systems.^[50] This force field was parametrized against a large dataset of DFT-D2 structures, including small fullerenes and different amorphous carbon phases. The original force field^[51] was fitted to experimental data such as heats of formation and structural data for hydrocarbon compounds and crystals, and quantum chemical energy curves for bond dissociations. Recently, a reparametrization version of the ReaxFF force field^[52] was developed (ReaxFF-20).^[53] The force field parameters were optimized to adequately capture the flattened carbon nanotube structure and the corresponding energy values. Its parameters were trained against the PCFF-IFF force field data, which shows good agreement with DFT-D2 data.

In addition to empirical force fields, we use a machine-learning force field, which was developed recently based on the Gaussian approximation potential (GAP-20).^[54] In this type of force field, kernel-based machine learning is employed instead of a classic empirical functional form and the force field is fitted to reference data computed using the optB88-vdW DFT method. In this context it should be noted that we have recently showed that the GAP-20 force field shows good performance in the prediction of the 1811 isomerization energy of C₆₀ fullerenes.^[55] Furthermore, Qian et al.^[56] reported that this force field could more accurately predict the structure and formation energies for different carbon systems. They also showed that the GAP-20 is the most appropriate force field for studying thermal and mechanical properties of carbon materials.

The initial fullerene geometries (C₂₀, C₆₀, C₈₀, C₁₈₀, C₂₄₀ and C₅₄₀) and the two/three shell carbon onion structures that we used in the first part of our study (i.e., C₂₀@C₆₀, C₂₀@C₈₀, C₆₀@C₁₈₀, C₈₀@C₂₄₀, C₆₀@C₂₄₀, C₂₄₀@C₅₄₀, C₆₀@C₂₄₀@C₅₄₀, and C₈₀@C₂₄₀@C₅₄₀) were taken from Hashmi and Lein.^[33] The double- and triple-layered structures are shown in the right panel of Fig. 1. The initial geometries for the giant C_N (N = 60 × n² for 1 ≤ n ≤ 10) fullerenes were taken from Noël et al.^[57] To generate the giant C₆₀ core multi-shell onion structures, we followed the same procedure that used by Hashmi and Lein^[37] and stacking the structures of the individual fullerenes, where each fullerene is aligned along its respective principal axes of inertia in all three possible orientations.

Fig. 1.  Interaction energies (E_int) of doubly and triply layered carbon onions calculated with commonly used carbon force fields compared with DFT interaction energies taken from Hashmi and Lein.^[33] The lines are fitted to the data to guide the eye. The E_int values for all force fields and DFT calculation are listed in Table 1. The right panel shows the snapshots of the DFT optimized double- and triple-layered onion structures.

Table 1.  Interaction energies (E_int, in kcal mol⁻¹) of double- and triple-layered carbon onion structures; the DFT E_int values are taken from Hashmi and Lein^[33]

All calculations are performed with the Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) package.^[58] The geometry relaxation and the energy minimization were performed using the conjugate gradient scheme, and the convergence of energy and force are set to 1.0⁻⁶ eV and 1.0⁻⁸ eVÅ⁻¹, respectively. The periodic boundary conditions were used in all three directions for all simulations. Note that to avoid the self-interaction of the fullerene and onion structures, we performed the simulations in a cubic box with a side length of 30 Å for double- and triple-layered onions and a 120 × 120 × 120 ų box for a giant onion structure. The visualizations are performed using the Open Visualization Tool (OVITO) software.^[59] The interaction energy of the onion structures are calculated directly from the total energy of the carbon onion and the total energies of the inner structure(s) and outer fullerene which can be expressed as

For example, for the C₆₀@C₂₄₀@C₅₄₀ onion structure the interaction energy can be calculated as follows:

In the first part of our study we evaluate the performance of ten force field methods developed for carbon nanomaterials for the interaction energies between the layers of eight carbon onion structures. We compare the force field results with density functional theory (DFT) reference interaction energies obtained by Hashmi and Lein.^[33] In particular, the DFT reference values were obtained at the RI-PBE0-D3/def2-TZVP//RI-PBE0-D3/def2-SVP level of theory. The calculated interaction energies (E_int) between the layers of the carbon onions for different force fields along with the DFT reference values are plotted in Fig. 1 and listed in Table 1. Before we discuss the force field results, we note that Hashmi and Lein^[33] also obtained second-order Møller-Plesset perturbation theory (MP2) interaction energies for a subset of the smaller carbon onions and that the DFT interaction energies for C₂₀ and C₆₀ nested in larger fullerenes are in qualitative agreement with the MP2 interaction energies. Namely, the interaction energy decreases as the size of the carbon onion increases. We begin by examining the interaction energies of the double-layered carbon onion structures, C₂₀@C₆₀, C₂₀@C₈₀, C₆₀@C₁₈₀, C₆₀@C₂₄₀, C₈₀@C₂₄₀, C₂₄₀@C₅₄₀. As can be seen from Fig. 1, the general trends of most of the force fields results (except the Tersoff and REBO-II force fields) agree with the DFT trends.

For the smallest system, C₂₀@C₆₀, most force fields tend to significantly underestimate the repulsive interaction energy between the nested fullerenes calculated at the DFT level. The LCBOP-I (see Table 1) is the only force field which significantly overestimates the interaction energy by ~25 %. Interestingly, the GAP-20 force field is the only method that is able to reproduce the DFT interaction energy within 10 kcal mol⁻¹. In particular, GAP-20 overestimates the DFT result by 9.3 kcal mol⁻¹. The ReaxFF-20 force field is a close second, and underestimates the DFT result by 25.2 kcal mol⁻¹. However, for the other force fields, i.e., Tersoff, Tersoff-S, REBO-II, REBO-II-S, AIREBO, AIREBO-M and ReaxFF-15, large underestimations are obtained (Table 1).

For the C₂₀@C₈₀ carbon onion, the ReaxFF-20 force field perfectly estimates the DFT interaction energy with a deviation of merely 2.8 kcal mol⁻¹. Three force fields (i.e., Tersoff-S, AIREBO-M, and ReaxFF-15) result in reasonably good performance with deviations from the DFT interaction energies ranging between 17.7 (ReaxFF-15) and 29.7 kcal mol⁻¹ (AIREBO-M); note that the deviations are given in absolute value. In contrast to the C₂₀@C₆₀ carbon onion for which the GAP-20 force field showed good performance, for the C₂₀@C₈₀ it results in a large deviation of 151.3 kcal mol⁻¹. The remaining force fields show performance that is on par or worse than that of the GAP-20 force field. For the C₂₀@C₁₈₀ carbon onion, REBO-II-S force field reproduces the DFT interaction energy with a small deviation of 4.2 kcal mol⁻¹. The second best force field is LCBOP-I which overestimates the DFT result by 19.1 kcal mol⁻¹. In contrast to previous cases, the predicted interaction energy by the ReaxFF force fields are underestimated with the deviation of 76.7 (ReaxFF-15) and 60.8 kcal mol⁻¹ (ReaxFF-20). In addition, similar to the C₂₀@C₈₀ carbon onion case, the GAP-20 force field results in a large deviation of 101.5 kcal mol⁻¹.

Moving on to the C₆₀@C₂₄₀ carbon onion, the ReaxFF-20 force field shows the best performance with an estimated interaction energy which is 11.5 kcal mol⁻¹ away from the DFT value. The AIREBO and ReaxFF-15 force fields are close seconds with deviation from the DFT result by 25.7 (AIREBO) and –33.6 kcal mol⁻¹ (ReaxFF-15). The calculated E_int values of the other force fields are considerably overestimated. It should be pointed out that the interactions energy estimated by ReaxFF force fields not only agrees with the DFT calculation by Hashmi and Lein^[33] but also agrees well with Stasyuk et al.^[32] and Casella et al.^[31] where they used different DFT functionals. Finally, we note, with the exception of the original Tersoff and REBO-II force fields, all the considered force fields show similar trends in the interaction energies to the DFT results.

The estimated interaction energy for the C₈₀@C₂₄₀ carbon onion by the ReaxFF-20 force field is close to the DFT results with a deviation of 5.3 kcal mol⁻¹. GAP-20 is the second-best force field and underestimates the DFT result by 13.0 kcal mol⁻¹. However, the other force fields show poor performance with large deviations ranging between 78.9 (LCBOP-I) and 144.2 kcal mol⁻¹ (AIREBO).

The C₆₀@C₂₄₀ and C₈₀@C₂₄₀ carbon onions are the smallest systems characterized by attractive interactions between the inner and outer fullerenes at the DFT level, where the attractive interaction energy is larger in the C₈₀@C₂₄₀ carbon onion. The ReaxFF-20 force field is able to nicely capture these trends in interaction energies both qualitatively and quantitatively. The GAP-20, ReaxFF-15, and LCBOP-I force fields also result in attractive interaction energies; however, they are unable to capture their relative magnitudes.

Finally, for the largest double-layered carbon onion (C₂₄₀@C₅₄₀) the ReaxFF-15 force field shows good performance with a deviation of just 12.6 kcal mol⁻¹ from the DFT E_int value. However, the other force fields tend to severely underestimate the interaction energies.

Now we examine the interaction energies of triple-layered C₆₀@C₂₄₀@C₅₄₀ and C₈₀@C₂₄₀@C₅₄₀ carbon onions. As can be seen from Fig. 1, these interaction energies pose a challenging target for all the force fields. The best performing force field is the ReaxFF-15 force field with deviations of over 100 kcal mol⁻¹ from the DFT interaction energies, whilst the other force fields result in significantly larger deviations (Table 1). Nevertheless, with few exceptions (e.g. Tersoff, REBO-II, REBO-II-S, and GAP-20) the considered force fields predict considerable attractive interactions for both triple-layered carbon onions. We note that the surprising and notable deterioration in performance of the GAP-20 force field when moving from the double-layered carbon onions to the triple-layered and the giant C₂₄₀@C₅₄₀ carbon onions may be attributed to the fact that the machine-learning training set included much smaller fullerenes.

Overall, the above results show that the new generation of the ReaxFF force field (ReaxFF-20) is the best force field for the double-layered carbon onions, with the older version (ReaxFF-15) being a close second-best performer (with a notable exception for the C₈₀@C₂₄₀ carbon onion). Whilst for the triple-layered carbon onions the ReaxFF-15 emerges as the best performing force field. However, the other force fields can at best predict qualitative trends in the DFT interaction energies of the carbon onions. Based on these results we choose the ReaxFF-15 force field to investigate the interaction energies involved in a giant C₆₀ core multi-shell onion structure. In order to produce a multi-layered onion structures we are using an sequence of

where n is the number of layers and is equal and greater than 3 (n ≥ 3), and . For example, a fourth-layered onion structure is , where ; hence, . Table 2 contains the nine multi-layered onion structures, number of carbon atoms and calculated interaction energy for each structures.

Table 2.  Interaction energy (E_int) and ΔE_int of nine different C₆₀ core carbon onion structures
Number of carbon atoms (N_atoms) in each onion structures are also provided. All the energy values are in kcal mol⁻¹

Fig. 2a illustrates the cross-sectional view (1 nm slice) of the giant C₆₀ core ten-layered ‘rainbow’ carbon onion structure that we investigate in this part of our study. Each colour represents different fullerene, and the red core structure is C₆₀ fullerene, and the violent outer layers is C₆₀₀₀ fullerene. Fig. 2b shows the experimental evidence of giant carbon onion produced by electron irradiation of a polyhedral particle at a high temperature (700°C). This panel is reprinted with permission from Banhart and Ajayan.^[60] During the past few decades, advancing the technology to synthesise a giant multi-layered carbon onion simulating these materials is an important area of computational nanotechnology. In this regard, in this part of our study, we investigate the interaction energy of the proposed giant multi-layered carbon onion.

Fig. 2.  Giant multi-layered carbon onion structure. (a) The cross-sectional view (10 Å slice) of the simulated C₆₀ core ten-layered carbon onion system. The geometry of individual fullerenes structures are taken from Noël et al.^[57] (b) Carbon onion generated by electron irradiation of a polyhedral particle at high temperature. Reprinted with permission from Banhart and Ajayan.^[60] © 1991 Wiley-VCH.

Table 2 gives the interaction energies for all the nine multi-layered onion structures, which are calculated using the ReaxFF-15 force field. The results for all nine multi-layered carbon onions reveal that the interaction energy increases rapidly with the number of layers and subsequent size of the onion structure. For example, the interaction energy for the three-layered C₆₀@C₂₄₀@C₅₄₀ onion amounts to –516.5 kcal mol⁻¹ and it increaces up to –8117.3 kcal mol⁻¹ for the largest carbon onion with 10 nested fullerenes. Inspection of Table 2 reveals that the increase in interaction energy with each consecutive shell (i.e., ) correlates linearly with the number of shells and following size of the multi-layered structure (with a squared correlation coefficient of R² = 0.9996). This linear correlation is illustrated in Fig. 3. This linear increase in the stabilization energy with each consecutive layer may have important thermodynamic consequences for describing the formation and growth of large carbon onions.

Fig. 3.  The ΔE_int of nine different C₆₀ core carbon onion structures as a function of number of layers. The average diameters of the carbon onion structures are also given in Å.

We evaluate the performance of ten well-known carbon force fields for the interaction energy of six different double-layered carbon onions (C₂₀@C₆₀, C₂₀@C₈₀, C₆₀@C₁₈₀, C₈₀@C₂₄₀, C₆₀@C₂₄₀ and C₂₄₀@C₅₄₀) and two triple-layered onions (C₆₀@C₂₄₀@C₅₄₀ and C₈₀@C₂₄₀@C₅₄₀). The reference interaction energy values have been obtained at the PBE0-D3/def2-TZVP level of theory. The ReaxFF-20 force field provides the best overall performance for the double-layered carbon onions, whilst the ReaxFF-15 force field provides the best overall performance for the triple-layered carbon onions. Interestingly, both force fields outperform the recently developed carbon machine-learning force field (GAP-20) for both double- and triple-layered carbon onions. We proceed to use the ReaxFF-15 force field to explore the interaction energies in a giant ten-layered carbon onion with a C₆₀ core and show that the interaction energy between the outer layer and the inner layers increases linearly with the number of layers in the carbon onion (with a squared correlation coefficient of R² = 0.9996). This linear increase in the stabilization energy with each consecutive layer may have important thermodynamic consequences for describing the formation and growth of large carbon onions.

ReaxFF-15 optimized geometries for the double-layered, triple-layered, and giant carbon onions considered in the present work are available on the Journal’s website.

The data that support this study are available in the article and accompanying online supplementary material.

Amir Karton is an Associate Editor of the Australian Journal of Chemistry but was blinded from the peer-review process for this paper.

This research was funded by the Australian Research Council (Project no. FT170100373).