# Unraveling the Crystallization Kinetics of the Ge2Sb2Te5 Phase Change
Compound with a Machine-Learned Interatomic Potential

Omar Abou El Kheir  Department of Materials Science, University of Milano-
Bicocca, Via R. Cozzi 55, I-20125 Milano, Italy  Luigi Bonati  Atomistic
Simulations, Italian Institute of Technology, Via Enrico Melen 83, 16152,
Genova, Italy  Michele Parrinello  Atomistic Simulations, Italian Institute of
Technology, Via Enrico Melen 83, 16152, Genova, Italy  Marco Bernasconi
Department of Materials Science, University of Milano-Bicocca, Via R. Cozzi
55, I-20125 Milano, Italy

###### Abstract

The phase change compound Ge2Sb2Te5 (GST225) is exploited in advanced non-
volatile electronic memories and in neuromorphic devices which both rely on a
fast and reversible transition between the crystalline and amorphous phases
induced by Joule heating. The crystallization kinetics of GST225 is a key
functional feature for the operation of these devices. We report here on the
development of a machine-learned interatomic potential for GST225 that allowed
us to perform large scale molecular dynamics simulations (over 10000 atoms for
over 100 ns) to uncover the details of the crystallization kinetics in a wide
range of temperatures of interest for the programming of the devices. The
potential is obtained by fitting with a deep neural network (NN) scheme a
large quantum-mechanical database generated within Density Functional Theory.
The availability of a highly efficient and yet highly accurate NN potential
opens the possibility to simulate phase change materials at the length and
time scales of the real devices.

This version of the article has been accepted for publication, after peer
review but is not the Version of Record and does not reflect post-acceptance
improvements, or any corrections. The Version of Record is available online
at: https://doi.org/10.1038/s41524-024-01217-6

Phase change materials, machine learning potentials, crystallization,
electronic memories, neural networks

##  I INTRODUCTION

In the last decades, the rise of the demand for data processing and storage
has stimulated a strong effort in the search of new computing architectures
and memory devices. Chalcogenide phase change materials Wuttig and Yamada
(2007) are at the heart of some of the most mature technologies suitable to
respond to these needs. Indeed, these materials are exploited in both emerging
non-volatile electronic memories, named phase change memories (PCM) Wuttig and
Yamada (2007); Noé  _et al._ (2017); Fantini (2020), and in neuromorphic and
in-memory computing devices Kuzum  _et al._ (2012); Tuma  _et al._ (2016);
Sebastian  _et al._ (2020). PCMs rely on a rapid (down to tens of ns) and
reversible transformation induced by Joule heating between the crystalline and
amorphous phases of the prototypical phase change compound Ge2Sb2Te5 (GST225)
Wuttig and Yamada (2007). Read out of the memory consists of the measurement
of the resistance of GST225 which differs by about three orders of magnitude
in the two phases Wuttig and Yamada (2007). Ge-rich GeSbTe alloys with
crystallization temperatures much higher than that of GST225 have also been
investigated for memories embedded in microcontrollers for automotive
applications Cappelletti  _et al._ (2020); Zuliani  _et al._ (2015). Moreover,
partial crystallization of the amorphous phase leads to different levels of
resistivity which is exploited in the realization of artificial synapses for
neuromorphic and in-memory computing Sebastian  _et al._ (2020); Zhang  _et
al._ (2019).
A key functional property for all these applications is the crystallization
kinetics of the amorphous phase between the glass transition (Tg) and the
melting (Tm) temperatures. This feature is, however, difficult to be
investigated experimentally because of the very high nucleation rates and
crystal growth velocities (a few m/s). Indeed, ultrafast differential scanning
calorimetry (DSC) was needed to measure the crystal growth velocity at the
high temperatures of interest for the operation of the devices Orava  _et al._
(2012). Information on the crystallization kinetics was, however, inferred
from DSC under several assumptions on the crystallization mechanisms based on
classical nucleation theory which require further validations.
On the other hand, atomistic insights on the early stage of the
crystallization process have been provided by molecular dynamics (MD)
simulations based on density functional theory (DFT) Zhang  _et al._ (2019);
Hegedüs and Elliott (2008); Lee and Elliott (2011a); Kalikka  _et al._ (2016);
Ronneberger  _et al._ (2015, 2018); Xu  _et al._ (2022); Rao  _et al._ (2017);
Lee and Elliott (2017). Several works on GST225 revealed that the high
nucleation rate can be ascribed to the presence in the supercooled liquid of
four-membered rings, which are the same building blocks of the cubic rocksalt
crystal, that act as seeds for the formation of crystalline nuclei Zhang  _et
al._ (2019); Hegedüs and Elliott (2008). Moreover, the high crystal growth
velocity is ascribed to the high fragility of the supercooled liquid which can
sustain high atomic mobility down to temperatures close to Tg, where the
thermodynamical driving force for crystallization is also high Zhang  _et al._
(2019). Although they provided crucial information on the early stage of
crystallization, DFT-MD methods suffer from limitations in system size and in
simulation time that prevent to address some important issues for the
operation of the memory devices. A DFT-MD study of the crystallization
kinetics on an extended temperature range to test the applicability of
classical nucleation, for instance, is still lacking for GST225.
In the last few years, the development of interatomic potentials based on the
fitting of a large DFT database by machine learning techniques emerged as a
viable approach to overcome these limitations of DFT-MD and enlarge the scope
of DFT methods Behler and Parrinello (2007); Behler (2016); Deringer  _et al._
(2019); Wang  _et al._ (2018). Concerning phase change materials, machine
learning schemes based on Neural Network (NN) methods have been exploited to
study the crystallization of the phase change materials GeTe Sosso  _et al._
(2012a, 2013); Lee  _et al._ (2020) and Sb Dragoni  _et al._ (2021); Shi  _et
al._ (2021). NN simulations of GeTe also allowed addressing the study of the
aging of the amorphous phase Gabardi  _et al._ (2015) that leads to an
increase of the electrical resistance with time (drift) which is detrimental
for the operation of the memory devices Raty  _et al._ (2015). For GST225, a
machine learning interatomic potential was recently developed with the
Gaussian Approximations Potential (GAP) method Mocanu  _et al._ (2018) which,
however, suffers from some inaccuracies in reproducing the DFT results on the
structural properties of the amorphous phase such as the fraction of homopolar
bonds which are believed to play a crucial role in the aging process Gabardi
_et al._ (2015); Raty  _et al._ (2015). Improvements of this potential have
been, however, very recently achieved Zhou  _et al._ (2023).
In this paper, we report on the development of an interatomic potential for
GST225 within the NN scheme implemented in the DeePMD code Wang  _et al._
(2018); Zhang  _et al._ (2018). We first validated the potential on the
structural, dynamical and thermodynamical properties of the liquid, amorphous
and crystalline phases. Then, we employed the NN potential to study the
crystallization kinetics over the wide temperature range of interest for the
operation of the devices, aiming at assessing the applicability of classical
nucleation theory. The simulations reveal that crystallization kinetics in the
temperature range 500-650 K is diffusion-controlled with an activation energy
corresponding to that of the self-diffusion coefficient. We also show that a
modified form Jackson (2004) of the phenomenological Wilson-Frenkel formula
Wilson (1900); Frenkel (1932) is suitable to fit the data on a wider
temperature range.

We remark that the efficient implementation of the DeePMD scheme allows
simulating tens of thousands of atoms for tens of nanoseconds at an affordable
computational cost. We envisage that this feature could be further exploited
to address some important issues for the operation of ultrascaled devices such
as the effect of confinement and nanostructuring on the crystallization
kinetics Kooi and Wuttig (2020), or the possible existence of a strong-to-
fragile transition in the supercooled liquid close to Tg which is particularly
relevant for the aging of the amorphous phase Pries  _et al._ (2022), just to
name a few.

##  II RESULTS

###  II.1 Generation and Validation of the Neural Network Potential

The NN potential was obtained by fitting DFT (see Methods) energies, forces
and the stress tensor of a database containing about 180000 supercell models
(configurations) of GST225 in the liquid, amorphous, cubic and hexagonal
phases by using the DeePMD-kit open-source package Wang  _et al._ (2018);
Zhang  _et al._ (2018).
Initially, the training set consisted of about 5000 configurations. Then, the
NN potential was refined by expanding the database in an iterative process
with atomic configurations generated from DFT-MD trajectories to enhance the
description of specific properties and from NN-MD trajectories whose energy
was badly described by the intermediate versions of the potential. In the
final database of about 180000 configurations, we covered a wide range of
thermodynamical conditions, i.e., different temperatures up to 2000 K and
several densities in the range 0.026-0.036 atom/Å3. Details of configurations
in the database at different conditions are reported in Tab. 1 while
information on the NN scheme are given in section on Methods. Each
configuration refers to the DFT energy, forces and stress of 108-atom cells
for the liquid, amorphous and hexagonal crystalline phases and of a 57- or
98-atom cell for the cubic crystal. The configurations of the supercooled
liquid were extracted from simulations of several tens of ps at fixed
temperature in the range 500-900 K, quenching refers to simulations in which
the system was cooled very rapidly from 990 to 300 K and metadynamics refers
to biased simulations Laio and Parrinello (202) to enhance the sampling in a
specific region of the phase space. The metadynamics method is based on the
identification of appropriate order parameters or collective variables (CVs)
that describe the slow modes of the process. An external potential is then
added to enhance the fluctuations of the selected CVs. The method allows
sampling the free energy surface by overcoming activation barriers much larger
than the thermal energy in a short time span. As a result, it allows to
collect a more heterogeneous set of configurations than standard MD
simulations, obtaining robust potentials that can also describe phase
transitions Bonati and Parrinello (2018); Niu  _et al._ (2020). In particular,
we performed metadynamics simulations in the supercooled liquid by using the
coordination numbers as CVs. The addition of such configurations to the
training set was crucial to reproduce the fraction of homopolar bonds in the
amorphous phase.

Table 1: Details of the database used for the training of the NN potential. phase/simulations | number of configurations
---|---
liquid | 4489
supercooled liquid | 64440
quenching | 70463
amorphous | 3620
cubic crystal | 8764
hexagonal crystal | 7935
metadynamics | 18524

