Unraveling the Crystallization Kinetics of the Ge 2 Sb 2 Te 5 Phase Change Compound with a Machine-Learned Interatomic Potential

The phase change compound Ge 2 Sb 2 Te 5 (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.


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 [1] 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) [1][2][3], and in neuromorphic and in-memory computing devices [4][5][6].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 Ge 2 Sb 2 Te 5 (GST225) [1].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 [1].Gerich GeSbTe alloys with crystallization temperatures much higher than that of GST225 have also been investigated for memories embedded in microcontrollers for automotive applications [7,8].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 [6,9].A key functional property for all these applications is the crystallization kinetics of the amorphous phase between the glass transition (T g ) and the melting (T m ) temperatures.This feature is, however, difficult to be investigated experimentally because of the very high nucleation rates and crystal growth veloc-ities (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 [10].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) [9,[11][12][13][14][15][16][17][18].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 [9,11].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 T g , where the thermodynamical driving force for crystallization is also high [9].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 [19][20][21][22].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 [23][24][25] and Sb [26,27].NN simulations of GeTe also allowed addressing the study of the aging of the amorphous phase [28] that leads to an increase of the electrical resistance with time (drift) which is detrimental for the operation of the memory devices [29].For GST225, a machine learning interatomic potential was recently developed with the Gaussian Approximations Potential (GAP) method [30] 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 [28,29].Improvements of this potential have been, however, very recently achieved [31].In this paper, we report on the development of an interatomic potential for GST225 within the NN scheme implemented in the DeePMD code [22,32].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 [33] of the phenomenological Wilson-Frenkel formula [34,35] 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 [36], or the possible existence of a strong-to-fragile transition in the supercooled liquid close to T g which is particularly relevant for the aging of the amorphous phase [37], just to name a few.

A. 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 [22,32].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.036atom/ Å3 .Details of configurations in the database at different conditions are reported in Tab.I 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 [38] 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 [39,40].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.
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/ Å [41][42][43][44].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 [45] and a Nosé-Hoover thermostat [46,47].
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 [48] which is very close to that of the liquid at 893 K (0.0307 atoms/ Å3 [49]).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.
Regarding the dynamical properties, we computed the self-diffusion coefficient D from NN-MD simulations at several temperatures above T m and below T m 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 [51] 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=6Dt 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 E a = 0.267 eV and D 0 =1.03x10 −3 cm 2 /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 [53] 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.3K and T 0 =330.6K (see Supplementary Figure 1).We have chosen the CG formula because it was used in the experimental work in Ref. [10] 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.
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 [54].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 crystalliquid 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 T m 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 T m 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 [49].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 T m 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 [55].We then repeated the same calculations in the temperature range 925-940 K by adding van der Waals interactions (vdW) according to Grimme [56] (D2) which prevents the coalescence of nanovoids that form by decreasing the density as it was also reported for GeTe [55].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   PBESol [52] and vdW-DF2 [49]) (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.different temperature we estimated T m =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 ∆H m =163 meV/atom which is close to the experimental values of 120 or 173 meV/atom [57,58].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 [59] 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 com-puted 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 ∆E m = 154 meV/atom.The value of ∆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 T m 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 T m 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 C v 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 T g =473 K reported in Ref. [60], where it was also proposed that crystallization in DSC occurs below T g .The decrease of C v with temperature in the supercooled liquid above T g is another feature typical of fragile liquids [61].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 T g .
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.
The structure of amorphous GST225 is similar to that emerged in previous DFT works [62,63].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 90 o is due to pyramidal and defective octahedral configurations, while the weak peak at about 170 o is due to axial bonds in defective octahedra.The shoulder at about 109 o is due to tetrahedra which are favored by homopolar Ge-Ge bonds [64] 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. [65].It is defined as , where the sum runs over the pairs of atoms bonded to a central atom j and forming a bonding angle θ 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 [66].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% [62,63].We remark, however, that a very recent paper [67] 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) [68] 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 extendedx-ray-fine-structure measurements [67].A similarly larger fraction of tetrahedra have been obtained in Ref. [69].Information on the medium range order is provided by the distribution of rings length which in GST225 is dominated by the four-membered rings [11,62].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 [11].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) [70] which highlights the presence of a short and a long axial bond in the defective octahedral configurations [49].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 [28,29].Concerning the dynamical properties, the NN potential reproduces well the DFT phonon density of states (DOS) of amorphous GST225 as shown in Fig. 5a.
The Crystalline Phases: At normal conditions, GST225 crystallizes in a hexagonal phase (space group P 3m1) [73][74][75].However, the amorphous phase crys- [71] (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 [72].
tallizes 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. [50]) 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. [76].Three different models of hexagonal GST225 have been proposed in literature differing in the distribution of Sb/Ge atoms in the cation sublattices [73][74][75].Here, we considered the Kooi stacking [74] 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 energyvolume points by the Birch-Murnaghan formula (see Supplementary Figures 7-8).The resulting parameters at equilibrium are compared with DFT data in Tab.III.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) [56] 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 [75].The NN poten-tial 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 [71] 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) [71].Hence, we calculated the phonon dispersion relations with the NN+D2 potential by employing the Phononpy code [72].The resulting phonon DOS are compared with previous DFT+D2 results [71] in Fig. 5c, while the phonon dispersion relations are compared with DFT+D2 results [71] 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.As a first application of the NN potential for GST225, we studied the kinetics of the crystallization process in the liquid phase supercooled below T m and in amorphous phase overheated above T g 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 dot 4 [77] 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 (t) = N (t)/2Aρ cubic N N , where N (t) is the number of atoms in the crystalline slab, ρ cubic N N is the theoretical equilibrium density of the cubic phase, A=6.15 × 6.15 nm 2 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) [34,35] v g = u kin (1 − exp(−∆µ/k B T )) where u kin is a kinetic prefactor and ∆µ 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 /λ 2 , where λ 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 [78].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 = λ and f =1 the kinetic prefactor has the form u kin = 6D λ that we used for instance in a previous simulation of the crystallization of GeTe from a crystal-liquid interface (heterogeneous crystallization) [79].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 (κ + − κ − ), where κ + /κ − are the rate of attachment/detachment of an atom to/from the crystalline surface.In turn κ + = νexp(−∆G * /k B T ) and κ − = νexp(−∆G * /k B T −∆µ/k B T ), where ∆G * is the activation energy for the attachment to the surface of an atom from the liquid with an attempt frequency ν.By assuming that κ + is equal to the rate of a jump in the diffusion process of a single atom in the liquid, it can be written as κ + = 6D/λ 2 .This approximation leads to v g = 6Ddf /λ 2 (1 − exp(−∆µ/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, λ 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 = 6Dd λ 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 ∆µ the expression given by Thompson and Spaepen ∆µ(T ) = ∆Hm(Tm−T )

2T
Tm+T ) [80].We set ∆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 ∆µ from the integration of the specific heat as ∆µ , where ∆C p is the difference in C p between the amorphous and the crystalline phases.We approximated ∆C p with ∆C v which was in turn computed from the caloric curve at fixed density as discussed in Sec. 2. The resulting ∆µ 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 λ 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 λ= 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 selfdiffusion 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 ef f that might embody a more short-ranged atomic motion.In glasses close to T g , the secondary β-relaxation is known to be the dominant source of atomic dynamics while the slower α relaxation controls the long range atomic diffusivity [81].At high temperature far from T g , the α and β relaxations actually coincide [81].The possibility that β-relaxation might enhance the crystallization kinetics of phase change materials close to T g has been put forward very recently [82,83].Indeed, the presence of a β-relaxation process in phase change materials has been recently identified from the so called β-wing in the temperature dependence of the loss modulus measured by dynamical mechanical spectroscopy (DSM) [82].Evidences of a link between the crystallization speed and the presence of β relaxation have been provided very recently for the eutectic alloy Ge 15 Sb 85 [83].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 β-relaxation dynamics due to the confinement [26].It is therefore of interest to investigate whether β-relaxation might also be of relevance for crystal growth in GST225.The β-relaxation can be detected in MD simulations by looking at the intermediate scattering function [81], as we did for instance for the GeTe phase change compound in our previous work [84].There, we also used four-point correlation functions [85] and isoconfigurational analysis [86] 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 β-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 T g is typical of a two steps relaxation dynamics with a faster β-relaxation and a slower α-relaxation that controls the long range atomic diffusivity after the plateau [81].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 βrelaxation dynamics starts to appear only very close to T g which is consistent with the presence of a β-wing in the experimental DSM data only below 440 K [82].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 [18].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 ∆µ which goes to zero at T m .This misfit could therefore be somehow reduced by changing ∆H m in the range 0.12-0.24eV/atom and T m 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 ∆H m and T m 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 ∆H m and T m 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 [78,87,88].This discrepancy has been ascribed to changes in the mobility of the supercooled liquid in the proximity of the crystal interface [78].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 >=4Dt 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. [78].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 [33] who multiplied κ + and κ − by a factor e −∆S/k B T where ∆S is the (positive) entropy difference between the liquid and the crystal.This expression, which is very seldom used [89][90][91], was justified by considering that e −∆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 [92].
At low T, ∆S is small and then the correction has a minor effect on v g , while ∆S ≈ 2.1k B at T m which leads to a large reduction of the crystal growth velocity (see Supplementary Figure 23).∆S as a function of temperature was obtained as ∆S = ∆H m /T m − Tm T ∆Cv T dT , with ∆C v given by the caloric curves in Fig. 3c.The modified WF formula v g = 6Dd/λ 2 e −∆S/k B T (1 − exp(−∆µ/k B T )) turns out to fit reasonably well the data at all temperatures as shown in Fig. 6a with λ=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 [93] proposed the different expression for the attachment/detachment rates κ ))e ∆µ/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 [94] does not work neither for reasonable values of the crystal/liquid interface energy in the range 0.05-1 J/m 2 [95] which enters as a parameter in the formula (see for instance Eq. 5 in Ref. [96]), which is consistent with the mostly continuous growth behavior observed for GST225.
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 dif-ferent 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 [79].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. [10] 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 (t) = 3N (t)/4πρ cubic N N 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.
Notice that the heterogeneous crystal growth velocity is higher than the homogeneous one as it was the case for GeTe [24,79].This difference can par- tially 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π)) 1/3 , where vol site is the volume associated with an adsorption site on the crystalline nucleus [97].If we take vol site = 4π/3(λ/2) 3 , the WF formula reads v g = 4D/λ(1−exp(−∆µ/k B T )) that we used for GeTe with λ= 3.0 Å [24].For GST225, we attempted to fit the data below 650 K with the more general formula v g = 8(vol site 3/(4π)) 1/3 D/λ 2 (1 − exp(−∆µ/k B T )), where vol 1/3 site is about half the lattice parameter of the cubic cell, i.e 3 Å.As expected, the resulting curve with λ= 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 −∆S/k B T , similarly to the heterogeneous case, as shown in Fig. 6c for λ= 2.05 Å.The difference in the values of λ 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 T g 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 22 ns at 680K rapidly with temperature above T g 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. [10].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 [10].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 [10].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 [11][12][13][14][15][98][99][100].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 [22].The NN potential is highly accurate as it reproduces the DFT results for a wide range of structural, dynami-cal 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 [10].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 β-relaxation seems to be present only at temperatures very close to T g .A modified form [33] of the WF formula turns out to be suitable to describe the crystal growth velocity on a wider temperature range from T g to T m .
We remark that the computational cost of the NN potential scales linearly with the system size [101] 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 Gerich GeSbTe alloys would also allow uncovering the mechanism of phase separation into Ge and less Gerich 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 [7].

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 [22,32].
DFT calculations were performed by using the Perdew-Burke-Ernzherof (PBE) exchange and correlation functional [102] and Goedecker-Teter-Hutter (GTH) norm conserving pseudopotentials with s and p valence electrons [103,104] as implemented in the CP2k package [105].Kohn-Sham orbitals were expanded in Gaussian-type orbitals of a triple-zetavalence 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 [105].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 [50,62,71,[106][107][108]. 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 Γ-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. [101,109]) 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.984eV 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. [110].The hyperparameters which control the learning process according to Ref. [22] are reported in the Supplementary Table 4.The NN-MD simulations were performed with the LAMMPS code [45] exploiting GPU acceleration with a timestep of 2 fs and a Nosé-Hoover thermostat [46,47].Finally, we exploited the Ovito [111] tool for the visualisation and the generation of all atomic snapshots of this manuscript.
Table I: First maximum and minimum of the pair correlation functions in the liquid phase of GST225 at 990 K from NN-MD and DFT-MD (in parenthesis) simulations.order parameter for the liquid and cubic phase of GST225 at 600 K.
The order parameter for atom i is defined by 11 m=−4 q 4m,i q * 4m,j where j runs over the N i neighbors up to a given cutoff (3.6 Å in our case) and q 4m,i = 1 Ni Ni j=1 Y 4m (r ij ) where Y 4m is the fourth order spherical harmonic with order m, r ij is the vector connecting the two atoms.We chose a threshold of Q dot 4 =0.87 to identify a crystalline atom.Figure 25: Crystal growth velocities (v g ) extracted from the motion of the crystal-liquid interface in the MD simulations along the [111] direction of the cubic phase (blue dots).The values of (v g ) for the (100) surface (see Fig. 6a in the article) are reported for the sake of comparison (red dots).The simulations for the (111) surface are performed with a 13689-atom supercell and a cross-section of A=2809 Å2 orthogonal to the growth direction.In the simulations at and above 675 K no crystal nucleation is observed.The data at these temperatures refer to an overcritical nucleus formed at a lower temperature which was then brought to the higher target temperature to monitor its growth.Table IV: Details of the hyperparameters used for the learning rate and the Loss Function (LF) employed for the training of the NN potential as defined in Ref. [14].

