From diamond to BC8 to simple cubic and back: kinetic pathways to post-diamond carbon phases from metadynamics

TL;DR

This study uses advanced simulations to reveal the atomistic mechanisms and pathways for transforming diamond into post-diamond carbon phases like BC8 and simple cubic under extreme pressures, guiding future experiments.

Contribution

It introduces a novel metadynamics approach with machine-learned potentials to uncover transformation pathways and identifies conditions favorable for synthesizing specific post-diamond phases.

Findings

Nucleation of liquid droplets initiates phase transformation.

BC8 phase is favored near 1.8 TPa and 3500-5000 K.

A new metastable P222 structure was discovered.

Abstract

The experimental observation of elusive post-diamond carbon phases at extreme pressures remains a major challenge in high-pressure science. Using metadynamics with coordination-number-based collective variables and SNAP machine-learned interatomic potential, we uncover atomistic mechanisms governing the transformation of cubic and hexagonal diamond into post-diamond phases above 1.5 TPa. The transition initiates via homogeneous nucleation of nanoscale liquid droplets, which rapidly crystallize into either BC8 (below 1.8 TPa) or simple cubic phases (above 2.1 TPa), once the liquid nucleus surpasses a critical size. Favorable conditions for synthesizing BC8 are identified near 1.8 TPa and 3500--5000 K. Decompression pathways from simple cubic and BC8 phases were also simulated to study possible experimental recovery of post-diamond carbon allotropes at ambient conditions. We also find a new metastable low-enthalpy structure with four-coordinated carbon atoms and space group P222. Our insights provide a theoretical foundation for experimental discovery of ultra-dense carbon phases under extreme conditions.

Full text

The diamond phase of carbon is renowned not only for its exceptional physical properties, such as extreme hardness, outstanding thermal conductivity, and optical transparency (Field, 1992, 2012), but also for its remarkable thermodynamic stability across an exceptionally wide range of pressures and temperatures (Lazicki et al., 2021; Smith et al., 2014; Nguyen-Cong et al., 2024). Post-diamond phases of carbon, first predicted by Yin and Cohen (Yin and Cohen, 1983; Yin, 1984), include the 8-atom body-centered (BC8) structure, stable between 1 and 3 TPa, and the simple cubic (SC) phase, stable above 3 TPa at $T=0$ (Fig. 1) Martinez-Canales et al. (2012). However, experimental validation of these predictions remains extremely challenging Smith et al. (2014); Lazicki et al. (2021).

Both experiments (Lazicki et al., 2021; Smith et al., 2014) and simulations (Nguyen-Cong et al., 2024) reveal the extreme metastability of cubic diamond (CD) at pressures up to 2 TPa, well beyond its thermodynamic phase boundary at temperatures below 4000 K. The absence of a transition to the thermodynamically stable BC8 phase above 1 TPa is attributed to the very strong covalent bonding in diamond, which creates a large kinetic barrier between the crystalline phases. A double shock compression pathway to BC8 phase was proposed based on molecular dynamics (MD) simulations Shi et al. (2023). A review of studies of diamond at extreme conditions is given in Ref. Li et al. (2024).

Very recently, a novel diamond-to-BC8 transformation mechanism has been uncovered in large-scale MD simulations using a quantum-accurate Spectral Neighbor Analysis Potential (SNAP) (Willman et al., 2022) and an initial nanocrystalline diamond (NCD) sample Nguyen-Cong et al. (2024). It involves an appearance of a metastable liquid carbon intermediate, which allows to overcome a significant crystal-to-crystal transition barrier by melting diamond first, followed by nucleation and growth of BC8 from the supercooled carbon liquid. The presence of grain boundaries in the NCD sample significantly facilitated melting. This raises an important question: can such melting and subsequent BC8 nucleation be observed in single-crystal diamond samples, which lack extended defects to seed the transformation, and on much longer time scales beyond the reach of conventional MD?

