Uncovering the effects of interface-induced ordering of liquid on crystal growth using machine learning

ABSTRACT The process of crystallization is often understood in terms of the fundamental microstructural elements of the crystallite being formed, such as surface orientation or the presence

of defects. Considerably less is known about the role of the liquid structure on the kinetics of crystal growth. Here atomistic simulations and machine learning methods are employed together

to demonstrate that the liquid adjacent to solid-liquid interfaces presents significant structural ordering, which effectively reduces the mobility of atoms and slows down the

crystallization kinetics. Through detailed studies of silicon and copper we discover that the extent to which liquid mobility is affected by interface-induced ordering (IIO) varies greatly

with the degree of ordering and nature of the adjacent interface. Physical mechanisms behind the IIO anisotropy are explained and it is demonstrated that incorporation of this effect on a

physically-motivated crystal growth model enables the quantitative prediction of the growth rate temperature dependence. SIMILAR CONTENT BEING VIEWED BY OTHERS MICROSCOPIC MECHANISMS OF

PRESSURE-INDUCED AMORPHOUS-AMORPHOUS TRANSITIONS AND CRYSTALLISATION IN SILICON Article Open access 16 January 2024 AUTONOMOUSLY REVEALING HIDDEN LOCAL STRUCTURES IN SUPERCOOLED LIQUIDS

Article Open access 30 October 2020 FAST CRYSTAL GROWTH OF ICE VII OWING TO THE DECOUPLING OF TRANSLATIONAL AND ROTATIONAL ORDERING Article Open access 04 July 2023 INTRODUCTION

Crystallization from the melt (Fig. 1) is a pervasive process in industry, from metal casting for structural applications to the Czochralski process for semiconductor wafer growth for

electronics. It is important to control and understand the crystal growth process because it is at this stage that the material’s microstructure morphology is created, which in turn defines

the material’s properties. Consequently, a great deal of effort has been put into understanding the complex interplay between structure, thermodynamics, and kinetics that governs the process

of crystal growth1,2,3. This has led to a mechanism-based understanding of crystallization4,5,6 in terms of the microstructural elements of the crystallite being formed. For example, the

character of the solid surface in contact with the liquid is known to affect the growth rate, with atomically rough surfaces leading to faster growth rates than flat low-index surfaces and

their vicinals. Considerably less attention has been put in understanding the effects that the liquid adjacent to the solid–liquid interface has on the process of crystal growth. Atomic

events leading to crystal growth are thermally activated processes taking place in the free-energy landscape illustrated in Fig. 2a. The rate of crystallization is proportional to \(\exp

(-\beta \Delta {E}_{{\rm{a}}})\), while the melting rate is proportional to \(\exp [-\beta (\Delta {E}_{{\rm{a}}}+\Delta \mu )]\), where Δ_E_a is the activation energy for solidification,

Δ_μ_ is the difference in chemical potential between the liquid and solid phases, _β_−1 ≡ _k_B_T_, and _k_B is the Boltzmann constant. The balance of these two rates results in the following

equation for the overall growth rate: $$r(T)=k(T)\left\{1-\exp \left[-\beta \Delta \mu (T)\right]\right\},$$ (1) where \(k(T)\equiv {k}_{0}\exp (-\beta \Delta {E}_{{\rm{a}}})\) is known as

the kinetic factor. In this model, known as the Wilson–Frenkel7,8 (WF) model, the activation energy for solidification is taken as the energy barrier for diffusion in the liquid, Δ_E_a = 

Δ_E_d, because crystallizing atoms must undergo the same self-diffusion process that occurs in the associated liquid phase. It is often found that the WF method cannot quantitatively predict

results from simulations or experiments9,10. This notorious discrepancy, while largely unsolved, has been attributed to changes in mobility of the supercooled liquid in the vicinity of the

crystal interface that would cause Δ_E_a > Δ_E_d, but no physical mechanism has been demonstrated to explain the origin of this effect. Here we employ atomistic simulations and machine

learning (ML) together to show that the solid–liquid interface induces partial ordering of the nearby liquid during crystal growth. Our approach is successfully applied to two different

families of materials: semiconductors and metals. We find that the interface-induced ordering (IIO) affects the mobility of liquid atoms and thus slows down the crystal growth kinetics. The

physical mechanism behind the IIO is explained and we demonstrate that by accounting for this effect it is possible to derive predictive models for crystal growth. RESULTS CRYSTAL GROWTH

SIMULATIONS We performed molecular dynamics (MD) simulations of crystalline silicon growth from its melt employing a simulation geometry akin to laboratory experiments of crystal growth: a

crystalline seed is introduced in the liquid and its growth is monitored over the course of the simulation (see Fig. 1 and Supplementary Video 1). This setup allows the different

microstructural elements of the growing crystallite to interact naturally (see Supplementary Videos 2 and 3), as they would in a crystal growth experiment. For this geometry Δ_μ_ = 

Δ_G_−_κ__γ_/_ρ_s, where Δ_G_ is the difference in free energy between the liquid and the crystal, and the second term is due to the Gibbs–Thomson effect, with _ρ_s being the density of the

solid, _γ_ the interfacial free energy, and \(\kappa =2/{{\mathcal{R}}}_{{\rm{eff}}}\) is a geometrical factor where \({{\mathcal{R}}}_{{\rm{eff}}}\) is the effective crystallite radius. All

the above parameters of the WF model were computed (Fig. 2b and c) in order to compare the model predictions against simulation results (for calculation details see “Methods” section and

Supplementary Notes 5 and 6). The comparison between model and simulations is shown in Fig. 3a, where it is evident that the WF model does not predict the growth rate for temperatures it was

not fitted to. MACHINE LEARNING ENCODING OF CRYSTALLIZATION EVENTS Historically, simpler simulation geometries have been favored as a way to isolate certain microstructural elements, which

are then probed separately10,11,12. Our use of the geometry shown in Fig. 1 makes the simulation more physically relevant at the expense of greatly diminishing the amount of information that

can be inferred due to the lack of a crystal growth model that accounts for all microstructural elements present and their respective interactions. Moreover, it also becomes challenging to

decipher the atomic events at play due to the sheer complexity of the environment that atoms are embedded in. Here these obstacles are overcome by employing ML algorithms to systematically

encode and classify the structure surrounding liquid atoms during crystallization events. Our approach builds on recently proposed ML strategies for the construction of a structural quantity

(namely softness \({\mathbb{S}}\)) that captures the propensity for atomic rearrangements to occur in disordered atomic environments, such as in glasses13,14 and inside grain boundaries15.

The structural characterization of local atomic environments is realized by assigning to each atom _i_ a local-structure fingerprint \({{\mathbf{x}}}_{i}\) constructed from a set of 21

radial structure functions13,16 \({\mathcal{G}}(r)\), as illustrated in Fig. 4a. Furthermore, atoms are labeled into three possible categories according to their first-neighbor’s

arrangement: liquid and crystal atoms have arrangement patterns statistically identical to the bulk liquid and bulk crystal, respectively. Meanwhile, crystallizing atoms have arrangement

patterns intermediary between the other two labels (see “Methods” section and Supplementary Note 2 for more details). It is possible to observe how these three groups of atoms are spread in

