In this study, we developed a machine learning-based neuroevolution potential (NEP) and employed MD simulations to investigate the annealing of amorphous carbon under varying temperature and density conditions. A total of nine representative amorphous carbon structures were identified, and their microscopic topological features were systematically analyzed through calculations of XRD patterns and structure factors. More importantly, we constructed a temperature-density topological phase diagram of amorphous carbon, uncovering the underlying principles governing structural formation under different thermodynamic conditions. A clear linear relationship was observed between critical temperature and density. Furthermore, we thoroughly examined the structural transitions of amorphous carbon under high-temperature and high-pressure conditions, providing valuable theoretical support for the synthesis, property tuning, and broader application of amorphous carbon materials.
Machine learning potentials (MLPs) is an emerging and powerful approach for accurately simulating material properties that are consistent with those observed in experiments. The core idea behind this method is to train machine learning (ML) models on large datasets of atomic configurations, which are typically labeled with energies, atomic forces, and virials calculated using first-principles methods such as density functional theory (DFT). Once trained, MLPs can perform MD simulations with computational efficiency comparable to empirical potentials while achieving the predictive accuracy of DFT. However, constructing high-quality and highly diverse training datasets remains one of the major challenges in MLP development. We constructed a dataset for carbon annealing simulations, comprising a total of 7,889 structures. The distribution of the dataset after dimensionality reduction using t-distributed stochastic neighbor embedding (t-SNE) is shown in Fig. 1a, with detailed construction procedures provided in Supplementary Table. 1. Each point represents an individual atomic structure, whose position is determined by structural descriptors of its local atomic environment. The Fig. 1a clearly shows that structurally similar data points exhibit obvious clustering, while dissimilar structures are well-separated, indicating that the dataset possesses good diversity and high-quality structural descriptors.
The NEP framework was employed to train the dataset and develop a MLP specifically designed for simulating the annealing process of amorphous carbon. The accuracy of the NEP framework has been well demonstrated in previous studies on complex systems such as short-range order in Pd-Cu-Ni-P alloys, lattice thermal conductivity, thermodynamic properties of water and ice, and structural stability of gold nanoclusters confined in zeolites. The training results are shown in Fig. 1b. The final root mean square errors (RMSE) for energy and forces on the training set reached 49.02 meV·atom and 565.5 meV·Å, respectively, while the corresponding values on the test set were 46.82 meV·atom and 561.32 meV·Å. These errors are comparable to, or even lower than, those reported for previous carbon-related MLPs, indicating that the NEP achieves a level of accuracy highly consistent with DFT calculations. The dataset also covers a wide range of energies and forces, ensuring the stability of NEP-based MD simulations under diverse conditions. Moreover, the achieved fitting accuracy is comparable to existing studies, validating the choice of hyperparameters used in training (see Supplementary Table. 2). The evolution of the loss functions over training steps and the virial fitting performance are shown in Supplementary Fig. 1, respectively.
To further assess the accuracy of the trained NEP, we performed a 15 ps relaxation simulation of an amorphous carbon system with a density of 3 g·cm at 4000 K, and calculated the radial distribution function (RDF), as shown in Fig. 1c. The NEP-simulated RDF closely matches the DFT results, confirming that NEP-driven MD simulations can reliably reproduce structural features consistent with DFT. In addition, we conducted structural optimizations and static energy calculations for diamond structures with a range of lattice constants, generating energy-volume (E-V) curves (Fig. 1d). The results show that NEP and DFT exhibit excellent agreement in describing the strain behavior of diamond structures. Furthermore, Fig. 1e presents a comparative analysis of the sp fraction of amorphous carbon across a wide density range (0.5-3.5 g·cm), as predicted by the NEP potential at 2000 K. These results are benchmarked against both experimental data and theoretical predictions from DFT and classical empirical potentials such as EDIP, Tersoff, and Brenner. The NEP predictions exhibit excellent agreement with experimental measurements and DFT calculations, accurately capturing the increase in sp content with density. In contrast, in addition to the EDIP potential, other traditional empirical potentials substantially underestimate the sp fraction, especially in the high-density regime. Although the EDIP potential exhibits limitations in describing long-range vdWs interactions and isolated sp bonds, its capability in representing sp bonding is comparable to that obtained from experiments and DFT calculations, which is consistent with recent literature report. This comparison highlights the accuracy of the NEP potential in simulating the structural evolution of amorphous carbon during annealing.
Using the trained NEP, we performed systematic large-scale annealing MD simulations of amorphous carbon by adjusting temperature and density. Figure 2 presents nine representative amorphous carbon structures obtained from different conditions. In Fig. 2a, a widely distributed DGN forms a three-dimensional, disordered, and porous carbon framework. Graphene sheets are connected via folding, stacking, or covalent bonding to form a non-periodic pore network, exhibiting a high degree of topological randomness. Previous studies have shown that the pore size of this structure decreases monotonically with increasing density. Due to the folding of graphene sheets, a large number of voids are formed experimentally, with porosity reaching up to 99.7% and a specific surface area as high as 850 m·g . We performed a statistical analysis of the ring structures in several typical DGN models (see supplementary Fig. 2). The results indicate that the carbon rings predominantly consist of 5-, 6-, 7-, and 8-membered rings, with 6-membered rings being the most abundant. Similar features have been reported by Wang et al., who showed that increasing density in nanoporous carbon leads to smaller, more numerous pores and the formation of elongated voids between stacked graphene fragments. Their ring statistics also revealed a six-membered ring-dominated topology, consistent with the highly graphitized sp network observed in our DGN structures. Martin et al. have pointed out that the topology of DGNs is mainly composed of saddle-shaped graphene nanoflakes with negative curvature, rich in out-of-plane topological defects. Our simulations further show that all DGN structures have a sp bonding fraction exceeding 96% (see supplementary Fig. 2), which is in excellent agreement with experimental results obtained via electron energy loss spectroscopy (EELS).
With increasing density, the graphene sheets formed under high-temperature conditions become compressed, resulting in high-density DGN (HD DGN, Fig. 2b). These structures are characterized by graphene sheets folded and stacked along specific directions, partially exhibiting parallel alignment, but still dominated by large, continuous folded regions. At slightly lower temperatures, some graphene sheets appear as fragmented domains, which we define as discontinuous DGN (Fig. 2c). In this structure, each nanoscale fragment can be viewed as a structural unit composed of several carbon rings, and the local microstructure is significantly distorted due to sp-sp hybridization. The formation of discontinuous DGNs is attributed to the inability of graphene sheets in multiple orientations to overcome energy barriers and rearrange into a unified low-energy structure at lower temperatures. In supplementary Fig. 3, we show the specific distribution of carbon ring types in this structure, which exhibits a significantly higher proportion of non-hexagonal rings compared to the 6-membered-ring-dominated configuration in conventional DGN. These short and partially aligned graphene sheets can be regarded as graphitic nuclei in the formation process of glassy carbon, making them an ideal structural model for this material. At even higher temperatures and densities, a large number of more or less parallel graphene sheets begin to form within the system. Although some overlapping still exists, the structure as a whole lacks long-range three-dimensional order. We define this structure as amorphous graphite (a-G, Fig. 2d). The formation mechanism of a-G originates from the system's ability to overcome energy barriers at elevated temperatures, enabling the development of locally larger graphitic domains.
When the system density increases to approximately 2.9 g·cm, the proportion of sp bonds rises significantly. At lower temperatures, an amorphous Diaphite (a-DG) structure is formed (Fig. 2e), while at higher temperatures, crystalline Diaphite is generated (Fig. 2f). In previous simulation studies, Zhu et al. discovered a-DG structures via MD simulations, revealing that they consist of disordered diamond-like regions mixed with graphene layers, where the local graphitic structure imparts medium-range order. They also identified a-DG as an intermediate state during the transformation from compressed multilayer graphene to diamond. Garvie et al. observed Diaphite structures in meteorites using HRTEM, showing they consist of few-layer graphene domains embedded within crystalline diamond, with topological correspondence at the interfaces. Nemeth et al. further proposed a classification framework for this new class of hybrid materials. Both a-DG and Diaphite are hybrid structures consisting of layered sp and bonded sp carbon domains, distinguished primarily by whether the sp-hybridized atoms form well-ordered crystalline diamond. Traditional empirical potentials such as EDIP have recently been successfully applied to describe mixed-hybrid microstructures and complex phase interfaces involving Diaphite, as demonstrated in recent studies. The NEP developed in this study also successfully reproduces this experimentally observed Diaphite structure. This result confirms the capability of the present NEP to accurately describe this complex carbon phase. Diaphite exhibits excellent fracture toughness and ductility due to the coexistence of sp and sp bonds.
As the density further increases to the range of 3.3-3.5 g·cm, the sp content rises sharply under extreme compression, and three types of diamond-like structures emerge in different temperature regimes: p-D (Fig. 2g), amorphous diamond (a-D, Fig. 2h), and nano-polycrystalline diamond (NPD, Fig. 2i). At sufficiently high temperatures and densities, diamond nuclei with various orientations continue to grow and eventually evolve into NPD. It is worth noting that p-D and NPD exhibit clear structural differences. Structurally, p-D consists of small-grained containing both cubic and hexagonal motifs, which are randomly embedded within an amorphous matrix. In contrast, NPD is composed of larger crystalline grains that exhibit well-defined long-range structural coherence. Recent experiments have successfully synthesized atomically a-D materials with densities exceeding 3 g·cm, and p-D has been obtained from C precursors within specific temperature windows under high-pressure conditions. Unlike a-D or NPD, p-D exhibits pronounced medium-range order, extending across several atomic shells. The p-D domains exhibit significant lattice distortion relative to ideal diamond structures, with the degree of deformation comparable to that observed in grain boundary atoms of ultrafine NPD.
To reveal the subtle differences among the nine amorphous carbon structures obtained, we calculated their XRD patterns and structure factors. Figure 3a shows the XRD results for the different structures. Although DGN, HD DGN, Discontinuous DGN, and a-G all exhibit graphitic characteristics, their diffraction curves display significant variations. DGN shows no distinct peaks, consistent with previously reported simulation results. HD DGN presents a clear peak around 1.7 Å, closely matching the experimental XRD pattern of glassy carbon. Discontinuous DGN and a-G exhibit characteristic peaks at 2.07 Å and 2.16 Å, respectively. It is noteworthy that the theoretical density of standard graphite is 2.26 g·cm, with a corresponding peak at 1.88 Å, whereas both Discontinuous DGN and a-G in our simulations reach densities of 2.7 g·cm, significantly higher than graphite. This indicates that the material undergoes volumetric compression, resulting in reduced interlayer spacing between graphene sheets (see supplementary Fig. 5). Moreover, the broadened peak at 2.07 Å in the Discontinuous DGN indicates a higher degree of disorder, similar to the a-DG structure, which exhibits extreme disorder due to the absence of diamond nuclei, showing only a weak diamond-related peak near 3.1 Å. Both Diaphite and NPD display characteristic diffraction peaks of crystalline diamond, including the (111), (220), and (311) planes. In NPD, the (220) and (311) peaks shift due to its higher density (3.5 g·cm); this shift is absent in NPD formed at a lower density of 3.3 g·cm, as shown in supplementary Fig. 6. In the XRD pattern of Diaphite, a distinct peak splitting appears around 3.0 Å, attributed to the presence of large hexagonal diamond domains accompanied by twinning and defect-induced nucleation behavior.
Furthermore, we performed Fourier transforms of the RDF to obtain structure factors, allowing for a deeper comparison of short-range order among the different structures (Fig. 3b). In the upper panel of Fig. 3b, DGN, HD DGN, and Discontinuous DGN exhibit clear differences in the low-Q region (Q < 4 Å), indicating fundamental distinctions in their degrees of short-range order. The sharper peaks observed for HD DGN and Discontinuous DGN suggest more locally ordered structures. The middle panel of Fig. 3b shows the structure factor of Diaphite, which features significantly enhanced peaks at positions I and I, reflecting the formation of diamond-like structures with a high proportion of sp bonding. In contrast, a-G and a-DG lack underlying crystalline frameworks, and their structure factor curves flatten beyond I without distinct peaks, exhibiting typical amorphous characteristics. NPD displays pronounced features at both I and I, indicative of diamond nucleation. In addition to calculating XRD curves and structural factors, recent studies have also been able to distinguish between p-D and NPD by calculating the free energy surface. Importantly, p-D and a-D have historically been difficult to distinguish due to their subtle structural differences. However, in this study, the high-precision MLP enables clear differentiation: p-D shows a slightly higher I peak than I, whereas a-D exhibits the opposite trend (lower panel of Fig. 3b). This observation is in excellent agreement with experimental results, further validating the accuracy of the potential in identifying microstructural features. These different amorphous carbon structures suggest that it may be possible in the future to regulate their microstructure for use in semiconductor information devices, silicon-based integrated circuits, gas sensors, and other applications.
To comprehensively illustrate the structural evolution of amorphous carbon under different temperature-density conditions, we constructed a phase diagram based on large-scale simulation results (Fig. 4a). Each symbol and its corresponding shaded area in the diagram represent a specific atomic topological structure. Overall, the phase diagram can be divided into three main regions: (I) graphite region, (II) transition region, and (III) diamond region. The simulation results indicate that the primary topological boundaries lie between ρ ≈ 2.8 g·cm and ρ ≈ 3.2 g·cm, which is in excellent agreement with experimentally reported. Within the range of ρ = 0.5-1.8 g·cm and T = 1500-4200 K, the phase diagram is dominated by the DGN structure; lower densities favor the folding of carbon atoms into graphene networks. These structures are characterized by abundant sp bonds and large carbon rings, resulting in an overall porous morphology. When ρ < 2.7 g·cm and T < 2800 K, the annealed structures are primarily Discontinuous DGN, corresponding to the trapezoidal region in the diagram. As temperature increases within the range of ρ = 2.0-2.9 g·cm, the structural evolution follows a sequence from Discontinuous DGN → HD DGN → a-G. The preferential formation of the a-G structure under our simulation conditions arises from the combined effects of thermodynamics, kinetics, and density constraints. Thermodynamically, sp-hybridized graphite is the most stable carbon phase at ambient pressure, with a lower formation enthalpy and higher entropy than diamond. Kinetically, high-temperature conditions significantly enhance the entropy-driven effect of the sp network, enabling the sp network dominated by hexagonal carbon rings to more easily overcome local energy barriers and expand into larger ordered structural domains. Additionally, the applied density range (2.0-2.9 g·cm) provides favorable interlayer spacing for vdWs interactions while suppressing diamond nucleation. These factors drive the system toward forming a-G with graphitic features.
Therefore, two graphitization pathways (blue dashed line) can be identified in the phase diagram. Graphitization path 1 represents the direct transition from DGN to a-G, highlighting the density-driven effect: when T ≥ 2800 K and ρ > 1.8 g·cm, the structure evolves into a-G without intermediate states. Graphitization path 2 occurs in the ρ ≈ 2.0-2.7 g·cm region, where, with increasing temperature, the system undergoes a sequential ordering process: Discontinuous DGN → HD DGN → a-G, revealing a temperature-induced stepwise transformation. When ρ lies between 2.8 and 3.1 g·cm (Region II), the system exhibits a mixed graphitic/diamond-like topology. At lower temperatures ( ~ 3000 K), a-DG forms, and upon heating, the amorphous diamond domains crystallize into diamonds, resulting in Diaphite. Further increasing the density to 2.9-3.5 g·cm (Region III), a-D is formed at low temperatures ( ~ 3000 K); continued heating leads to the emergence of a p-D intermediate state, which eventually transforms into NPD. There are also two diamondization pathways in the diamond region (black dotted line): diamondization path 1 involves a direct Diaphite → NPD transition driven by increasing density at high temperatures, a path commonly observed in high-pressure experiments. Diamondization path 2 occurs at constant density through temperature elevation, following the sequence a-D → p-D → NPD, highlighting the role of p-D as a transitional phase, this crystallization pathway is also commonly seen in previous studies. Figure 4b shows the temperature-density distribution of the sp/sp ratio: in the ρ < 2.8 g·cm region, deep blue indicates a very low sp fraction, while for ρ > 2.8 g·cm, the sharp increase in sp bonding results in a gradual shift toward deep red, consistent with the structural analysis. It is expected that more precise physical parameters will be proposed in future studies to distinguish between the nine types of amorphous carbon structures and to achieve a more accurate description of the phase diagram. Possible strategies include the development of machine-learned structural order parameters trained on local atomic environments, the use of localized electronic structure descriptors, and topological invariants derived from ring statistics or Voronoi-based tessellation. These methods may provide higher sensitivity to subtle local differences and thus enhance structural classification accuracy.
This phase diagram reveals in detail correspondence between temperature-density conditions and the topological structures of amorphous carbon, prompting deeper reflection on the underlying physical principles. First, within specific density ranges -- for example, at a density of 2.7 g·cm -- Discontinuous DGN forms at 2800 K, while a-G forms at 4000 K. This indicates that, within our simulation, temperature acts as a topological order parameter largely independent of density. This model-specific finding contrasts with conventional density-driven phase transition mechanisms and underscores the significant role of thermal effects in the local structural evolution of amorphous carbon. Second, we observed that at each given density, the critical temperature separating two phases exhibits an approximately linear increase with increasing density (indicated by the red dashed lines in the Fig. 4a). Moreover, the slopes of these critical temperature curves appear to be similar. Since the slope reflects the topological reconstruction energy barrier, this suggests that higher densities require higher temperatures to induce phase transitions. This interpretation is supported by our calculations of atomic root mean square displacement (RMSD) for samples annealed at 3000 K across the density range of 0.5-3.5 g·cm. As shown in Fig. 4c, at the same temperature, lower-density structures exhibit larger atomic RMSD due to more available free volume, whereas higher-density structures restrict atomic motion. Therefore, it can be inferred that higher temperatures are required at higher densities to drive atomic rearrangement in confined spaces. Additionally, the graphitization and diamondization paths appear to share structural similarities. In both cases, path 1 is driven by density, while path 2 is governed by temperature. Both path 2 transformations involve intermediate phases -- HD DGN and p-D, respectively. Finally, at 3.3 g·cm and 3200 K, p-D forms at a unique triple point where a-DG, Diaphite, and NPD coexist. This special condition marks a topological three-phase equilibrium. Previous metadynamics simulations have also suggested that p-D formation occurs within a temperature window of 3000-4000 K, aligning closely with experimentally observed values when interpreted within the Arrhenius framework; temperatures outside this range tend to yield either a-D or NPD. By decoding the phase transformation mechanisms of amorphous carbon in the temperature-density space, this study opens up further questions about how amorphous carbon structures dynamically evolve under realistic pressure-driven conditions.
As shown in Fig. 5, we achieved phase transitions between different amorphous carbon structures by tuning the applied pressure. The initial configuration was a DGN structure annealed at 3000 K with a density of 0.5 g·cm. Figure 5a presents the evolution of sp and sp atomic fractions, as well as crystalline atoms, under two target pressures: 40 GPa and 60 GPa. Since the evolution trends under both conditions are largely consistent, our analysis focuses on the results at 40 GPa. Based on the simulation settings, we can divide it into five stages.
In stage I, under an initial high temperature of 4000 K, no significant changes in physical parameters were observed, indicating that the DGN structure exhibits excellent thermal stability under high-temperature conditions -- consistent with experimental findings that carbon materials possess high-temperature resistance. Entering stage II, as the pressure gradually increases, the system undergoes complex structural evolution. Before 700 ps, the atomic fractions of sp and sp remain nearly unchanged (Fig. 5a, upper panel), with sp bonds dominating. Although the pressure at this point has reached approximately 20 GPa and the density about 2.5 g·cm (Fig. 5b), the structure primarily transforms into HD DGN (Fig. 5c①). As the pressure continues to rise, the sp fraction decreases sharply while the sp fraction increases rapidly. By 882 ps, the system evolves into a mixed structure of graphitic fragments and amorphous tetrahedral carbon -- a-DG (Fig. 5c②). The graphitic domains gradually collapse under further volumetric compression, and by 1000 ps, the system transitions into an a-D state, predominantly composed of tetrahedral structures (Fig. 5c③), with the sp content reaching approximately 63%. In stage III, crystalline atoms begin to appear, marking the onset of diamond nucleation. A p-D emerges at 1154 ps (Fig. 5c④), though it exists only briefly before rapidly crystallizing under thermal activation, eventually forming a stable NPD (Fig. 5c⑤). The system remains in this stable state for an extended period until entering stage IV. In stage IV, as the pressure is reduced while maintaining high temperature, the previously stable NPD nuclei begin to dissociate. The structure transforms into layered graphitic configurations, with a transient Diaphite phase appearing at 3004 ps (Fig. 5c⑥). This intermediate phase also exists for a short time before evolving into a Discontinuous DGN structure at 3013 ps (Fig. 5c⑦). Detailed descriptions of this transition process are provided in supplementary Fig. 7. Finally, in stage V, the system stabilizes into a DGN structure with a density of approximately 1.5 g·cm at 3750 ps (Fig. 5c⑧), which is slightly higher than the initial density of 0.5 g·cm. The densities of all typical structures found in Fig. 5c are consistent with those shown in the phase diagram (Fig. 4a).
Throughout the entire phase transition process, all of the previously discussed eight structures were observed, except for a-G. This absence is mainly attributed to the initial structure being DGN: although some partially folded graphene sheets may align into parallel layered structures under pressure, their lifetime is too short to evolve into large-scale a-G domains. By tuning the applied pressure, we achieved a controllable phase transition from DGN to diamond-dominated structures (p-D and NPD), highlighting the critical role of pressure in the topological reconstruction of amorphous carbon. A comparison between the 40 GPa and 60 GPa results reveals that during stage II, sp atoms reach their peak concentration more quickly under 60 GPa (Fig. 5a, lower panel). However, the crystalline atom content in later stages is actually lower than that at 40 GPa -- approximately 52% versus 63%. This suggests that high pressure alone is insufficient to overcome the nucleation barrier of amorphous carbon; additional thermal activation is required. Notably, although the applied pressure is symmetric, the processes of diamond nucleation and dissociation exhibit significant asymmetry. The system first undergo a fully a-D before gradually crystallizing into diamond structures under the influence of temperature -- a process corresponding to diamondization path 2 in Fig. 4a. Along this path, the system's density remains nearly unchanged, which aligns well with the "temperature-dominated" nucleation mechanism discussed earlier in Fig. 4a. Previous studies have indicated that the transformation from graphite to diamond requires the assistance of nanoscale coherent interfaces to consume remaining graphite and complete the transition(i.e., Diamondization Path 1). This pathway was not observed in our simulations, likely due to the absence of a-G in the system. Under high temperatures, the short-timescale compression renders the DGNs unstable, making it prone to directly transform into an a-D structure. As a result, the phase transition pathway a-D → p-D → NPD emerges as the optimal crystallization route for the system. Upon pressure release and volume expansion, diamond tends to dissociate into graphitic layered structures under continued high-temperature conditions. This behavior is attributed to the lower free energy of graphite compared to diamond.