A powerful technique to study the fundamental mechanisms of solid–solid phase transformations is metadynamics (Laio and Parrinello, 2002) – an enhanced sampling method that efficiently explores free energy landscapes and helps overcome barriers in MD simulations of rare events (for a review see Ref. Bussi and Laio (2020)). This approach is particularly useful for studying first-order structural transitions in crystals (Martoňák et al., 2003, 2005, 2006; Oganov, 2010; Hromadová et al., 2013; Plašienka and Martoňák, 2015; Badin and Martoňák, 2021), where activated processes are inaccessible to standard MD, as the transitions occur on timescales too long to be sampled. Metadynamics adds to the interaction potential of the system a history-dependent bias potential that acts on a selected set of collective variables (CVs) – low-dimensional descriptors of the system’s important degrees of freedom relevant to the structural transition. Early studies employed supercell vectors Martoňák et al. (2003, 2005, 2006) or average coordination number Zipoli et al. (2004); Hu et al. (2015) as CVs while recently more complex CVs were constructed Piaggi and Parrinello (2017, 2018); Piaggi et al. (2017), including neural-network-based path collective variables Rogal et al. (2019). In Ref.Badin and Martoňák (2021) it was shown that the average coordination number is a simple, generic, and efficient CV that not only captures displacive mechanisms of pressure-induced structural transformations but also enables access to nucleation processes in sufficiently large systems. Very recently, coordination number based CV was also used to study complex structural transformations in coesite(Vrba and Martoňák, 2025).

An early metadynamics study Sun et al. (2009) using ab initio MD and supercell vectors as CVs Martoňák et al. (2006) did not observe the diamond-to-BC8 transformations. At 2 TPa and 4000 K, CD transformed to the SC phase rather than to BC8. The BC8 phase was observed only upon subsequent decompression of the SC phase to 1 TPa at 5000 K. The limitation of that study was the very small system size (64 atoms), which significantly restricted the exploration of phase space.

In this work, we take the next step forward by employing the recently developed quantum-accurate SNAP model for carbon (Willman et al., 2022), which enables metadynamics simulations of much larger systems and longer timescales Nguyen-Cong et al. (2021); Willman et al. (2024). By focusing on phase transformations in compressed and heated single-crystal diamond, our study complements direct SNAP MD simulations of phase transformations in NCD (Nguyen-Cong et al., 2024) and extends the prior simulations of Sun et al. (Sun et al., 2009) to include homogeneous nucleation mechanisms of new carbon phases. The goal is to determine the pressure–temperature conditions under which cubic (CD, cF8) and closely related metastable hexagonal diamond (HD, hP4, lonsdaleite) transform into BC8. In addition, we explore potential transformations induced by decompression of the SC and BC8 phases and address the important question of whether BC8 can remain metastable upon decompression to lower pressures and potentially be recoverable at ambient conditions.

The metadynamics investigation of diamond-to-post-diamond transformations was performed using two coordination numbers as CVs – focusing at the first and second coordination shells in carbon phases (for more detailed information on the CVs, see description and Fig. S1 in Supplemental Material SM (see also Refs. Thompson et al. (2022); Tribello et al. (2014); Stukowski (2009); Erhard et al. (2024); Perdew et al. (1996); Kresse and Hafner (1993); Kresse and Furthmüller (1996, 1996); Kresse and Joubert (1999); Monkhorst and Pack (1976); Momma and Izumi (2011) therein)). The MD and metadynamics simulations were performed with the LAMMPS Thompson et al. (2022) and PLUMED Tribello et al. (2014) packages. The use of the SNAP potential (Willman et al., 2022) enables simulations of relatively large atomic systems up to 4,096 atoms across a wide range of pressures (1.5–2.1 TPa) and temperatures (2000–5000 K). As shown in Ref. (Badin and Martoňák, 2021), system size plays a crucial role: its increase results in a change in the fundamental mechanism of the phase transition, from martensitic, collective atomic rearrangements to nucleation and growth. We now present in Figs. 2-7 (prepared by the OVITO package(Stukowski, 2009)) typical kinetic pathways observed upon compression and decompression at specific pressure-temperature points, shown at the phase diagram in Fig. 1. The time evolution of CVs, fractions of the relevant phases, as well as the exploration of the space of CVs in all presented simulations are provided in Figs. S2-S7 of Supplemental Material SM . The identification of structural phases was performed using the Machine Learning Structure Identifier (MLSI) algorithm described in Ref. Erhard et al. (2024) (see also Supplemental Material SM ).