the \({{\mathbb{R}}}^{21}\)-space of local-structure fingerprints \({{\mathbf{x}}}_{i}\) with the help of an algorithm known as principal component analysis (PCA). With this method, a

dimensionality reduction transformation is performed to create a \({{\mathbb{R}}}^{2}\) representation of the \({{\mathbb{R}}}^{21}\) data, as shown in Fig. 4b. Superimposed in this figure

is also the trajectory of an atom that undergoes crystallization over the course of the simulation. Atoms assume varied local-structure fingerprints \({{\mathbf{x}}}_{i}\) depending on both

the surrounding liquid structure and the nearby interface morphology. In order to quantify these variations in microstructure we proceed as follows. First, an ML algorithm known as support

vector machine17,18,19 is employed to find the hyperplane that optimally separates the crystallizing atoms from the liquid atoms in the \({{\mathbb{R}}}^{21}\)-space of

\({{\mathbf{x}}}_{i}\). Then, the distance of each atom _i_ from the hyperplane (\({{\mathbb{S}}}_{i}\), known as softness13,14,15) is measured: atoms with \({{\mathbb{S}}}_{i}> 0\) lie

on the crystallizing side of the hyperplane, while \({{\mathbb{S}}}_{i}<0\) atoms lie on the liquid side. This approach is found to correctly classify liquid and crystallizing atoms with

an accuracy of 96%. It is important to realize that \({\mathbb{S}}\) is not an order parameter because it was not designed to track the change from the liquid to the solid phase. Instead,

\({\mathbb{S}}\) measures the propensity of an atom in the liquid phase to undergo the process of crystallization. Shown in Fig. 4c is a simulation snapshot with atoms colored according to

their softness value (see also Supplementary Videos 4 and 5). In this figure \({\mathbb{S}}\) is seen to capture the structural signs of dynamical heterogeneity in the supercooled liquid far

from the crystal, with clear indications of strong spatial correlations. These fluctuating heterogeneities have recently been shown to be preferential sites for crystal nucleation20,21.

Thus, Fig. 4c establishes that \({\mathbb{S}}\) is indeed capable of capturing subtle signs of structural ordering in liquids. LOCAL-STRUCTURE DEPENDENT (LSD) CRYSTAL GROWTH MODEL It is

possible now to decompose the crystal growth rate of Fig. 3a as a function of the local structure using \({\mathbb{S}}\). For example, in Fig. 3b it is shown that the total growth rate at

_T_ ≈ 1233 K varies by almost four orders of magnitude as the local structure changes. For this reason, we propose to address the limitations of the WF model by taking into account the local

structure surrounding the crystallizing atoms through an explicit dependence on \({\mathbb{S}}\): $$r(T,{\mathbb{S}})=k(T,{\mathbb{S}})\left\{1-\exp \left[-\beta \Delta \mu

(T,{\mathbb{S}})\right]\right\},$$ (2) with \(\Delta \mu (T,{\mathbb{S}})=\Delta G(T)-\kappa \gamma ({\mathbb{S}})/{\rho }_{{\rm{s}}}\). Indeed, accounting for the information about the

local structure contained in plots such as Fig. 3b results in a crystal growth model with predictive capabilities, as shown in Fig. 3a (see Supplementary Notes 3 and 4 for details on the

model calculation). Notice how the LSD model is capable of predicting the growth rate for a wide range of temperatures (i.e. _T_ < 1388 K) not included in the model parametrization. In

particular, the experimentally measured growth rate and its slope show much better agreement with our LSD model, Eq. (2), than with the WF model, Eq. (1). In the Supplementary Note 3 we show

that the variables introduced by the dependence on \({\mathbb{S}}\) are not independent parameters. Thus, the improved reproduction and prediction of simulation results cannot be attributed

simply to Eq. (2) exhibiting higher capacity or flexibility in modeling complex relationships when compared to Eq. (1). We now turn to investigate the ramifications of the LSD model

\(r(T,{\mathbb{S}})\) and uncover the source of its predictive capabilities. The kinetic factor \(k(T,{\mathbb{S}})\), shown in Fig. 5a, is observed to be a strong function of the local

structure, varying by as much as three orders of magnitude with \({\mathbb{S}}\). For each value of \({\mathbb{S}}\) the kinetic factor shows an Arrhenius-like temperature dependence

$$k(T,{\mathbb{S}})={k}_{0}({\mathbb{S}})\exp [-\beta \Delta {E}_{{\rm{a}}}({\mathbb{S}})].$$ This striking outcome suggests that each value of \({\mathbb{S}}\) corresponds to a thermally

activated and independent crystallization channel with well-characterized energy scale. Such a picture is reminiscent of our traditional understanding of crystallization in terms of the

solid–liquid interface morphology, with different values of \({\mathbb{S}}\) encoding the influence of different microstructural elements. But here \({\mathbb{S}}\) encodes more than just

the crystal local microstructure: it also encodes the variation in the structure of the liquid. The variation of liquid properties with local structure is reflected in the dependence of the

activation energy barrier of these crystallization channels with \({\mathbb{S}}\) shown in Fig. 5b, where it can be seen that \(\Delta {E}_{{\rm{a}}}({\mathbb{S}})\) varies over 1 eV with

\({\mathbb{S}}\). Additionally, the activation energy decreases monotonically with \({\mathbb{S}}\) and seems to approach the energy barrier for diffusion Δ_E_d (Fig. 2b) asymptotically.

Hence, the mobility of liquid atoms close to the solid–liquid interface seems to vary greatly, from a negligible change (Δ_E_a ≈ Δ_E_d) compared to bulk liquid to a dramatic reduction in

mobility due to the increase in Δ_E_a. This change in the liquid structure due to the presence of the solid–liquid interface is known as IIO, Supplementary Fig. 5b shows that the structural

change does indeed lead to local ordering. Varying \({\mathbb{S}}\) also has a pronounced effect on the Arrhenius prefactor \({k}_{0}({\mathbb{S}})\), as indicated in Fig. 5c, which

decreases by three orders of magnitude with \({\mathbb{S}}\). Because \({\mathrm{ln}}\,[{k}_{0}({\mathbb{S}})]\) can be interpreted as the product of an entropic contribution15 to the free

energy barrier and a term involving the population of crystallizing atoms with softness \({\mathbb{S}}\), the observed decrease in the prefactor indicates that there are less rearrangement

pathways leading liquid atoms to the activated state (i.e. to crystallization) as \({\mathbb{S}}\) increases. Hence, Fig. 5b and c together indicate that from all observed local-structure

arrangements surrounding crystallizing atoms, only very few lead to low-energy barriers. Additionally, Fig. 3b indicates that these few channels with low-energy barriers are the ones

contributing the most to the overall growth rate. Next, we examine how the free energy of the solid–liquid interface to which atoms attach varies with \({\mathbb{S}}\), which should give us

a glimpse of the microstructure at the crystallite surface. Figure 6a shows that \(\gamma ({\mathbb{S}})\) decreases monotonically with \({\mathbb{S}}\), starting at large

values—corresponding to high-index interfaces—and reaching interfacial free energy values characteristic of low-index interfaces in silicon. This finding implies that the decrease in

