Development of a machine learning interatomic potential for exploring pressure-dependent kinetics of phase transitions in germanium Open Access

In recent years, there has been a remarkable transformation in how we approach atomistic simulations, and at the forefront of this change is the increase in the popularity of machine learning interatomic potentials (MLIPs). This innovative approach offers a way to significantly reduce the computational cost of ab initio calculations, typically based on density functional theory (DFT), with a great accuracy/efficiency trade-off allowing us to study larger systems and/or longer time-scale processes. The remarkable efficiency and precision of MLIPs have made them essential tools, gaining widespread recognition and adoption for exploring various atomistic systems.^1–11 

Our objective in this work was the development of an interatomic potential suitable for investigating crystal phase transitions, taking germanium as a particularly interesting system due to its potential applications.¹² The emphasis lies in showcasing and testing the developed potential, rather than conducting a comprehensive study of phase transitions in Ge. It is important to note that this endeavor necessitates the creation of an MLIP capable of accurately representing pressure-dependent kinetic barriers. Merely training the MLIP relying solely on configurations close to equilibrium and snapshots from molecular dynamics near metastable crystal structures would not provide adequate insight into transition mechanisms. This limitation arises because crucial saddle points defining kinetic barriers would be under-represented or excluded from the training dataset, leading to unreliable predictions. To address this, we incorporated configurations into our dataset along minimum energy paths obtained through solid-state nudged elastic band (ssNEB) calculations under varying stress conditions, to provide insight into transition barriers.

Empirical potentials, including Stillinger-Weber,^13,14 Tersoff,¹⁵ and modified embedded-atom method (MEAM),^16,17 have been developed to significantly reduce computational costs compared to the DFT calculations, but may be unreliable when dealing with metastable phases¹⁸ and, more generally, with kinetics.

The interesting electronic properties of germanium, such as high intrinsic electron mobility, make it useful for novel high-speed electronic and photonic devices. In addition, its narrow bandgap enables efficient infrared detection in optoelectronic devices, rendering Ge a valuable material for advancing modern technology. Ge is a well-known semiconductor that typically exhibits a cubic diamond (CD) structure (space group Fd$3̄$m1) at ambient pressure and temperature. Through anvil press and nanoindentation experiments, CD transforms into a metallic phase with the β-Sn structure (I4₁/amd) upon applying non-hydrostatic pressure at around 10 GPa,^19–21 showing remarkable stability over a wide range of pressures. Upon increasing pressure, β-Sn undergoes a series of transitions into other high-density metallic phases, including simple hexagonal (SH) and HCP.^22,23 Upon unloading from β-Sn, CD is not recovered and different metastable allotropes are reached instead. Experiments showed the presence of metastable BC8 (Ia$3̄$⁠), ST12 (P4₃ $2̄12̄$⁠), and R8 (R$3̄$⁠) phases and led to believe that the formation of such allotropes depends on stress conditions and rate of decompression from β-Sn.^24–28 Sustained interest in ST12 stems from its optical properties,²⁹ with recent advancements leading to the synthesis of a large pure bulk sample,³⁰ the production of ST12 nanoparticles,^31,32 and also the fabrication of nanowires.³³ BC8 and R8 are observed to transform to a hexagonal diamond (HD) phase under ambient pressure and temperature.^34,35 This is of great technological interest given HD’s predicted narrow direct bandgap, rendering it well-suited for optoelectronic applications compatible with silicon.¹² It is also noteworthy that amorphous Ge (a-Ge) has been extensively documented in various experiments: observed as a byproduct during the unloading from β-Sn, as a pathway to achieve β-Sn through loading and also through observed transitions from a low-density to a high-density form of a-Ge.^36–38 