The energy root-mean-square error (RMSE) between the DFT values and those
predicted by the NN potential is 8.4 and 8.6 meV/atom for the training and
test sets. The RMSE on forces is 159 meV/Å for both sets. The distributions of
NN errors on energies and forces are given in Fig. 1. The typical average
error obtained with DeePMD for highly disordered phases of multicomponent
systems like ours (i.e liquid and/or amorphous phases) are in the range 2-7
meV/atom and 90-145 meV/Å Urata (2022); Zhai and Wang (2023); Liu  _et al._
(2023); Huang  _et al._ (2023). The NN potential has been validated on the
properties of the liquid, amorphous, cubic and hexagonal phases as described
in the separate sections below. All the NN-MD simulations were performed with
the LAMMPS code Thompson  _et al._ (2022) and a Nosé-Hoover thermostat Nosé
(1984); Hoover (1985).

Figure 1:  Accuracy of the neural network potential. Cumulative fraction of
the absolute errors of the NN potential in training and test data sets for a)
the energies per atom ($\Delta$E=$|$EDFT-ENN$|$) and b) forces
($\Delta$F=$|$FDFT-FNN$|$). c) The distribution of the absolute error as a
function of DFT energy.

The Liquid Phase: The structural properties of the liquid phase of GST225 have
been obtained from NN-MD simulation with a 999-atom supercell at 990 K and
compared with those obtained from DFT-MD simulation with a 216-atom supercell
at the same temperature. In both models, we used the experimental density of
the amorphous phase of 0.0309 atoms/Å3 Njoroge  _et al._ (2002) which is very
close to that of the liquid at 893 K (0.0307 atoms/Å3 Schumacher  _et al._
(2016)). To generate the liquid models we first performed a 7 ps long MD
simulation at 2000 K to properly randomize the atoms in the box. Then, we
performed a second equilibration at 990 K for 10 ps. Finally, we evaluated the
structural properties over trajectories 10 ps long.

The pair correlation functions, the bond angles distribution function and the
distribution of coordination numbers are compared with DFT results in Fig. 2.
The structural properties obtained from NN-MD simulations are in excellent
agreement with those obtained from DFT-MD which suggests that the NN potential
reliably describes the liquid phase.

Figure 2:  Structural and dynamical properties of the liquid phase. Structural
and dynamical properties of GST225 in the liquid phase from DFT (blue curves)
and NN (red curves) simulations with a 216-atom or 999-atom models,
respectively. a) Total and partial radial distribution functions at 990 K. The
position of the first maximum and minimum of each partial correlation function
are given in the Supplementary Table 1. b) Angle distribution function
resolved per central atomic species at 990 K. The data were normalized to the
number of triplets in each model. c) Distribution of coordination numbers
resolved per chemical species at 990 K computed by integrating the partial
pair correlation functions up to a bonding cutoff which corresponds to 3.2 Å
for all pairs except for Sb-Te for which we use 3.4 Å as was done in Ref.
Caravati  _et al._ (2009). d) The Arrhenius plot of the diffusion coefficient
$D$ as a function of temperature compared with previous DFT results with the
different functionals PBE Rizzi  _et al._ (2014), PBESol Flores-Ruiz and
Micoulaut (2018) and vdW-DF2 Schumacher  _et al._ (2016)) (see Supplementary
Table 2 for the corresponding Arrhenius parameters). The NN data of $D$ below
700 K are shown with a black border as they have been excluded from the
Arrhenius fit.

Regarding the dynamical properties, we computed the self-diffusion coefficient
$D$ from NN-MD simulations at several temperatures above Tm and below Tm in
the supercooled liquid spanning the range 500-1200 K. The density was fixed to
the value of the experimental amorphous phase (0.0309 atoms/Å3) as it was done
in the previous DFT-PBE work Rizzi  _et al._ (2014) that we take as a
reference for the validation. The self-diffusion coefficients was obtained
from the mean square displacement (MSD) and the Einstein relation MSD=6$D$t
from equilibrated trajectories at constant energy over time intervals from 40
ps at high temperatures to 300 ps at low temperatures. At temperatures above
700 K the data can be fitted with the Arrhenius function
$D$=$D_{0}$exp$({-E_{a}}/{k_{B}T})$ (see Fig. 2d) with Ea= 0.267 eV and
$D_{0}$=1.03x10-3 cm2/s which are similar to other DFT values reported in the
literature (see Supplementary Table 2). Below 700 K deviations from the
Arrhenius law are present due to the fragility of the system. For a fragile
liquid, the self-diffusion coefficient can be fitted in a wider range of
temperatures with the Cohen-Grest (CG) formula Cohen and Grest (1979) as
$\log_{10}(D(T))=A-{2B}/(T-T_{0}+[(T-T_{0})^{2}+4CT]^{1/2})$, which for GST225
yields A=-2.45, B=602 K, C=17.3 K and $T_{0}$=330.6 K (see Supplementary
Figure 1). We have chosen the CG formula because it was used in the
experimental work in Ref. Orava  _et al._ (2012) to fit the kinetic prefactor
in the crystal growth velocity inferred from DSC data on which we will come
back later in the discussion of the crystallization kinetics. The diffusion
coefficient as a function of temperature resolved for the different atomic
species is also reported in the Supplementary Figure 1.

Figure 3:  Melting and glass transition temperatures. a) The potential energy
as a function of time at different temperatures of the cubic-liquid interface
model. Simulations are performed at constant volume at a density of 0.0309
atom/Å3 which corresponds to the theoretical equilibrium density of the cubic
crystal (see Sec. 2.3). The simulation cell with edges
6.15$\times$6.15$\times$11.08 nm3 was initially prepared with an interface
lying on the $xy$ plane separating a 2 nm thick slab in the cubic phase from a
liquid slab 9 nm thick. b) Snapshots along the trajectories of panel a)
showing the initial configuration and the movement of the interface between
the two phases. The melting temperature is estimated to be in the range
920-930 K. c) Energy (potential plus kinetic) of the supercooled liquid as a
function of temperature in simulations at constant volume. The lines are a
quadratic fit below and above Tg. The resulting heat capacity at constant
volume is given in the inset. The data at each temperature are obtained by
averaging over 40 ps simulations of a 999-atom cell, initially equilibrated at
1000 K and then cooled down to the target temperature in 20-100 ps.

As a further step in the validation of the NN potential, we estimated the
melting temperature of the crystalline cubic phase of GST225 by means of the
phase coexistence method Allen and Tildesley (2017). To this aim, we prepared
a 12960-atom model of the cubic-liquid interface initially set at 900 K and at
the theoretical equilibrium density of the cubic phase. The crystal-liquid
interface corresponds to the (001) surface of the cubic crystal. Then we
carried out several independent simulations at constant volume and at
different temperatures in the range 880-950 K to monitor the potential energy
as a function of time (Fig. 3a), starting from the initial configuration in
Fig. 3b. At temperatures above Tm we expect the crystalline region to melt as
it is the case at 950 and 940 K (see snapshot II in Fig. 3b ) while for
temperatures below Tm we expect the crystalline region to grow, as it is
indeed the case below 920 K. This set of simulations suggests that the melting
point of the NN potential is within the range 920-930 K which is very close
the experimental value of 900 K Schumacher  _et al._ (2016). As stated above,
these results refer to simulations with the density fixed to that of the cubic
phase, also for the liquid. A better estimation of Tm should, however, be
obtained from NPT simulations to describe the density change across melting
and the thermal expansion of the crystal. NPT simulations with the PBE
functional are, however, problematic as they underestimate significantly the
equilibrium density due to the coalescence of nanovoids as it was observed for
GeTe Sosso  _et al._ (2012b). We then repeated the same calculations in the
temperature range 925-940 K by adding van der Waals interactions (vdW)
according to Grimme Grimme (2006) (D2) which prevents the coalescence of
nanovoids that form by decreasing the density as it was also reported for GeTe
Sosso  _et al._ (2012b). Starting from the same initial configuration, we
first carried NPT-MD simulation allowing the cell edges to change at fixed
angles. The equilibrium density is reached on a time scale of 30 ps which is
much shorter than that required for crystallization. From several NPT
simulations at different temperature we estimated Tm=940 K which is close to
the previous NVT results with no vdW interactions (see Supplementary Figure
2). At this temperature, the latent heat of melting is $\Delta H_{m}$=163
meV/atom which is close to the experimental values of 120 or 173 meV/atom
Senkader and Wright (2004); Zhao  _et al._ (2022). To calculate the latent
heat at 940 K, we performed NPT simulation for a 999-atom model of the liquid
and a 900-atom special quasi-random-structure model Zunger  _et al._ (1990)
for the cubic crystalline phase.

To assess the error in the latent heat introduced in case one does not take
into account the density change upon melting (i.e. NVT simulations), we also
computed the energy of the liquid phase at the equilibrium density of the
cubic crystalline phase at 940 K. The resulting heat of melting at constant
volume is $\Delta E_{m}$ = 154 meV/atom. The value of $\Delta E_{m}$ computed
at constant volume without vdW interactions is instead 166 meV/atom. We expect
this latter value to differ from the latent heat computed in the NPT ensemble
by a similar error of about 10 meV/atom as found for the simulations with vdW
interactions. We also estimated the change in Tm due to the choice of NVT
conditions by using the Clausius-Clapeyron equation with the calculated latent
heat at 940 K, and the theoretical equilibrium density of the two phases
within NN+D2 simulations. This yields a change in Tm of about 10 K.

All the results reported in the following still refer to simulations without
vdW interactions as we want to validate the NN potential over DFT-PBE data.
vdW interactions might be later added in simulations with the NN potential
even by using different schemes.