Arrhenius prefactor \({k}_{0}({\mathbb{S}})\) with softness (Fig. 5c) leads to fewer rearrangement pathways because crystallization events with large positive values of \({\mathbb{S}}\)

happen at low-index surfaces and their vicinals, which naturally offer less crystallization sites than high-index interfaces. Despite the scarcity of crystallization sites offered by

low-index interfaces and their vicinals, Fig. 3b shows that they contribute the most to the overall growth rate, with 70% of all atoms attaching to interfaces with \({\mathbb{S}}\ge 0.75\).

This observation is confirmed by direct measurement of the distribution of crystal surfaces to which crystallizing atoms attach: Fig. 6b reveals strong preferential attachment to a wide

variety of (111) vicinals. The high-intensity spot around (435) corresponds to step–step separation distances from 15 to 24 Å (Fig. 6c), indicating that the majority of the crystallization

events take place on vicinal surfaces with well-separated steps, which is exactly what is expected for silicon22. Figure 6b also shows a smaller amount of events occurring at high-index

interfaces, further validating the above observations. Notice in Supplementary Video 1 that the crystallite also exhibits signs of rough interfaces, thus it is possible that a small fraction

of the identified high-index are actually rough. APPLICABILITY TO A DIFFERENT FAMILY OF MATERIALS In order to verify that our approach in creating LSD predictive models of crystal growth is

not particular to silicon (or semiconductors) we apply it in the development of a crystal growth model for an elemental metal, namely copper (see “Methods” section and Supplementary Note 7

for simulation details). The resulting model is shown in Fig. 7a, where it can be seen that the LSD model of copper also correctly predicts the growth rate at temperatures at which it was

not parametrized on (i.e. it is a predictive model), while the WF model is not capable of reproducing simulation results at temperatures to which it was not fitted, similarly to what was

observed for silicon in Fig. 3a. Moreover, analysis of the parameters of the LSD model of copper (Fig. 7b) shows that all parameters present the same trend with \({\mathbb{S}}\) as observed

for silicon, including the Arrhenius behavior for the kinetic factor (shown in Supplementary Fig. 19c). The major difference observed between the LSD models of silicon and copper is that

\({k}_{0}({\mathbb{S}})\) and \(\Delta {E}_{{\rm{a}}}({\mathbb{S}})\) assume much larger values for copper. We attribute this to the predominance of rough interfaces in metallic systems. In

contrast to semiconductors, interfaces in metallic systems typically do not advance by the lateral motion of steps. Instead, metal interfaces often advance by atomic attachment directly on

top of them, leading to growth normal to the interface itself. This growth mechanism is reflected in Fig. 7a (inset), where it is seen that atomic attachments occur directly on (001) and

(111) interfaces—leading to normal growth—instead of vicinals of (111) as observed for silicon in Fig. 6b (compare also Fig. 1 to Supplementary Fig. 14 and Supplementary Video  1 to

Supplementary Video 7). Normal growth leads to the formation of atomically rough interfaces that offer a much larger amount of atomic disorder than well-structured high-index interfaces.

Hence, rough interfaces present a larger availability of sites for liquid atoms to attach, leading to much higher values for \({k}_{0}({\mathbb{S}})\) due to the numerous atomic pathways

leading to crystallization. The predominance of rough interfaces also explains the larger values of \(\Delta {E}_{{\rm{a}}}({\mathbb{S}})\) observed for copper, but this explanation will be

postponed until the “Discussion” section, where the connection will be discussed in the light of the effects of IIO on crystal growth. DISCUSSION Solid–liquid interfaces in equilibrium are

known to affect the structure of the nearby liquid by imparting some amount of order on it23,24,25,26,27,28,29,30,31. Here, we have established that the IIO of the liquid also occurs during

the process of crystal growth—a dynamic situation in which the solid–liquid interface is not in equilibrium. The observed IIO seems to decrease the mobility of liquid atoms through changes

in the activation barrier for crystallization Δ_E_a (Fig. 5b), effectively slowing down the crystallization kinetics. Comparison of Figs. 5b and 6a reveals that the IIO of the liquid is

anisotropic, i.e. it depends on the interface orientation and microscopic levels of roughness. The trend (illustrated in Fig. 8) is such that low-index surfaces and their vicinals

(corresponding to large positive \({\mathbb{S}}\) in Fig. 6a) cause weak ordering, resulting in smaller activation energies (i.e. \(\Delta {E}_{{\rm{a}}}({\mathbb{S}})\) close to the energy

barrier for atomic diffusion in the liquid bulk Δ_E_d). Meanwhile, high-index interfaces (negative \({\mathbb{S}}\) in Fig. 6a) cause strong ordering of the liquid, which becomes rigid and

results in activation energies much larger than Δ_E_d. However, even in the case of strong IIO the activation energy (1.75 eV) is still much smaller than the ≈4.6 eV32 barrier for

vacancy-mediated self-diffusion in crystalline silicon. This indicates that the structural order of the liquid affected by IIO is nowhere as substantial as crystalline order. The physical

cause of the IIO anisotropy is that the interaction between the crystal surface and the liquid is mediated by the amount of dangling bonds on the crystal surface. Thus, strong liquid

ordering (and slower mobility) is observed at high-index interfaces because these interfaces interact more strongly with the liquid, since they present more dangling bonds when compared to

low-index interfaces and its vicinals. This mechanism is illustrated in Fig. 8a and b, while its effect on the free-energy landscape of the system is illustrated schematically in Fig. 8c.

Notice that this mechanism also explains why copper has much larger values of \(\Delta {E}_{{\rm{a}}}({\mathbb{S}})\) (Fig. 7b) while having similar energy barrier for diffusion in the

liquid Δ_E_d: rough interfaces are predominant in copper and these interfaces have much stronger interactions with the liquid when compared to low-index and vicinal surfaces, which are

predominant in silicon. Dynamical heterogeneities present in the liquid (Fig. 4c) also affect the coordination of atoms20,33,34. For this reason, it is reasonable to expect that they

contribute to the  ≈1 eV dispersion in \(\Delta {E}_{{\rm{a}}}({\mathbb{S}})\) observed in Fig. 5b. Nonetheless, there is no evident reason to believe that their effect is anisotropic since

dynamical heterogeneities have origin in random thermal fluctuations. In conclusion, we have discovered that the IIO of liquids strongly affects the process of crystal growth in metals and

semiconductors. It is found that the modified structure of the liquid nearby solid–liquid interfaces reduces the mobility of liquid atoms, an effect shown to be essential in order to build a

predictive model of the growth rate temperature dependence. Indeed, the construction of such predictive model was only possible by identifying and incorporating in the model the family of

all thermally activated events— each with its own energy scale—leading liquid atoms to the crystal phase. Our work elevates the liquid structure to the same level of importance as the

crystal surface morphology in understanding crystallization, a knowledge that can enable material advances through the incorporation of liquid-structure engineering as a novel pathway for

synthesis. Our results were only made possible by employing atomistic simulations and ML together. The strength of this combined approach is that one can perform complex simulations and yet

glean physical insight from notoriously haphazard atomic environments. This innovative application of ML in materials science blends conventional scientific methods with data science tools