The Vienna Ab Initio Simulation Package (VASP)^39,40 was used to perform all the DFT calculations. The Ge 3d 4s 4p electrons were treated as valence electrons by the projector augmented-wave (PAW) method.⁴¹ Perdew–Burke–Ernzerhof (PBE)⁴² was used as the exchange–correlation functional as it is a popular choice given its reasonable accuracy over a wide range of systems. Since metallic phases of Ge are demanding in terms of convergence, we conducted convergence tests to find the optimal balance between accuracy and efficiency. The truncation energy of the plane wave basis was set to 500 eV, and the k-point grid spacing was fixed at 0.14 Å⁻¹. Formation energies were computed by taking the energy per atom difference between a given structurally optimized crystal phase and the most stable one, the CD phase. Specifically, using the selected k-point grid spacing, we observed that the formation energies for the 4H-HD, 2H-HD, BC8, ST12, and R8 phases have already reached convergence, showing variations of less than 1 meV/atom upon further parameter refinement. However, for denser metal structures, such as β-Sn, Fmmm, SH, HCP, BCC, and FCC, there are still oscillations, with some reaching up to ∼5 meV/atom in their computed energies. We still consider this result as satisfactory since different DFT functionals would predict formation energies with even larger discrepancies. Since we are interested in pressure-dependent phase transitions, we also monitored the convergence in terms of pressure. Using the selected k-point grid spacing of 0.14 Å⁻¹, we observed that semiconducting phases have already reached a convergence of about 0.01 GPa. Some denser metallic structures, more demanding in terms of convergence, still present oscillations, with some limited examples reaching $≈1GPa$⁠.

Molecular dynamics (MD) simulations were carried out in LAMMPS⁴³ using the DeePMD extension.⁴⁴ To sample data, both the canonical (NVT) and the isothermal–isobaric (NPT) ensemble were used with an integrator time step of 1 fs, while the temperature and pressure were set using a Nose–Hoover thermostat and barostat⁴⁵ with a relaxation time of 0.5 and 5.0 ps, respectively.

To perform local structural optimizations, the Fast Inertial Relaxation Engine (FIRE) optimizer⁴⁶ within the Atomic Simulation Environment (ASE)⁴⁷ was employed by setting up the trained MLIP as a calculator to obtain potential energies, forces, and stresses.

When investigating transitions, a critical parameter of focus is the activation energy barrier associated with the reaction mechanism. This is significant because following transition state theory,⁴⁸ the transition rate is directly proportional to the exponential of the free energy barrier. Many different methods have been developed to accurately locate such transition states,^49–55 but only a few take into account variable cell transformations. In this work, a generalized solid-state nudged elastic band (ssNEB) method⁵⁶ with the climbing-image^57,58 was employed for determining reaction pathways of crystal phase transitions involving both atomic and cell degrees of freedom. We utilized this well-established technique by means of its implementation for the ASE package in the TSASE code.⁵⁹ 

In order to build the MLIP for Ge, a reference dataset for the training process is generated using an iterative training procedure, exploiting the concept of active learning.^60–63 

To thoroughly explore the kinetics of solid–solid phase transitions, we integrated configurations explored by performing ssNEB calculations into our dataset. This inclusion is pivotal, as solely relying on sampling configurations near metastable crystal structures would not be sufficient for studying transition barriers.

The outlined procedure, schematically shown in Fig. 1, can be summarized as follows: we selected many different Ge-crystalline phases from the Materials Project (MP) database.⁶⁴ These include CD (mp-32), 4H-HD (mp-1091415), 2H-HD (mp-1007760), BC8 (mp-1080106), ST12 (mp-137), R8 (mp-128), β-Sn (mp-78), SH (mp-1224349), Imma (mp-1061054), Fmmm (mp-148), HCP (mp-1008733), BCC (mp-998883), and FCC (mp-12093). After replicating these cells, the lattice parameters and the atomic positions are perturbed to obtain many different configurations. A starting incomplete dataset is built based on these structures that are labeled with energies, forces, and stresses computed through DFT calculations. An initial, coarse model trained on this dataset is used to explore more configurations. This is done by sampling snapshots from NVT and NPT MD, and intermediate images from ssNEBs. The energies, forces, and stresses of these structures are computed again through single-point DFT calculations and added to the existing starting dataset.