Figure 1 :
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 (∆E=|E DF T -E N N |) and b) forces (∆F=|F DF T -F N N |).c) The distribution of the absolute error as a function of DFT energy.

Figure 2 :
Figure2: 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 999atom 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 Table1.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.[50].d) The Arrhenius plot of the diffusion coefficient D as a function of temperature compared with previous DFT results with the different functionals PBE[51], PBESol[52] and vdW-DF2[49]) (see Supplementary Table2for 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.

Figure 3 :
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×6.15×11.08nm 3 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 T g .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.

Figure 4 :
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.

Figure 5 :
Figure 5: Phonons.a) Phonon density of states (DOS) of the amorphous phase obtained from Γ-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 Γ-point of a 270-atom supercell.Phonon frequencies in panel a) andb) 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.[71] (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[72].

Figure 6 :
Figure 6: Crystal growth velocity.a) Crystal growth velocities (v g ) extracted from the motion of the crystalliquid 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 T g (see text).The light blue line refers to the results of Ref.[10]  inferred from DSC data available below 650 K.

Figure 7 :
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.

Figure 1 :Figure 2 :
Figure 1: a) Cohen-Grest fit of the self-diffusion coefficient as a function of temperature (green curve) compared with the Arrhenius fit (blue curve, see Fig. 2d in the article).b) The self-diffusion coefficient as a function of temperature resolved per chemical species.