As a final result on the supercooled liquid, we report in Fig. 3c the energy
as a function of temperature computed at a constant density of 0.0309
atoms/Å3. The energy displays clearly two slopes as a function of temperature
which give rise to a jump in the specific heat shown in the inset of Fig. 3c.
The jump in Cv at about 500 K can be identified with the glass transition
temperature which turns out to be very similar to the latest experimental
value of Tg=473 K reported in Ref. Pries  _et al._ (2019), where it was also
proposed that crystallization in DSC occurs below Tg. The decrease of Cv with
temperature in the supercooled liquid above Tg is another feature typical of
fragile liquids Angell (1995). We also repeated the same simulations at
constant (zero) pressure with the NN+D2 potential, the resulting enthalpy and
$C_{p}$ as a function of temperature, reported in the Supplementary Figure 3,
yield a very similar estimate of Tg.

The Amorphous Phase: The NN potential was then validated on the structural
properties of the amorphous phase. A 999-atom model of amorphous GST225 at the
experimental density (0.0309 atom/Å3) was generated by quenching from the melt
at 990 K to 300 K in 100 ps. The pair correlation functions, the bond angles
distribution functions and the distribution of the coordination numbers are
compared in Fig.4 with DFT results from four 216-atom models generated with
the same quenching protocol and at the same density of the NN model. The
resulting average partial coordination numbers are compared in Table 2. The
comparison of the structural properties of two NN models of different sizes
containing 999 or 7992 atoms (see Supplementary Figure 4) shows that the
structural properties are well converged in the 999-atom cell. The agreement
between NN and DFT simulations is overall excellent, including small but very
important details such as the fraction of homopolar Ge-Ge and Sb-Sb bonds.

Figure 4:  Structural properties of the amorphous phase. Structural properties of amorphous GST225 at 300 K from NN (red curves) and DFT (blue curves) simulations. DFT data are averaged over four independent 216-atom models while NN data are obtained from a 999-atom model. a) Total and partial radial distribution functions. The position of the first maximum and minimum of the different functions is reported in the Supplementary Table 3. b) Angular distribution function resolved per central atomic species. The data are normalized to the number of triplets in each model. c) Distribution of coordination numbers resolved per chemical species, obtained by integrating the pair correlation functions up to a bonding cutoff corresponding to 3.2 Å for all pairs except for Sb-Te for which a longer cutoff value of 3.4 Å was used. The spread over the four DFT models are indicated by error bars. d) Distribution of the local order parameter $q$ for tetrahedricity for four-coordinated Ge atoms for the NN model and for each of our four DFT models. Table 2: Average coordination numbers for different pairs of atoms in amorphous GST225 at 300 K generated from DFT (in parenthesis) and NN simulations. Error bars are obtained from the analysis of four models. - | Ge | Sb | Te
---|---|---|---
Total |  $4.19\pm$0.01 (4.17$\pm$0.1) | 4.46$\pm$0.02 (4.40$\pm$0.05) | 3.16$\pm$0.02 (3.16$\pm$0.05)
With Ge | 0.33$\pm$0.04 (0.34$\pm$0.05) | 0.33$\pm$0.02 (0.26$\pm$0.05) | 1.41$\pm$0.02 (1.42$\pm$0.02)
With Sb | 0.33$\pm$0.02(0.26$\pm$0.05) | 0.51$\pm$0.03 (0.51$\pm$0.1) | 1.45$\pm$0.02 (1.45$\pm$0.05)
With Te | 3.53$\pm$0.06 (3.57$\pm$0.05) | 3.64$\pm$0.05 (3.63$\pm$0.1) | 0.29$\pm$0.02 (0.28$\pm$0.02)

The structure of amorphous GST225 is similar to that emerged in previous DFT
works Caravati  _et al._ (2007); Akola and Jones (2007). Te atoms are mostly
three-fold coordinated in a pyramidal geometry, Sb atoms are both three-fold
coordinated in a pyramidal geometry (three bonding angles of 90∘) and four- or
five-fold coordinated in a defective octahedral environment (octahedral
bonding angles but coordination lower than six), most of the Ge atoms are in
pyramidal or defective octahedral geometry while a minority fraction of Ge
atoms are in tetrahedral geometries. The bonding geometry is revealed by the
coordination numbers and by the angle distribution function where the peak at
about 90o is due to pyramidal and defective octahedral configurations, while
the weak peak at about 170o is due to axial bonds in defective octahedra. The
shoulder at about 109o is due to tetrahedra which are favored by homopolar Ge-
Ge bonds Deringer  _et al._ (2014) that are present in the liquid and survive
in the amorphous phase due to fast quenching. A quantitative measure of the
fraction of tetrahedral environments can be obtained from the local order
parameter $q$ introduced in Ref. Errington and Debenedetti (2001). It is
defined as $q=1-\frac{3}{8}\sum_{i>k}(\frac{1}{3}+\cos\theta_{ijk})^{2}$,
where the sum runs over the pairs of atoms bonded to a central atom $j$ and
forming a bonding angle $\theta_{ijk}$.The order parameter evaluates to $q$=1
for the ideal tetrahedral geometry, to $q$=0 for the 6-fold coordinated
octahedral site, to $q$=5/8 for a 4-fold coordinated defective octahedral
site, and $q$=7/8 for a pyramidal geometry. The distribution of the local
order parameter $q$ for tetrahedricity is reported in Fig. 4d for four-
coordinated Ge atoms. The bimodal shape corresponds to tetrahedral and
defective octahedral geometries. We quantified the fraction of Ge atoms in a
tetrahedral environment by integrating the $q$ parameter between 0.8 and 1 as
discussed in previous works Spreafico  _et al._ (2011). In the NN model, 30%
of Ge atoms are in the tetrahedral geometry to be compared with an average
value of 23% in our DFT models. DFT calculations in literature report a
fraction of tetrahedral Ge in the range 27-35% Caravati  _et al._ (2007);
Akola and Jones (2007). We remark, however, that a very recent paper Guerboub
_et al._ (2023) reports on the comparison of two 504-atom models of amorphous
GST225 generated with two different pseudopotentials for Ge, based either on
the Troullier-Martins (TM)Troullier and Martins (1991) or GTH schemes.
Although the two models reproduce similarly well the experimental neutron
diffraction data, the fraction of tetrahedral Ge is 36 $\%$ in the GTH model
while it goes up to 65 $\%$ in the TM model which better reproduces the Ge-Ge
and Ge-Te bond lengths inferred from extended-x-ray-fine-structure
measurements Guerboub  _et al._ (2023). A similarly larger fraction of
tetrahedra have been obtained in Ref. Micoulaut and Flores-Ruiz (2021).
Information on the medium range order is provided by the distribution of rings
length which in GST225 is dominated by the four-membered rings Hegedüs and
Elliott (2008); Caravati  _et al._ (2007). This feature has been regarded as
the precursor for rapid crystallization as the four-membered ring is also the
structural unit of the rocksalt crystalline phase Hegedüs and Elliott (2008).
The distribution of rings length reported in the Supplementary Figure 5 for
the NN and DFT models shows that the NN potential is able to reproduce also
this crucial feature. Another useful descriptor of the structure of the
amorphous phase is given by the angular-limited three-body correlation
function (ALTBC) Bichara  _et al._ (1992) which highlights the presence of a
short and a long axial bond in the defective octahedral configurations
Schumacher  _et al._ (2016). The ALTBC functions for NN and DFT simulations of
amorphous GST225 are in good agreement as well (see Supplementary Figure 6).
The NN potential is thus able to reproduce very well the structural properties
of the amorphous phase including very crucial details such as the fraction of
Ge atoms in tetrahedral configurations and the fraction of homopolar Ge-Ge
bonds which might rule the aging of the amorphous phase and the consequent
increase (drift) of the electrical resistance with time, as it occurs in the
amorphous phase of the parent compound GeTe Gabardi  _et al._ (2015); Raty
_et al._ (2015). Concerning the dynamical properties, the NN potential
reproduces well the DFT phonon density of states (DOS) of amorphous GST225 as
shown in Fig. 5a.

Figure 5:  Phonons. a) Phonon density of states (DOS) of the amorphous phase
obtained from $\Gamma$-point frequencies of a 216- and 999-atom supercell with
the NN potential and of a 216-atom supercell with DFT. b) DFT and NN phonon
DOS of the cubic phase of GST225. The phonons were computed at the $\Gamma$\-
point of a 270-atom supercell. Phonon frequencies in panel a) and b) are
obtained by diagonalizing the dynamical matrix from finite atomic
displacements. Each frequency is then broadened by a Gaussian function with a
width of 1 cm-1. c) Phonon density of states of the crystalline hexagonal
phase of GST225 computed with the NN+D2 potential and by DFT+D2 from Ref.
Campi  _et al._ (2017) (Kooi phase, see text). The NN+D2 DOS for the hexagonal
phase are obtained from the the force constant matrix of a 12x12x4 supercell
and the Phononpy code Togo and Tanaka (2015).