Multiple neural network (NN) models are then trained using identical training sets but with different random seeds. Consequently, these models would lead to diverse predictions for the same configuration. If predictions are consistent among the different models, it would suggest that the configuration aligns closely with the training set; conversely, discrepancies indicate unreliability in predictions. By employing this approach, we can select configurations exhibiting unreliable predicted energies, forces, and stresses; compute these observables again through single-point DFT calculations; and then, iteratively expand the training set by incorporating such data deemed most relevant.

The sequential addition of inaccurately estimated configurations into the training dataset during iterations is a crucial aspect of constructing the MLIP.

The final dataset, upon which the actual model is trained, is composed of roughly 2700 structures making up a total of roughly 112 thousand atomic environments, and it is split into 87% training and 13% validation sets. The split was done randomly, except for certain structures that we specifically included in the training set: phases at their minimum energy configuration and with lattice parameter perturbations close to the values corresponding to this minimum. Since many crystal-phase transitions are pressure-induced, structures covering a wide range of pressure have been sampled, reaching almost 100 GPa for some configurations. Most of the sampling was conducted within a pressure range of 5 GPa tensile to 30 GPa compressive strain. This range was chosen because the most technologically relevant crystal phase transformations in Ge, such as the formation of hexagonal Ge and ST12, occur after loading of CD to β-Sn and subsequent unloading; these experiments are typically performed below 20 GPa. Structures containing interstitial atoms and vacancies were obtained by sampling MD starting from manually crafted configurations. Disordered structures were obtained by heating and then cooling some given crystal phases to get highly defective and amorphous-like configurations. These highly varied structures allow for a more comprehensive sampling of possible atomic environments, leading to a more robust interatomic potential. Transition structures were obtained by sampling intermediate images from ssNEB calculations between different crystal phases over a wide range of pressures. In addition, after replicating these intermediate images, the lattice parameters and the atomic positions are perturbed to sample more different configurations close to the transition pathway found. This simple yet effective procedure is an alternative to other methods, such as transition path sampling,⁶⁵ metadynamics,⁶⁶ umbrella sampling,⁶⁷ or many more recent ones,^68–72 employed for an extensive sampling of the configuration space close to transition states and along a reaction coordinate.

The dataset resulting from our iterative procedure is presented in Table I.

Utilizing our dataset, the NN potential was trained using the Deep Potential Molecular Dynamic package (DeePMD-kit).^73,74 The deep potential (DP) model decomposes the total energy of the system into atomic contributions by employing a local geometric description of the atomic environment. In this approach, a multilayer NN is used as a regression model to describe the relationship between atomic configurations and energy. For each local system, an embedding network converts the atomic coordinate information into a descriptor matrix, which is fed into a fitting network mapping the matrix to the local atomic energy. An in-depth explanation of the NN architecture is given in the original paper.⁷³ 

To accurately capture the structural information of many Ge crystal phases and different configurations, we used two types of descriptors constructed from all the information (both angular and radial) in the so-called DeepPot-SE framework:⁷⁴ a, “se_e3” type embedding that takes angles between two neighboring atoms as input and an “se_e2_a” type embedding that takes the distance between atoms as input. For the former descriptor, we used a cutoff value of r_c = 3.30 Å; a smoothing cutoff r_s = 2.00 Å; and a three-layer embedding network containing 8, 16, and 32 neurons. For the latter, r_c = 6.60 Å and r_s = 3.30 Å were used along with a three-layer embedding network containing 16, 32, and 64 neurons with 32 axis neurons. For the fitting net, three hidden layers with 160, 120, and 80 neurons in each layer were employed. We trained and tested multiple models with various architectures and hyperparameters to empirically identify optimal values that balance accuracy and efficiency.

Figure 2 shows the regression plots for energies, forces, and stresses concerning both the training and validation datasets. Our model demonstrates a satisfactory performance, with root mean square errors (RMSEs) of ∼5 meV for energies and 0.1 eV/Å for forces. Stress prediction is also robust, exhibiting an RMSE of less than 0.5 GPa over a broad spectrum spanning ∼100 GPa.