to produce physically consistent predictive models and novel conceptual knowledge. METHODS SILICON CRYSTAL GROWTH SIMULATIONS The MD simulations were performed using the Large-Scale

Atomic/Molecular Massively Parallel Simulator (LAMMPS35) software, with the interactions between silicon atoms described by the Stillinger–Weber36 interatomic potential. The timestep was

selected as ~1/56th of the period of the highest-frequency phonon mode of this system, or Δ_t_ = 1 fs. The crystal growth simulations contained 500,000 atoms and were initialized with a

spherical crystalline seed of ~3000 atoms in the diamond cubic structure. The lattice parameter for the atoms in the crystal seed was chosen taking into account dilation due to thermal

expansion, then the remainder of the simulation cell was filled with randomly distributed atoms at the equilibrium liquid density for that temperature at zero pressure. The system was

equilibrated by first relaxing the liquid atoms using a Conjugate Gradient37 algorithm for 200 steps. Next, the liquid was equilibrated at finite temperature using the

Bussi–Donadio–Parrinello38 (BDP) thermostat for 3 ps with a damping parameter of 0.1 ps. Finally, liquid atoms were equilibrated for 2 ps at zero pressure and finite temperature using the

same thermostat just described and a chain Nosé–Hoover barostat39,40,41,42,43 with damping parameter of 1 ps and a chain length of three, allowing only for isotropic dilation/contraction of

the system. During the entirety of this equilibration process the crystalline seed atoms were kept frozen at their equilibrium crystal structure with fixed lattice parameters (i.e. they did

not dilate/contract with the liquid atoms). After equilibration the BDP thermostat and chain Nosé–Hoover barostat were applied to the entire system, both with damping parameter of 1.0 ps, to

maintain the system at finite temperature and zero pressure for a total of 3 ns during which snapshots were recorded every 1 ps. The crystal growth process can be seen in the Supplementary

Video 1. Snapshots saved from the MD simulations were subsequently relaxed using 20 steps of the Steepest-Descent37 algorithm. The crystal growth simulations were performed at temperatures

ranging from 1125 to 1500 K in intervals of 25 K. The damping parameter for the thermostat was selected conservatively such that the liquid diffusivity was not affected by the presence of

the thermostat, i.e. it had the same value within the statistical uncertainty as the diffusivity computed without a thermostat. Thus, the thermostating of the crystal growth simulation was

performed gently as to not affect the kinetics of the system. See Supplementary Note 6 for more details on the diffusivity calculations. PHASE IDENTIFICATION In order to identify to which

phase (liquid or crystal) each particle belongs, we used the order parameter introduced by Rein ten Wolde et al.44. The complete description and analysis of the construction of this order

parameter can be found in the Supplementary Note 2. Ultimately, this method provides us with a parameter _α__i_(_t_) for each atom _i_ at the time _t_ of each MD snapshot. The physical

interpretation of this parameter is that _α__i_ is the fraction of bonds that atom _i_ makes that resemble bonds in a perfect crystal structure. As shown in the Supplementary Note 2, the

parameter _α__i_ correctly identifies atoms in the perfect crystal or bulk liquid with accuracy of 100% within the statistical uncertainty. It is important to notice that although _α__i_ can

discern between liquid and crystal atoms, it does not differentiate between crystalline structures. We confirm that the silicon atoms are indeed crystallizing in the diamond cubic structure

by performing the polyhedral template-matching45,46 analysis. ENCODING ATOMIC DYNAMICS (ML LABELING) The dynamics of each atom was encoded using the time evolution of the _α__i_(_t_) order

parameter. A representative plot47 of _α__i_(_t_) is shown in Supplementary Fig. 4b and c. Notice that due to thermal fluctuations the instantaneous value of _α__i_(_t_) for atoms in the

liquid and crystal phases might differ from their perfect values of 0.0 and 1.0, respectively, even after the short Steepest-Descent relaxation. Hence, we perform a moving-window average of

_α__i_(_t_) with window length of 20 ps and use the window-averaged \({\bar{\alpha }}_{i}(t)\) to label the atomic dynamics as illustrated in Supplementary Fig. 4a. Atoms with \({\bar{\alpha

}}_{i}(t)=0\) for _t_ ∈ [_t_0−_τ__ℓ_, _t_0 + _τ__ℓ_] receive label _y__i_ = −1 at time _t_0. These are atoms deep in the liquid phase that will not be transitioning to the crystal state in

the near future, neither have tried to transition in the near past. From the analysis of curves such as in Supplementary Fig. 4c we choose _τ__ℓ_ = 15 ps as a reasonable value. Next we

identify atoms that have just started to move out of the bulk liquid (i.e. crystallizing atoms) as those within a 20 ps window from the point where \({\bar{\alpha }}_{i}=0.25\), i.e. _y__i_ 

= 1 for _t_ ∈ [_t_0−_τ_, _t_0 + _τ_] where \({\bar{\alpha }}_{i}({t}_{0})=0.25\) and _τ_ = 10 ps. See Supplementary Note 2 for more details on the labeling process. LOCAL-STRUCTURE

FINGERPRINT (ML FEATURES) The local structure surrounding each atom was characterized using a set of radial structure functions13,16: $${{\mathcal{G}}}_{i}(r)=\mathop{\sum

}\limits_{j=1}^{n(i)}\exp \left[-{({r}_{ij}-r)}^{2}/2{\sigma }^{2}\right],$$ where _i_ is the atom whose local structure is being described, _n_(_i_) is the number neighbors of _i_ within a

cutoff radius _r_cut, _r__i__j_ is the distance between atom _i_ and one of its neighbors _j_, _r_ and _σ_ are two parameters that define the radial structure function. These smoothly

varying functions of _r_ count the number of neighbors of _i_ at a distance _r_. In this interpretation, parameter _r_ represents the radial distance from _i_ at which we are counting the

number of neighbors, while _σ_ adjusts how smoothly the function varies as atoms move in and out of the distance _r_ vicinity. We have used a grid-search48 algorithm to perform the

hyperparameters tuning (see Supplementary Note 1 for more details), resulting in _σ_ = 0.5 Å, _r_cut = 10.8 Å, and _r__n_ = (2.0 + 0.4_n_) Å, with _n_ = 0, 1, 2, …, 20. With this set of 21

radial structure functions—one for each value of _r_—the local-structure fingerprint of each atom _i_ was built as a vector:

$${{\mathbf{x}}}_{i}=\left[{{\mathcal{G}}}_{i}({r}_{1}),{{\mathcal{G}}}_{i}({r}_{2}),\ldots ,{{\mathcal{G}}}_{i}({r}_{21})\right].$$ SOFTNESS CALCULATION The data was assembled by pairing

the dynamic labels _y__i_ with their corresponding structural fingerprint \({{\mathbf{x}}}_{i}\). Then, 10,000 \(({y}_{i},{{\mathbf{x}}}_{i})\) pairs were randomly selected and equally

divided between the _y_ = −1 and _y_ = 1 classes to train a support vector machine17,18,19 (SVM) classifier. The SVM algorithm finds the hyperplane of the form \({\rm{w}}\cdot

{\mathbf{x}}-b=0\) that optimally separates the two classes, where W and _b_ are the parameters that define this hyperplane. Before training the SVM classifier the elements of the