The Crystalline Phases: At normal conditions, GST225 crystallizes in a
hexagonal phase (space group $P{\bar{3}}m1$)Petrov  _et al._ (1968); Kooi and
De Hosson (2002); Matsunaga  _et al._ (2004). However, the amorphous phase
crystallizes in a metastable cubic rocksalt crystal which is the phase of
interest for the operation of the memory devices. The metastable cubic phase
consists of a NaCl structure with the anionic sublattice occupied by Te and
the cationic one occupied by Ge, Sb and 20% of vacancies. In this benchmark,
the cubic phase was modeled by a 300-site supercell (the same as in Ref.
Caravati  _et al._ (2009)) with 30 vacancies and 270 atoms at the
stoichiometric GST225 composition. The hexagonal phase contains instead 9
atoms in the primitive unit cell arranged along the $c$ direction with a
ABCABC stacking. Each formula unit forms a lamella separated from the others
by a so-called vdW gap, although the interlamella interaction is not just a
vdW contact as discussed in Ref. Wuttig  _et al._ (2018). Three different
models of hexagonal GST225 have been proposed in literature differing in the
distribution of Sb/Ge atoms in the cation sublattices Petrov  _et al._ (1968);
Kooi and De Hosson (2002); Matsunaga  _et al._ (2004). Here, we considered the
Kooi stacking Kooi and De Hosson (2002) with Sb atoms occupying the cation
planes close to the vdW gap. We computed the equation of state at zero
temperature of the cubic and hexagonal phase by fitting the energy-volume
points by the Birch-Murnaghan formula (see Supplementary Figures 7-8). The
resulting parameters at equilibrium are compared with DFT data in Tab. 3. The
theoretical NN (and DFT-PBE) equilibrium density of the cubic phase (0.0309
atom/Å3) turns out to be equal to the experimental density of the amorphous
phase. By adding vdW interactions (D2)Grimme (2006) the equilibrium density of
the cubic phase raises to 0.0332 atom/Å3 which is closer to the experimental
value of 0.0328 atom/Å3 Matsunaga  _et al._ (2004). The NN potential
reproduces well also the phonon DOS obtained by DFT as shown in Fig. 5b,c. In
the hexagonal phase, GST225 features phonon instabilities when the DFT-PBE
scheme is applied Campi  _et al._ (2017) which, as expected, are also present
in NN calculations as it should be. These instabilities at the PBE level are
removed by including the semiempirical vdW correction due to Grimme (D2) Campi
_et al._ (2017). Hence, we calculated the phonon dispersion relations with the
NN+D2 potential by employing the Phononpy code Togo and Tanaka (2015). The
resulting phonon DOS are compared with previous DFT+D2 results Campi  _et al._
(2017) in Fig. 5c, while the phonon dispersion relations are compared with
DFT+D2 results Campi  _et al._ (2017) in the Supplementary Figure 9. Although
reasonable, the agreement with DFT results is less satisfactorily for the
phonon dispersion relations than for all the other properties analyzed so far.
In fact, it is usually rather difficult to reproduce phonon dispersion
relations by a NN potential not explicitly devised for this property.

In summary, although the RMSE for the energies and forces are not very small
(8 meV/atom and 159 meV/Å), albeit similar to other NN potentials in
literature for disordered multicomponent materials, the validation of the
potential over the properties of liquid, amorphous and crystalline phases is
excellent. Overall, we judge that our potential is sufficiently accurate to
address the study of many properties of this system including the
crystallization kinetics which is the subject of the next section.

Table 3: Fitting parameters of the Birch-Murnaghan equation of state of cubic
and hexagonal GST225 from NN and DFT-PBE calculations at zero temperature. The
energy (E0), volume (V0), bulk modulus (B) and derivative of B (B′) at
equilibrium are reported.

| E0 (eV/atom) | V0 (Å3/atom) | B (GPa) | B′
---|---|---|---|---
DFT cubic | -180.7424 | 32.42 | 25.1 | 7.53
NN cubic | -180.7408 | 32.35 | 22.5 | 9.41
DFT hexagonal | -180.7883 | 31.13 | 19.8 | 21.75
NN hexagonal | -180.7890 | 32.13 | 18.7 | 23.38

###  II.2 Simulation of the Crystallization Process

As a first application of the NN potential for GST225, we studied the kinetics
of the crystallization process in the liquid phase supercooled below Tm and in
amorphous phase overheated above Tg by evaluating the crystal growth velocity
($v_{g}$) as a function of temperature in the range 500-940 K.

The NN potential allowed us to perform simulations with several thousands of
atoms for overall 100 ns that provided the crystal growth velocity $v_{g}$ as
a function of temperature with greater details than reported previously by DFT
simulations. To this end, we first extended to a wider range of temperatures
the simulations of the 12960-atom model of the liquid-crystal interface
discussed in previous sections. The model was first quenched in 40-80 ps from
900 K to each target temperature to monitor the evolution of the crystalline
slab. The number of crystalline atoms is quantified by using the local order
parameter for crystallinity Q${}_{4}^{dot}$ Steinhardt  _et al._ (1983) that
is suitable to distinguish atoms in the crystalline phase from atoms in
liquid/amorphous environments as shown in the Supplementary Figure 10. Then,
the evolution of the crystalline interface was monitored to estimate the
crystal growth velocity as $v_{g}=dL(t)/dt$ where $L(t)$ is the effective
(half) thickness of the crystalline slab given by
$L\left(t\right)=N(t)/{2A\rho_{NN}^{cubic}}$, where $N(t)$ is the number of
atoms in the crystalline slab, $\rho_{NN}^{cubic}$ is the theoretical
equilibrium density of the cubic phase, $A$=6.15 $\times$ 6.15 nm2 is the
cross-section of the cell orthogonal to the growth direction, and the factor
two at the denominator accounts for the presence of two growing surfaces. The
evolution of $L(t)$ as a function of time at several temperatures is reported
in the Supplementary Figure 11. Since the crystal-liquid interface lays on the
(001) plane of the cubic phase, the crystal growth velocity corresponds to the
growth along the [001] direction of the cubic crystal. Before analyzing the
results, we verified that our thermostat was effective in getting rid of the
latent heat of crystallization released during crystal growth that in the real
system diffuses away very fast due to electronic thermal conductivity of the
liquid. To this aim, we considered slices (bins) 10 Å wide at different
distances from the liquid-crystal interface as shown in the Supplementary
Figure 12. The local temperature and the fraction of crystalline atoms in the
different slices are shown in the Supplementary Figures 13-17 for different
average temperatures. The local temperature is indeed rather uniform across
the liquid-like slab at different distances from the surface. The crystal
growth velocity as a function of temperature is often described by the
phenomenological Wilson–Frenkel formula (WF) Wilson (1900); Frenkel (1932)
$v_{g}=u_{kin}(1-exp(-\Delta\mu/k_{B}T))$ where $u_{kin}$ is a kinetic
prefactor and $\Delta\mu$ is the free energy difference between the
crystalline and supercooled liquid phases.

For a diffusion-controlled growth, the kinetic prefactor is given by
$u_{kin}=6Ddf/\lambda^{2}$, where $\lambda$ is the typical jump distance of
atoms in the elementary diffusion process, $D$ is the diffusion coefficient,
$d$ is the interlayer spacing along the growth direction and $f$ represents
the fraction of surface sites where a new atom can be incorporated Kerrache
_et al._ (2008). The WF formula is typically adequate to describe a continuous
growth of a rough surface as it seems to be the case here; a snapshot of the
growing surface is given in the Supplementary Figure 18.

By setting $d=\lambda$ and $f$=1 the kinetic prefactor has the form
$u_{kin}=\frac{6D}{\lambda}$ that we used for instance in a previous
simulation of the crystallization of GeTe from a crystal-liquid interface
(heterogeneous crystallization) Sosso  _et al._ (2015). We spend a few words
on the justification of the WF formula in view of possible different choices
for the kinetic prefactor that we will discuss later on. The crystal growth
velocity can be expressed as $v_{g}=df(\kappa^{+}-\kappa^{-}$), where
$\kappa^{+}/\kappa^{-}$ are the rate of attachment/detachment of an atom
to/from the crystalline surface. In turn $\kappa^{+}=\nu exp(-\Delta
G^{*}/k_{B}T)$ and $\kappa^{-}=\nu exp(-\Delta
G^{*}/k_{B}T-\Delta\mu/k_{B}T)$, where $\Delta G^{*}$ is the activation energy
for the attachment to the surface of an atom from the liquid with an attempt
frequency $\nu$. By assuming that $\kappa^{+}$ is equal to the rate of a jump
in the diffusion process of a single atom in the liquid, it can be written as
$\kappa^{+}=6D/\lambda^{2}$. This approximation leads to
$v_{g}=6Ddf/\lambda^{2}(1-exp(-\Delta\mu/k_{B}T))$. We remark, however, that
$f$ should also include possible other corrections to the sticking. Due to the
uncertainties in the form of the kinetic prefactor, $\lambda$ is typically
considered as a fitting parameter.