To assess the reliability of the model, a diverse test set consisting of snapshots from MD at different temperatures and pressures was also prepared. In addition, the main results concerning transition states are thoroughly discussed in the dedicated Sec. III B. In detail, the test set comprises 125 diamond structures (CD, 4H-HD, and 2H-HD), corresponding to 7888 atomic environments, 74 low-density semiconductor structures (BC8, ST12, and R8) (3984 atomic environments), 97 high-density metal structures (6144 atomic environments), and 146 disordered or defective structures, encompassing amorphous configurations as well as those with interstitial atoms or vacancies (13 401 atomic environments). Even on this test set, our model exhibits strong predictive capability, showing minimal deviations in the regression plots shown in Fig. 2.

Further evaluation focused on the prediction of formation energies. To determine the DFT formation energies, we computed the energy difference per atom between a specific structurally optimized crystalline configuration and the cubic diamond reference. Similarly, for estimating the DP formation energies, we followed the same procedure, but utilizing structures optimized for the DP. The comparison of the results obtained from such calculations is presented in Fig. 3 and Table II.

The model performs remarkably well, with errors smaller than 5 meV for all the tested structures, except for the BCC phase. As expected, the predictive ability of the model on the 4H-HD, 2H-HD, BC8, ST12, and R8 phases is slightly better than the other higher energy ones. Particularly noteworthy is the accurate prediction of the 6H-HD phase, despite its absence in the training set. As an intermediate phase positioned between CD and 4H-HD and characterized by distinct stacking faults, the inclusion of CD, 4H-HD, and 2H-HD in the training set seems to be enough to adequately encapsulate the local atomic environments of the 6H-HD structure.

In addition, we analyzed the accuracy of energy vs volume curves, as shown in Fig. 4, which provide insights into the volumetric compression–expansion behavior of different phases. While the model performs admirably for diamond phases and the other semiconductor structures BC8, ST12, and R8, it exhibits slightly lower accuracy for metallic phases, especially in the volumetric expansion region. This is because apart from the aforementioned issue concerning convergence, the training set emphasized the sampling of these phases under pressure conditions rather than volumetric dilation, and it is generally understood that NN models are interpolative, providing reliable results only within their training domain. This decision was made because these phases are only competitive under high pressures, rendering volumetric dilation less pertinent from a physical standpoint.

Three additional phases, which were not part of the training set, are also tested: Pbam, P4₂/ncm, and P4₁2₁2. These structures, recently predicted theoretically, display significantly lower energies compared to the majority of metastable polymorphs observed in Ge. For this reason, there is speculation that these phases might emerge as metastable forms in experiments utilizing indentation or anvil cells.^76,77 The comparison presented in Table II, between formation energies computed via DFT and those predicted by the DP, reveals a larger error for such crystal phases absent from the training set. To address this, any investigation aimed at these structures would necessitate their inclusion in the training set as well.

Crystal phase transitions in germanium were extensively investigated utilizing ssNEB calculations and NPT MD. Our research revealed several pathways between different crystal phases, helping understand the phase transition mechanism and kinetics. As our objective in this paper is to showcase and assess the developed potential, rather than conduct an exhaustive study of phase transitions in Ge, in the following, we delve into a few illustrative examples.

First, a collective transformation involving a 16-atom cell was observed, connecting the CD and BC8 phases. Multiple cell geometries have been selected as the starting and ending points for the ssNEB calculation. We chose the geometries that required the least amount of atomic displacement and bond breaking. The minimum energy path for this transition, shown in Fig. 5(a), reveals a relatively high predicted activation energy of ∼0.28 eV/atom. It is indeed commonly acknowledged that the reaction barrier during solid phase transition can often be significant, as numerous chemical bonds are simultaneously formed or broken as the crystal transitions from one phase to another.