Figure 3 :
Figure 3: Enthalpy of the supercooled liquid as a function of temperature in NPT simulations.The lines are quadratic fit above and below T g .The heat capacity at constant pressure as a function of temperature obtained from the enthalpy fit is shown in the inset.

Figure 4 :
Figure 4: Structural properties of two NN models of amorphous GST225 at 300 K with 999 (blue curves) or 7992 (red curves) atoms.a) Total and partial pair correlation functions.b) Bond angle distribution functions 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.d) Distribution of the local order parameter q for tetrahedricity for four-coordinated Ge atoms.

Figure 5 :Figure 6 : 1 )Figure 7 :Figure 8 :FrequencyFigure 9 :Figure 10 :
Figure 5: Distribution of rings length in amorphous GST225 at 300 K in the NN 999-atom model and averaged over four DFT models, compared with previous DFT results from Refs. 7 and 8.The distribution is computed by using the algorithm in Ref. 9.

Figure 12 :Figure 13 :Figure 14 :Figure 15 :Figure 16 :Figure 17 :
Figure 12: Snapshot of the slices (bins) 10 Å wide at different distances from the interface with the crystalline slab over which we computed the local temperature in Figure 13-17 and the local diffusion coefficient in Figure 22.

Figure 18 :
Figure 18: A snapshot of the (100) growing surface at 650 K from two different views.The surface contour is computed with the code Ovito [12].Only crystalline atoms are shown.Atoms on the liquid-crystal interface are depicted by violet.