fingerprints were standardized48 to have zero mean and standard deviation of one. The optimal hyperplane found, denoted by the parameters W* and _b_*, correctly separates the two classes

with an accuracy of 96%. See Supplementary Note 1 for more details about how these parameters are found and an in-depth analysis of the quality of the classifier computed. All results shown

here were obtained using data from the crystal growth simulation at _T_ = 1500 K to train the SVM classifier. However, the results can be reproduced within the statistical uncertainty when

training at any other temperature, as shown in the Supplementary Note 3. Once the SVM classifier has been trained it was applied to the entire data set, composed of 27.5 million data points

(excluding the data used for training, hyperparameter tuning, and cross validation). The value of softness14 for each data point (or atom) is the signed distance from the hyperplane, or

\({{\mathbb{S}}}_{i}={{\rm{w}}}^{* }\cdot {{\mathbf{x}}}_{i}-{b}^{* }\) for each atom _i_. PARAMETER ESTIMATION AND TEMPERATURE EXTRAPOLATION In order to test how predictive the LSD and WF

models are we performed the model parameterization of both models using only the data collected for low undercooling (i.e. _T_ ≥ 1388 K) and observed how well the model predicts the

temperature dependence for higher undercoolings (i.e. temperatures as low as 1128 K). As shown in Fig. 3a the LSD model is capable of predicting the growth rate at temperatures it was not

parametrized on, while the WF model only reproduces the simulation results at temperatures it was fitted to. We attribute this to the fact that the LSD model accounts for the Arrhenius

family of thermally activated atomic events leading to crystal growth, as labeled by \({\mathbb{S}}\). This is fundamental physical information that is not incorporated in the WF model. For

Figs. 5 and 6a only, the \({\mathbb{S}}\) dependence of the parameters of the LSD model (i.e. Δ_E_a, _k_0, and _γ_) was measured after reparametrizing this model using the data from all

simulations with _T_ ≥ 1206 K. The reparametrization was necessary only in order to reduce the statistical uncertainty of the parameters measured. This range of temperatures was chosen

because below _T_ = 1206 K the kinetic factor \(k(T,{\mathbb{S}})\) showed signs of non-Arrhenius behavior for some values of \({\mathbb{S}}\) due to the approaching glass transition

temperature. Notice that below _T_ = 1206 K the bulk liquid diffusivity also shows signs of departure from the Arrhenius behavior. Thus, there are no reasons to expect the Arrhenius behavior

for \(k(T,{\mathbb{S}})\) to hold below _T_ = 1206 K because the mobility of liquid atoms close to the crystal surface is smaller than in the bulk due to IIO effects. PRINCIPAL COMPONENT

ANALYSIS  Figure 4b was obtained by applying the PCA48 to a data set containing equal amounts of crystallizing, liquid, and crystal, for a total of 60,00 data points. Crystal atoms were

defined as those for which \({\bar{\alpha }}_{i}(t)=1.0\) (bulk crystal atoms) or \({\bar{\alpha }}_{i}(t)=0.75\) (stacking fault atoms) for _t_ ∈ [_t_0−_τ__x_, _t_0 + _τ__x_] with _τ__x_ = 

15 ps. From the PCA we obtained the components of each data point along the two eigenvectors with largest eigenvalues, which are used to plot Fig. 4b. The atom trajectory was obtained by

applying the same PCA transformation along a single 3 ns trajectory of an atom in a simulation at 1500 K. GROWTH RATE DETERMINATION The number of atoms in the crystallite _N_(_t_) at any

given time _t_ was determined as the number of atoms with \({\bar{\alpha }}_{i}(t)> 0.25\). From this information the effective crystallite radius \({{\mathcal{R}}}_{{\rm{eff}}}(t)\)

(shown in Supplementary Fig. 10a) was computed assuming a spherical shape (which results in \(\kappa =2/{{\mathcal{R}}}_{{\rm{eff}}}\), where _κ_ is the geometrical factor in the

Gibbs–Thomson term). The growth rate for each temperature, shown in Fig. 3a, was determined by a linear fit of _N_(_t_) over the time interval for which \({{\mathcal{R}}}_{{\rm{eff}}}\in

[80\mathring{\rm{A}} ,100\mathring{\rm{A}} ]\). This interval is such that the crystallite is small enough to not be affected by finite-size effects, but large enough to give the system time

to equilibrate into a steady-state growth condition. See Supplementary Note 4 for a more detailed analysis. Error bars in Fig. 3a represent the 95% confidence interval as computed using the

bootstrap method with 1000 samples of the same size as the original distribution. INTERFACE TEMPERATURE When studying crystal growth, it is important to differentiate between the

temperature of the supercooled liquid surrounding the crystal (but far from the interface) from the solid–liquid interface temperature. Under steady-growth conditions these two temperatures

will differ because of the latent heat released at the interface and the finite rate of heat transport. Here, the interface temperature was computed by considering only the kinetic energy of

atoms with \({\bar{\alpha }}_{i}\in (0.15,0.75)\). This interval of \(\bar{\alpha }\) was selected because it includes both, interfacial liquid and interfacial crystal atoms. Supplementary

Fig. 10b shows the interface temperature under steady-state growth as a function of the surrounding liquid bath temperature. See Supplementary Note 4 for more details. CRYSTAL SURFACE

ANALYSIS  Figure 6b was obtained by constructing a polyhedral surface mesh around the crystallite (i.e. atoms with \({\bar{\alpha }}_{i}(t)\,> \, 0.25\)) using the algorithms by

Stukowski49 (as implemented in Ovito46) with a probe-sphere radius of 3.0 Å and a smoothing level of 10. From this mesh the surface directions were inferred and averaged over the time

interval for which the crystal growth occurs in a steady state. The data for constructing Fig. 6b was obtained by finding the orientation of the closest surface to each crystallizing atom.

The only atomically smooth surfaces in silicon are {111} surfaces22, consequently steps can only exist in these surfaces. Hence, the step–step separation distance shown in Fig. 6c was

computed assuming that (111) faceting occurs at all surface orientations. SOLID AND LIQUID FREE ENERGIES The accurate calculation of the solid and liquid free energies is crucial in crystal

growth studies. As shown in the Supplementary Note 5, employing approximations such as the quasi-harmonic approximation results in the underestimation of the predicted growth rates by as

much as 36%. For this reason, we performed the solid and liquid free energies calculations using state-of-the-art nonequilibrium thermodynamic integration methods that make no approximating

assumptions on the physical characteristics of the system. The crystal free energy was determined using the nonequilibrium Frenkel–Ladd50,51,52 (FL) and the reversible scaling52,53 (RS)

methods, following closely the approach described by Freitas52. For both methods a system of 21,952 silicon atoms in the diamond cubic structure was employed. The thermodynamic switching was

performed in 200 ps for each direction, before which the system was equilibrated for 20 ps. The FL switching was realized for temperatures ranging from 100 to 2000 K in intervals of 100 K.

For each temperature the switching was repeated in five independent simulations to estimate the statistical uncertainty. Similarly, the RS switching was also repeated five times. The