Similarly, another pathway connecting the CD and β-Sn phases was identified, showing a remarkable agreement between ab initio and NN prediction, as shown in Fig. 5(b). This transformation is quite easy to identify because the relative positions of atoms in the cell remain constant during the transformation, only the lattice parameters change. This pathway, previously reported for Silicon,⁷⁸ occurs similarly in Ge. The activation barrier, estimated to be ∼0.32 eV/atom, presents a significant obstacle. However, the transformation is not hindered because stress can be applied to modify the MEP and the transition energy barrier, which directly affects the rate of phase transition. When subjected to external stress, the energy landscape is transformed into an enthalpy landscape, accounting for the work done by external stresses. To explore this phenomenon, we conducted multiple calculations using the ssNEB method with various levels of uniaxial stress. Through this procedure, we determined that at ∼7 GPa of compressive uniaxial stress, the transition is expected to be barrierless. To our knowledge, such an agreement between ab initio and NN on the calculation of pressure-dependent energy barriers is unprecedented.

NPT MD simulations were employed to investigate such pressure-induced transformation. To examine the structural transformation, we utilized a simulation cell measuring ∼56 × 56 × 76 Å, containing 10 368 atoms. The temperature was set at 300 K (room temperature), with uniaxial compressive stress along the z direction at 7 GPa. Snapshots from such simulated transition are shown in Fig. 6(a).

To ensure confidence in the model’s accuracy, we calculated the deviation in the predicted force for each atom at each step using an ensemble of three neural networks. We then computed the mean deviation at each step and recorded the maximum deviation value among all atoms. As shown in Fig. 6(b), force predictions are considered satisfactory, with mean and maximum deviations below 0.05 and 0.15 eV/Å, respectively, even at the transition midpoint where deviations are higher.

We then identified a pathway connecting the ST12 and β-Sn phases, shown by its corresponding MEP in Fig. 5(c). This transformation pathway, involving a 12-atom cell, aligns with a previously proposed atomistic mechanism for such transition.⁸⁰ The activation energy barrier here is ∼0.18 eV/atom, but it can be reduced by applying stress. In particular, we observed that the application of compressive planar stress, in this case, effectively reduces the barrier, leading to a transition becoming barrierless at ∼12 GPa of applied planar stress.

In addition, the impact of hydrostatic pressure on the R8 to β-Sn transition was investigated. The results indicate a gradual reduction in the activation barrier with increasing pressure, as shown in Fig. 5(d). Interestingly, a transition pathway via the BCC phase emerges above 20 GPa. In this case, more intermediate images have been added to the pathway, which has been broken up into separate NEB calculations between all local minima.

It is noteworthy that, as shown in Fig. 5, a comparison of the ssNEB calculations performed using the MLIP with those conducted using DFT reveals that the prediction of the energy of intermediate images in the transition from R8 to β-Sn is quite poor, with errors reaching nearly 30 meV. This discrepancy arises because this transition process was examined as a test after the development of the MLIP, thus lacking intermediate structures in the training set.

NPT MD simulations were once again utilized to explore this pressure-induced transformation. A rhomboidal prism-shaped simulation cell housing 12 288 atoms was employed, with the temperature maintained at 200 K and hydrostatic pressure at 23 GPa. The pressure valued was chosen by conducting ss-NEB calculation on a simulation cell like the one used for the results shown in Fig. 5(d). We estimated the kinetic barrier at this pressure value to be about 8 meV/atom. Such a low barrier would allow us to see the transition directly during MD in a short simulation.

Snapshots capturing during MD simulation are shown in Fig. 7(a). The nucleation of the BCC structure can be observed emerging from the R8 bulk and subsequently transform into polycrystalline β-Sn with defects. As shown in Fig. 7(b), the mean deviations are ∼0.05 eV/Å and decrease below 0.02 eV/Å as the transition progresses and more β-Sn forms.

The transition from β-Sn to R8, in reverse order with respect to the one investigated here, has been experimentally observed to occur in nanoindentation upon unloading at specific rates and stress conditions. However, the R8 to β-Sn transition under loading investigated in this work has been less studied, primarily because obtaining R8 from β-Sn is a prerequisite. Therefore, for such a transition, there are no existing values to compare our estimation to. Still, our findings could be valuable for better understanding cyclic nanoindentation experiments, such as those described in Ref. 28.

