Carbon materials and their unique properties have been extensively studied by molecular dynamics, thanks to the wide range of available carbon bond order potentials (CBOPs). Recently, with the increase in popularity of machine learning (ML), potentials such as Gaussian approximation potential (GAP), trained using ML, can accurately predict results for carbon. However, selecting the right potential is crucial as each performs differently for different carbon allotropes, and these differences can lead to inaccurate results. This work compares the widely used CBOPs and the GAP-20 ML potential with density functional theory results, including lattice constants, cohesive energies, defect formation energies, van der Waals interactions, thermal stabilities, and mechanical properties for different carbon allotropes. We find that GAP-20 can more accurately predict the structure, defect properties, and formation energies for a variety of crystalline phase carbon compared to CBOPs. Importantly, GAP-20 can simulate the thermal stability of C60 and the fracture of carbon nanotubes and graphene accurately, where CBOPs struggle. However, similar to CBOPs, GAP-20 is unable to accurately account for van der Waals interactions. Despite this, we find that GAP-20 outperforms all CBOPs assessed here and is at present the most suitable potential for studying thermal and mechanical properties for pristine and defective carbon.

Carbon is one of the most widely studied and important elements across many branches of science. It has many allotropes each with remarkable properties, such as three-dimensional bulk materials e.g., graphite, diamond, and amorphous carbon and low-dimensional materials e.g., graphene, carbon nanotubes (CNTs), and fullerenes.1–6 With the development of research methods comes the ability to understand materials on very small spatial and temporal scales, though not all can be verified on an experimental level. High-performance computing has enabled theoretical studies of the structure,7,8 growth mechanism,9–12 electrical13,14 and thermal transport,15–17 and mechanical properties18–24 of a vast range of carbon allotropes.

Density functional theory (DFT) models the electronic structure of materials and has the advantage of high precision, but it is too computationally expensive for large-scale molecular dynamics (MD) simulations. Conversely, carbon bond order potentials (CBOPs) can simulate large-scale systems over time, by using empirical potentials, but with lower accuracy. Therefore, significant research has been devoted to developing empirical potentials for carbon, which can simulate large-scale systems while maintaining high accuracy.25–34 In 1988, Tersoff was the first to propose a multi-body potential model for carbon materials employing a bond order formalism, namely, the Tersoff potential,25 which gives a highly accurate description of the structure and energy of amorphous carbon. Brenner et al. developed reactive empirical bond order potentials Brenner-I (REBO-I) and Brenner-II (REBO-II) based on the Tersoff potential with the addition of torsional interactions.26,30 Both Tersoff and REBO potentials are good at representing covalent bonds but ignore long-range interactions, such as van der Waals interaction. Later, Stuart et al. added a Lennard-Jones (LJ) term to the REBO potential to make up for the missing long-range interaction, creating the adaptive intermolecular reactive bond order (AIREBO) potential.28 However, AIREBO exhibits an unphysically divergent power-law repulsion when under pressure, which was solved by O’Connor et al., who used quantum chemistry to fit a Morse potential instead of a LJ potential.33 Furthermore, Lindsay and Broido optimized the Tersoff and REBO parameter sets to improve the lattice dynamics and phonon dispersion accuracy, considerably improving the calculated thermal transport properties for graphene and CNTs compared to experiments.32 Los and Fasolino built a long-range bond-order potential (LCBOP) for carbon,31 which gives a good description of the diamond to the graphite phase transformation. van Duin et al. developed the ReaxFF potential to study large-scale reactive chemical systems for hydrocarbons,29 and in 2015, Srinivasan et al. reparameterized the original ReaxFF potential to create ReaxFFC2013, which more accurately describes condensed carbon phases.34 Justo et al. developed the environment-dependent interatomic potential (EDIP) for silicon as an improvement over the Tersoff potential.35 Marks and Jiang et al. extended EDIP for carbon and silicon–carbon systems.27,36 CBOPs have been used to explore the mechanical properties of carbon materials ever since they were developed. However, some CBOPs have been shown to cause a non-physical force enhancement in tensile test simulations, and a 1.92–2.0 Å C–C cutoff was suggested to avoid this phenomenon.37–40 He et al.38 carried out MD simulations on the Stone–Wales (SW) defect in graphene and revealed that defect orientation would dominate the mechanical properties. Song et al.39 explored the fracture behavior of polycrystalline graphene via MD simulation, showing that the fracture first occurs at the grain boundary. Xu et al.40 performed MD simulations on defective graphene and found that graphene with a large number of vacancy defects exhibits super-ductility.

In recent years, machine learning (ML) has been applied to the development of interatomic potentials.41–52 Unlike the CBOPs mentioned above, machine learning interatomic potentials (ML-IAPs) do not rely on fixed mathematical expressions but instead learn the mathematical representation of the potential energy surface (PES) via training, which allows for more accurate predictions of energy, forces, and so on. For example, Deringer and Csányi built a Gaussian approximation potential (GAP) for liquid and amorphous carbon that showed close-to-DFT accuracy.43 Wen and Tadmor employed a neural network to model short-range interactions and a theoretically motivated analytical term to model long-range dispersion in their hybrid neural network potential (hNN-Grx) for multilayer graphene, which can be used to study thermal conductivity in monolayer graphene and friction in bilayer graphene.49 Wen and Tadmor also developed a set of Dropout Uncertainty Neural Network (DUNN) potentials for carbon to predict uncertainty for static (e.g., energy) and dynamical (e.g., phonon dispersion) properties.51 Hedman et al. used deep learning to train neural network potentials on the nine-carbon allotrope dataset (CA-9), which reproduced ab initio results with high accuracy.52 Finally, Rowe et al. developed the GAP-17 potential for graphene using the GAP model46 and, most recently, the GAP-20 potential for various crystalline phase carbon and amorphous carbon,50 including dispersion corrections and long-range interactions. Although current ML-IAPs have shown excellent accuracy in predicting static and dynamic properties of carbon allotropes, there is still a clear lack of testing related to the prediction of mechanical properties.

Although previous work has compared classic potential and DFT, they are limited to amorphous carbon and elastic properties of graphene.53,54 Thus, in this work, we compare GAP-20 and widely used CBOPs (Tersoff, mo-Tersoff, AIREBO, mo-AIREBO, ReaxFFC2013, EDIP, and LCBOP) with DFT in terms of lattice constants, cohesive energies, defect properties, van der Waals (vdW) interaction, thermal stability, and mechanical properties for the carbon allotropes diamond, graphite, graphene, and CNTs. As a result of our comprehensive testing, we propose the most suitable potential for simulating crystalline phases of carbon.

The crystalline phases of carbon studied in this work are shown in Fig. 1. We employ periodic boundary conditions (PBCs) for all systems, with 15 Å of vacuum separating the graphene sheets in the z-direction and at least 10 Å of vacuum perpendicular to the (10,0) and (5,5) CNT axis to prevent self-interaction. C60 is placed in a cubic simulation box with a side length of 20 Å.

FIG. 1.

The structures of crystalline phase carbon in the present work are (a) graphite, (b) graphene, (c) C60, (d) (10,0) CNT, (e) (5,5) CNT, and (f) diamond.

FIG. 1.

The structures of crystalline phase carbon in the present work are (a) graphite, (b) graphene, (c) C60, (d) (10,0) CNT, (e) (5,5) CNT, and (f) diamond.

Close modal

We select several widely used and representative potentials for testing, summarized in Table I. mo-AIREBO is the modified AIREBO where the C–C bond cutoff is set to 2.0 Å, and this form is mostly used to simulate the fracture behavior of carbon materials as mentioned above.

TABLE I.

Assessed classical potentials for carbon. The column Elements denotes the elements for which the potential is parameterized, and vdW denotes the inclusion of van der Waals interactions.

PotentialElementsvdWReferencesPotentialElementsvdWReferences
Tersoff No 25  ReaxFFC2013 Yes 34  
AIREBO C/H Yes 28  EDIP C/Si No 36  
LCBOP Yes 31  mo-AIREBO C/H Yes 37  
mo-Tersoff No 32  GAP-20 Yes 50  
PotentialElementsvdWReferencesPotentialElementsvdWReferences
Tersoff No 25  ReaxFFC2013 Yes 34  
AIREBO C/H Yes 28  EDIP C/Si No 36  
LCBOP Yes 31  mo-AIREBO C/H Yes 37  
mo-Tersoff No 32  GAP-20 Yes 50  

All simulations with the potentials listed in Table I were performed using the large-scale atomic/molecular massively parallel simulator (LAMMPS).55 Lattice relaxation was performed on the PBC direction of all models with the conjugate gradient (CG) algorithm, and the convergence of energy and force is set to 10−10 eV and 10−8 eV/Å, respectively. In the tensile test, tensile strain is applied in the axial direction of the CNT and the armchair and zigzag directions of graphene, and the cell stress is obtained after relaxation. The fully relaxed structure is then taken as the initial structure for the next strain step. Here, the effective thickness of graphene and the CNT wall is considered to be 3.34 Å to get effective area Seff = l0 × 3.34 Å2 and effective stress σeff=σbulk×S0Seff, where l0 represents the graphene edge length and circumference of carbon nanotubes and S0 is the area of the simulation box perpendicular to the stretching direction.18 In the thermal stability test, the Newton equations of motion were integrated using the Verlet algorithm with a time step of 0.5 fs. The temperature was set and maintained at 100–5000 K using the Nosé–Hoover thermostat in the NVT ensemble.56,57 The entire system is simulated for 10 ps where the potential energy for the last 5 ps is collected and analyzed.

All DFT calculations were carried out using the Vienna ab initio simulation package (VASP).58,59 The generalized gradient approximation in the Perdew–Burke–Ernzerhof functional is used for the exchange-correlation function.60 DFT-D3 was also introduced to ensure that the van der Waals was correctly described.61 The plane-wave energy cutoff is 400 and 600 eV for structure optimization and van der Waals interaction calculation, respectively (the detailed plane-wave energy cutoff test is summarized in Table S1). The SCF convergence criterion was set to 10−4 eV and relaxation was performed until the forces acting on the atoms was smaller than 10−2 eV/Å. The Brillouin zone was represented by a Γ centered k-point mesh with a grid spacing of 2π × 0.04 Å−1. Furthermore, a 10 ps long ab initio molecular dynamics simulation, with a time step of 1.0 fs, in the NVT ensemble56,62 was employed to test the thermal stability of C60.