_S_-shaped function was employed in the FL switching, while the RS switching was performed with _T_i = 100 K and _T_f = 2000 K under the constant d_T_/d_t_ constraint. The system’s

center-of-mass was kept fixed for the FL and RS simulations, while a Langevin54,55 thermostat with damping parameter of 0.1 ps was applied. For the RS method a chain Nosé–Hoover barostat

with damping parameter of 1 ps and chain length of three was used to keep zero pressure. The absolute free energies and a comparison with the harmonic and quasi-harmonic approximations56 can

be seen in Supplementary Fig. 11. Liquid free energies were computed using the Uhlenbeck–Ford57,58 (UF) and RS53,58 methods, following closely the approach described by Leite58. The liquid

free energy calculations had the same number of atoms, switching time, equilibration time, and thermostat as the crystal free-energy calculations. The liquid density was the equilibrium

density at zero pressure, with the thermal expansion dilation taken into account. For the UF method we used _p_ = 50, _σ_ = 1.5 Å, and a cutoff radius of _r_c = 5_σ_. The UF switching was

performed linearly with time, while the RS switching had the same time dependence as the crystal with _T__i_ = 2000 K and _T_f = 1100 K (the lower final temperature _T_f was chosen to avoid

the liquid vitrification at low temperatures). For both methods—UF and RS—the switchings were repeated in five independent simulations to estimate the statistical uncertainty. COPPER All

results for copper were obtained from simulations that followed the exact same specifications as the simulations described above for silicon. The only modifications performed are described

in this section. The interaction between copper atoms was described using the embedded-atom method59 interatomic potential of Foiles et al.60. The timestep was selected as approximately

1/66th of the period of the highest-frequency phonon mode of this system, or Δ_t_ = 2 fs. The crystal growth simulations contained 1,000,000 atoms (notice that this is twice the size of the

simulations for silicon) and were initialized with a spherical crystalline seed of ~24,000 atoms (eight times larger than for silicon) in the face-centered cubic structure. The difference in

system size allowed us to explore much lower undercoolings for the crystal growth simulations, which ran for a total of 2 ns per temperature. In contrast to the simulations for silicon, the

snapshots saved from MD simulations were not subsequently relaxed before computing structural parameters because it has been shown that energy minimizations lead to significant

crystallization in metallic systems61. The crystal growth simulations were carried out at temperatures ranging from 900 to 1200 K in intervals of 25 K. The growth rate for each temperature,

shown in Fig. 7a, was determined by a linear fit of _N_(_t_) over the time interval for which \({{\mathcal{R}}}_{{\rm{eff}}}\in [100{{\AA}} ,120{{\AA}} ]\), as detailed in Supplementary Note

 4. The local-structure fingerprint \({{\mathbf{x}}}_{i}\) for the copper atoms was composed of a set of 37 radial structure functions defined by the following parameters: _σ_ = 0.3 Å,

_r_cut = 20.6 Å, and _r__n_ = (2.0 + 0.5_n_) Å, with _n_ = 0, 1, 2, …, 37. These parameters were obtained through the same hyperparameter optimization process applied to silicon, as

described in the Supplementary Note 1. The SVM classifier trained with the data collected for copper had accuracy of 97%. The results presented here were obtained using data from the

simulation at _T_ = 1200 K to train the SVM classifier. The total size of the data set for copper was 77.9 million data points. All results for the LSD model for copper (including Fig. 7 and

Supplementary Fig. 19c) were obtained using data from simulations at _T_ ≥ 1136 K. Below this temperature the kinetic factor showed signs of non-Arrhenius behavior for some values of

\({\mathbb{S}}\) due to the approaching glass transition temperature, similarly to what was observed for silicon. Notice that below _T_ = 1136 K the bulk liquid diffusivity also shows signs

of departure from the Arrhenius behavior. Thus, there are no reasons to expect the Arrhenius behavior for \(k(T,{\mathbb{S}})\) to hold below _T_ = 1136 K because the mobility of liquid

atoms close to the crystal surface is smaller than in the bulk due to IIO effects. The solid and liquid free energies were computed for systems containing 19,652 copper atoms. The FL and RS

methods were applied with a switching time of 400 ps for each direction, preceded by an equilibration time of 40 ps. The FL switching was realized for temperatures ranging from 100 to 1300 K

in intervals of 100 K. The RS switching for the solid was performed with _T_i = 100 K and _T_f = 1300 K, while for the liquid we used _T_i = 2000 K and _T_f = 900 K. For the UF free-energy

calculations we used _p_ = 75, _σ_ = 1.3 Å, and a cutoff radius of _r_c = 5_σ_. DATA AVAILABILITY The data that support the findings of this study are available from the corresponding author

upon reasonable request. CODE AVAILABILITY The code and scripts used to generate the results in this paper can be downloaded from the following repository:

https://github.com/freitas-rodrigo/CrystallizationMechanismsFromML. Any custom code that is not currently available in the repository can be subsequently added to the repository upon request

to the corresponding author. REFERENCES * Asta, M. et al. Solidification microstructures and solid-state parallels: Recent developments, future directions. _Acta Mater._ 57, 941 (2009). CAS

  Google Scholar  * Libbrecht, K. G. Physical dynamics of ice crystal growth. _Annu. Rev. Mater. Res._ 47, 271 (2017). CAS  Google Scholar  * Burton, W.-K., Cabrera, N. & Frank, F. The

growth of crystals and the equilibrium structure of their surfaces. _Philos. Trans. R. Soc. Lond. Ser. A_ 243, 299 (1951). ADS  MathSciNet  MATH  Google Scholar  * Chernov, A. A. _Modern

Crystallography III: Crystal Growth_, Vol. 36 (Springer Science & Business Media, 2012). * Pimpinelli, A. & Villain, J. _Physics of Crystal Growth_ (Cambridge University Press,

1999). * Saito, Y. _Statistical Physics of Crystal Growth_ (World Scientific, 1996). * Wilson, H. W. On the velocity of solidification and viscosity of super-cooled liquids. _Lond. Edinb.

Dublin Philos. Mag. J. Sci._ 50, 238 (1900). MATH  Google Scholar  * Frenkel, J. _Kinetic Theory of Liquids_ (The Clarendon Press, Oxford, 1946). * Stolk, P., Polman, A. & Sinke, W.

Experimental test of kinetic theories for heterogeneous freezing in silicon. _Phys. Rev. B_ 47, 5 (1993). ADS  CAS  Google Scholar  * Sun, G., Xu, J. & Harrowell, P. The mechanism of the

ultrafast crystal growth of pure metals from their melts. _Nat. Mater._ 17, 881 (2018). ADS  CAS  PubMed  Google Scholar  * Hoyt, J. & Asta, M. Atomistic computation of liquid

diffusivity, solid-liquid interfacial free energy, and kinetic coeffcient in Au and Ag. _Phys. Rev. B_ 65, 214106 (2002). ADS  Google Scholar  * Freitas, R., Frolov, T. & Asta, M. Step

free energies at faceted solid surfaces: Theory and atomistic calculations for steps on the Cu (111) surface. _Phys. Rev. B_ 95, 155444 (2017). ADS  Google Scholar  * Cubuk, E. D. et al.