We attempted to reproduce the crystal growth velocity extracted from MD
simulations with the WF relation and $u_{kin}=\frac{6Dd}{\lambda^{2}}$ (i.e.
$f$=1) with d=3.0 Å for the growth along the [001] direction of the cubic
phase, by using the theoretical $D$ obtained from the CG fit of the MD data
and by using for $\Delta\mu$ the expression given by Thompson and Spaepen
$\Delta\mu(T)=\frac{\Delta H_{m}(T_{m}-T)}{T_{m}}\frac{2T}{T_{m}+T)}$ Thompson
and Spaepen (1979). We set $\Delta H_{m}$ = 166 meV/atom, $T_{m}$ = 940 K as
estimated from our MD simulation (see previous sections). We also checked that
the Thompson-Spaepen formula is fairly accurate in our case by computing
$\Delta\mu$ from the integration of the specific heat as $\Delta\mu(T)=\Delta
H_{m}(1-\frac{T}{T_{m}})-\int^{T_{m}}_{T}\Delta
C_{p}dT+T\int^{T_{m}}_{T}\frac{\Delta C_{p}}{T}dT$, where $\Delta C_{p}$ is
the difference in $C_{p}$ between the amorphous and the crystalline phases. We
approximated $\Delta C_{p}$ with $\Delta C_{v}$ which was in turn computed
from the caloric curve at fixed density as discussed in Sec. 2. The resulting
$\Delta\mu$ in the temperature range of interest is nearly indistinguishable
from that obtained from the Thompson-Spaepan approximation as shown in the
Supplementary Figure 19. Actually, we have not been able to reasonably fit the
data over all temperatures by using just $\lambda$ as a free parameter.
Therefore we restricted the fit to the lower temperatures (below 650 K) which
yields the WF curve shown in Fig. 6a with a physically reasonable value of
$\lambda$= 2.42 Å. The WF curve reproduces the crystal growth velocity at low
temperatures but largely overestimates $v_{g}$ at high temperatures. We first
discuss the behavior at low temperatures below 650 K. The good fitting with
the WF formula means that $u_{kin}$ is thermally activated with an activation
energy close to that of the self-diffusion coefficient. We remark that the
self-diffusion coefficient describes the long-scale atomic diffusion process
of individual atom while the term that enters in the kinetic prefactor
$u_{kin}$ for crystallization is actually an effective diffusion coefficient
$D_{eff}$ that might embody a more short-ranged atomic motion. In glasses
close to Tg, the secondary $\beta$-relaxation is known to be the dominant
source of atomic dynamics while the slower $\alpha$ relaxation controls the
long range atomic diffusivity Cavagna (2009). At high temperature far from Tg,
the $\alpha$ and $\beta$ relaxations actually coincide Cavagna (2009). The
possibility that $\beta$-relaxation might enhance the crystallization kinetics
of phase change materials close to Tg has been put forward very recently Peng
_et al._ (2020); Cheng  _et al._ (2022). Indeed, the presence of a
$\beta$-relaxation process in phase change materials has been recently
identified from the so called $\beta$-wing in the temperature dependence of
the loss modulus measured by dynamical mechanical spectroscopy (DSM) Peng  _et
al._ (2020). Evidences of a link between the crystallization speed and the
presence of $\beta$ relaxation have been provided very recently for the
eutectic alloy Ge15Sb85 Cheng  _et al._ (2022). Molecular dynamics simulations
have also revealed that the stabilization of the amorphous phase of Sb in
ultrathin films is possibly due to the reduction of the $\beta$-relaxation
dynamics due to the confinement Dragoni  _et al._ (2021). It is therefore of
interest to investigate whether $\beta$-relaxation might also be of relevance
for crystal growth in GST225. The $\beta$-relaxation can be detected in MD
simulations by looking at the intermediate scattering function Cavagna (2009),
as we did for instance for the GeTe phase change compound in our previous work
Sosso  _et al._ (2014). There, we also used four-point correlation functions
Donati  _et al._ (2002) and isoconfigurational analysis Widmer-Cooper and
Harrowell (2007) to investigate dynamical heterogeneities in the supercooled
liquid phase which is associated with the fragility and the breakdown of the
Stokes-Einstein relation between $D$ and the viscosity. The availability of a
NN potential for GST225 now allows for an extension of our previous analysis
on GeTe to the ternary compound. However, we leave this new chapter on the
properties of the supercooled liquid phase for a future work as the present
one is already very dense of information. Therefore, here the
$\beta$-relaxation is addressed in a simpler manner by looking at the MSD as a
function of time and at different temperatures. A plateau in the MSD plotted
in a log-log scale close to Tg is typical of a two steps relaxation dynamics
with a faster $\beta$-relaxation and a slower $\alpha$-relaxation that
controls the long range atomic diffusivity after the plateau Cavagna (2009).
In GST225, a clear plateau is present at 400 K, although an inflection starts
to appear in the MSD at 500 K which is the lowest temperature at which we have
investigated the crystallization kinetics (see Supplementary Figure 20). This
suggests that the $\beta$-relaxation dynamics starts to appear only very close
to Tg which is consistent with the presence of a $\beta$-wing in the
experimental DSM data only below 440 K Peng  _et al._ (2020). We mention that
previous DFT molecular dynamics simulations suggested that the fast
crystallization in GST225 is characterized by concerted atomic motions favored
by the presence of flexible axial bonds Lee and Elliott (2017). Our results on
the crystal growth velocity suggest that these atomic motions would still
feature an activated behavior with an activation energy close to that of the
self-diffusion coefficient.

Turning now to the crystal growth velocity above 650 K, it is clear from Fig.
6a that the WF formula is unable to reproduce the data from the simulations.
The behaviour at high temperature is particularly sensitive to the value of
$\Delta\mu$ which goes to zero at Tm. This misfit could therefore be somehow
reduced by changing $\Delta H_{m}$ in the range 0.12-0.24 eV/atom and Tm in
the range 860-940 K due to the uncertainties in these figures discussed in
previous sections. The results for the crystal growth velocities obtained from
the WF formula and different values of $\Delta H_{m}$ and Tm are reported in
the Supplementary Figure 21. Still a sizable disagreement between the WF
formula and the crystal growth velocities extracted from the simulations is
present even for the best choice of $\Delta H_{m}$ and Tm in the ranges given
above.

There are actually different examples in metals and semiconductors with
diffusion-controlled crystallization kinetics in which the WF formula does not
quantitatively predict results from simulations or experiments Kerrache  _et
al._ (2008); Freitas and Reed (2020); Stolk  _et al._ (1993). This discrepancy
has been ascribed to changes in the mobility of the supercooled liquid in the
proximity of the crystal interface Kerrache  _et al._ (2008). To assess the
origin of this discrepancy in GST225, we have thus analyzed the local atomic
mobility as a function of the distance from the crystal-liquid interface. The
local 2D diffusion coefficient in obtained from the Einstein relation
$<x^{2}>+<y^{2}>$=4$D$t in NVT simulations at different average temperatures
where the WF formula fails. The $x$ and $y$ directions lye in the interface
plane. The local $D$ is computed in the slices at different distances from the
interface shown in the Supplementary Figure 12. The local $D$ reported in the
Supplementary Figure 22 is indeed lower at the interface than deep in the
liquid as also discussed in Ref. Kerrache  _et al._ (2008). However, the
reduction in the mobility closer to the interface is not sufficiently large to
justify the overestimation of $v_{g}$ given by the WF formula once the bulk
value of $D$ is used.

In the attempt to fit the $v_{g}$ data at high temperatures, we reconsidered
the form of $u_{kin}$ following the prescription by Jackson Jackson (2004) who
multiplied $\kappa^{+}$ and $\kappa^{-}$ by a factor $e^{-\Delta S/k_{B}T}$
where $\Delta S$ is the (positive) entropy difference between the liquid and
the crystal. This expression, which is very seldom used Herlach  _et al._
(2018); Buta  _et al._ (2007); Burke  _et al._ (1988), was justified by
considering that $e^{-\Delta S/k_{B}T}$ is the ratio of the number
configurations in the crystal and in the liquid and that the rate at which an
atom can be incorporated in the crystal depends on the rate at which an
ergodic sampling of configurations in the liquid would find a crystalline
configuration Ediger  _et al._ (2008).

At low T, $\Delta S$ is small and then the correction has a minor effect on
$v_{g}$, while $\Delta S\approx 2.1k_{B}$ at Tm which leads to a large
reduction of the crystal growth velocity (see Supplementary Figure 23).
$\Delta S$ as a function of temperature was obtained as $\Delta S=\Delta
H_{m}/T_{m}-\int^{T_{m}}_{T}\frac{\Delta C_{v}}{T}dT$, with $\Delta C_{v}$
given by the caloric curves in Fig. 3c. The modified WF formula
$v_{g}=6Dd/\lambda^{2}e^{-\Delta S/k_{B}T}(1-exp(-\Delta\mu/k_{B}T))$ turns
out to fit reasonably well the data at all temperatures as shown in Fig. 6a
with $\lambda$=1.65 Å which is still a reasonable number for the typical
jumping distance. The overall small misfit at high temperatures could now be
accounted for by the slightly lower mobility close to the surface that we
discussed above (see Supplementary Figure 22). As a final remark, we also
mention that Kelton and Greer Kelton and Greer (1986) proposed the different
expression for the attachment/detachment rates $\kappa^{+}=\nu exp(-\Delta
G^{*}/k_{B}T+\Delta\mu/2k_{B}T)$, $\kappa^{-}=\nu exp(-\Delta
G^{*}/k_{B}T-\Delta\mu/2k_{B}T)$ which leads to
$v_{g}=6Ddf/\lambda^{2}(1-exp(-\Delta\mu/k_{B}T))e^{\Delta\mu/2k_{B}T}$. The
fit of the data with this latter formula is, however, less satisfactorily, as
shown in the Supplementary Figure 24.

We also mention that the formula for $v_{g}$ given by the two-dimensional
surface nucleation growth model Jackson  _et al._ (1967) does not work neither
for reasonable values of the crystal/liquid interface energy in the range
0.05-1 J/m2 Kohary and Wright (2013) which enters as a parameter in the
formula (see for instance Eq. 5 in Ref. Ferreira Nascimento and Zanotto
(2010)), which is consistent with the mostly continuous growth behavior
observed for GST225.

Figure 6: Crystal growth velocity. a) Crystal growth velocities ($v_{g}$)
extracted from the motion of the crystal-liquid interface in the MD
simulations (red dots). The crystal growth is along the [001] direction of the
cubic phase. The dark blue and green lines correspond to WF and modified-WF
fits (see text). b) Potential energy as a function of time in simulations of
the homogeneous crystallization of the supercooled liquid phase at different
temperatures. The inset shows a snapshot of the formation of several crystal
nuclei (with different colors) at 500 K, liquid-like atoms are not shown. c)
Crystal growth velocities ($v_{g}$) for the homogeneous (red dots)
crystallization. The dark blue and green lines correspond to WF and modified-
WF fits (see text). d) Crystal growth velocity (red dots) of a model of the
amorphous phase overheated at different temperatures above Tg (see text). The
light blue line refers to the results of Ref. Orava  _et al._ (2012) inferred
from DSC data available below 650 K.

We have also repeated the simulation of the crystallization from the crystal-
liquid interface along the [111] direction of the cubic phase. The results at
different temperatures are compared with those of the (100) surface in the
Supplementary Figure 25. The crystal growth velocity at the (100) and (111)
surfaces are essentially the same, which is somehow surprising given that in
simulations of GeTe the crystal growth velocity was significantly lower for
the (111) than for the (100) surface Sosso  _et al._ (2015). This behavior in
GST225 might be ascribed to the high roughness of the growing surfaces (see
Supplementary Figure 18).