We begin by employing the potentials summarized in Table I to predict the lattice constants of CNTs, graphene, graphite, and diamond. The calculated lattice constants are presented in Table II. To quantify the potentials relative to the DFT, we calculate the relative error and the accuracy grade, PLC.

TABLE II.

Calculated lattice constants of crystalline phase carbon in the present work.

Lattice constants (Å) (error, %)
Potential(5,5) CNT(10,0) CNTGraphene-ZZGraphene-ACGraphite-cDiamondPLC (%)
Tersoff 2.527 (2.78) 4.360 (2.17) 2.530 (2.71) 4.382 (2.62) ⋯ 3.566 (0.50) 97.84 
mo-Tersoff 2.506 (1.93) 4.317 (1.18) 2.492 (1.18) 4.316 (1.09) ⋯ 3.645 (2.73) 98.38 
AIREBO 2.432 (1.09) 4.184 (1.96) 2.419 (1.78) 4.190 (1.87) 6.716 (3.16) 3.557 (0.24) 98.32 
mo-AIREBO 2.432 (1.09) 4.184 (1.96) 2.419 (1.78) 4.190 (1.87) 6.716 (3.16) 3.557 (0.24) 99.32 
LCBOP 2.461 (0.10) 4.246 (0.49) 2.459 (0.15) 4.260 (0.24) 6.692 (2.79) 3.565 (0.48) 99.29 
EDIP 2.582 (5.02) 4.458 (4.48) 2.595 (5.34) 4.494 (5.25) ⋯ 3.561 (0.37) 95.91 
ReaxFFC2013 2.469 (0.42) 4.247 (0.48) 2.462 (0.03) 4.265 (0.13) 6.520 (0.14) 3.578 (3.58) 99.66 
GAP-20 2.461 (0.09) 4.251 (0.37) 2.465 (0.07) 4.269 (0.02) 6.706 (3.0) 3.593 (1.26) 99.20 
PBE-D3 2.459 4.267 2.463 4.270 6.510 3.548 100.0 
Experiment ⋯ ⋯ 2.460 (Ref. 634.260 (Ref. 636.672 (Ref. 643.566 (Ref. 6 
Lattice constants (Å) (error, %)
Potential(5,5) CNT(10,0) CNTGraphene-ZZGraphene-ACGraphite-cDiamondPLC (%)
Tersoff 2.527 (2.78) 4.360 (2.17) 2.530 (2.71) 4.382 (2.62) ⋯ 3.566 (0.50) 97.84 
mo-Tersoff 2.506 (1.93) 4.317 (1.18) 2.492 (1.18) 4.316 (1.09) ⋯ 3.645 (2.73) 98.38 
AIREBO 2.432 (1.09) 4.184 (1.96) 2.419 (1.78) 4.190 (1.87) 6.716 (3.16) 3.557 (0.24) 98.32 
mo-AIREBO 2.432 (1.09) 4.184 (1.96) 2.419 (1.78) 4.190 (1.87) 6.716 (3.16) 3.557 (0.24) 99.32 
LCBOP 2.461 (0.10) 4.246 (0.49) 2.459 (0.15) 4.260 (0.24) 6.692 (2.79) 3.565 (0.48) 99.29 
EDIP 2.582 (5.02) 4.458 (4.48) 2.595 (5.34) 4.494 (5.25) ⋯ 3.561 (0.37) 95.91 
ReaxFFC2013 2.469 (0.42) 4.247 (0.48) 2.462 (0.03) 4.265 (0.13) 6.520 (0.14) 3.578 (3.58) 99.66 
GAP-20 2.461 (0.09) 4.251 (0.37) 2.465 (0.07) 4.269 (0.02) 6.706 (3.0) 3.593 (1.26) 99.20 
PBE-D3 2.459 4.267 2.463 4.270 6.510 3.548 100.0 
Experiment ⋯ ⋯ 2.460 (Ref. 634.260 (Ref. 636.672 (Ref. 643.566 (Ref. 6 

As per the results in Table II, most of the potentials tested can accurately determine lattice constants of the carbon allotropes tested. GAP-20, LCBOP, and ReaxFFC2013 optimized CNTs and graphene with an error of less than 0.5% relative to DFT. However, for the diamond lattice, GAP-20 and ReaxFFC2013 are not as accurate as the Tersoff, AIREBO, LCBOP, and EDIP potentials. The lattice constant in the c direction of graphite, dc, represents the double interlayer distance in bulk graphite and is mainly affected by long-range interactions. Here, relaxation with ReaxFFC2013 results in a value close to that of PBE-D3, indicating that the equilibrium position of the long-range interaction is similar, while GAP-20, mo-AIREBO, and AIREBO agree within ∼3%. Compared to the other lattice constants, the choice of exchange-correlation functional plays a larger role, using optB88-vdW50 results in a dc of 6.65 Å, higher than 6.51 Å of PBE-D3 and closer to the values obtained by GAP-20 and the AIREBO potentials rather than ReaxFFC2013. Absent values for the Tersoff potentials and EDIP are due to them not containing van der Waals (vdW) interactions; hence, they should not be used alone to optimize graphite. Many previous works have used additional terms such as Lennard-Jones or Kolmogorov–Crespi (K–C) to solve the problem of vdW interactions.65 In general, all potentials can predict the lattice constant of crystalline carbon materials very well. GAP-20 and ReaxFFC2013 are the best of the ones tested here with overall accuracies of 99.20% and 99.66%, respectively.

The cohesive energy describes the aggregation of atoms to form a crystal. It is also a key parameter for phase transition simulations, where its magnitude affects the difficulty of a phase transition. In this work, the cohesive energy per atom, Ec, is calculated by the following equation:

Ec=EnEatom/n,
(1)

where E is the total energy of crystalline phase carbon, Eatom is the energy of a single carbon atom in vacuum, and n is the number of carbon atoms in the crystalline phase. The cohesive energy per atom calculated by Eq. (1) for all potentials listed in Table I is summarized in Fig. 2(a). The experimentally measured cohesive energy of graphite and diamond is −7.374 and −7.346 eV, respectively,66 and for most of the carbon allotropes tested here, the cohesive energy calculated using the CBOPs is around −7.4 eV. ReaxFFC2013 predicts, however, cohesive energies of around −5.0 eV due to the high single atom energy of Eatom = −2.06 eV, which is much larger than that of DFT. Compared to the CBOPs, the cohesive energies obtained using GAP-20 are much closer to the DFT results, giving it a high accuracy of about 95%.

FIG. 2.

(a) Cohesive energy of carbon allotropes for each potential. (b) Formation energy of carbon allotropes relative to graphite.

FIG. 2.

(a) Cohesive energy of carbon allotropes for each potential. (b) Formation energy of carbon allotropes relative to graphite.

Close modal

To compare the accuracy in energy for the different potentials, we calculated the formation energy of carbon allotropes relative to graphite,

Ef=EnEgraphite/n,
(2)

where E is the total energy of the carbon allotrope, Egraphite is the energy of a single atom in graphite, and n is the number of carbon atoms in the allotrope. As shown in Fig. 2(b), the formation energy, as calculated by Eq. (2), has the following order for most potentials: graphite < graphene < diamond < (10,0) CNT < (5,5) CNT < C60, in agreement with previous work.50 

Graphene and graphite have similar thermodynamically stable phases, where the formation energy of graphene (0.05 eV per atom) is lowest compared to diamond, C60, and the CNTs. Due to the lack of long-range interaction, CBOPs, such as Tersoff and EDIP, cannot correctly distinguish between graphite and graphene. For diamond, the GAP-20 result (0.12 eV) is closest to DFT (0.09 eV), while all other potentials (except mo-Tersoff) underestimate the formation energy of diamond. mo-Tersoff shows an abnormally high formation energy of 1.4 eV per atom.

For the (5,5) and (10,0) CNTs, the GAP-20 and ReaxFFC2013 results agree very well with DFT, both around 0.18 and 0.23 eV, respectively. The other potentials underestimate the value compared to DFT. EDIP calculates Ef to be almost 0 eV for both CNTs, which is most likely because it was parameterized to describe Si–C systems, not one-dimensional carbon. For C60, the formation energy calculated by LCBOP, ReaxFFC2013, and GAP-20 is around 0.4 eV, close to the DFT value of 0.43 eV. EDIP also underestimates the formation energy of C60 by 0.2 eV per atom, while Tersoff, mo-Tersoff, and AIREBO overestimate it.

In general, the CBOPs tested here overestimate the cohesive energy of carbon allotropes, while GAP-20 is in good agreement with DFT. Although some CBOPs, such as LCBOP (for C60) and ReaxFFC2013 (for CNTs and graphene), can predict formation energy accurately, they exhibit an inability to generalize as they are specifically designed for a particular system or set of systems.

Carbon materials, especially low-dimensional ones, often have defects as a result of their synthesis process.67,68 With the advancement of simulation methods, it is possible to explore defect formation and healing mechanisms,69–71 but it is important to have an accurate potential for the application. Therefore, in this part, we compare the defect formation energy and geometry of different defects using the potentials in Table I.

The defect formation energy, Edf, is obtained using the following equation:

Edf=EdefectmEatom+Eperfect,
(3)

where Edefect is the energy of the defective structure, Eatom is the energy of a single atom in perfect graphene, Eperfect is the total energy of the defect-free structure, and m is the number of the introduced (positive) or removed (negative) atoms.

A 10 × 6 × 1 supercell of graphene is used to eliminate the defect self-interaction, and for each perfect structure, the lattice constants and atomic degrees of freedom are fully relaxed until the structure optimization reaches convergence. For the defective structures, only the atomic degrees of freedom are relaxed until convergence is reached.

Single vacancy (SV) is the most common point defect in graphene, where one atom in the graphene lattice is missing. It has been observed in experiments via STM72 and TEM.73 The structure and formation energy of the SV defect, relaxed with different potentials, are shown in Fig. 3 and Fig. S1. The formation energy for one SV defect obtained by DFT is 7.74 eV, which agrees well with the 7.7 eV obtained using the optB88-vdW functional report by Rowe et al.50 AIREBO and LCBOP result in 7.65 and 7.60 eV, respectively, and the GAP-20 and EDIP are also within 1 eV of DFT, all lower than the reference value. Surprisingly, mo-Tersoff significantly underestimates the SV formation energy, resulting in 0.53 eV, and ReaxFFC2013 overestimated it by 2.3 eV, resulting in 10.2 eV. Thus, ReaxFFC2013 and mo-Tersoff are not suitable for describing the SV defects in graphene.

FIG. 3.

The defect geometry and formation energy of a single vacancy defect relaxed by (a) PBE, (b) GAP-20, (c) ReaxFFC2013, and (d) LCBOP.

FIG. 3.

The defect geometry and formation energy of a single vacancy defect relaxed by (a) PBE, (b) GAP-20, (c) ReaxFFC2013, and (d) LCBOP.

Close modal

The geometry and formation energy are important properties for defects. To quantify the difference in structure, consider the symmetry of the SV defect, and we select three atom distances (d1, d2, and d3) shown in Fig. 3(a) to characterize the defect geometry. The results are summarized in Table III and show that d3 (2.46 Å for the defect-free structure) varies between 2.5 and 2.6 Å for most potentials, indicating that the three adjacent atoms do not undergo large displacement after one atom is removed. In contrast, d3 for ReaxFFC2013 and mo-Tersoff is larger than 2.9 Å. This noticeable difference causes the vacancy by the missing atom to become more triangular.

TABLE III.

Defect formation energy and atomic distances around the single vacancy. Here, the three atomic distances d1, d2, and d3 characterize the defect geometry.

Atomic distance (Å)
PotentialsEdf (eV)d1d2d3PSV (%)
Tersoff 6.89 1.456 1.445 2.613 93.00 
mo-Tersoff 0.53 1.407 1.351 3.034 49.57 
AIREBO 7.65 1.384 1.395 2.557 98.92 
mo-AIREBO 7.65 1.384 1.395 2.557 98.92 
LCBOP 7.60 1.407 1.397 2.671 98.10 
EDIP 6.88 1.505 1.548 2.367 90.42 
ReaxFFC2013 10.23 1.426 1.346 2.936 80.68 
GAP-20 6.99 1.426 1.390 2.492 94.67 
PBE-D3 7.74 1.418 1.394 2.544 100.0 
Atomic distance (Å)
PotentialsEdf (eV)d1d2d3PSV (%)
Tersoff 6.89 1.456 1.445 2.613 93.00 
mo-Tersoff 0.53 1.407 1.351 3.034 49.57 
AIREBO 7.65 1.384 1.395 2.557 98.92 
mo-AIREBO 7.65 1.384 1.395 2.557 98.92 
LCBOP 7.60 1.407 1.397 2.671 98.10 
EDIP 6.88 1.505 1.548 2.367 90.42 
ReaxFFC2013 10.23 1.426 1.346 2.936 80.68 
GAP-20 6.99 1.426 1.390 2.492 94.67 
PBE-D3 7.74 1.418 1.394 2.544 100.0 

To measure the accuracy of the potentials compared to DFT, we combine the formation energy and geometry to estimate the accuracy of a single vacancy defect relative to DFT, PSV. In Table III, the formation energy and geometry are weighted equally at 50%. Here, AIREBO and LCBOP show very high PSV of 98.92% and 98.10%, respectively, followed by GAP-20, Tersoff, and EDIP around 90%, while ReaxFFC2013 and mo-Tersoff have a lower PSV.

Double vacancy (DV) defect is the removal of two adjacent atoms. As for the SV defect, we show the defect geometry in Fig. 4 and Fig. S2 and summarize the formation energy and the atom distance around the defect in Table IV.

FIG. 4.

The defect geometry and formation energy of the double vacancy obtained by (a) PBE, (b) GAP-20, (c) ReaxFFC2013, and (d) LCBOP.

FIG. 4.

The defect geometry and formation energy of the double vacancy obtained by (a) PBE, (b) GAP-20, (c) ReaxFFC2013, and (d) LCBOP.

Close modal
TABLE IV.

Defect formation energy and atomic distance around the double vacancy. Here, the three atomic distances d1, d2, and d3 characterize the defect geometry.

Atomic distance (Å)
PotentialEdf (eV)d1d2d3PDV (%)
Tersoff 9.19 1.454 1.444 2.611 82.68 
mo-Tersoff 0.76 1.390 1.348 3.010 43.30 
AIREBO 10.24 1.376 1.392 2.551 74.83 
mo-AIREBO 10.24 1.376 1.392 2.551 74.83 
LCBOP 10.01 1.396 1.395 2.653 75.72 
EDIP 9.19 1.509 1.550 2.372 83.19 
ReaxFFC2013 13.40 1.496 1.308 2.903 50.19 
GAP-20 7.99 1.440 1.476 1.897 96.70 
PBE-D3 7.55 1.451 1.455 1.895 100.00 
Atomic distance (Å)
PotentialEdf (eV)d1d2d3PDV (%)
Tersoff 9.19 1.454 1.444 2.611 82.68 
mo-Tersoff 0.76 1.390 1.348 3.010 43.30 
AIREBO 10.24 1.376 1.392 2.551 74.83 
mo-AIREBO 10.24 1.376 1.392 2.551 74.83 
LCBOP 10.01 1.396 1.395 2.653 75.72 
EDIP 9.19 1.509 1.550 2.372 83.19 
ReaxFFC2013 13.40 1.496 1.308 2.903 50.19 
GAP-20 7.99 1.440 1.476 1.897 96.70 
PBE-D3 7.55 1.451 1.455 1.895 100.00 

As shown in Figs. 4(a) and 4(b), after optimization by DFT and GAP-20, the DV defect will form a 5-8-5 structure where the length of d3 in the pentagon is close to 1.89 Å. For the other potentials, the length of d3 is greater than 2.3 Å, which exceeds the C–C bond length. Except for d3, the bond lengths of d1 and d2 for the CBOPs are close to the bond length obtained via DFT. For the DV defect, the formation energy obtained by DFT is 7.55 eV, which is lower than for the SV defects. However, the formation energy obtained by most potentials, except mo-Tersoff, is higher than the DFT reference value.

The accuracy of the potentials for the DV defect, PDV (evaluated with the same criteria as for the SV defect), is presented in Table IV. Here, AIREBO and LCBOP, which perform best for the SV defect, cannot accurately model the DV defect nor can any other potential apart from GAP-20. Achieving a high accuracy of 96.70% as it reproduces the 5-8-5 structure obtained with DFT, the key relationship between energy and structure implies that only GAP-20 is sufficient among the potentials assessed here to model DVs in graphene.

Stone–Wales defects, unlike SV and DV defects, have no missing atoms. Instead, the graphene lattice has been reconstructed such that one C–C bond is rotated by 90°, which results in four hexagons transformed into two pentagons and two heptagons (5-7-7-5), as shown in Fig. 5 and Fig. S3. The formation energy Edf, bond lengths d1 and d2, and angle θ, respectively, as well as the accuracy relative to DFT PSW are presented in Table V.

FIG. 5.

Atomic structures and formation energy of the SW defect optimized by (a) PBE, (b) GAP-20, (c) AIREBO, and (d) mo-AIREBO. Interestingly, mo-AIREBO cannot reproduce the 5-7-7-5 structure, resulting in high formation energy.

FIG. 5.

Atomic structures and formation energy of the SW defect optimized by (a) PBE, (b) GAP-20, (c) AIREBO, and (d) mo-AIREBO. Interestingly, mo-AIREBO cannot reproduce the 5-7-7-5 structure, resulting in high formation energy.

Close modal
TABLE V.

Defect formation energy, atomic distance, and C–C–C bond angle for the Stone–Wales defect heptagon. Here, the atomic distances d1 and d2 and the bond angle θ characterize the defect geometry and PSW is the accuracy of a Stone–Wales defect.

AtomicC–C–C
distance (Å)angle (deg)
PotentialEdf (eV)d1d2θPSW (%)
Tersoff 11.8 1.522 1.329 120.0 48.82 
mo-Tersoff 6.15 1.521 1.362 118.815 86.63 
AIREBO 6.23 1.459 1.332 115.396 86.40 
mo-AIREBO 27.26 1.732 1.487 103.115 48.09 
LCBOP 5.29 1.444 1.312 121.066 96.68 
EDIP 5.57 1.583 1.369 117.551 91.47 
ReaxFFC2013 4.07 1.446 1.320 119.672 90.20 
GAP-20 5.23 1.471 1.305 121.365 97.47 
PBE-D3 4.99 1.461 1.305 121.266 100.00 
AtomicC–C–C
distance (Å)angle (deg)
PotentialEdf (eV)d1d2θPSW (%)
Tersoff 11.8 1.522 1.329 120.0 48.82 
mo-Tersoff 6.15 1.521 1.362 118.815 86.63 
AIREBO 6.23 1.459 1.332 115.396 86.40 
mo-AIREBO 27.26 1.732 1.487 103.115 48.09 
LCBOP 5.29 1.444 1.312 121.066 96.68 
EDIP 5.57 1.583 1.369 117.551 91.47 
ReaxFFC2013 4.07 1.446 1.320 119.672 90.20 
GAP-20 5.23 1.471 1.305 121.365 97.47 
PBE-D3 4.99 1.461 1.305 121.266 100.00 

The pentagons and heptagons predicted by DFT are like regular polygons, and the SW formation energy is 4.99 eV. Potentials such as GAP-20 and EDIP are close to DFT both in terms of structure and formation energy. However, mo-AIREBO fails to optimize the 5-7-7-5 structure and returns a highly overestimated formation energy, Ef = 27.94 eV, despite mo-AIREBO and AIREBO performing very similarly for the SV and DV defects. mo-AIREBO is unable to describe the SW defect due to the modification of the C–C bond cutoff, making it unsuitable for modeling significant changes in the atomic environment.

Among the potentials compared to DFT, GAP-20 and LCBOP perform the best, while EDIP and ReaxFFC2013 maintain accuracies above 90%. Although the Tersoff potential accurately predicts the SW defect structure, it overestimates the formation energy.

The hexagonal structure of graphene is very stable, and the insertion of additional carbon atoms is energetically expensive. The formation energy and structural information of a single ad-atom defect are presented in Table VI. A single carbon atom may be embedded at the bridge or top position, as shown in Figs. 6(a) and 6(b). All potentials optimize the ad-atom to be embedded at the bridge position, except for AIREBO. The formation energy obtained via DFT is 6.54 eV, while the other potentials range from 4.21 to 9.79 eV. Surprisingly, mo-AIREBO achieved a high accuracy of 92.03%.

TABLE VI.

Formation energy and atomic distances around the single ad-atom defect. Here, the three atomic distances d1, d2, and d3 characterize the defect geometry and Padd-1 is the accuracy of a single atom defect.

Atom distance (Å)
PotentialEdf (eV)d1d2d3Padd-1 (%)
Tersoff 7.99 1.740 1.740 1.700 82.35 
mo-Tersoff 9.79 1.679 1.679 1.716 69.76 
AIREBO 4.23 1.421 2.434 1.526 70.85 
mo-AIREBO 6.18 1.713 1.713 1.629 92.03 
LCBOP 7.32 1.390 1.390 2.329 83.07 
EDIP 4.21 1.606 1.568 1.730 78.67 
ReaxFFC2013 5.1 1.506 1.506 2.248 81.45 
GAP-20 5.25 1.390 1.390 2.007 82.62 
PBE-D3 6.54 1.512 1.512 1.556 100.00 
Atom distance (Å)
PotentialEdf (eV)d1d2d3Padd-1 (%)
Tersoff 7.99 1.740 1.740 1.700 82.35 
mo-Tersoff 9.79 1.679 1.679 1.716 69.76 
AIREBO 4.23 1.421 2.434 1.526 70.85 
mo-AIREBO 6.18 1.713 1.713 1.629 92.03 
LCBOP 7.32 1.390 1.390 2.329 83.07 
EDIP 4.21 1.606 1.568 1.730 78.67 
ReaxFFC2013 5.1 1.506 1.506 2.248 81.45 
GAP-20 5.25 1.390 1.390 2.007 82.62 
PBE-D3 6.54 1.512 1.512 1.556 100.00 
FIG. 6.

Atomic structure and formation energy of the ad-atom optimized by (a) DFT and (b) ARIEBO. Ad-dimer structure and formation energy obtained by (d) DFT and (e) mo-AIREBO. mo-AIREBO cannot correctly predict the change in curvature. Structure characteristics of the (c) ad-atom defect and (f) ad-dimer defect geometry. Here, the red atoms represent the added atoms.

FIG. 6.

Atomic structure and formation energy of the ad-atom optimized by (a) DFT and (b) ARIEBO. Ad-dimer structure and formation energy obtained by (d) DFT and (e) mo-AIREBO. mo-AIREBO cannot correctly predict the change in curvature. Structure characteristics of the (c) ad-atom defect and (f) ad-dimer defect geometry. Here, the red atoms represent the added atoms.

Close modal

Adding two carbon atoms to graphene is more complex since they will first form a dimer and then be embedded in a hexagon to form a new 7-5-5-7 structure, with a change in the graphene curvature, as shown in Figs. 6(d) and 6(e) and Fig. S4. After optimization with DFT, an ad-dimer formation energy of 6.18 eV is obtained, with GAP-20 resulting in a value of 5.91 eV. GAP-20 also reproduces the correct geometry and achieves an accuracy of 96.82%. Similar to the performance of each potential for the single ad-atom, the Tersoff and LCBOP potentials overestimate the formation energy of the ad-dimer, while EDIP and ReaxFFC2013 underestimate it. Interestingly, AIREBO reproduces the structure of the ad-atom and ad-dimer defects but underestimates the formation energy. While mo-AIREBO cannot reproduce the structures but more accurately matches the energies calculated by DFT. The formation energy and structure of the ad-dimer defects are summarized in Table VII.

TABLE VII.

Defect formation energy, atomic distances, and C–C–C bond angles around the ad-dimer defect. Here, the three atomic distances d1, d2, and d3 and the two bond angles θ1 and θ2 characterize the defect geometry and Padd-2 is the accuracy of an ad-dimer defect.

Atom distance (Å)Angle (deg)
PotentialEdf (eV)d1d2d3θ1θ2Padd-2 (%)
Tersoff 7.99 1.640 1.478 1.460 122.46 128.05 82.27 
mo-Tersoff 9.8 1.487 1.704 1.503 60.15 136.80 20.71 
AIREBO 4.23 1.517 1.419 1.419 115.50 131.67 83.37 
mo-AIREBO 6.17 1.686 1.573 1.282 74.72 111.67 58.11 
LCBOP 7.33 1.495 1.417 1.447 119.68 131.64 89.66 
EDIP 3.01 1.558 1.504 1.552 116.76 118.25 24.35 
ReaxFFC2013 4.86 1.406 1.417 1.423 120.13 131.77 88.74 
GAP-20 5.91 1.478 1.392 1.418 121.51 131.17 96.82 
PBE-D3 6.18 1.428 1.410 1.419 116.74 132.53 100.00 
Atom distance (Å)Angle (deg)
PotentialEdf (eV)d1d2d3θ1θ2Padd-2 (%)
Tersoff 7.99 1.640 1.478 1.460 122.46 128.05 82.27 
mo-Tersoff 9.8 1.487 1.704 1.503 60.15 136.80 20.71 
AIREBO 4.23 1.517 1.419 1.419 115.50 131.67 83.37 
mo-AIREBO 6.17 1.686 1.573 1.282 74.72 111.67 58.11 
LCBOP 7.33 1.495 1.417 1.447 119.68 131.64 89.66 
EDIP 3.01 1.558 1.504 1.552 116.76 118.25 24.35 
ReaxFFC2013 4.86 1.406 1.417 1.423 120.13 131.77 88.74 
GAP-20 5.91 1.478 1.392 1.418 121.51 131.17 96.82 
PBE-D3 6.18 1.428 1.410 1.419 116.74 132.53 100.00 

We have shown that CBOPs can describe some but not all defects well. For example, Tersoff hasa 93% accuracy in SV defects, but it performs poorly in describing the larger structural changes associated with SW, resulting in only 48.82% accuracy. The reason for this is that CBOPs often only consider crystal properties at the beginning of their development, with parameters obtained for defect-free structures. Thus, there will be large errors in describing the properties of defects. On the contrary, the training set for GAP-20 contains defect-free and defective structures; therefore, it can more accurately describe the structure and formation energy of defects.

The edges of graphene and CNTs are of great significance as different edges bring about different properties, for example, different edges will affect the growth of graphene.74–76 Here, we define the edge energy, Eedge, as the energy density of the edge given by the following equation:

Eedge=EcutEperfect/2l,
(4)

where Ecut is the energy of the graphene or CNTs with the exposed edge, Eperfect is the energy of the periodic structure, and l is the length (perimeter) of the exposed edge. All edge energies for graphene and CNTs calculated using Eq. (4) are summarized in Table VIII. The armchair edge energy, EedgeAC = 1.01 eV/Å, for graphene is lower than the zigzag edge energy, EedgeZZ = 1.33 eV/Å, as predicted by DFT and consistent with the results obtained by Baran et al., 77 Li et al.,78 and Hedman and Larsson79,80 Among the potentials, only GAP-20 and ReaxFFC2013 successfully predict this trend (EedgeAC < EedgeZZ). EDIP, LCBOP, ReacFFC2013, and GAP-20 are within 5% of the DFT value for the armchair edge energy, but only GAP-20 is within 10% of the DFT calculated zigzag edge energy, showing that a high accuracy in energy is needed to correctly describe the relationship between armchair and zigzag edges. The edges of (10,0) and (5,5) CNTs are zigzag and armchair, respectively, where Eedge(5,5) < Eedge(10,0) as predicted by DFT and consistent with the trend of edge energies in graphene. Only GAP-20 and AIREBO can accurately reproduce this trend for the (5,5) and (10,0) edge energies. For the (5,5) edge energy, EDIP predicts the same value as DFT, while most potentials slightly overestimate it, except for Tersoff and ReaxFFC2013 (which give higher values) and mo-Tersoff (which severely underestimates it). For the (10,0) edge, all potentials underestimate the edge energy with AIREBO, mo-AIREBO, and LCBOP showing almost no difference to the (5,5) edge. Only GAP-20 comes within 20% of the DFT calculated (10,0) edge energy.

TABLE VIII.

Armchair and zigzag edge energies of graphene and (10,0) and (5,5) CNTs, and Ped is the accuracy of the edge energy.

Graphene (eV/Å)CNT (eV/Å)
PotentialArmchairZigzag(5,5)(10,0)Ped (%)
Tersoff 1.07 0.92 1.02 0.88 80.81 
mo-Tersoff 0.72 0.19 0.33 0.15 33.06 
AIREBO 1.18 1.10 0.98 0.99 84.67 
mo-AIREBO 1.18 1.10 0.98 0.99 84.67 
LCBOP 1.02 1.02 0.98 0.98 86.89 
EDIP 1.06 0.90 0.94 0.68 78.83 
ReaxFFC2013 1.24 1.49 1.16 0.99 79.78 
GAP-20 1.11 1.24 0.99 1.05 89.83 
PBE-D3 1.01 1.33 0.94 1.29 100.0 
Graphene (eV/Å)CNT (eV/Å)
PotentialArmchairZigzag(5,5)(10,0)Ped (%)
Tersoff 1.07 0.92 1.02 0.88 80.81 
mo-Tersoff 0.72 0.19 0.33 0.15 33.06 
AIREBO 1.18 1.10 0.98 0.99 84.67 
mo-AIREBO 1.18 1.10 0.98 0.99 84.67 
LCBOP 1.02 1.02 0.98 0.98 86.89 
EDIP 1.06 0.90 0.94 0.68 78.83 
ReaxFFC2013 1.24 1.49 1.16 0.99 79.78 
GAP-20 1.11 1.24 0.99 1.05 89.83 
PBE-D3 1.01 1.33 0.94 1.29 100.0 

For all graphene and CNT edges examined, GAP-20 has an accuracy of 89.83%, while Tersoff, AIREBO, and LCBOP also perform impressively, slightly worse than GAP-20 but still with more than 80% accurate.

The van der Waals interaction dominates the mobility of stacked graphene layers where the direction and difficulty of slip depend on the potential energy surface (PES). Here, the vdW interaction for bi-layer graphene, EvdW, is calculated with the following equation:

EvdW=Ebi2E0/N,
(5)

where Ebi is the energy of bi-layer graphene, E0 is the energy of a single graphene sheet, and N is the number of atoms in the bi-layer system. Both AA stacking and AB stacking are considered, and the vdW interaction for different inter-layer distances, d, is plotted in Fig. 7(a). In the AIREBO potential, the vdW interaction is modeled by a Lennard-Jones term, which leads to an exaggerated inter-layer interaction in the repulsive region and an underestimation in the attractive region. LCBOP also underestimates the attractive part but describes the repulsive part more accurate than AIREBO, especially for the AB stacking. Among the CBOPs, the curve of ReaxFFC2013 in Fig. 7(a) is always below the DFT curve, indicating that it underestimates the repulsive interaction and overestimates the attractive interaction. Notably, the GAP-20 potential has two local minima for both the AA and AB stacking, as shown in Fig. 7(a), which is far from the optB88-vdw results used in the GAP-20 training set, indicating that the potential may have been overfitted during training. Figure 7(b) shows the energy difference between AA stacking and AB stacking, EvdWAAAB = EvdWAAEvdWAB, where EvdWAAAB exhibits a continuous decrease for DFT and CBOPs, while GAP-20 shows a non-continuous decrease due to the two-local minima. In addition, the EAA−AB obtained by the CBOPs drops to almost 0 meV at 3.2 Å. This small energy difference indicates that the CBOPs can no longer distinguish between AA stacking and AB stacking in bilayer graphene after this distance compared to PBE-D3 and optB88-vdw, which fall to 0 meV after 4 Å.

FIG. 7.

(a) van der Waals interaction for AA stacking and AB stacking as a function of inter-layer distance, d. (b) Difference between the van der Waals interaction for AA and AB stacking. (c)–(e) PES for bilayer graphene obtained via PBE+D3, AIREBO, and GAP-20 at a fixed 3.4 Å interlayer distance. The maximum region (red color) and minimal region (blue color) represent the AA stacking and AB stacking, respectively. Note the difference in scale for (c)–(e).

FIG. 7.

(a) van der Waals interaction for AA stacking and AB stacking as a function of inter-layer distance, d. (b) Difference between the van der Waals interaction for AA and AB stacking. (c)–(e) PES for bilayer graphene obtained via PBE+D3, AIREBO, and GAP-20 at a fixed 3.4 Å interlayer distance. The maximum region (red color) and minimal region (blue color) represent the AA stacking and AB stacking, respectively. Note the difference in scale for (c)–(e).

Close modal

The calculated PES for bilayer graphene is reported in Figs. 7(c)7(e) and Fig. S5 for a fixed interlayer distance of 3.4 Å. It is worth noting that when using the current DFT-D3 method for PES scanning, the global minimum is distributed around AB stacking. The shape of the PES is reproduced with the CBOPs, but with an order of magnitude lower energy barrier. A lower energy barrier makes the layers in graphite more prone to sliding. Although GAP-20 produces a barrier height similar to DFT, the shape of the PES is very different compared to DFT and the CBOPs. The sliding path between AB stacking and AA stacking has a higher barrier, which causes the layers in graphite to be locked. The CBOPs tested here benefit from an added long-range interaction, but they are still relatively inaccurate and cannot distinguish between AA stacking and AB stacking in bilayer graphene beyond 3.2 Å separation. Although the training set used for the GAP-20 potential was created using optB88-vdw, the trained potential cannot accurately reproduce the PES of bi-layer graphene present here.

MD simulations with temperature ramping were performed on C60 to compare the temperature responses between the potentials and DFT. The potential energy per atom has a two-stage relationship with temperature, as shown in Fig. 8(a). First, the potential energy increases linearly with temperature and then abruptly rises, indicating the breaking of C–C bonds and the simultaneous release of energy. The temperature at which this occurs is considered as the melting point and is summarized in Fig. 8(b). As seen in Fig. 8(a), mo-Tersoff exhibits an abnormal relationship where the potential energy abruptly falls instead of rising, indicating that energy is gained when C–C bonds are broken, a clearly unphysical behavior, making it unsuitable for thermal studies. DFT predicts that C60 remains intact at 4500 K and decomposes at 5000 K, with a slightly deformed structure at 4500 K, as shown in Fig. 8(c). Most potentials assessed here predict C60 melting at 4000 K, with mo-AIREBO and mo-Tersoff being unsuitable for thermal studies due to changes in their original parameters causing premature melting. From a structural point of view, mo-AIREBO and mo-Tersoff are very similar, showing C60 decomposing to amorphous carbon at 2500 and 2000 K, respectively. As the temperature rises, the structure turns to a long carbon chain at 3000 K. GAP-20 results in C60 melting and forming carbon chains above 4500 K. The atomic structures for the rest of the CBOPs are given in Fig. S6.

FIG. 8.

(a) Potential energy of C60 as a function of temperature. (b) Melting point of C60 as predicted by the potentials and DFT. (c) Atomic structures of C60 at 10 ps for different temperatures with DFT, GAP-20, mo-Tersoff, and mo-AIREBO.

FIG. 8.

(a) Potential energy of C60 as a function of temperature. (b) Melting point of C60 as predicted by the potentials and DFT. (c) Atomic structures of C60 at 10 ps for different temperatures with DFT, GAP-20, mo-Tersoff, and mo-AIREBO.

Close modal

Here, we employ GAP-20 and CBOPs to perform uniaxial tensile tests on graphene and (5,5) and (10,0) CNTs to screen the most suitable potential for modeling mechanical properties compared to DFT. We consider the strain energy, calculated via

ΔE=EsE0/n,
(6)

where Es and E0 are the potential energy with and without the strain and n is the number of atoms in the structure.

Tensile strain is applied along the axial direction of (5,5) and (10,0) CNTs, as shown in Figs. 9(a) and 9(d). From the DFT obtained stress–strain curve of the (5,5) CNT in Fig. 9(b), a fracture strength of 90 GPa at 33% strain is observed. GAP-20 most closely reproduces this behavior, making the (5,5) CNT completely fracture at 26% strain with a strength of 120 GPa, close to the theoretical and experimental results reported by Zhang et al.81 and Lourie et al.82 The AIREBO potential exhibits two stages in the (5,5) tensile test; before 30% strain, the stress–strain curve is similar to DFT. Then, the stress unphysically increases until the tube fractures and the same phenomenon can be seen for the other CBOPs (dashed lines). Changing the C–C bond cutoff like in the mo-AIREBO results in the normal fracture behavior. As seen in Fig. 9(e), the mechanical properties of the (10,0) tube are like that of the (5,5) tube where the facture strain predicted by DFT and GAP-20 is around 20%. For both the (5,5) and (10,0) tubes, mo-AIREBO shows a slow decrease in the stress before fracture, indicating the breaking of individual C–C bonds. At the point of fracture, the strain energy reaches a maximum and then abruptly decreases, as shown in Figs. 9(c) and 9(f), indicating that the structure of the CNT has been destroyed. Snapshots of the CNT after fracture are plotted in Fig. 10, showing that the CBOPs cannot accurately reproduce the DFT fracture results. Unlike the CBOPs, GAP-20 can accurately describe the C–C bond break and fracture of the CNT into halves without overstretching the CNT, as observed for CBOPs.

FIG. 9.

Schematic diagrams for the uniaxial tensile test of (a) (5,5) and (d) (10,0) CNTs. Stress–strain (b) and (c) and strain energy (e) and (f) curves for the (5,5) and (10,0) CNTs under the tensile load. More detailed information for some potentials can be found in Figs. S7 and S8.

FIG. 9.

Schematic diagrams for the uniaxial tensile test of (a) (5,5) and (d) (10,0) CNTs. Stress–strain (b) and (c) and strain energy (e) and (f) curves for the (5,5) and (10,0) CNTs under the tensile load. More detailed information for some potentials can be found in Figs. S7 and S8.

Close modal
FIG. 10.

Snapshots of (a) (5,5) and (b) (10,0) CNTs after fracture and the initial structure before the tensile test with different potentials. The CNTs show a clear fracture for DFT and GAP-20, while the CBOPs give unphysical results.

FIG. 10.

Snapshots of (a) (5,5) and (b) (10,0) CNTs after fracture and the initial structure before the tensile test with different potentials. The CNTs show a clear fracture for DFT and GAP-20, while the CBOPs give unphysical results.

Close modal

The tensile test procedure for graphene is slightly different than for the CNTs. Here, the zigzag (armchair) direction will be relaxed when the strain is acting along on the armchair (zigzag) direction to ensure uniaxial load, as shown in Fig. 11(a). When the strain acts in the armchair direction, Figs. 11(b) and 11(c), DFT predicts fracture at 25% strain with a strength of 115 GPa, which agrees well with previous published results (130 and 121 GPa for experimental and DFT results, respectively).18,83 Similar to the CNT tensile tests, all CBOPs aside from mo-AIREBO show an unphysical increase in stress before fracture and cannot accurately reproduce the mechanical behavior of graphene. GAP-20 agrees well with DFT, with a minor overestimation of the fracture strength, while mo-AIREBO overestimates it and predicts fracture to occur at almost 40% strain along the armchair direction rather than 25%. Figures 11(e) and 11(f) show similar results for the zigzag direction where the facture strain predicted by GAP-20 is 20%, very close to the DFT-predicted value of 23%. GAP-20 slightly overestimates the fracture strength, 135.5 GPa compared to 101.1 GPa for DFT, while mo-AIREBO overestimates it by 15.2%. Overall, GAP-20 agrees well with the DFT strain energy curves in Figs. 11(c) and 11(f).

FIG. 11.

Schematic diagrams of the uniaxial tensile test acting along the (a) armchair and (d) zigzag direction. Stress–strain (b) and (c) and strain energy (e) and (f) curves for graphene under the tensile load along the armchair and zigzag direction. More detailed information for some potentials can be found in Figs. S9 and S10.

FIG. 11.

Schematic diagrams of the uniaxial tensile test acting along the (a) armchair and (d) zigzag direction. Stress–strain (b) and (c) and strain energy (e) and (f) curves for graphene under the tensile load along the armchair and zigzag direction. More detailed information for some potentials can be found in Figs. S9 and S10.

Close modal

The structures after fracture are plotted in Figs. 12(a) and 12(b), revealing that the CBOPs cause graphene to be over-stretched, thus overestimating the ductility and incorrectly preventing fracture. Tersoff and mo-Tersoff result in a structural transformation and deformation, while ReaxFFC2013 sees the formation of multiple defects. In contrast, GAP-20 caused the graphene to undergo brittle fracture at a reasonable strain, consistent with DFT, while mo-AIREBO only reproduced a DFT-like behavior for the uniaxial tensile test in the armchair direction of graphene.

FIG. 12.

Snapshots of graphene after fracture when the tensile load is applied along the (a) zigzag and (b) armchair direction of graphene with DFT and different potentials.

FIG. 12.

Snapshots of graphene after fracture when the tensile load is applied along the (a) zigzag and (b) armchair direction of graphene with DFT and different potentials.

Close modal

To evaluate the accuracy of the predicted mechanical properties for the various potentials, the fracture strain, fracture strength, critical strain energy, and Young’s modulus are summarized in Tables IX and X. One of the reasons why CBOPs cannot model the fracture behavior correctly is due to the C–C bond cutoff, which defines the extent of the covalent bond. This causes C–C bonds that should be broken during the tensile test to still maintain a covalent interaction. We demonstrate that CBOPs cannot correctly describe the breaking of C–C bonds, which leads to unphysical results. However, they still predict mechanical properties such as Young’s modulus, within the elastic deformation region, close to previous experimental results (1.0 TPa) and theoretical predictions (1.04 TPa) for graphene.83–87 The GAP-20 training set contains structures placed under varying degrees of stress, which is sufficient for the trained potential to learn the bond-breaking process. Furthermore, mo-AIREBO and GAP-20 produce results similar to DFT in the elastic and plastic regions, making them suitable for simulating mechanical properties, especially the fracture behavior. However, in combination with the defects and thermal stability test results shown here, GAP-20 is the more suitable potential for simulating mechanical properties of the defect system at high temperatures.

TABLE IX.

Mechanical properties of CNTs, including failure strain, failure strength, critical strain energy, and Young’s modulus.

Fracture strainFracture strengthStrain energyYoung’s
(%)(GPa)(eV)modulus (TPa)
Potential(5,5)(10,0)(5,5)(10,0)(5,5)(10,0)(5,5)(10,0)PCNT (%)
Tersoff 45 55 293.4 299.0 3.27 4.11 0.984 1.14 28.37 
mo-Tersoff 29 56 202.0 377.7 1.44 3.32 0.852 0.960 48.05 
AIREBO 43 51 369.5 391.4 2.71 2.71 0.757 1.010 32.91 
mo-AIREBO 34 19 135.3 112.2 1.59 1.39 0.820 1.010 69.34 
LCBOP 55 55 249.1 407.1 3.12 3.88 0.734 0.883 26.66 
EDIP 45 46 230.8 252.7 2.54 2.49 0.852 0.988 34.97 
ReaxFFC2013 33 39 159.3 183.8 1.59 1.79 0.764 0.825 53.81 
GAP-20 26 20 120.0 129.45 1.08 0.76 0.732 0.951 82.54 
PBE-D3 33 22 90.0 101.71 1.42 0.86 0.839 0.955 100.00 
Fracture strainFracture strengthStrain energyYoung’s
(%)(GPa)(eV)modulus (TPa)
Potential(5,5)(10,0)(5,5)(10,0)(5,5)(10,0)(5,5)(10,0)PCNT (%)
Tersoff 45 55 293.4 299.0 3.27 4.11 0.984 1.14 28.37 
mo-Tersoff 29 56 202.0 377.7 1.44 3.32 0.852 0.960 48.05 
AIREBO 43 51 369.5 391.4 2.71 2.71 0.757 1.010 32.91 
mo-AIREBO 34 19 135.3 112.2 1.59 1.39 0.820 1.010 69.34 
LCBOP 55 55 249.1 407.1 3.12 3.88 0.734 0.883 26.66 
EDIP 45 46 230.8 252.7 2.54 2.49 0.852 0.988 34.97 
ReaxFFC2013 33 39 159.3 183.8 1.59 1.79 0.764 0.825 53.81 
GAP-20 26 20 120.0 129.45 1.08 0.76 0.732 0.951 82.54 
PBE-D3 33 22 90.0 101.71 1.42 0.86 0.839 0.955 100.00 
TABLE X.

Mechanical properties of graphene, including failure strain, failure strength, critical strain energy, and Young’s modulus.

Fracture strainFracture strengthStrain energyYoung’s
(%)(GPa)(eV)modulus (TPa)
PotentialZZACZZACZZACZZACPGRA (%)
Tersoff 56 45 292.1 342.4 4.25 3.51 1.09 1.10 18.55 
mo-Tersoff 50 45 318.5 354.8 2.81 1.91 0.99 0.48 18.82 
AIREBO 51 47 415.3 454.2 2.78 3.04 0.98 0.85 22.09 
mo-AIREBO 20 38 116.4 155.1 0.74 3.28 0.88 0.85 68.33 
LCBOP 57 58 420.8 292.7 4.17 0.85 0.95 0.85 28.53 
EDIP 47 39 243.0 261.2 2.67 0.85 0.95 0.92 40.37 
ReaxFFC2013 39 45 158.6 212.7 1.56 0.92 0.65 0.92 49.59 
GAP-20 20 24 135.5 144.9 0.77 0.90 0.92 0.90 86.37 
PBE-D3 23 25 101.0 115.2 0.91 0.96 0.96 0.96 100.00 
Fracture strainFracture strengthStrain energyYoung’s
(%)(GPa)(eV)modulus (TPa)
PotentialZZACZZACZZACZZACPGRA (%)
Tersoff 56 45 292.1 342.4 4.25 3.51 1.09 1.10 18.55 
mo-Tersoff 50 45 318.5 354.8 2.81 1.91 0.99 0.48 18.82 
AIREBO 51 47 415.3 454.2 2.78 3.04 0.98 0.85 22.09 
mo-AIREBO 20 38 116.4 155.1 0.74 3.28 0.88 0.85 68.33 
LCBOP 57 58 420.8 292.7 4.17 0.85 0.95 0.85 28.53 
EDIP 47 39 243.0 261.2 2.67 0.85 0.95 0.92 40.37 
ReaxFFC2013 39 45 158.6 212.7 1.56 0.92 0.65 0.92 49.59 
GAP-20 20 24 135.5 144.9 0.77 0.90 0.92 0.90 86.37 
PBE-D3 23 25 101.0 115.2 0.91 0.96 0.96 0.96 100.00 

For each property tested in this work, we plot the accuracy of each CBOP and GAP-20 in Fig. 13 and compare their overall accuracy in Fig. 14. Clearly, GAP-20 performs better than all CBOPs studied here, especially in the simulation of cohesive energy, defects, thermal stability, and mechanical properties. Among them, the accurate prediction of mechanical properties is lacking for CBOPs, and GAP-20 makes up for this deficiency with high accuracy.

FIG. 13.

Accuracy of CBOPs and GAP-20 in cohesive energy, point defect, ad-atom defect, edge energy, thermal stability, and mechanical properties.

FIG. 13.

Accuracy of CBOPs and GAP-20 in cohesive energy, point defect, ad-atom defect, edge energy, thermal stability, and mechanical properties.

Close modal
FIG. 14.

The overall accuracy of the tested potentials indicates that the ML-IAP GAP-20 performs significantly better than the other potentials.

FIG. 14.

The overall accuracy of the tested potentials indicates that the ML-IAP GAP-20 performs significantly better than the other potentials.

Close modal

In summary, CBOPs (Tersoff, mo-Tersoff, AIREBO, mo-AIREBO, LCBOP, EDIP, and ReaxFFC2013) and GAP-20 ML-IAP are tested with respect to lattice constants, cohesive energy, defect properties, thermal stability, van der Waals interaction, and mechanical properties for crystalline phase carbon materials.

From our tests, we find that all potentials can predict the crystal structure, but only GAP-20 is accurate within 95% of the DFT-calculated cohesive energy. In reproducing the formation energy and optimizing the structures of common graphene point defects (SV, DV, and SW), only GAP-20 and LCBOP reach 90% accuracy compared to DFT. Notably, mo-AIREBO cannot accurately obtain the 5-7-7-5 structure and significantly overestimates the formation energy of SW defects. For the ad-atom defect, the accuracy of GAP-20, LCBOP, and ReaxFFC2013 is higher than 85%. Although mo-AIREBO most accurately deals with the ad-atom defect, it cannot reproduce the induced curvature graphene from the ad-dimer defect such as EDIP and mo-Tersoff. GAP-20 can accurately predict the edge energies of graphene and the (5,5) CNT and agrees that EedgeAC < EedgeZZ, but like the rest of the CBOPs, it underestimates the edge energy of the (10,0) CNT. In the thermal stability test, mo-Tersoff and mo-AIREBO underestimate the melting point of C60 significantly to 2000 K, while most of the potentials observe C60 withstand temperatures above 3000 K, with GAP-20 best reproducing the potential energy increase with temperature. We find that neither the CBOPs nor GAP-20 can describe van der Waals interactions correctly. CBOPs underestimate the interaction energy by an order of magnitude compared to DFT, and GAP-20 cannot reproduce the shape of the potential energy surface obtained with DFT, although it produces comparable energy barriers. Furthermore, we find that CBOPs cannot accurately reproduce the fracture behavior of graphene and CNTs. In addition, while reducing the C–C cutoff can remedy this, it then renders them unable to predict defects and thermal properties. In contrast, GAP-20 can correctly describe bond breakage, defects, and structures at high temperatures, making it suitable to accurately model the mechanical properties of carbon materials with defects and under high temperatures. Our results show that machine learning is a promising method for constructing universally applicable potentials for simulating and predicting material properties with much higher accuracy compared to traditional potentials that are restricted to specific applications due to their more limited parameterization.

See the supplementary material for the complete atomic structures of defective graphene, PES of bilayer graphene, and stress–strain curves of CNT and graphene when using different CBOPs.

The authors acknowledge the support from the Institute for Basic Science (Grant No. IBS-R019-D1), South Korea. The computational resources from CMCM, IBS, are also acknowledged.

The data that support the findings of this study are available from the corresponding author upon reasonable request.

1.
S.
Iijima
, “
Helical microtubules of graphitic carbon
,”
Nature
354
(
6348
),
56
58
(
1991
).
2.
H. W.
Kroto
,
J. R.
Heath
,
S. C.
O’Brien
,
R. F.
Curl
, and
R. E.
Smalley
, “
C60: Buckminsterfullerene
,”
Nature
318
(
6042
),
162
163
(
1985
).
3.
K. S.
Novoselov
, “
Electric field effect in atomically thin carbon films
,”
Science
306
(
5696
),
666
669
(
2004
).
4.
D. D. L.
Chung
, “
Review graphite
,”
J. Mater. Sci.
37
(
8
),
1475
1489
(
2002
).
5.
J.
Robertson
, “
Amorphous-carbon
,”
Adv. Phys.
35
(
4
),
317
374
(
1986
).
6.
B. V.
Spitsyn
,
L. L.
Bouilov
, and
B. V.
Derjaguin
, “
Vapor growth of diamond on diamond and other surfaces
,”
J. Cryst. Growth
52
,
219
226
(
1981
).
7.
R. H.
Baughman
,
H.
Eckhardt
, and
M.
Kertesz
, “
Structure‐property predictions for new planar forms of carbon: Layered phases containing sp2 and sp atoms
,”
J. Chem. Phys.
87
(
11
),
6687
6699
(
1987
).
8.
B.
Wu
,
X.
Jia
,
Y.
Wang
,
J.
Hu
,
E.
Gao
, and
Z.
Liu
, “
Superflexible C68-graphyne as a promising anode material for lithium-ion batteries
,”
J. Mater. Chem. A
7
(
29
),
17357
17365
(
2019
).
9.
T.
Wu
,
X.
Zhang
,
Q.
Yuan
,
J.
Xue
,
G.
Lu
,
Z.
Liu
,
H.
Wang
,
H.
Wang
,
F.
Ding
,
Q.
Yu
,
X.
Xie
, and
M.
Jiang
, “
Fast growth of inch-sized single-crystalline graphene from a controlled single nucleus on Cu–Ni alloys
,”
Nat. Mater.
15
(
1
),
43
47
(
2016
).
10.
L.
Qiu
,
F.
Ding
,
L.
Qiu
, and
F.
Ding
, “
Contact-induced phase separation of alloy catalyst to promote carbon nanotube growth
,”
Phys. Rev. Lett.
123
(
25
),
256101
(
2019
).
11.
J.
Dong
,
L.
Zhang
,
X.
Dai
, and
F.
Ding
, “
The epitaxy of 2D materials growth
,”
Nat. Commun.
11
(
1
),
5862
(
2020
).
12.
X.
Kong
,
J.
Zhuang
,
L.
Zhu
, and
F.
Ding
, “
The complementary graphene growth and etching revealed by large-scale kinetic Monte Carlo simulation
,”
npj Comput. Mater.
7
(
1
),
14
(
2021
).
13.
D.
Coello-Fiallos
,
T.
Tene
,
J. L.
Guayllas
,
D.
Haro
,
A.
Haro
, and
C.
Vacacela Gomez
, “
DFT comparison of structural and electronic properties of graphene and germanene: Monolayer and bilayer systems
,”
Mater. Today
4
(
7
),
6835
6841
(
2017
).
14.
P.
Lucignano
,
D.
Alfè
,
V.
Cataudella
,
D.
Ninno
, and
G.
Cantele
, “
Crucial role of atomic corrugation on the flat bands and energy gaps of twisted bilayer graphene at the magic angle θ ∼ 1.08°
,”
Phys. Rev. B
99
(
19
),
195419
(
2019
).
15.
F.
Hao
,
D.
Fang
, and
Z.
Xu
, “
Mechanical and thermal transport properties of graphene with defects
,”
Appl. Phys. Lett.
99
(
4
),
041901
(
2011
).
16.
A.
Bagri
,
S.-P.
Kim
,
R. S.
Ruoff
, and
V. B.
Shenoy
, “
Thermal transport across twin grain boundaries in polycrystalline graphene from nonequilibrium molecular dynamics simulations
,”
Nano Lett.
11
(
9
),
3917
3921
(
2011
).
17.
L.
Lindsay
,
W.
Li
,
J.
Carrete
,
N.
Mingo
,
D. A.
Broido
, and
T. L.
Reinecke
, “
Phonon thermal transport in strained and unstrained graphene from first principles
,”
Phys. Rev. B
89
(
15
),
155426
(
2014
).
18.
F.
Liu
,
P.
Ming
, and
J.
Li
, “
Ab initio calculation of ideal strength and phonon instability of graphene under tension
,”
Phys. Rev. B
76
(
6
),
064120
(
2007
).
19.
M. A. N.
Dewapriya
,
S. A.
Meguid
, and
R. K. N. D.
Rajapakse
, “
Atomistic modelling of crack-inclusion interaction in graphene
,”
Eng. Fract. Mech.
195
,
92
103
(
2018
).
20.
T.
Belytschko
,
S. P.
Xiao
,
G. C.
Schatz
, and
R. S.
Ruoff
, “
Atomistic simulations of nanotube fracture
,”
Phys. Rev. B
65
(
23
),
235430
(
2002
).
21.
M.
Topsakal
and
S.
Ciraci
, “
Elastic and plastic deformation of graphene, silicene, and boron nitride honeycomb nanoribbons under uniaxial tension: A first-principles density-functional theory study
,”
Phys. Rev. B
81
(
2
),
024107
(
2010
).
22.
K.
Cao
,
S.
Feng
,
Y.
Han
,
L.
Gao
,
T.
Hue Ly
,
Z.
Xu
, and
Y.
Lu
, “
Elastic straining of free-standing monolayer graphene
,”
Nat. Commun.
11
(
1
),
284
(
2020
).
23.
T.
Zhang
,
X.
Li
, and
H.
Gao
, “
Fracture of graphene: A review
,”
Int. J. Fract.
196
(
1-2
),
1
31
(
2015
).
24.
Y.
Wei
and
R.
Yang
, “
Nanomechanics of graphene
,”
Natl. Sci. Rev.
6
(
2
),
324
348
(
2019
).
25.
J.
Tersoff
, “
Empirical interatomic potential for carbon, with applications to amorphous carbon
,”
Phys. Rev. Lett.
61
(
25
),
2879
2882
(
1988
).
26.
D. W.
Brenner
, “
Empirical potential for hydrocarbons for use in simulating the chemical vapor deposition of diamond films
,”
Phys. Rev. B
42
(
15
),
9458
9471
(
1990
).
27.
N. A.
Marks
, “
Generalizing the environment-dependent interaction potential for carbon
,”
Phys. Rev. B
63
(
3
),
035401
(
2000
).
28.
S. J.
Stuart
,
A. B.
Tutein
, and
J. A.
Harrison
, “
A reactive potential for hydrocarbons with intermolecular interactions
,”
J. Chem. Phys.
112
(
14
),
6472
6486
(
2000
).
29.
A. C. T.
van Duin
,
S.
Dasgupta
,
F.
Lorant
, and
W. A.
Goddard
, “
ReaxFF: A reactive force field for hydrocarbons
,”
J. Phys. Chem. A
105
(
41
),
9396
9409
(
2001
).
30.
D. W.
Brenner
,
O. A.
Shenderova
,
J. A.
Harrison
,
S. J.
Stuart
,
B.
Ni
, and
S. B.
Sinnott
, “
A second-generation reactive empirical bond order (REBO) potential energy expression for hydrocarbons
,”
J. Condens. Matter Phys.
14
(
4
),
783
802
(
2002
).
31.
J. H.
Los
and
A.
Fasolino
, “
Intrinsic long-range bond-order potential for carbon: Performance in Monte Carlo simulations of graphitization
,”
Phys. Rev. B
68
(
2
),
024107
(
2003
).
32.
L.
Lindsay
and
D. A.
Broido
, “
Optimized Tersoff and Brenner empirical potential parameters for lattice dynamics and phonon thermal transport in carbon nanotubes and graphene
,”
Phys. Rev. B
81
(
20
),
205441
(
2010
).
33.
T. C.
O’Connor
,
J.
Andzelm
, and
M. O.
Robbins
, “
AIREBO-M: A reactive model for hydrocarbons at extreme pressures
,”
J. Chem. Phys.
142
(
2
),
024903
(
2015
).
34.
S. G.
Srinivasan
,
A. C. T.
van Duin
, and
P.
Ganesh
, “
Development of a ReaxFF potential for carbon condensed phases and its application to the thermal fragmentation of a large fullerene
,”
J. Phys. Chem. A
119
(
4
),
571
580
(
2015
).
35.
J. F.
Justo
,
M. Z.
Bazant
,
E.
Kaxiras
,
V. V.
Bulatov
, and
S.
Yip
, “
Interatomic potential for silicon defects and disordered phases
,”
Phys. Rev. B
58
(
5
),
2539
2550
(
1998
).
36.
C.
Jiang
,
D.
Morgan
, and
I.
Szlufarska
, “
Carbon tri-interstitial defect: A model for the DII center
,”
Phys. Rev. B
86
(
14
),
144118
(
2012
).
37.
O. A.
Shenderova
,
D. W.
Brenner
,
A.
Omeltchenko
,
X.
Su
, and
L. H.
Yang
, “
Atomistic modeling of the fracture of polycrystalline diamond
,”
Phys. Rev. B
61
(
6
),
3877
3888
(
2000
).
38.
L.
He
,
S.
Guo
,
J.
Lei
,
Z.
Sha
, and
Z.
Liu
, “
The effect of Stone–Thrower–Wales defects on mechanical properties of graphene sheets—A molecular dynamics study
,”
Carbon
75
,
124
132
(
2014
).
39.
Z.
Song
,
V. I.
Artyukhov
,
B. I.
Yakobson
, and
Z.
Xu
, “
Pseudo Hall–Petch strength reduction in polycrystalline graphene
,”
Nano Lett.
13
(
4
),
1829
1833
(
2013
).
40.
L. Q.
Xu
,
N.
Wei
, and
Y. P.
Zheng
, “
Mechanical properties of highly defective graphene: From brittle rupture to ductile fracture
,”
Nanotechnology
24
(
50
),
505703
(
2013
).
41.
A. P.
Bartók
,
M. C.
Payne
,
R.
Kondor
, and
G.
Csányi
, “
Gaussian approximation potentials: The accuracy of quantum mechanics, without the electrons
,”
Phys. Rev. Lett.
104
(
13
),
136403
(
2010
).
42.
J.
Behler
, “
Perspective: Machine learning potentials for atomistic simulations
,”
J. Chem. Phys.
145
(
17
),
170901
(
2016
).
43.
V. L.
Deringer
and
G.
Csányi
, “
Machine learning based interatomic potential for amorphous carbon
,”
Phys. Rev. B
95
(
9
),
094203
(
2017
).
44.
A.
Kamath
,
R. A.
Vargas-Hernández
,
R. V.
Krems
,
T.
Carrington
, and
S.
Manzhos
, “
Neural networks vs Gaussian process regression for representing potential energy surfaces: A comparative study of fit quality and vibrational spectrum accuracy
,”
J. Chem. Phys.
148
(
24
),
241702
(
2018
).
45.
K. T.
Schütt
,
H. E.
Sauceda
,
P.-J.
Kindermans
,
A.
Tkatchenko
, and
K.-R.
Müller
, “
SchNet—A deep learning architecture for molecules and materials
,”
J. Chem. Phys.
148
(
24
),
241722
(
2018
).
46.
P.
Rowe
,
G.
Csányi
,
D.
Alfè
, and
A.
Michaelides
, “
Development of a machine learning potential for graphene
,”
Phys. Rev. B
97
(
5
),
054303
(
2018
).
47.
L.
Zhang
,
J.
Han
,
H.
Wang
,
W. A.
Saidi
,
R.
Car
, and
W.
E
, “
End-to-end symmetry preserving inter-atomic potential energy model for finite and extended systems
,”
Adv.Neural Inf. Process.
31
,
4436
4446
(
2018
).
48.
O. T.
Unke
and
M.
Meuwly
, “
PhysNet: A neural network for predicting energies, forces, dipole moments, and partial charges
,”
J. Chem. Theory Comput.
15
(
6
),
3678
3693
(
2019
).
49.
M.
Wen
and
E. B.
Tadmor
, “
Hybrid neural network potential for multilayer graphene
,”
Phys. Rev. B
100
(
19
),
195419
(
2019
).
50.
P.
Rowe
,
V. L.
Deringer
,
P.
Gasparotto
,
G.
Csányi
, and
A.
Michaelides
, “
An accurate and transferable machine learning potential for carbon
,”
J. Chem. Phys.
153
(
3
),
034702
(
2020
).
51.
M.
Wen
and
E. B.
Tadmor
, “
Uncertainty quantification in molecular simulations with dropout neural network potentials
,”
npj Comput. Mater.
6
(
1
),
124
(
2020
).
52.
D.
Hedman
,
T.
Rothe
,
G.
Johansson
,
F.
Sandin
,
J. A.
Larsson
, and
Y.
Miyamoto
, “
Impact of training and validation data on the performance of neural network potentials: A case study on carbon using the CA-9 dataset
,”
Carbon Trends
3
,
100027
(
2021
).
53.
I. V.
Lebedeva
,
A. S.
Minkin
,
A. M.
Popov
, and
A. A.
Knizhnik
, “
Elastic constants of graphene: Comparison of empirical potentials and DFT calculations
,”
Physica E
108
,
326
338
(
2019
).
54.
C.
de Tomas
,
I.
Suarez-Martinez
, and
N. A.
Marks
, “
Graphitization of amorphous carbons: A comparative study of interatomic potentials
,”
Carbon
109
,
681
693
(
2016
).
55.
S.
Plimpton
, “
Fast parallel algorithms for short-range molecular dynamics
,”
J. Comput. Phys.
117
(
1
),
1
19
(
1995
).
56.
S.
Nosé
, “
Unified formulation of the constant temperature molecular dynamics methods
,”
J. Chem. Phys.
81
(
1
),
511
519
(
1984
).
57.
W. G.
Hoover
, “
Canonical dynamics: Equilibrium phase-space distributions
,”
Phys. Rev. A
31
(
3
),
1695
1697
(
1985
).
58.
G.
Kresse
and
J.
Hafner
, “
Ab initio molecular dynamics for open-shell transition metals
,”
Phys. Rev. B
48
(
17
),
13115
13118
(
1993
).
59.
G.
Kresse
and
J.
Furthmüller
, “
Efficiency of ab-initio total energy calculations for metals and semiconductors using a plane-wave basis set
,”
Comput. Mater. Sci.
6
(
1
),
15
50
(
1996
).
60.
J. P.
Perdew
,
K.
Burke
, and
M.
Ernzerhof
, “
Generalized gradient approximation made simple
,”
Phys. Rev. Lett.
77
(
18
),
3865
3868
(
1996
).
61.
S.
Grimme
,
J.
Antony
,
S.
Ehrlich
, and
H.
Krieg
, “
A consistent and accurate ab initio parametrization of density functional dispersion correction (DFT-D) for the 94 elements H-Pu
,”
J. Chem. Phys.
132
(
15
),
154104
(
2010
).
62.
S.
Nosé
, “
Constant temperature molecular dynamics methods
,”
Prog. Theor. Phys. Suppl.
103
,
1
46
(
1991
).
63.
C.
Chen
,
J.
Avila
,
H.
Arezki
,
V. L.
Nguyen
,
J.
Shen
,
M.
Mucha-Kruczyński
,
F.
Yao
,
M.
Boutchich
,
Y.
Chen
,
Y. H.
Lee
, and
M. C.
Asensio
, “
Large local lattice expansion in graphene adlayers grown on copper
,”
Nat. Mater.
17
(
5
),
450
455
(
2018
).
64.
Y.
Baskin
and
L.
Meyer
, “
Lattice constants of graphite at low temperatures
,”
Phys. Rev.
100
(
2
),
544
(
1955
).
65.
A. N.
Kolmogorov
and
V. H.
Crespi
, “
Registry-dependent interlayer potential for graphitic systems
,”
Phys. Rev. B
71
(
23
),
235415
(
2005
).
66.
H.
Shin
,
S.
Kang
,
J.
Koo
,
H.
Lee
,
J.
Kim
, and
Y.
Kwon
, “
Cohesion energetics of carbon allotropes: Quantum Monte Carlo study
,”
J. Chem. Phys.
140
(
11
),
114702
(
2014
).
67.
T.
Niu
,
M.
Zhou
,
J.
Zhang
,
Y.
Feng
, and
W.
Chen
, “
Growth intermediates for CVD graphene on Cu(111): Carbon clusters and defective graphene
,”
J. Am. Chem. Soc.
135
(
22
),
8409
8414
(
2013
).
68.
L.
Colombo
,
X.
Li
,
B.
Han
,
C.
Magnuson
,
W.
Cai
,
Y.
Zhu
, and
R. S.
Ruoff
, “
Growth kinetics and defects of CVD graphene on Cu
,”
ECS Trans.
28
(
5
),
109
114
(
2019
).
69.
Q.
Yuan
,
Z.
Xu
,
B. I.
Yakobson
, and
F.
Ding
, “
Efficient defect healing in catalytic carbon nanotube growth
,”
Phys. Rev. Lett.
108
(
24
),
245505
(
2012
).
70.
S.
Karoui
,
H.
Amara
,
C.
Bichara
, and
F.
Ducastelle
, “
Nickel-assisted healing of defective graphene
,”
ACS Nano
4
(
10
),
6114
6120
(
2010
).
71.
A. J.
Page
,
Y.
Ohta
,
Y.
Okamoto
,
S.
Irle
, and
K.
Morokuma
, “
Defect healing during single-walled carbon nanotube growth: A density-functional tight-binding molecular dynamics investigation
,”
J. Phys. Chem. C
113
(
47
),
20198
20207
(
2009
).
72.
M. M.
Ugeda
,
I.
Brihuega
,
F.
Guinea
, and
J. M.
Gómez-Rodríguez
, “
Missing atom as a source of carbon magnetism
,”
Phys. Rev. Lett.
104
(
9
),
096804
(
2010
).
73.
M. H.
Gass
,
U.
Bangert
,
A. L.
Bleloch
,
P.
Wang
,
R. R.
Nair
, and
A. K.
Geim
, “
Free-standing graphene at atomic resolution
,”
Nat. Nanotechnol.
3
(
11
),
676
681
(
2008
).
74.
H.
Shu
,
X.
Chen
,
X.
Tao
, and
F.
Ding
, “
Edge structural stability and kinetics of graphene chemical vapor deposition growth
,”
ACS Nano
6
(
4
),
3243
3250
(
2012
).
75.
Q.
Yuan
,
B. I.
Yakobson
, and
F.
Ding
, “
Edge-catalyst wetting and orientation control of graphene growth by chemical vapor deposition growth
,”
J. Phys. Chem. Lett.
5
(
18
),
3093
3099
(
2014
).
76.
C. O.
Girit
,
J. C.
Meyer
,
R.
Erni
,
M. D.
Rossell
,
C.
Kisielowski
,
L.
Yang
,
C.-H.
Park
,
M. F.
Crommie
,
M. L.
Cohen
,
S. G.
Louie
, and
A.
Zettl
, “
Graphene at the edge: Stability and dynamics
,”
Science
323
(
5922
),
1705
1708
(
2009
).
77.
J. D.
Baran
,
W.
Kołodziejczyk
,
P.
Larsson
,
R.
Ahuja
, and
J. A.
Larsson
, “
On the stability of single-walled carbon nanotubes and their binding strengths
,”
Theor. Chem. Acc.
131
(
9
),
1270
(
2012
).
78.
Y.
Li
,
R.
Ahuja
, and
J. A.
Larsson
, “
Communication: Origin of the difference between carbon nanotube armchair and zigzag ends
,”
J. Chem. Phys.
140
(
9
),
091102
(
2014
).
79.
D.
Hedman
and
J. A.
Larsson
, “
Length dependent stability of single-walled carbon nanotubes and how it affects their growth
,”
Carbon
116
,
443
447
(
2017
).
80.
D.
Hedman
and
J. A.
Larsson
, “
Analytical modelling of single-walled carbon nanotube energies: The impact of curvature, length and temperature
,”
SN Appl. Sci.
2
(
3
),
367
(
2020
).
81.
P.
Zhang
,
P. E.
Lammert
, and
V. H.
Crespi
, “
Plastic deformations of carbon nanotubes
,”
Phys. Rev. Lett.
81
(
24
),
5346
5349
(
1998
).
82.
O.
Lourie
,
D. M.
Cox
, and
H. D.
Wagner
, “
Buckling and collapse of embedded carbon nanotubes
,”
Phys. Rev. Lett.
81
(
8
),
1638
1641
(
1998
).
83.
C.
Lee
,
X.
Wei
,
J. W.
Kysar
, and
J.
Hone
, “
Measurement of the elastic properties and intrinsic strength of monolayer graphene
,”
Science
321
(
5887
),
385
388
(
2008
).
84.
M. M.
Shokrieh
and
R.
Rafiee
, “
Prediction of Young’s modulus of graphene sheets and carbon nanotubes using nanoscale continuum mechanics approach
,”
Mater. Des.
31
(
2
),
790
795
(
2010
).
85.
B. G.
Demczyk
,
Y. M.
Wang
,
J.
Cumings
,
M.
Hetman
,
W.
Han
,
A.
Zettl
, and
R. O.
Ritchie
, “
Direct mechanical measurement of the tensile strength and elastic modulus of multiwalled carbon nanotubes
,”
Mater. Sci. Eng.: A
334
(
1-2
),
173
178
(
2002
).
86.
J.-W.
Jiang
,
J.-S.
Wang
, and
B.
Li
, “
Young’s modulus of graphene: A molecular dynamics study
,”
Phys. Rev. B
80
(
11
),
113405
(
2009
).
87.
F.
Memarian
,
A.
Fereidoon
, and
M.
Darvish Ganji
, “
Graphene Young’s modulus: Molecular mechanics and DFT treatments
,”
Superlattices Microstruct.
85
,
348
356
(
2015
).

Supplementary Material