Identifying structural ow defects in disordered solids using machine-learning methods. _Phys. Rev. Lett._ 114, 108001 (2015). ADS  CAS  PubMed  Google Scholar  * Schoenholz, S. S., Cubuk, E.

D., Sussman, D. M., Kaxiras, E. & Liu, A. J. A structural approach to relaxation in glassy liquids. _Nat. Phys._ 12, 469 (2016). CAS  Google Scholar  * Sharp, T. A. et al. Machine

learning determination of atomic dynamics at grain boundaries. _Proc. Natl Acad. Sci. USA_ 115, 10943 (2018). ADS  CAS  PubMed  Google Scholar  * Behler, J. & Parrinello, M. Generalized

neural-network representation of high-dimensional potential-energy surfaces. _Phys. Rev. Lett._ 98, 146401 (2007). ADS  PubMed  Google Scholar  * Cortes, C. & Vapnik, V. Support-vector

networks. _Mach. Learn._ 20, 273 (1995). MATH  Google Scholar  * Pedregosa, F. et al. Scikit-learn: Machine learning in Python. _J. Mach. Learn. Res._ 12, 2825 (2011). MathSciNet  MATH 

Google Scholar  * Chang, C. -C., & Lin, C. -J. LIBSVM: A library for support vector machines. _ACM Trans. Intell. Syst. Technol._ 2, 27:1–27:27 (2011). * Fitzner, M., Sosso, G. C., Cox,

S. J. & Michaelides, A. Ice is born in low-mobility regions of supercooled liquid water. _Proc. Natl Acad. Sci. USA_ 116, 2009 (2019). ADS  CAS  PubMed  Google Scholar  * Moore, E. B.

& Molinero, V. Structural transformation in supercooled water controls the crystallization rate of ice. _Nature_ 479, 506 (2011). ADS  CAS  PubMed  Google Scholar  * Fujiwara, K. Crystal

growth behaviors of silicon during melt growth processes. _Int. J. Photoenergy_ 2012, (2012). * Oh, S. H., Kauffmann, Y., Scheu, C., Kaplan, W. D. & Rühle, M. Crystal growth behaviors

of silicon during melt growth processes. _Science_ 310, 661 (2005). ADS  CAS  PubMed  Google Scholar  * Oh, S. H. et al. Oscillatory mass transport in vapor-liquid-solid growth of sapphire

nanowires. _Science_ 330, 489 (2010). ADS  CAS  PubMed  Google Scholar  * Panciera, F. et al. Surface crystallization of liquid Au-Si and its impact on catalysis. _Adv. Mater._ 31, 1806544

(2019). Google Scholar  * Kaplan, W. D. & Kauffmann, Y. Structural order in liquids induced by interfaces with crystals. _Annu. Rev. Mater. Res._ 36, 1 (2006). ADS  CAS  Google Scholar 

* de Jong, A. et al. Complex geometric structure of a simple solid-liquid interface: Gan (0001)-ga. _Phys. Rev. Lett._ 124, 086101 (2020). ADS  PubMed  Google Scholar  * Huisman, W. J. et

al. Layering of a liquid metal in contact with a hard wall. _Nature_ 390, 379 (1997). ADS  CAS  Google Scholar  * Reedijk, M., Arsic, J., Hollander, F., De Vries, S. & Vlieg, E. Liquid

order at the interface of kdp crystals with water: Evidence for icelike layers. _Phys. Rev. Lett._ 90, 066103 (2003). ADS  CAS  PubMed  Google Scholar  * Teich, E. G., van Anders, G.,

Klotsa, D., Dshemuchadse, J. & Glotzer, S. C. Clusters of polyhedra in spherical confinement. _Proc. Natl Acad. Sci. USA_ 113, E669 (2016). CAS  PubMed  Google Scholar  * Spaepen, F. A

structural model for the solid-liquid interface in monatomic systems. _Acta Metall._ 23, 729 (1975). CAS  Google Scholar  * Bracht, H., Silvestri, H., Sharp, I. & Haller, E. Selfand

foreign-atom diffusion in semiconductor isotope heterostructures. II. Experimental results for silicon. _Phys. Rev. B_ 75, 035211 (2007). ADS  Google Scholar  * Sciortino, F., Geiger, A.

& Stanley, H. E. Effect of defects on molecular mobility in liquid water. _Nature_ 354, 218 (1991). ADS  CAS  Google Scholar  * Sastry, S. & Angell, C. A. Liquid-liquid phase

transition in supercooled silicon. _Nat. Mater._ 2, 739 (2003). ADS  CAS  PubMed  Google Scholar  * Plimpton, S. Fast parallel algorithms for short-range molecular dynamics. _J. Comput.

Phys._ 117, 1 (1995). ADS  CAS  MATH  Google Scholar  * Stillinger, F. H. & Weber, T. A. Computer simulation of local order in condensed phases of silicon. _Phys. Rev. B_ 31, 5262

(1985). ADS  CAS  Google Scholar  * Bulatov, V. & Cai, W. _Computer Simulations of Dislocations_, Vol. 3 (Oxford University Press on Demand, 2006). * Bussi, G., Donadio, D. &

Parrinello, M. Canonical sampling through velocity rescaling. _J. Chem. Phys._ 126, 014101 (2007). ADS  PubMed  Google Scholar  * Martyna, G. J., Tobias, D. J. & Klein, M. L. Constant

pressure molecular dynamics algorithms. _J. Chem. Phys._ 101, 4177 (1994). ADS  CAS  Google Scholar  * Parrinello, M. & Rahman, A. Polymorphic transitions in single crystals: A new

molecular dynamics method. _J. Appl. Phys._ 52, 7182 (1981). ADS  CAS  Google Scholar  * Tuckerman, M. E., Alejandre, J., López-Rendón, R., Jochim, A. L. & Martyna, G. J. A