In this work, an accurate and efficient MLIP tailored for investigating pressure-dependent crystal phase transitions in germanium has been developed. The resulting DP is several orders of magnitude faster than the DFT calculations and scales linearly with the number of atoms. Regression plots comparing NN predicted energies, forces, and stresses with the DFT calculations performed on a suitable test set, together with the shown accurate prediction of formation energies and energy–volume curves of many relevant crystal phases, demonstrate the robustness of the model in describing such metastable allotropes of Ge. Moreover, the iterative refinement of the model through insertion in the training set of transition structures sampled from ssNEB calculations performed with the model itself has been shown as a simple yet effective practice to get the model to predict accurate activation barriers.

Our MLIP, which can be downloaded together with the full database,⁸¹ will serve as a valuable tool facilitating further advancements in the study of Ge crystal phase transitions and their applications.

A.F. expresses gratitude to Naman Katyal and Jiyoung Lee for their insightful feedback. A.F., F.M., and E.S. acknowledge financial support from ICSC–Centro Nazionale di Ricerca in High-Performance Computing, Big Data and Quantum Computing, funded by the European Union–NextGenerationEU. A.F., F.R., F.M., and E.S. acknowledge the CINECA consortium under the ISCRA initiative for the availability of high-performance computing resources and support. The work in Austin was supported by the National Science Foundation (Grant No. CHE-2102317) and the Texas Advanced Computing Center. E.S. acknowledges financial support from the National Recovery and Resilience Plan (NRRP), Mission 4, Component 2, Investment 1.1, Call for tender No. 104 published on February 2, 2022, by the Italian Ministry of University and Research (MUR), funded by the European Union–NextGenerationEU–Project Title “SiGe Hexagonal Diamond Phase by nanoIndenTation (HD- PIT)”–CUP H53D23000780001-Grant Assignment Decree No. 957 adopted on June 30, 2023 by the Italian Ministry of Ministry of University and Research (MUR).

The authors have no conflicts to disclose.

A. Fantasia: Data curation (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Software (equal); Writing – original draft (equal). F. Rovaris: Data curation (supporting); Formal analysis (supporting); Investigation (equal); Methodology (equal); Software (equal); Supervision (supporting); Writing – review & editing (supporting). O. Abou El Kheir: Conceptualization (supporting); Methodology (equal); Software (equal); Writing – review & editing (supporting). A. Marzegalli: Conceptualization (supporting); Data curation (supporting); Investigation (supporting); Validation (supporting). D. Lanzoni: Formal analysis (equal); Methodology (supporting); Software (supporting); Supervision (supporting); Writing – review & editing (supporting). L. Pessina: Data curation (supporting); Investigation (supporting); Visualization (supporting). P. Xiao: Conceptualization (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Supervision (equal); Writing – review & editing (supporting). C. Zhou: Data curation (equal); Methodology (equal); Supervision (supporting); Validation (supporting). L. Li: Conceptualization (equal); Methodology (equal); Software (equal); Supervision (supporting); Writing – review & editing (supporting). G. Henkelman: Conceptualization (equal); Funding acquisition (equal); Methodology (equal); Software (equal); Supervision (equal); Writing – review & editing (supporting). E. Scalise: Conceptualization (equal); Funding acquisition (equal); Investigation (equal); Methodology (equal); Project administration (equal); Resources (equal); Software (supporting); Supervision (equal); Validation (equal); Writing – review & editing (supporting). F. Montalenti: Conceptualization (equal); Funding acquisition (equal); Investigation (equal); Methodology (equal); Project administration (equal); Resources (equal); Supervision (equal); Validation (equal); Writing – original draft (equal); Writing – review & editing (equal).

The data that support the findings of this study are available from the corresponding author upon reasonable request. Also, the full database and potential file can be downloaded from the Repository quoted in Ref. 81.