At 1.5 TPa and temperatures between 3500 and 5000 K, CD transforms into BC8 via nucleation of an intermediate disordered region. In a representative case at 4000 K (Fig. 2, Fig. S2), a localized disordered zone repeatedly forms and recrystallizes back to diamond until it exceeds a critical size of $\sim$ 350 atoms at $2.7\,\mathrm{ns}$, when BC8 rapidly nucleates and grows within the droplet. This transient disordered state acts as a precursor, lowering the barrier for BC8 emergence. Complete conversion to a BC8 nucleus occurs by $2.9\,\mathrm{ns}$, followed by formation of an interface with the parent diamond lattice. The nucleus then expands within the matrix, fully consuming the initial phase by $3.2\,\mathrm{ns}$ (Fig. 2). The estimated free-energy barrier of $\sim$ 500 eV (see Supplemental Material SM for information about the barrier calculation) is substantial but consistent with expectations for homogeneous nucleation in a defect-free and strongly bonded crystal. At a lower temperature of 3000 K, however, no phase transition is observed even over 7 ns of metadynamics simulation. Although transient disorder and spontaneous crystallization into simple motifs such as SC or hexagonal occur, the system ultimately reverts to diamond, indicating that nucleation of the liquid precursor is thermodynamically suppressed at lower temperatures.

At a higher pressure of 1.8 TPa and temperatures from 4000 to 5000 K, a similar transformation to BC8 is observed. However, at lower temperatures (down to 2500 K), a metastable SC structure becomes more prevalent. For example, at 3000 K, a disordered droplet is also formed, but it first crystallizes into an SC crystal, followed by melting and recrystallization into BC8 (Fig. 3 and Fig. S3). Growth then proceeds rapidly, fully consuming the parent phase. The transient appearance of the SC phase can be attributed to its lower structural complexity (one atom per primitive cell vs. eight in BC8) and increased energetic accessibility under elevated pressure, making it a kinetically favorable intermediate. At 2500 K we observe a similar sequence of events, while at 2000 K, the diamond phase remains stable and no transformation is observed.

Upon increasing the pressure to 2.1 TPa, while sampling temperatures between 4500 and 5000 K, the BC8 phase forms through a disordered droplet. However, at lower temperatures between 2500 and 4000 K the transformation instead favors the SC phase over BC8. This behavior at 4000 K is illustrated in Fig. 4 and Fig. S4. Initially, a disordered region emerges and crystallizes into an SC cluster containing up to 450 atoms, which subsequently reverts to diamond. In a second nucleation event, however, a SC nucleus forms within a liquid-like droplet and irreversibly grows into a slab of SC crystal, eventually consuming the entire simulation volume.

We also performed a similar series of simulations starting from metastable HD, since, to our knowledge, this transition pathway has not yet been explored. At 1.8 TPa and 4000 K, we again observed the initial nucleation of a liquid-like droplet. In this case, however, crystallization into the BC8 phase occurred rapidly, within $\sim$ 8 ps. The BC8 phase then grew quickly to consume the entire hexagonal lattice within 72 ps (Fig.5 and Fig. S5). In this case, the fast crystallization from the droplet preserved the orientational relationship between the parent and product crystal structures: the [111] direction of BC8 is parallel to the hexagonal axis of HD, but different orientations were also observed. The observed fast crystallization may indicate that BC8 could be more readily synthesized experimentally from HD than from CD. At a higher pressure of 2.1 TPa, as in the case of CD, HD transformed into the SC phase.