We have also studied the crystallization in the bulk (homogeneous
crystallization) of the supercooled liquid to directly compare our results
with the experimental DSC data in Ref. Orava  _et al._ (2012) that refer to
these conditions. To this aim, we generated a 7992-atom cubic model of the
liquid phase at 990 K as discussed previously. Then, we quenched the model
down to different temperatures in 80 ps MD simulations in the NVT ensemble and
at the experimental density of the amorphous phase (0.0309 atom/Å3). Finally,
we performed long simulations up to 20 ns to crystallize the supercooled
liquid and to extract the crystal growth velocity. The potential energy as a
function of time for simulations at different temperatures shown in Fig. 6b
reveals the onset of the crystallization with an incubation time that
increases with temperature. Overcritical nucleus/nuclei form on a time scale
of 0.2 to 3 ns in the range 500-600 K, and after about 7 ns at 650 K.
Nucleation was not observed at and above 700 K in simulations lasting over 20
ns. For temperatures where nucleation was not observed after a few tens of
nanoseconds, the crystal growth velocities were estimated by heating at the
target temperature a configuration with an overcritical nucleus generated at a
lower temperature. Below 600 K the number of overcritical nuclei increases as
temperature decreases (see inset of Fig. 6b). For temperatures above 650 K, we
observe a single crystallite which gives rise to a uniform single crystal
filling the simulation box. Crystalline atoms are assigned to a crystalline
nucleus when they fall within a 3.6 Å from the outermost atoms of the nucleus.
The critical nucleus increases in size with temperature and it contains about
40 atoms at 500 K and 50 atoms at 650 K.

The evolution of the crystalline nuclei was monitored to evaluate $v_{g}$
which for a spherical nucleus with radius $R$ is given by $v_{g}={dR(t)/dt}$
with
$R\left(t\right)=\left({3N(t)}/{4\pi\rho_{NN}^{cubic}}\right)^{\frac{1}{3}}$,
where $N$ is the number of atoms in the nucleus. This assumption is valid only
in the early stage of crystallization when the nuclei do not interact with
each other or with their periodic image. The evolution of $R(t)$ of the nuclei
at several temperatures is reported in the Supplementary Figure 26. The
resulting $v_{g}$, averaged over different nuclei, are shown in Fig. 6c.

Figure 7: Coarsening. Simulation of crystallites coarsening in a 27000-atom
supercell at 680 K. Several nuclei, shown by different colors in the left
panel, coalesce after 22 ns into only two crystallites shown in the right
panel.

Notice that the heterogeneous crystal growth velocity is higher than the
homogeneous one as it was the case for GeTe Sosso  _et al._ (2013, 2015). This
difference can partially be accounted for by a geometric factor. The
interlayer distance $d$ entering in the WF formula for the crystal growth
velocity at the crystal/liquid interface turns for a spherical nucleus into
the factor $4/3(vol_{site}3/(4\pi))^{1/3}$, where $vol_{site}$ is the volume
associated with an adsorption site on the crystalline nucleus Raoux  _et al._
(2010). If we take $vol_{site}=4\pi/3(\lambda/2)^{3}$, the WF formula reads
$v_{g}=4D/\lambda(1-exp(-\Delta\mu/k_{B}T))$ that we used for GeTe with
$\lambda$= 3.0 ÅSosso  _et al._ (2013). For GST225, we attempted to fit the
data below 650 K with the more general formula
$v_{g}=8(vol_{site}3/(4\pi))^{1/3}D/\lambda^{2}(1-exp(-\Delta\mu/k_{B}T))$,
where $vol_{site}^{1/3}$ is about half the lattice parameter of the cubic
cell, i.e 3 Å. As expected, the resulting curve with $\lambda$= 3.02 Å, shown
in Fig. 6c, overestimates the crystal growth velocity at high temperature. A
better fit is obtained by adding the entropic factor $e^{-\Delta S/k_{B}T}$,
similarly to the heterogeneous case, as shown in Fig.6c for $\lambda$= 2.05 Å.
The difference in the values of $\lambda$ obtained from the best fit for the
homogeneous and heterogeneous case is due to the fact that the geometric
factor discussed above is in both case an approximated form because the
surface is rough and the nucleus is not spherical. Moreover, the sticking is
not expected to be the same for an extended surface and for a small
nanoparticle.

Finally, we repeated the simulations of the homogeneous crystallization for
the amorphous phase overheated above Tg in a 7992-atom cubic cell. The
amorphous model, first equilibrated at 300 K, was heated at different target
temperatures in the range 500-800 K with a heating rate of 5 K/ps. Note that
the fraction of Ge atoms in tetrahedral configurations decreases very rapidly
with temperature above Tg and before crystal nucleation sets in, as shown by
the distribution of $q$ order parameter reported in the Supplementary Figure
27. Crystal nucleation was observed only below 650 K in the simulation time of
6 ns. The evolution in time of the radius of the nuclei at different
temperatures is reported in in the Supplementary Figure 28. The resulting
$v_{g}$ as a function of temperature for the overheated amorphous phase is
compared in Fig. 6d with experimental DSC data from Ref. Orava  _et al._
(2012). The crystal growth velocities for the overheated amorphous phase are
very similar to those obtained for the supercooled liquid at the same
temperature and they are very close to the experimental data from ultrafast
DSC Orava  _et al._ (2012). At 600 K, for instance, the crystal growth
velocity from NN-MD is 4.0 m/s in the homogeneous crystallization of the
supercooled liquid, to be compared with the experimental value of about 2.5
m/s Orava  _et al._ (2012). Other values obtained from DFT-MD simulations in
literature at the same temperature range from 5.5 m/s in small models, to
0.3-1 m/s for intermediate models and 0.5 m/s for the largest models Hegedüs
and Elliott (2008); Lee and Elliott (2011a); Kalikka  _et al._ (2016);
Ronneberger  _et al._ (2015, 2018); Lee and Elliott (2011b); Kalikka  _et al._
(2012, 2014).
Just to provide a last example of the capability of the NN potential, we show
in Fig. 7 the evolution in time of a 27.000-atom model at 680 K in which we
observed the coarsening of crystalline nuclei on the time scale of 20 ns.

##  III DISCUSSION

In summary, we have devised a machine learning neural network interatomic
potential for GST225 based on the DeePMD scheme which employs deep neural
networks to fit a large DFT database Wang  _et al._ (2018). The NN potential
is highly accurate as it reproduces the DFT results for a wide range of
structural, dynamical and thermodynamical properties of the crystalline,
amorphous and liquid phases. The NN potential has been exploited in large
scale (12000 atoms for over 100 ns) MD simulations of the supercooled liquid
and overheated amorphous phases which yielded the crystal growth velocities in
a wide temperature range of interest for the operation of the memory devices
in good agreement with experimental data from ultrafast DSC Orava  _et al._
(2012). In the temperature range 500-650 K, the crystal growth velocity
extracted from the simulations show an activated behavior controlled by the
self-diffusion coefficient. The analysis of the MSD showed that the
$\beta$-relaxation seems to be present only at temperatures very close to Tg.
A modified form Jackson (2004) of the WF formula turns out to be suitable to
describe the crystal growth velocity on a wider temperature range from Tg to
Tm.

We remark that the computational cost of the NN potential scales linearly with
the system size Lu  _et al._ (2021) which would allow simulating the entire
volume of the active material of ultrascaled memories (linear scale of 10 nm)
on the timescale of the device operation. This would allow addressing several
issues such as crystallites coarsening and evolution of grain boundaries which
are particularly relevant for multiscale programming exploited in the
neuromorphing computing, just to name a few others besides those already
mentioned in the introduction. Extension of the potential to Ge-rich GeSbTe
alloys would also allow uncovering the mechanism of phase separation into Ge
and less Ge-rich ternary alloys which is believed to be responsible for the
raise of the crystallization temperature upon Ge enrichment of interest for
applications in memories embedded in microcontrollers Cappelletti  _et al._
(2020).

##  IV METHODS

The NN potential was obtained by fitting DFT energies, forces and the stress
tensor of a database containing about 180000 supercell models (configurations)
of GST225 in the liquid, amorphous, cubic and hexagonal phases by using the
DeePMD-kit open-source package Wang  _et al._ (2018); Zhang  _et al._ (2018).

DFT calculations were performed by using the Perdew-Burke-Ernzherof (PBE)
exchange and correlation functional Perdew  _et al._ (1996) and Goedecker-
Teter-Hutter (GTH) norm conserving pseudopotentials with $s$ and $p$ valence
electrons Goedecker and Teter (1996); Krack (2005) as implemented in the CP2k
package VandeVondele  _et al._ (2005). Kohn-Sham orbitals were expanded in
Gaussian-type orbitals of a triple-zeta-valence plus polarization basis set,
while a basis set of plane waves up to a kinetic energy cutoff of 100 Ry is
used to represent the charge density as implemented in the Quickstep scheme
VandeVondele  _et al._ (2005). This DFT-PBE framework has been shown in
previous works to reproduce well the experimental structural and dynamical
properties of GST225 in the liquid, amorphous and crystalline phases Caravati
_et al._ (2007, 2009); Sosso  _et al._ (2011); Campi  _et al._ (2017);
Baratella  _et al._ (2022); Cobelli  _et al._ (2021).
The liquid and amorphous phases were modeled by a 108-atom cubic supercell.
The cubic metastable phase was modeled by a cubic 57-atom and by an
ortorhombic 98-atom supercells. The hexagonal phase was modeled by a 108-atom
supercell. To generate the database configurations, DFT-MD simulations were
performed in the Born-Oppheneimer approximation with a timestep of 2 fs and by
restricting the Brillouin Zone (BZ) integration to the supercell
$\Gamma$-point. Then, for a subset of the atomic configurations extracted from
the MD trajectories, energy and forces were computed by properly integrating
the BZ with a higher cutoff of 400 Ry and then were used as a training set for
the NN. Integration of the Brillouin Zone was performed over a 3x3x3 k-point
mesh for 108-atom cell or a 4x4x4 mesh for smaller cells.
In the construction of the NN potential, we employed two-body and three-body
embedding descriptors (see Ref. Lu  _et al._ (2021); Wang  _et al._ (2022)) by
including the third or the first coordination shell, respectively. To this
end, we set the distance cutoff to 7 Å for the two-body descriptor and to 3.8
Å for the three-body descriptor.