Liouville-operator derived measure-preserving integrator for molecular dynamics simulations in the isothermal{isobaric ensemble. _J. Phys. A_ 39, 5629 (2006). ADS  MathSciNet  CAS  MATH 

Google Scholar  * Shinoda, W., Shiga, M. & Mikami, M. Rapid estimation of elastic constants by molecular dynamics simulation under constant stress. _Phys. Rev. B_ 69, 134103 (2004). ADS

  Google Scholar  * Tuckerman, M. _Statistical Mechanics: Theory and Molecular Simulation_ (Oxford University Press, 2010) * Rein ten Wolde, P., Ruiz-Montero, M. J. & Frenkel, D.

Numerical calculation of the rate of crystal nucleation in a Lennard-Jones system at moderate undercooling. _J. Chem. Phys._ 104, 9932 (1996). ADS  Google Scholar  * Larsen, P. M., Schmidt,

S. & Schiøtz, J. Robust structural identification via polyhedral template matching. _Model. Simul. Mater. Sci. Eng._ 24, 055007 (2016). ADS  Google Scholar  * Stukowski, A. Visualization

and analysis of atomistic simulation data with OVITO-the Open Visualization Tool. _Model. Simul. Mater. Sci. Eng._ 18, 015012 (2009). ADS  Google Scholar  * Hunter, J. D. Matplotlib: A 2D

Graphics Environment. _Comput. Sci. Eng._ 9, 90 (2007). Google Scholar  * Goodfellow, I., Bengio, Y. & Courville, A. _Deep Learning_ (MIT Press, 2016) . * Stukowski, A. Computational

analysis methods in atomistic modeling of crystals. _J. Miner. Met. Mater. Soc._ 66, 399 (2014). CAS  Google Scholar  * Frenkel, D. & Ladd, A. J. New Monte Carlo method to compute the

free energy of arbitrary solids. Application to the fcc and hcp phases of hard spheres. _J. Chem. Phys._ 81, 3188 (1984). ADS  CAS  Google Scholar  * de Koning, M. & Antonelli, A.

Einstein crystal as a reference system in free energy estimation using adiabatic switching. _Phys. Rev. E_ 53, 465 (1996). ADS  Google Scholar  * Freitas, R., Asta, M. & de Koning, M.

Nonequilibrium free-energy calculation of solids using LAMMPS_Comput. Mater. Sci._ 112, 333 (2016). CAS  Google Scholar  * de Koning, M., Antonelli, A. & Yip, S. Optimized freeenergy

evaluation using a single reversible-scaling simulation. _Phys. Rev. Lett._ 83, 3973 (1999). ADS  Google Scholar  * Schneider, T. & Stoll, E. Molecular-dynamics study of a

three-dimensional one-component model for distortive phase transitions. _Phys. Rev. B_ 17, 1302 (1978). ADS  CAS  Google Scholar  * Allen, M. P. & Tildesley, D. J. _Computer Simulation

of Liquids_ (Oxford university press, 2017). * Fultz, B. _Phase Transitions in Materials_ (Cambridge University Press, 2014). * Paula Leite, R., Freitas, R., Azevedo, R. & de Koning, M.

The uhlenbeck-ford model: Exact virial coeffcients and application as a reference system in uid-phase freeenergy calculations. _J. Chem. Phys._ 145, 194101 (2016). ADS  PubMed  Google

Scholar  * Leite, R. P. & de Koning, M. Nonequilibrium freeenergy calculations of uids using LAMMPS. _Comput. Mater. Sci._ 159, 316 (2019). Google Scholar  * Daw, M. S. & Baskes, M.

I. Embedded-atom method: Derivation and application to impurities, surfaces, and other defects in metals. _Phys. Rev. B_ 29, 6443 (1984). ADS  CAS  Google Scholar  * Foiles, S., Baskes, M.

& Daw, M. S. Embedded-atommethod functions for the fcc metals cu, ag, au, ni, pd, pt, and their alloys. _Phys. Rev. B_ 33, 7983 (1986). ADS  CAS  Google Scholar  * Sun, G., Xu, J. &

Harrowell, P. The mechanism of the ultrafast crystal growth of pure metals from their melts. _Nat. Mater._ 17, 881 (2018). ADS  CAS  PubMed  Google Scholar  * Stukowski, A. & Albe, K.

Extracting dislocations and on-dislocation crystal defects from atomistic simulation data. _Model. Simul. Mater. Sci. Eng._ 18, 085001 (2010). ADS  Google Scholar  * Stukowski, A., Bulatov,

V. V. & Arsenlis, A. Automated identification and indexing of dislocations in crystal interfaces. _Model. Simul. Mater. Sci. Eng._ 20, 085007 (2012). ADS  Google Scholar  * Galvin, G.,

Mayer, J. & Peercy, P. Solidification kinetics of pulsed laser melted silicon based on thermodynamic considerations. _Appl. Phys. Lett._ 46, 644 (1985). ADS  CAS  Google Scholar  *

Hedler, A., Klaumünzer, S. L. & Wesch, W. Amorphous silicon exhibits a glass transition. _Nat. Mater._ 3, 804 (2004). ADS  CAS  PubMed  Google Scholar  * Shou, W. & Pan, H.

Silicon-wall interfacial free energy via thermodynamics integration. _J. Chem. Phys._ 145, 184702 (2016). ADS  PubMed  Google Scholar  Download references ACKNOWLEDGEMENTS The authors would

like to thank Dr. Vasily Bulatov for providing the script for projecting crystal directions into the standard triangle, Prof. Suneel Kodambaka for bringing to our attention the literature on

interface-induced ordering of liquids, and Dr. Pablo Damasceno for the clarifications regarding the entropic origins of layering in weakly bonded liquids. The authors also acknowledge the

fruitful discussions with Prof. Qian Yang, Gowoon Cheon, Evan Antoniuk, and Yanbing Zhu. This work was supported by the Department of Energy National Nuclear Security Administration under

Award Number DE-NA0002007, National Science Foundation grants DMREF-1922312, and CAREER-1455050. AUTHOR INFORMATION AUTHORS AND AFFILIATIONS * Department of Materials Science and

Engineering, Stanford University, Stanford, CA, 94305, USA Rodrigo Freitas & Evan J. Reed Authors * Rodrigo Freitas View author publications You can also search for this author inPubMed 

Google Scholar * Evan J. Reed View author publications You can also search for this author inPubMed Google Scholar CONTRIBUTIONS R.F. performed all calculations and data analysis. R.F. and

E.J.R. worked jointly on the interpretation of the data, project design, and manuscript preparation. CORRESPONDING AUTHOR Correspondence to Rodrigo Freitas. ETHICS DECLARATIONS COMPETING

INTERESTS The authors declare no competing interests. ADDITIONAL INFORMATION PEER REVIEW INFORMATION _Nature Communications_ thanks the anonymous reviewer(s) for their contribution to the

peer review of this work. PUBLISHER’S NOTE Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations. SUPPLEMENTARY INFORMATION

SUPPLEMENTARY INFORMATION DESCRIPTION OF ADDITIONAL SUPPLEMENTARY FILES SUPPLEMENTARY MOVIE 1 SUPPLEMENTARY MOVIE 2 SUPPLEMENTARY MOVIE 3 SUPPLEMENTARY MOVIE 4 SUPPLEMENTARY MOVIE 5

SUPPLEMENTARY MOVIE 6 SUPPLEMENTARY MOVIE 7 SUPPLEMENTARY MOVIE 8 RIGHTS AND PERMISSIONS OPEN ACCESS This article is licensed under a Creative Commons Attribution 4.0 International License,

which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link

to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons license, unless

indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons license and your intended use is not permitted by statutory regulation or

exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/. Reprints

and permissions ABOUT THIS ARTICLE CITE THIS ARTICLE Freitas, R., Reed, E.J. Uncovering the effects of interface-induced ordering of liquid on crystal growth using machine learning. _Nat

Commun_ 11, 3260 (2020). https://doi.org/10.1038/s41467-020-16892-4 Download citation * Received: 17 April 2020 * Accepted: 26 May 2020 * Published: 26 June 2020 * DOI:

https://doi.org/10.1038/s41467-020-16892-4 SHARE THIS ARTICLE Anyone you share the following link with will be able to read this content: Get shareable link Sorry, a shareable link is not

currently available for this article. Copy to clipboard Provided by the Springer Nature SharedIt content-sharing initiative