Figure 19 :Figure 20 : 1 )Figure 21 :Figure 22 :
Figure 19: The free energy difference ∆µ between the cubic crystal and the liquid as a function of temperature from the Thompson-Spaepen formula (TS) and from the integration of the heat capacity C p (see article).

Figure 23 :Figure 24 :
Figure 23: The difference in entropy between the liquid and the crystal ∆S as a function of temperature obtained as ∆S = ∆H m /T m − Tm T ∆Cv T dT , with ∆C v given by the caloric curves in Fig. 3c in the article.

Figure 26 :
Figure 26:  The evolution of the radius of the crystalline nuclei as a function of time in simulations at different temperatures in the supercooled liquid and their linear fit in the range highlighted by horizontal bars which yields the crystal growth velocity.In the simulations at and above 675 K no crystal nucleation is observed.The data at these temperatures refer to an overcritical nucleus formed at a lower temperature which was then brought to the higher target temperature to monitor its growth.

Figure 27 :
Figure27: Distribution of the local order parameter q for tetrahedricity (see article) for four-coordinated Ge atoms in the amorphous phase at different temperature below and above T g .The data are obtained by averaging over a time span always shorter than the incubation time for crystallization.

Figure 28 :
Figure 28: The evolution of the radius of the crystalline nuclei as a function of time in simulations at different temperatures in the overheated amorphous phase and their linear fit in the range highlighted by horizontal bars yielding the crystal growth velocity.
Learning rate initial value = 0.001 final value = 3.51 10 −8 Decay of the learning rate prefactor = 0.95 steps = 7275 Energy prefactor in the LF initial value = 1 fin al value = 1000 Forces prefactor in the LF initial value = 8000 final value = 800 Virial prefactor in the LF initial value = 0.2 final value = 100 Batch size (configurations) training = 16 test = 128 Training epochs 200

Table I :
Details of the database used for the training of the NN potential.

Table II :
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.

Table III :
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 (E 0 ), volume (V 0 ), bulk modulus (B) and derivative of B (B ′ ) at equilibrium are reported.

Table II :
Activation energy for the self-diffusion coefficient in liquid GST225 from this work (NN-MD simulations) and previous works in literature.NN+D2 refer to simulations with the semiempirical correction due to Grimme 1 (D2) to include van der Waals interactions.

Table III :
Position ( Å) of the first maximum and minimum of the pair correlation functions in the amorphous phase of GST225 at 300 K from DFT and NN simulations.DFT data are given in parenthesis.The Ge-Ge homopolar bond in the four DFT models generated for this benchmark is slightly longer (2.66±0.05Å)thanthatreported in Ref.5(2.62 Å) while it is shorter than that reported in Ref.6(2.84Å).