We set the maximum number of neighbors to 30, 30 and 40 for Ge, Sb and Te for
the two-body descriptor and to 10 for the three-body descriptors for all
species. The embedding network has 3 hidden layers with 20, 40 and 80 neurons
for the two-body descriptors, and 3 hidden layers with 3, 6 and 12 neurons for
the three-body descriptors. Finally, the network for the fitting of energy and
forces consists of 4 hidden layers with 120, 60, 30 and 15 neurons with atomic
reference for the chemical species of -102.297, -146.521 and -218.984 eV for
Ge, Sb and Te atoms, obtained from isolated atoms calculations with CP2k. In
the embedding and fitting network, we have used the hyperbolic tangent as an
activation function. We have also exploited the residual neurons as discussed
in Ref. He  _et al._ (2016). The hyperparameters which control the learning
process according to Ref. Wang  _et al._ (2018) are reported in the
Supplementary Table 4. The NN-MD simulations were performed with the LAMMPS
code Thompson  _et al._ (2022) exploiting GPU acceleration with a timestep of
2 fs and a Nosé-Hoover thermostat Nosé (1984); Hoover (1985). Finally, we
exploited the OvitoStukowski (2010) tool for the visualisation and the
generation of all atomic snapshots of this manuscript.

Data availability

The NN potential, the training DFT database and atomic trajectories of the
crystallization process are available at
https://doi.org/10.24435/materialscloud:a8-45

Acknowledgements

The project has received funding from European Union NextGenerationEU through
the Italian Ministry of University and Research under PNRR M4C2I1.4 ICSC
Centro Nazionale di Ricerca in High Performance Computing, Big Data and
Quantum Computing (Grant No. CN00000013)

Code availability

LAMMPS, and DeePMD are free and open source codes available at https://
lammps.sandia.gov and http://www.deepmd.org, respectively.

Competing interests

The authors declare no competing interests.

Authors contributions

O.A.E.K. and M.B. designed the research, O.A.E.K. and M.B. performed the
research and analyzed data with the help of L.B. and M.P. for the generation
of the NN potential. O.A.E.K. and M.B. wrote the paper and all the authors
edited the manuscript before submission.