Our findings therefore suggest that the pressure range in which BC8 can form directly from diamond is relatively narrow. Under excessive overpressurization beyond 2 TPa, the system preferentially stabilizes in the SC phase, bypassing BC8. To further explore the diversity of structural transformations, we examined the potential reverse transformations of SC phase upon its decompression to 1.5 TPa at various temperatures. The evolution at 3500 K is shown in Fig. 6 and Fig. S6. Surprisingly, we found that the sample first transforms into a novel crystalline intermediate, which subsequently converts into BC8. Upon quenching to $T=0$, this intermediate is identified as a metastable crystal with a complex four-coordinated structure with 16 atoms in the unit cell and P222 symmetry (see Supplemental Material SM for details of the structure (Fig. S8), including enthalpy differences, and electronic density of states (Fig. S9)). These results indicate that even if diamond is driven to the SC phase by extreme compression, BC8 can still emerge upon subsequent decompression, consistent with prior observations Sun et al. (2009).

Finally, we investigated the key question of BC8 metastability upon quenching to ambient conditions. Early theoretical work Mailhiot and McMahan (1991) predicted that, at low pressure, BC8 may preferentially transform to HD rather than CD. To test this, we performed metadynamics simulations using a 1024-atom BC8 supercell at 3000 K and pressures of 500, 300, and 200 GPa. In all cases, the system had to overcome a substantial barrier before transforming into a mixture of CD and HD phases. The transformation pathway observed at 200 GPa and 3000 K is shown in Fig. 7 and Fig. S7. The pressure dependence of the barrier (Fig. S10) suggests that it decreases rapidly below 200 GPa. Accordingly, we carried out a standard NPT MD simulation at 100 GPa and 3000 K, where we observed a spontaneous transformation of BC8 into a CD/HD mixture that initiated immediately and completed within  750 ps. To probe the quenchability of BC8 at room temperature, we further conducted a 5-ns NPT MD simulation at 300 K and zero pressure. In this case, BC8 was found to be dynamically stable, suggesting the possibility of recovering this phase under ambient pressure provided the temperature remains sufficiently low.

To conclude, our metadynamics simulations revealed multiple kinetic pathways for the synthesis of the post-diamond BC8 phase in carbon at terapascal pressures (shown on the phase diagram in Fig. 1). We found that pressures around 1.8 TPa and temperatures between 3500 and 5000 K are optimal for direct BC8 formation, starting from both CD and HD diamond. In the latter case, the BC8 crystallization from the disordered intermediate droplet appears to proceed faster. In addition, alternative transformation routes between post-diamond phases at these extreme conditions were explored. BC8 is formed upon decompression of the SC phase, and might be quenchable to ambient pressure at room temperature. From a methodological standpoint, this study demonstrates that metadynamics employing simple generic collective variables such as coordination numbers, combined with quantum-accurate nachine learning potentials is a powerful approach for uncovering novel phase transformation mechanisms in complex materials. It enables the resolution of atomistic details of nontrivial nucleation processes under high-temperature conditions, where entropic contributions are substantial and cannot be neglected.

The work at USF is supported by Lawrence Livermore National Laboratory’s Academic Collaborative Team award, DOE/NNSA (Award Nos. DE-NA-0004089 and DE-NA0004234), DOE/FES (Award No. DE-SC0024640), NSF (Award No. 2421937), and NASA (Award No. 80NSSC25K7172). The computational work at USF is performed using leadership-class HPC systems: OLCF Frontier at Oak Ridge National Laboratory (ALCC and INCITE Awards Nos. MAT198 and MAT261), ALCF Aurora at Argonne National Laboratory (ALCC Award “FusAblator”), Perlmutter at National Energy Research Scientific Computing Center (NERSC award FES-ERCAPm3993 and ERCAP0034777) and TACC Frontera at the University of Texas at Austin (LRAC Award No. DMR21006).