## References

  * Wuttig and Yamada (2007) M. Wuttig and N. Yamada, Nat. Mater. 6, 824 (2007).
  * Noé  _et al._ (2017) P. Noé, C. Vallée, F. Hippert, F. Fillot, and J.-Y. Raty, Semicond. Sci. Technol. 33, 013002 (2017).
  * Fantini (2020) P. Fantini, J. Phys. D: Appl. Phys. 53, 283002 (2020).
  * Kuzum  _et al._ (2012) D. Kuzum, R. G. Jeyasingh, B. Lee, and H.-S. P. Wong, Nano Lett. 12, 2179 (2012).
  * Tuma  _et al._ (2016) T. Tuma, A. Pantazi, M. Le Gallo, A. Sebastian, and E. Eleftheriou, Nat. Nanotechnol. 11, 693 (2016).
  * Sebastian  _et al._ (2020) A. Sebastian, M. Le Gallo, R. Khaddam-Aljameh, and E. Eleftheriou, Nat. Nanotechnol. 15, 529 (2020).
  * Cappelletti  _et al._ (2020) P. Cappelletti, R. Annunziata, F. Arnaud, F. Disegni, A. Maurelli, and P. Zuliani, J. Phys. D 53, 193002 (2020).
  * Zuliani  _et al._ (2015) P. Zuliani, E. Palumbo, M. Borghi, G. Dalla Libera, and R. Annunziata, Solid State Electron. 111, 27 (2015).
  * Zhang  _et al._ (2019) W. Zhang, R. Mazzarello, M. Wuttig, and E. Ma, Nat. Rev. Mater. 4, 150 (2019).
  * Orava  _et al._ (2012) J. Orava, A. Greer, B. Gholipour, D. Hewak, and C. Smith, Nat. Mater. 11, 279 (2012).
  * Hegedüs and Elliott (2008) J. Hegedüs and S. Elliott, Nat. Mater. 7, 399 (2008).
  * Lee and Elliott (2011a) T. H. Lee and S. R. Elliott, Phys. Rev. Lett. 107, 145702 (2011a).
  * Kalikka  _et al._ (2016) J. Kalikka, J. Akola, and R. Jones, Phys. Rev. B 94, 134105 (2016).
  * Ronneberger  _et al._ (2015) I. Ronneberger, W. Zhang, H. Eshet, and R. Mazzarello, Adv. Funct. Mater. 25, 6407 (2015).
  * Ronneberger  _et al._ (2018) I. Ronneberger, W. Zhang, and R. Mazzarello, MRS Commun. 8, 1018 (2018).
  * Xu  _et al._ (2022) Y. Xu, Y. Zhou, X.-D. Wang, W. Zhang, E. Ma, V. L. Deringer, and R. Mazzarello, Adv. Mater. 34, 2109139 (2022).
  * Rao  _et al._ (2017) F. Rao, K. Ding, Y. Zhou, Y. Zheng, M. Xia, S. Lv, Z. Song, S. Feng, I. Ronneberger, R. Mazzarello, _et al._ , Science 358, 1423 (2017).
  * Lee and Elliott (2017) T. H. Lee and S. R. Elliott, Adv. Mater. 29, 1700814 (2017).
  * Behler and Parrinello (2007) J. Behler and M. Parrinello, Phys. Rev. Lett. 98, 146401 (2007).
  * Behler (2016) J. Behler, J Chem. Phys. 145, 170901 (2016).
  * Deringer  _et al._ (2019) V. L. Deringer, M. A. Caro, and G. Csányi, Adv. Mater. 31, 1902765 (2019).
  * Wang  _et al._ (2018) H. Wang, L. Zhang, J. Han, and W. E, Comput. Phys. Commun. 228, 178 (2018).
  * Sosso  _et al._ (2012a) G. C. Sosso, G. Miceli, S. Caravati, J. Behler, and M. Bernasconi, Phys. Rev. B 85, 174103 (2012a).
  * Sosso  _et al._ (2013) G. C. Sosso, G. Miceli, S. Caravati, F. Giberti, J. Behler, and M. Bernasconi, J. Phys. Chem. Lett. 4, 4241 (2013).
  * Lee  _et al._ (2020) D. Lee, K. Lee, D. Yoo, W. Jeong, and S. Han, Comput. Mater. Sci. 181, 109725 (2020).
  * Dragoni  _et al._ (2021) D. Dragoni, J. Behler, and M. Bernasconi, Nanoscale 13, 16146 (2021).
  * Shi  _et al._ (2021) M. Shi, J. Li, M. Tao, X. Zhang, and J. Liu, Mater. Sci. Semicond 136, 106146 (2021).
  * Gabardi  _et al._ (2015) S. Gabardi, S. Caravati, G. C. Sosso, J. Behler, and M. Bernasconi, Phys. Rev. B 92, 054201 (2015).
  * Raty  _et al._ (2015) J.-Y. Raty, W. Zhang, J. Luckas, C. Chen, R. Mazzarello, B. Bichara, and M. Wuttig, Nat. Commun. 6, 7467 (2015).
  * Mocanu  _et al._ (2018) F. C. Mocanu, K. Konstantinou, T. H. Lee, N. Bernstein, V. L. Deringer, G. Csányi, and S. R. Elliott, J. Phys. Chem. B 122, 8998 (2018).
  * Zhou  _et al._ (2023) Y. Zhou, W. Zhang, E. Ma, and V. L. Deringer, Nat. Electr. https://doi.org/10.1038/s41928-023-01030-x (2023).
  * Zhang  _et al._ (2018) L. Zhang, J. Han, H. Wang, R. Car, and W. E, Phys. Rev. Lett. 120, 143001 (2018).
  * Jackson (2004) K. A. Jackson,  _Kinetic Processes_ (Wiley-VCH, Weinheim, Germany, 2004).
  * Wilson (1900) H. W. Wilson, Lond. Edinb. Dublin Philos. Mag. J. Sci. 50, 238 (1900).
  * Frenkel (1932) J. Frenkel, Phisik. Zeit. Sowjetunion 1, 498 (1932).
  * Kooi and Wuttig (2020) B. J. Kooi and M. Wuttig, Adv. Mater. 32, 1908302 (2020).
  * Pries  _et al._ (2022) J. Pries, H. Weber, J. Benke-Jacob, I. Kaban, S. Wei, M. Wuttig, and P. Lucas, Adv. Func. Mater. 32, 2202714 (2022).
  * Laio and Parrinello (202) A. Laio and M. Parrinello, Proc. Natl. Acad. Sci. U.S.A. 99, 12562 (202).
  * Bonati and Parrinello (2018) L. Bonati and M. Parrinello, Phys. Rev. Lett. 121, 265701 (2018).
  * Niu  _et al._ (2020) H. Niu, L. Bonati, P. M. Piaggi, and M. Parrinello, Nat. Commun. 11, 2654 (2020).
  * Urata (2022) S. Urata, J. Phys. Chem. C 126, 21507 (2022).
  * Zhai and Wang (2023) B. Zhai and H. Wang, Comput. Mater. Sci. 216, 111843 (2023).
  * Liu  _et al._ (2023) J. Liu, R. Liu, Y. Cao, and M. Chen, Phys. Chem. Chem. Phys. 25, 983 (2023).
  * Huang  _et al._ (2023) X. Huang, K. Luo, Y. Shen, Y. Yue, and Q. An, Energy and AI 11, 100210 (2023).
  * Thompson  _et al._ (2022) A. P. Thompson, H. M. Aktulga, R. Berger, D. S. Bolintineanu, W. M. Brown, P. S. Crozier, P. J. in ’t Veld, A. Kohlmeyer, S. G. Moore, T. D. Nguyen, R. Shan, M. J. Stevens, J. Tranchida, C. Trott, and S. J. Plimpton, Comp. Phys. Comm. 271, 108171 (2022).
  * Nosé (1984) S. Nosé, J. Chem. Phys. 8, 511 (1984).
  * Hoover (1985) W. J. Hoover, Phys. Rev. A 31, 1695 (1985).
  * Njoroge  _et al._ (2002) W. K. Njoroge, H.-W. Wöltgens, and M. Wuttig, J. Vac. Sci. Technol. A 20, 230 (2002).
  * Schumacher  _et al._ (2016) M. Schumacher, H. Weber, P. Jóvári, Y. Tsuchiya, T. G. Youngs, I. Kaban, and R. Mazzarello, Sci. Rep. 6, 1 (2016).
  * Caravati  _et al._ (2009) S. Caravati, M. Bernasconi, T. Kühne, M. Krack, and M. Parrinello, J. Phys. Condens. Matter 21, 255501 (2009).
  * Rizzi  _et al._ (2014) M. Rizzi, N. Ciocchini, S. Caravati, M. Bernasconi, P. Fantini, and D. Ielmini, in  _2014 IEEE International Electron Devices Meeting_ (2014) pp. 29.6.1–29.6.4.
  * Flores-Ruiz and Micoulaut (2018) H. Flores-Ruiz and M. Micoulaut, J. Chem. Phys. 148, 034502 (2018).
  * Cohen and Grest (1979) M. H. Cohen and G. Grest, Phys. Rev. B 20, 1077 (1979).
  * Allen and Tildesley (2017) M. P. Allen and D. J. Tildesley,  _Computer simulation of liquids_ (Oxford University Press, Oxford (UK), 2017).
  * Sosso  _et al._ (2012b) G. C. Sosso, J. Behler, and M. Bernasconi, Phys. Status Solidi B 249, 1880 (2012b).
  * Grimme (2006) S. Grimme, J. Comput. Chem. 27, 1787 (2006).
  * Senkader and Wright (2004) S. Senkader and C. D. Wright, J. Appl. Phys. 95, 504 (2004).
  * Zhao  _et al._ (2022) J. Zhao, J. Hui, Z. Ye, T. Lai, M. Y. Efremov, H. Wang, and L. H. Allen, Adv. Mater. Inter. 9, 2200429 (2022).
  * Zunger  _et al._ (1990) A. Zunger, S.-H. Wei, L. G. Ferreira, and J. E. Bernard, Phys. Rev. Lett. 65, 353 (1990).
  * Pries  _et al._ (2019) J. Pries, S. Wei, M. Wuttig, and P. Lucas, Adv. Mater. 31, 1900784 (2019).
  * Angell (1995) C. A. Angell, Science 267, 1924 (1995).
  * Caravati  _et al._ (2007) S. Caravati, M. Bernasconi, T. Kühne, M. Krack, and M. Parrinello, Appl. Phys. Lett. 91, 171906 (2007).
  * Akola and Jones (2007) J. Akola and R. Jones, Phys. Rev. B 76, 235201 (2007).
  * Deringer  _et al._ (2014) V. L. Deringer, W. Zhang, M. Lumeij, S. Maintz, M. Wuttig, R. Mazzarello, and R. Dronskowski, Angew. Chem. Int. Ed. 53, 10817 (2014).
  * Errington and Debenedetti (2001) J. R. Errington and P. G. Debenedetti, Nature 409, 318 (2001).
  * Spreafico  _et al._ (2011) E. Spreafico, S. Caravati, and M. Bernasconi, Phys. Rev. B 84, 144205 (2011).
  * Guerboub  _et al._ (2023) M. Guerboub, S. D. W. Wendji, C. Massobrio, A. Bouzid, M. Boero, G. Ori, and E. Martin, J. Chem. Phys. 158, 084504 (2023).
  * Troullier and Martins (1991) N. Troullier and J. L. Martins, Phys. Rev. B 43, 8861 (1991).
  * Micoulaut and Flores-Ruiz (2021) M. Micoulaut and H. Flores-Ruiz, Phys. Rev. B 103, 134206 (2021).
  * Bichara  _et al._ (1992) C. Bichara, A. Pellegatti, and J.-P. Gaspard, Phys. Rev. B 47, 5002 (1992).
  * Campi  _et al._ (2017) D. Campi, L. Paulatto, G. Fugallo, F. Mauri, and M. Bernasconi, Phys. Rev. B 95, 024311 (2017).
  * Togo and Tanaka (2015) A. Togo and I. Tanaka, Scr. Mater. 108, 1 (2015).
  * Petrov  _et al._ (1968) I. Petrov, R. Imamov, and Z. Pinsker, Sov. Phys. Crystallogr. 13, 339 (1968).
  * Kooi and De Hosson (2002) B. Kooi and J. T. M. De Hosson, J. Appl. Phys. 92, 3584 (2002).
  * Matsunaga  _et al._ (2004) T. Matsunaga, N. Yamada, and Y. Kubota, Acta Crystallographica Section B: Structural Science 60, 685 (2004).
  * Wuttig  _et al._ (2018) M. Wuttig, V. L. Deringer, X. Gonze, C. Bichara, and J.-Y. Raty, Adv. Mater. 30, 1803777 (2018).
  * Steinhardt  _et al._ (1983) P. J. Steinhardt, D. R. Nelson, and M. Ronchetti, Phys. Rev. B 28, 784 (1983).
  * Kerrache  _et al._ (2008) A. Kerrache, J. Horbach, and K. Binder, EPL 81, 58001 (2008).
  * Sosso  _et al._ (2015) G. C. Sosso, M. Salvalaglio, J. Behler, M. Bernasconi, and M. Parrinello, J. Phys. Chem. C 119, 6428 (2015).
  * Thompson and Spaepen (1979) C. V. Thompson and F. Spaepen, Acta Metall. 27, 1855 (1979).
  * Cavagna (2009) A. Cavagna, Phys. Rep. 476, 51 (2009).
  * Peng  _et al._ (2020) S.-X. Peng, Y. Cheng, J. Pries, S. Wei, H.-B. Yu, and M. Wuttig, Sci. Advanc. 6, eaay6726 (2020).
  * Cheng  _et al._ (2022) Y. Cheng, Q. Yang, J. Wang, T. Dimitriadis, M. Schumacher, M. J. Zhang, H. Müller, N. Amini, F. Yang, A. Schoekel, J. Pries, R. Mazzarello, M. Wuttig, H.-B. Yu, and S. Wei, Nat. Commun. 13, 7352 (2022).
  * Sosso  _et al._ (2014) G. C. Sosso, J. Colombo, J. Behler, E. Del Gado, and M. Bernasconi, J. Phys. Chem. B 118, 13621 (2014).
  * Donati  _et al._ (2002) C. Donati, S. Franz, S. C. Glotzer, and G. Parisi, J. Non Cryst. Sol. 307, 215 (2002).
  * Widmer-Cooper and Harrowell (2007) A. Widmer-Cooper and P. Harrowell, J. Chem. Phys. 126, 154503 (2007).
  * Freitas and Reed (2020) R. Freitas and E. J. Reed, Nat. Commun. 11, 3260 (2020).
  * Stolk  _et al._ (1993) P. A. Stolk, A. Polman, and W. C. Sink, Phys. Rev. B 47, 5 (1993).
  * Herlach  _et al._ (2018) D. Herlach, D. Simons, and P.-Y. Pichon, Phil. Trans. R. Soc. A 376, 20170205 (2018).
  * Buta  _et al._ (2007) D. Buta, M. Asta, and J. J. Hoyt, J. Chem. Phys. 127, 074703 (2007).
  * Burke  _et al._ (1988) F. E. Burke, J. Q. Broughton, and G. H. Gilmer, J. Chem. Phys. 89, 1030 (1988).
  * Ediger  _et al._ (2008) M. D. Ediger, P. Harrowell, and L. Yu, J. Chem. Phys. 128, 034709 (2008).
  * Kelton and Greer (1986) K. F. Kelton and A. L. Greer, J. Non-Cryst. Solids 79, 295 (1986).
  * Jackson  _et al._ (1967) K. A. Jackson, D. R. Uhlmann, and J. D. Hunt, J. Cryst. Growth 1, 1 (1967).
  * Kohary and Wright (2013) K. Kohary and C. D. Wright, Phys. Status Solidi B 250, 944 (2013).
  * Ferreira Nascimento and Zanotto (2010) M. L. Ferreira Nascimento and E. D. Zanotto, J. Chem. Phys. 133, 174701 (2010).
  * Raoux  _et al._ (2010) S. Raoux, W. Welnic, and D. Ielmini, Chem. Rev. 110, 240 (2010).
  * Lee and Elliott (2011b) T. H. Lee and S. R. Elliott, Phys. Rev. B 84, 094124 (2011b).
  * Kalikka  _et al._ (2012) J. Kalikka, J. Akola, J. Larrucea, and R. O. Jones, Phys. Rev. B 86, 144113 (2012).
  * Kalikka  _et al._ (2014) J. Kalikka, J. Akola, and R. Jones, Phys. Rev. B 90, 184109 (2014).
  * Lu  _et al._ (2021) D. Lu, H. Wang, M. Chen, L. Lin, R. Car, W. E, W. Jia, and L. Zhang, Comput. Phys. Commun. 259, 107624 (2021).
  * Perdew  _et al._ (1996) J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 77, 3865 (1996).
  * Goedecker and Teter (1996) S. Goedecker and M. Teter, Phys. Rev. B 54, 1703 (1996).
  * Krack (2005) M. Krack, Theor. Chem. Acc. 114, 145 (2005).
  * VandeVondele  _et al._ (2005) J. VandeVondele, M. Krack, F. Mohamed, M. Parrinello, T. Chassaing, and J. Hutter, Comput. Phys. Commun. 167, 103 (2005).
  * Sosso  _et al._ (2011) G. C. Sosso, S. Caravati, R. Mazzarello, and M. Bernasconi, Phys. Rev. B 83, 134201 (2011).
  * Baratella  _et al._ (2022) D. Baratella, D. Dragoni, and M. Bernasconi, Phys. Status Solidi RRL 16, 2100470 (2022).
  * Cobelli  _et al._ (2021) M. Cobelli, D. Dragoni, S. Caravati, and M. Bernasconi, Phys. Rev. Mater. 5, 045004 (2021).
  * Wang  _et al._ (2022) X. Wang, Y. Wang, L. Zhang, F. Dai, and H. Wang, Nucl. Fusion 62, 126013 (2022).
  * He  _et al._ (2016) K. He, X. Zhang, S. Ren, and J. Sun, in  _2016 IEEE Conference on Computer Vision and Pattern Recognition (CVPR)_ (2016) p. 770.
  * Stukowski (2010) A. Stukowski, Modelling Simul. Mater. Sci. Eng. 18 (2010).