Neural surrogates for crystal growth dynamics with variable supersaturation: explicit vs. implicit conditioning Matteo Rigoni,1 Daniele Lanzoni,1, 2 Francesco Montalenti,1 and Roberto Bergamaschini1
arXiv:2604.21753v1 [cond-mat.mtrl-sci] 23 Apr 2026
1
Dept. of Materials Science, University of Milano-Bicocca, 20125 Milano, Italy 2 Dept. of Physics, University of Genova, 16146 Genova, Italy
Simulations of crystal growth are performed by using Convolutional Recurrent Neural Network surrogate models, trained on a dataset of time sequences computed by numerical integration of Allen-Cahn dynamics including faceting via kinetic anisotropy. Two network architectures are developed to take into account the effects of a variable supersaturation value. The first infers it implicitly by processing an input mini-sequence of a few evolution frames and then returns a consistent continuation of the evolution. The second takes the supersaturation parameter as an explicit input along with a single initial frame and predicts the entire sequence. The two models are systematically tested to establish strengths and weaknesses, comparing the prediction performance for models trained on datasets of different size and, in the first architecture, different lengths of input mini-sequence. The analysis of point-wise and mean absolute errors shows how the explicit parameter conditioning guarantees the best results, reproducing with high-fidelity the ground-truth profiles. Comparable results are achievable by the mini-sequence approach only when using larger training datasets. The trained models show strong conditioning by the supersaturation parameter, consistently reproducing its overall impact on growth rates as well as its local effect on the faceted morphology. Moreover, they are perfectly scalable even on 256 times larger domains and can be successfully extended to more than 10 times longer sequences with limited error accumulation. The analysis highlights the potential and limits of these approaches in view of their general exploitation for crystal growth simulations.
I.
INTRODUCTION
Machine Learning (ML) and particularly Deep Learning (DL) is currently revolutionizing materials science, both from an experimental and a computational standpoint [1–4]. In the latter context, one of the perhaps most intriguing possibilities brought by data-driven approaches is the significant reduction of the computational costs associated with simulations, with the possibility of screening new materials and processes faster than ever before [5–8], with potential disruptive impacts on the speed of technological and fundamental discoveries. One of the applications that could most benefit from the use of ML methods are continuum models for the temporal evolution of materials. This class of simulations allows the study of nano- to macro-scale processes, such as crystalline materials growth, microstructure formation and evolution, mechanical and electrical responses, etc. Among these applications, problems involving interface motion are particularly challenging, since they often require computationally expensive numerical procedures and the solution of challenging coupled PDEs. One of the traditional ways to tackle this kind of issue is through phase-field (PF) methods [9, 10], in which the description of the system geometry or microstructure is traced implicitly by the evolution of one or multiple order parameters. Still, PF approaches could be nonetheless expensive, due to the requirement of fine meshes and of the solution of stiff nonlinear equations. Indeed, several ML methods have been proposed in recent years to mitigate these issues. These have exploited many of the state-of-the-art methods in DL, such as autoencoders for the so-called ”latent space models” [11–14],
Convolutional Neural Networks (CNNs) [15, 16], neural operators [13, 17], and graph NN [18] to cite a few. A comprehensive review can be found at Ref. [19]. Among other approaches, the use of Convolutional Recurrent Neural Networks (CRNN) [20, 21] is particularly interesting, as they allow to jointly capture spatial (Convolutional aspect) and temporal (Recurrent aspect) correlations in data. Several studies have already shown how CRNN may be used to learn and reproduce at a smaller computational cost simulations of the evolution of technologically and theoretically relevant mesoscale phenomena [22–25], such as grain growth, solidification processes, and spinodal decomposition. Moreover, the fully-convolutional nature of the approach allows for generalizations to arbitrary domain size, increasing the range of applicability of these models. One of the peculiarities of Recurrent NNs is their possibility of both generating and processing temporal sequences. The second property has been explicitly used in some works exploiting CRNN to predict the microstructural and morphological evolution of materials based on short initial sequences [23, 26]. Perhaps more importantly, this possibility would be particularly interesting in experimental contexts and situations in which the dynamics is driven by hidden or unknown factors. Moreover, instead of explicitly providing the CRNN model with all the quantity determining the evolution of a system, one could let the NN automatically infer the correct evolution law from a small sequence instead of from a single initial condition. Indeed, this idea was exploited in Ref. [22], where the solidification latent heat variable was hidden from the CRNN and was implicitly inferred by the model from the provided 10 initial steps.
2 The opposite possibility of explicitly providing as input the parameters defining the dynamics and a single initial condition, which is closest to a traditional simulation scheme, has also been recently explored. Some works [16, 25] have indeed conditioned the NN prediction with external parameters. For example, in Ref. [25], the mismatch between lattice parameters for pure phases was passed to the CRNN to predict the spinodal decomposition in coherent alloys. A direct comparison of these two competing approaches using similar architectures, however, is not currently available in the literature. Our work addresses this question and quantitatively compares the prediction performances of CRNN models derived from Refs. [25, 27] for the case study of crystal growth. Indeed, the acceleration of simulations and composition/process screening offered by ML tools are particularly appealing in this field and interest in these applications has substantially increased in recent years [28–31]. In the case at hand, we consider the prototypical Allen-Cahn dynamics [10, 32] in two dimensions. In order to mimic crystal faceting, key for realistic and reliable simulations, we introduce an anisotropic kinetic coefficient producing hexagonal morphologies. Variable values of supersaturation, i.e., of the difference in chemical potential between the gas or liquid mother phase and the crystal itself, are admitted. Hence, since the supersaturation directly controls the actual growth rate of the crystalline phase and the finedetails of faceting, the evolution sequence from a given initial state is parametrically dependent on its value. Although very simple, this model contains the fundamental aspects that characterize the faceted growth of crystals while giving the possibility of constructing large datasets with little computational effort, thus enabling the extensive testing of the NN behaviour. The study finds that explicitly conditioning on the supersaturation value reliably increases the quality of predictions and yields more data-efficient NN models. In fact, we find that to obtain the same accuracy, a ≈ 15× larger dataset is required for CRNN exploiting small sequences for implicit inference with respect to explicitly conditioned ones. Based on our findings, this strategy should therefore be regarded as the better option, i.e., explicit information should always be preferred over the automatic correlation discovery capabilities of deep learning methods, whenever possible. The paper is organized as follows. In Sect. II we define the Allen-Cahn model with the phase-field approach and provide the technical description of the two NN architectures. Then in Sect. III we analyse in-depth the NN capabilities. First, we analyse the training results (Sect. III A). Then, we quantitatively estimate the accuracy of the different trained models on extensive test sets (Sect. III B), inspecting their sensitivity to different supersaturation values. Finally, in Sect. III D, we exploit the full potential of our CRNN architecture for largescale simulations, beyond the training domain size, and in Sect. III E we analyse the effect on prediction quality
=0.2
20 τ
40 τ
60 τ
80 τ
=0.8 FIG. 1. Example of evolution sequences from the same initial configuration for the two limiting values of ∆µ, i.e. 0.2 and 0.8 considered in our training set. Each time unit τ corresponds to 25 time steps.
due to the density of seeds in the initial conditions. II. A.
METHODS
Phase-field simulations
A simple two-dimensional Allen-Cahn model is considered as a prototypical description of crystal growth under isothermal conditions [9]. In particular, we consider an order parameter φ distinguishing the crystal phase (φ = 1) from its liquid or gaseous mother phase (φ = 0) and set a standard Ginzburg-Landau energy functional G[φ] having the two phases as minima, biased by a supersaturation parameter ∆µ: Z h i ϵ G[φ] = |∇φ|2 + w(φ) − p(φ)∆µ dx (1) Ω 2 with ϵ the width of the diffuse boundary between the two phases, Ω the integration domain, and a bulk energy set by the symmetric double-well potential w(φ) = (18/ϵ)φ2 (1 − φ)2 and a smooth, bias function set as p(φ) = φ3 (10 − 15φ + 6φ2 ). The growth process is then traced implicitly by the time evolution of φ as determined by the local chemical potential µ = δG/δφ as dφ δG = −k(α) = −k(α) −ϵ∇2 φ + w′ (φ) − p′ (φ)∆µ dt δφ (2) with k(α) a kinetic coefficient, eventually dependent on the local profile orientation α = arctan(∇φy /∇φx ) in order to account for anisotropic growth rates. In the present work, we consider k(α) = 1 + β cos(N α), with β = 0.8 and N = 6 to produce hexagonally faceted crystals. For the sake of simplicity, all seeds share the same crystallographic axis thus mimicking epitaxial growth. Eq. 2 is here solved numerically by exploiting a simple forward Euler integration scheme and finite differences on a square grid with periodic boundary conditions. A value of ϵ = 6 pixels is set and a time step of 0.05/ϵ is used.
3
B.
Neural Network approach
In this work, the Convolutional Recurrent Neural Network (CRNN) developed in [24, 27] is specialized to reproduce the Allen-Cahn dynamics and adapted to tackle a heterogeneous dataset of time evolution sequences as determined by randomly chosen supersaturation parameters ∆µ. The CRNN objective is to generate a sequence of subsequent stages during the crystalline material growth, provided one or more snapshots of previous morphologies. In this work, we focus on two variants of the architecture, whose main difference is whether the model is aware of the actual value of the supersaturation parameter. In practice, both models are implemented as a convolutional recurrent GRU [20, 27], composed of three recurrent stacked layers as sketched in Fig. 2(a). Kernel size used is 5 × 5, with circular padding encoding periodic boundary conditions by construction. Hidden states have 16 channels, for a total of about 140000 parameters for both models. The next state of the system φt+τ is obtained by summing the output of the NNs to the previous phase-field map φt , i.e., the models are required to learn the residual between subsequent states of the system. At variance with the work in Ref. [24], no mass conservation is required in the current case, hence the NN output does not need any additional processing. A schematic for both architectures is shown in Fig. 2(b) and Fig. 2(c). The first NN, which will be called NNseq in the following, does not receive ∆µ as an explicit parameter, but rather is required to infer it by processing as input a mini-sequence of s frames. This is made possible by the recurrent layers in the NN architecture, which act as a sort of ”memory” capable of implicitly extracting temporal information from input sequences. While this approach is particularly appealing for applications where the driving forces for the dynamics are not known or not completely specified (e.g., in experimental contexts), it also has the downside of always requiring multiple timeframes as input for the NN model. In this sense, the goal of NNseq could be understood as a sequence-completion task. On the other hand, the second model NNpar is explic-
ConvGRU 1
(a)
ConvGRU 2
CRNN
In Fig. 1 we report two evolution sequences starting from the same initial profile and setting a different supersaturation ∆µ. As made evident by the similarity between the profile obtained for ∆µ = 0.2 at time 80τ and the one for ∆µ = 0.8 at 20τ , the main effect of ∆µ is to directly control the overall growth rate, providing a relative rescaling of the time-scale. A closer inspection however shows substantial differences in the contours, more rounded for the lowest ∆µ and characterized by sharper corners for the highest value, as due to the different weights of the gradient-energy contribution and the double-well imbalance (∝ ∆µ) to the chemical potential in eq. 2.
ConvGRU 3 toScalarField
prediction
input (s frames)
(b) NNseq
CRNN ... CRNN
CRNN
input
(c) NNpar
CRNN
CRNN ...
prediction +
+
+
CRNN
CRNN
CRNN ...
FIG. 2. (a) Schematics of the used convolutional recurrent neural network (CRNN) internal architecture. (b,c) NNseq and NNpar schematics, showing the difference in the input handling. The internal architecture of the CRNN block is the one reported in panel (b).
itly provided the actual value of the supersaturation ∆µ. To this aim, the numerical value of the supersaturation is converted into a uniform spatial tensor and concatenated to the φ map for the current timestep, using the same strategy as in [25]. The input of NNpar model is therefore composed of a two-channel image, with the first channel containing the material morphology as described by the phase field φ and the second being a constant map equal to ∆µ. This simple but effective encoding scheme can, in principle, be straightforwardly extended to nonconstant driving forces. Once ∆µ and the initial condition are explicitly known, the full evolution of the system is completely determined, since the underlying equation of motion Eq. 2 is first-order in time. For this reason, NNpar is provided with a single snapshot to predict the full growth.
4 RESULTS AND DISCUSSION A.
Training and Validation
A dataset of 7500 PF simulation sequences is constructed for the training of both NN models. Each case consists of a series of 200 time-frames (taken every τ = 25 integration time-steps) reporting the φ field as evolving from a random distribution of elliptical seeds of different eccentricity and orientation toward complete coverage of the domain. A 128 × 128 domain is considered and each grid value φ(i, j) ∈ [0, 1] defines a pixel in the gray-scale frame image. The value of supersaturation ∆µ is uniformly sampled within the range [0.2, 0.8], corresponding to a ≈ 4× variation in the net growth rates (see Fig. 1) and encompassing both rounded and sharp facet corners. A set of Nset sequences is extracted from the full dataset for the training of the NN, with a 4 : 1 random partitioning between training and validation. For each sample, only a segment of T = 50 contiguous time frames is selected at random to expose the NN to diverse evolution stages. The NN training is performed using the standard implementation of the Adam optimizer [33]. The loss function L is defined as the common mean squared error (MSE) between the PF ground-truth (φPF ) and NNpredicted (φNN ) profiles, averaged on all time-frames t of each sequence and on all Nb samples in the training batch: L(ϑ) =
Nts X T X
2 1 NN ⟨ φPF n (t) − φn (t|ϑ) ⟩ Nb T n=1 t=1
of input frames while requesting the prediction of the increasing number of remaining sequence steps, down to the actual mini-sequence of length s for NNseq or single frame for NNpar decided as input. During the process, the loss is evaluated on the predicted frames. While resulting in an apparent raise of the loss, this gradual increment of the complexity of the NN task has been demonstrated to return more stable and efficient convergence of the model. For the present work, different models have been trained by varying both the size of the dataset Nset and, for NNseq , the length s of the input mini-sequence, as listed in Table I. All models showed convergence of the training procedure, with the training and validation losses steadily decreasing, albeit for fluctuations introduced by the stochastic Adam algorithm, with no sign of overfitting. In Fig. 3, we report the evolution of the training and validation losses obtained over 300 training epochs for the
10-2 10-3 10-4
with ϑ the set of NN parameters and ⟨ . ⟩ indicates the spatial average. The same loss function is also computed on the validation set. Mirror symmetry is taken into account as in Refs. [24, 27] via data augmentation. A mini-batch size of Nb = 3 has been used. Given the recurrent architecture, the curriculum learning technique [34] is implemented, starting from the prediction of only the last frame of a sequence in the first training epoch and then gradually reducing the number
10-5
NNseq
NNpar
TABLE I. List of NN models investigated in the present study for both NNseq and NNpar architectures, labelled according to dataset size Nset and length s of the input sequence.
0
50
100
150
200
250
300
Training epoch
(b) NNpar
100
training validation
10-1
Loss
ID Nset s NN1.5k 1500 5 5 NN5k 5000 5 5 NN7.5k 7500 5 5 5k NN3 5000 3 NN5k 5000 7 7 NN0.5k 500 1 ∆µ NN1.5k 1500 1 ∆µ 5k NN∆µ 5000 1
training validation
10-1
(3)
arch.
(a) NNseq
100
Loss
III.
10-2 10-3 10-4 10-5
0
50
100
150
200
250
300
Training epoch FIG. 3. Training and validation losses during training for the mini-sequence model NNseq (a) and the parameter model NNpar (b). The cases for smaller and larger dataset sizes are shown for both.
5
FIG. 4. Evolution sequence for a representative case (∆µ = 0.6) in the test set comparing the ground-truth PF solution (a) with the corresponding predictions from both NNseq (b,c) and NNpar (d,e). The pixel-by-pixel error δ is superimposed on the predicted images by colour bar. (f) Time evolution of the MAE between NN predictions and the true PF solution.
models with larger and smaller datasets for both NNseq (a) and NNpar (b). Reasonably, the models trained on the largest datasets exhibit a faster decay rate of the loss for each epoch, leading to a final value about one order of magnitude lower than those trained on the smallest. This trend is consistently observed also for the interme1.5k diate models NN5k 5 and NN∆µ (see Fig. S1 of Supplementary Material). The effect of considering a different mini-sequence length s is instead less evident, but slightly lower loss values for longer input sequences may be observed (see Fig. S2 of Supplementary Material). Notably, NN7.5k and NN5k 5 ∆µ return similar loss values despite the different architecture and dataset size, suggesting similar prediction performances. This expectation will, however, be confuted by the extensive testing in the following sections. After each training run, the model with the lowest validation loss within the last 50 epochs is selected as the best performing one and used for all the following analyses.
B.
Testing the NN prediction performances
In order to evaluate the predictive performances of the trained models, an independent test set of 2500 simulations is generated with the same criteria used for the training and validation sets. All models in Table I are then used to predict the full 200τ -long time evolution sequence for all these cases to assess their performance. Notice that this task exceeds by 4× the 50-frame-long sequences to which the NNs are exposed during training, thus serving as a generalization test to longer simulation times. As a first evaluation, predicted frames can be ”visually” compared with the corresponding ground-truth ones computed from the numerical integration of the PF eq. (2). In Fig. 4 and Fig. 5 we consider two representative evolutions for different values of ∆µ (0.6 and 0.3 respectively) and compare the predictions for both NNseq and NNpar , again focusing on the models trained on the smallest and largest datasets so to appreciate the best and worst performance for both classes (see Fig. S3
6
FIG. 5. Evolution sequence for a representative case (∆µ = 0.3) in the test set comparing the ground-truth PF solution (a) with the corresponding predictions from both NNseq (b,c) and NNpar (d,e). The pixel-by-pixel error δ is superimposed on the predicted images by colour bar. (f) Time evolution of the MAE between NN predictions and the true PF solution.
and Fig. S4 of Supplementary Material for the analyses on all other NN models). In panels (b-e) of both figures, the pixel-by-pixel signed difference between the NNpredicted φ field and the PF one, i.e. δ = φNN − φPF , is computed and superimposed by colour map to the blackand-white output images so to evidence the regions where the predictions are in excess (δ > 0) or defect (δ < 0) with respect to the true ones (a). For a more quantitative estimation of the prediction errors, we also monitored the Mean Absolute error (MAE) defined as the domain average of the absolute value of δ, i.e., MAE=⟨|δ|⟩, as a function of time. Since the φ field is bounded in the range [0, 1], the MAE value can be directly taken as the fraction of mispredicted pixels in the image. The time evolutions of MAE by the different NN for both evolutions in Fig. 4 and Fig. 5 (see also Fig. S3(f) and Fig. S4(f)) are reported in the respective panels (f). In the case of Fig. 4, we observe a relatively fast growth leading to full domain coverage after about 120τ . We can observe how the evolution sequences are well cap-
tured by all models. Most discrepancies are localized in the regions where crystalline units coalesce, i.e., where the dynamics abruptly change because of the change in topology, causing a slight delay or speed up that does not significantly impact the overall evolution. The largest errors, still leading to a MAE always less than 0.06 are found for the models trained on smaller datasets NN1.5k 5 and NN0.5k ∆µ , while the best performance is achieved by NN5k ∆µ , which is always below 0.01 MAE. It can be noticed how, as time advances, the MAE generally grows because of error accumulation, resulting in steeper peaks when coalesce events occur and returning to zero when the domain is fully covered. In the case of Fig. 5, characterized by a slower dynamics leading just to a partial filling of the domain, we instead observe a neat difference in the prediction performance with the NNseq models struggling to advance the profile at the right rate in contrast to the NNpar models which are still providing a good reproduction of the sequence but for local errors. While NN1.5k results are 5 completely unreliable as tracing a constantly faster dy-
7
NNseq
NNpar
max(MAE)
0.5
0.1 0.05
0.01 0.005
0.9 perc.
IQR
{ median
0.1 perc.
FIG. 6. Distribution of the maximum of MAE for the predictions of time evolutions in a test set of 2500 sequences conformal to the training dataset for the different NN models. The full distribution is shown in colour by violin plots. The boxes correspond to the interquartile range (IQR). The median value and the [0.1, 0.9] percentile range are also shown.
namics that lead to MAE accumulation beyond 0.2, the prediction given by NN7.5k , although more reasonable, 5 returns a MAE error that peaks two times larger than 5k 0.5k the one of NN0.5k ∆µ . Both NN∆µ and NN∆µ are in line with the MAE errors found in the previous case of Fig. 4. Notably, the NN0.5k ∆µ model reports a MAE error below 7.5k NN5 for a large part of the evolution. Although specific of the selected simulation cases, the analysis referred to Fig. 4 and Fig. 5 is consistent for all simulations composing the test set. Each case indeed results in an error accumulation in the early growth stages, reaching a maximum at intermediate coverage and then going back to zero at later times. This error accumulation becomes critical for NNseq models when considering low supersaturations. To better quantify the overall accuracy of the NN prediction across the entire test set, we now inspect the statistical distribution of the maximum MAE value for each evolution sequence as obtained from any NN model in Table I. It must be noted that this analysis provides the worst-case scenario for the prediction error, as within the time sequence, the discrepancy is generally lower than that. The results are reported in Fig. 6, showing the error distributions as violin plots along with the corresponding box-plot representation tracing the median and interquartile range as well as the [0.1, 0.9] percentile range. As already noted from the loss comparisons in Fig. 3 and by the sequence comparisons in Fig. 4 and Fig. 5, the size of the training dataset is crucial to achieve better accuracy. In the case of the NNseq architecture, only the largest dataset of 7500 sequences guarantees that more than 90% of the predictions yield max MAE errors below an acceptable 0.075 threshold, with a median value of 0.043. By comparing the NNseq models differing just in the mini-
sequence length s, we clearly see that the choice of s = 5 is the optimal one. Indeed, for shorter sequence NN5k 3 the NN performs poorly with a broad error distribution with median error almost double the one of NN5k 5 , probably due to the insufficient information to infer the proper evolution rate. On the other hand, the NN5k 7 model, taking a 2-frame longer initial sequence than NN5k 5 does not show any benefit from it, exhibiting almost the same error distribution at the price of more input frames. NNpar models perform significantly better than the NNseq ones, to the point that the worse of the firsts, NN0.5k ∆µ , is on par with the best among the seconds, i.e. NN7.5k , yielding 90% of the predictions beyond the same 5 0.075 error threshold and a median value of 0.048, despite being trained on a 15 times smaller dataset and corroborating the superior data efficiency of this approach. The best performer among our tests is NN5k ∆µ , yielding more than 90% of the predictions with errors below 0.018 and a median value of 0.011. The supremacy of NNpar was expected as the explicit conditioning by supplying the correct ∆µ value as input strongly simplifies the task with respect to inferring it implicitly from a short mini-sequence, as requested to NNseq 7.5k models. The similarity between NN0.5k in∆µ and NN5 dicates that the performance gap of NNseq can be closed only at the price of a more demanding training on significantly larger datasets. It is also worth to point out that the difference between NNpar and NNseq revealed by this testing was not apparent from the training and validation losses in Fig. 3, thus underlying the importance of independent testing.
C.
NN performances as a function of supersaturation
As suggested by the comparison of Fig. 4 and Fig. 5, the NN models, especially NNseq ones, could behave differently according to the actual value of ∆µ. We here inspect a [0.1, 1.0] range of ∆µ values, slightly exceeding the one of the training set, evaluating the NN extrapolation capabilities too. To this purpose, we prepared new test sets for each value of ∆µ, each composed of 100 sequences. To avoid biases due to the initial random configuration, the same 100 initial profiles are adopted for all test sets. Moreover, we do not fix the number of frames of the test sequence but we let the simulations run until reaching complete filling of the domain to determine the actual maximum of MAE on the complete dynamics. Domain filling is simply evaluated though the coverage θ defined as the average value of φ over the whole domain. The complete filling condition is defined by θ > 0.999. Following this procedure, the length of the test sequences changes from less than 100τ for the highest ∆µ values to more than 800τ for the lowest ∆µ, thus providing in this latter case a further test of time extrapolation. The results of this analysis, for the best-performing models in the two classes, i.e. NN7.5k and NN5k 5 ∆µ , are illustrated in
8
(a)
1
training range
0.5
max(MAE)
Fig. 7, reporting for each sampled ∆µ a box-plot showing the median value, interquartile range and [0.1, 0.9] percentile range of the distribution of MAE maxima (see also Fig. S5 of Supplementary Material for the same analysis on the original test set, for all trained models). 7.5k It can be seen that the NN5k for ∆µ outperforms NN5 any ∆µ, with maximum MAE errors which remains well below 0.04 for 90% of the test cases for all ∆µ values within the training range. A pronounced increase in the error is, however, observed when extrapolating for ∆µ ≤ 0.2 with predictions becoming unreliable right below that threshold. Extrapolation toward higher ∆µ beyond 0.8 is instead possible, even if the error increases sharply, thus limiting the accessible range to ∆µ ≤ 1.0. The behaviour of maximum MAE errors for NN7.5k in5 stead shows a more pronounced dependence on ∆µ also within the training range. The best results are found when approaching the higher bound of ∆µ = 0.8. At intermediate values, the error tends to stabilize below a threshold of 0.08 for 90% of the test cases, similarly to the overall values found from Fig. 6. However, when decreasing ∆µ below 0.3, the error rapidly increases so that a large fraction of the NN predictions within the lowest ∆µ portion of the training range are actually subject to large MAE beyond 0.1. The model also fails in extrapolating beyond the training range as errors quickly increase. By this analysis, we conclude that the NN7.5k 5 should only be trusted for predictions within the range of ∆µ ∈ [0.3, 0.8] despite the training including also lower values. A possible interpretation for this failure at low ∆µ may be related to the fact that such low values imply a very small profile change between the few frames in the mini-sequence, which limits the NN model’s capability of recognizing the actual evolution rate. As from Sect. II A it was noted how a change in ∆µ affects the smoothness of the growth front morphology, beyond a mere rescaling of the growth rate. Assessing the trained models’ ability to reproduce such fine details is the ultimate test of their robustness. In Fig. 7(b-d) we compare the profiles obtained by simulations started from the same initial configuration of Fig. 1 for three different values of ∆µ = 0.25, 0.50 and 0.75 at times 120τ , 60τ and 40τ scaled to return equivalent growth stages. As evident from the PF ground-truth profiles, the morphologies are similar, except for the aforementioned trend toward sharper edges with increasing ∆µ. Both NN7.5k 5 and NN5k ∆µ models return an overall satisfactory reproduction of such profiles. Significant local discrepancies are only evident for NN7.5k at low ∆µ, in line with the 5 previous analysis. At a close inspection it can also be appreciated how the predicted profiles are characterized by rounder corners for low ∆µ while sharper edges are obtained for larger values (see also Fig. S6 of Supplementary Material), indicating that the trained models fully learned the role of the ∆µ parameter in controlling also the finer details in the profile evolutions. It is particularly remarkable that MAE prediction errors for NN5k ∆µ and for NN7.5k at sufficiently large ∆µ are typically lower 5
0.9 perc.
IQR
{ median
0.1 perc.
0.1 0.05
0.01 0.005
0.2
0.4
Δμ
0.6
0.8
δ: -1
1.0 +1
(b) 120 τ =0.25
MAE=0.043
0.008
0.027
0.007
0.011
0.007
(c) 60 τ =0.50
(d) 40 τ =0.75
PF FIG. 7. (a) Dependency of the maximum of MAE distribution on the supersaturation value ∆µ for both NN7.5k and 5 NN5k ∆µ . Each box represents the interquartile range (IQR) of the error distribution on a dataset of 100 sequences starting from the same initial profiles and evolving according to the corresponding supersaturation ∆µ up to complete filling of the domain. The median value and the [0.1, 0.9] percentile range are also shown. The shaded areas represent ∆µ values outside the range of the training set. (b-d) Comparison between true PF and NN-predicted profiles, for an equivalent growth stage from the same initial configuration of Fig. 1 and different ∆µ. The pixel-by-pixel error δ is superimposed on the predicted images by colour bar along with the value of MAE.
than the discrepancies between the different PF profiles at same evolution stages (typical ⟨|δ|⟩ ≈ 0.02 − 0.04), thus providing a further proof that the trained models perform better than rescaling over time the same averaged features.
9
FIG. 8. Evolution sequence for crystal growth under the same conditions of training set but on a 2048 × 2048 domain (∆µ = 0.318) as obtained by models (a) NN7.5k and (b) NN5k 5 ∆µ . Pixel-by-pixel errors δ are superimposed on the evolution frames by colour map. Insets provide a magnification of the 128 × 128 left-bottom corner. (c) Time evolution of MAE or the sequences in (a) and (b). (d) Analysis of the distribution of MAE maxima across 50 evolution sequences for different domain sizes L by both violin plot and box plots showing the median, interquartile range (IQR), and [0.1, 0.9] percentile range.
D.
Generalization to large domains
Both NN approaches discussed here exploit a fullyconvolutional architecture so that the trained models can be, in principle, applied to any domain size, eventually much larger than the one considered for its training. This is particularly appealing in the perspective application of this class of models, since it would allow training on small-size, low-cost simulation domains, while promising the exploitation for larger domains at a fraction of the computational cost that using PF would imply. To prove this potential, we then analysed the NN performances in predicting the evolution sequence for 50 test cases defined with the same conditions of the training dataset but extending on L×L domain sizes of 256×256, 512×512, 1024×1024 and 2048×2048 collocation points. An example of NN predicted profiles on the largest
2048 × 2048 domain is reported in Fig. 8 for the NN7.5k 5 (a) and NN5k ∆µ (b) best models (see Fig. S7 for additional examples on other domain sizes). The simulation corresponds to a value of ∆µ = 0.318 requiring a time of about 300τ to reach complete filling of the domain, i.e., 6 times the duration of the training sequence. A close inspection of the pixel-by-pixel error δ, highlighted in the insets, clearly shows local discrepancies in the NN7.5k predicted 5 profiles, while for NN5k they are barely distinguishable. ∆µ The evaluation of the MAE during the time evolution, reported in Fig. 8(c), shows a progressive error accumulation reaching a maximum around t = 170τ , when most of the largest domains coalesce, and then decreasing to zero as the domain is fully covered. The MAE peak values are as small as 0.035 for NN7.5k and 0.010 for NN5k 5 ∆µ . While the reported case of Fig, 8 can be taken as representative of the test set simulations, a complete quantita-
10
FIG. 9. Evolution sequence for crystal growth starting from a low coverage configuration of θ0 = 0.023 on a 512 × 512 domain, with ∆µ = 0.5, as obtained by the models (a) NN7.5k and (b) NN5k 5 ∆µ . Pixel-wise errors δ are superimposed on the evolution frames by colour map. (c) Time evolution of MAE for the sequences in (a) and (b). (d) Dependence of the MAE maxima distribution as a function of the initial profile coverage. Box-plots show the median, interquartile range (IQR) and [0.1, 0.9] percentile range for 50 evolution sequences with ∆µ = 0.5. The gray shaded area corresponds to the training set coverage range.
tive analysis of the NN accuracy can only be achieved by inspecting the actual distribution of MAE maxima across the test case, as already done for the original 128×128 domain size in Fig. 6. In Fig. 8(d), we report by violin plots and box plots the resulting distributions for the tests on the different domain sizes, still comparing NN7.5k and 5 NN5k ∆µ . We can see that the errors are consistent across all sizes, with the NN7.5k capable of returning 90% of the 5 predictions with maximum MAE errors below 0.06 and the NN5k ∆µ remaining just around 0.01. We then conclude that the NN can reliably be extended to larger domain sizes, ideally approaching the experimental ones, with no loss in accuracy with respect to the original training conditions.
E.
Initial coverage effect
A last key aspect that has been neglected so far in our tests is the role of the different initial distribution of seeds. All the previous analysis indeed applied the same protocol to generate the initial profiles, returning fairly consistent seed distributions that correspond to an initial
coverage θ of 0.11 ± 0.015. However, when considering larger domains, it becomes natural to consider different seed densities, in particular enabling the analysis of more sparse configurations. In this last section, we inspect how the trained NN can possibly generalize with respect to variable initial coverages. In Fig. 9 we consider the time evolution in a 512 × 512 domain as an example, initiated from a seed distribution corresponding to a coverage of about θ0 = 0.023, i.e., about 5 times lower than the ones of training. The supersaturation condition is set to ∆µ = 0.5. The complete evolution up to full domain filling extends for about 500τ , so that the comparison extends 10 times beyond the training sequence duration. In Fig. 9(a), we have the evolution predicted by the NN7.5k model, while in 5 panel (b) the one from NN5k is reported. The pixel-wise ∆µ errors in the predictions are highlighted by the δ colour map, while the evolution of MAE is shown in Fig. 9(c). It can be observed that both models return quite consistent evolutions despite the unusual initial configuration with respect to training conditions. Consistent with the rest of the study, NN5k ∆µ performs better, returning a maximum MAE of just 0.02 while NN7.5k returns larger errors 5
11 at merging domains, reaching a still acceptable MAE of 0.08. To fully characterize the relationship between NN accuracy and initial surface coverage, we generated a series of test sets composed of 50 cases and initialized with prescribed coverage θ0 both below and above the typical value used in the NN training. A fixed ∆µ = 0.5 is here considered. For each θ0 test set, we then analyse the distribution of MAE maxima and report the corresponding median, interquartile range, and [0.1, 0.9] percentile range by box-plots in Fig. 9(d). Once again, the best performing models NN7.5k and NN5k 5 ∆µ are compared. Consistent with the previous analyses, NN5k ∆µ proves to be the most accurate, returning MAE maxima as low as 0.02 for more than 90% of the test cases above a θ0 = 0.03 threshold. Predictions are still reliable for θ0 = 0.02 but fail for lower coverage. Similarly, also NN7.5k performs 5 well for all coverages above the θ0 = 0.03 threshold but becomes unreliable below that. We find that in both cases the NN can reliably predict the evolution of profiles starting from initial coverage as low as 0.03. This holds true for different ∆µ values within the range of low-prediction error (i.e., [0.3, 0.9]; ∆µ = 0.3 and 0.7 are reported in Fig. S8 of Supplementary Material).
IV.
CONCLUSIONS
Two Convolutional Recurrent Neural Network architectures have been analysed in depth, demonstrating their ability to approximate the numerical solution of the Allen-Cahn equation for crystal growth. Notably, the NN surrogates are found to properly recognize the fine effects of variable supersaturation values on the faceted growth morphologies, beyond a mere rescaling of growth rates on averaged features, thus proving their suitability for applications to more realistic and detailed models of specific crystalline materials. As expected, explicit conditioning by supplying the critical parameter as input, i.e. supersaturation in the present case, ensures the best predictions, yielding the lowest mean absolute error in the cases analysed. Our best model (NN5k ∆µ ) is found capable of producing quantitative predictions of evolution sequences on the whole spectrum of tested conditions, with local errors limited to just a few pixels (more than 90% of the tested cases are predicted with less than 1.8% of wrong pixels). Notably, we find that even when training on 10 times smaller dataset (model NN0.5k ∆µ ) the predictions are still generally acceptable (the prediction error grows up to 7.5%), indicating that the approach can still be usable even if training data were scarcely available. On the other hand, the NNseq architecture, inferring the supersaturation implicitly from short sequences, gen-
erally requires larger training datasets to produce consistent levels of prediction accuracy. In particular, our analysis suggests that an increment in the dataset size of at least a factor 15 should be needed to match performances. Despite this higher training cost, this fully datadriven approach could still be the only choice whenever tackling a problem for which the leading parameters are unknown, e.g. dealing with experimental data. Fortunately, the model shows consistent performance already for a mini-sequence as short as 5 frames, limiting the effort needed to produce a suitable input. The trained models were also successfully applied to larger computational domains, with no loss of accuracy. Furthermore, they also provided good generalization capabilities for moderate variations in the initial crystal seed density, and for simulation more than 10 times longer than the training ones, thus making it possible to apply the approach even to more realistic scales. It is worth noting that the numerical solution of the present Allen-Cahn does not pose severe computational bottlenecks. Indeed, the execution times for the finitedifference algorithm are generally faster than the NN evaluation on the same CPU (here, Intel Core i7-12600k). A significant benefit of a factor ≈ 10× is achieved only if performing the NN evaluation on a GPU (NVIDIA RTX 3060). Indeed, in this work, the main goal was not to achieve acceleration, but took advantage of the inexpensive model to perform extensive testing. On the other hand, we expect substantial speed-ups offered by the NN surrogate whenever considering more complex, nonlinear dynamics eventually requiring more advanced integration schemes and more costly techniques, e.g., finite element method.
DATA AVAILABILITY STATEMENT
The datasets used to train and test the model are openly available in Materials Cloud Archive at https://doi.org/10.24435/materialscloud:yv-sy, reference number 2026.85. The code used to train the NN model is freely available on GitHub at https://github.com/dlanzo/CRANE.
ACKNOWLEDGMENTS
FM, RB, and DL acknowledge financial support from ICSC—Centro Nazionale di Ricerca in High-Performance Computing, Big Data and Quantum Computing, funded by the European Union—NextGenerationEU. D.L. acknowledges financial support from ICSC SPOKE 7 CNR, project INNOVATOR, CUP B93C22000620006, CN00000013.
12
[1] J. Wei, X. Chu, X.-Y. Sun, K. Xu, H.-X. Deng, J. Chen, Z. Wei, and M. Lei, Machine learning in materials science, InfoMat 1, 338 (2019). [2] K. Choudhary, B. DeCost, C. Chen, A. Jain, F. Tavazza, R. Cohn, C. W. Park, A. Choudhary, A. Agrawal, S. J. Billinge, et al., Recent advances and applications of deep learning methods in materials science, npj Computational Materials 8, 59 (2022). [3] J. Schmidt, M. R. Marques, S. Botti, and M. A. Marques, Recent advances and applications of machine learning in solid-state materials science, npj computational materials 5, 83 (2019). [4] I. Peivaste, S. Belouettar, F. Mercuri, N. Fantuzzi, H. Dehghani, R. Izadi, H. Ibrahim, J. Lengiewicz, M. Belouettar-Mathis, K. Bendine, et al., Artificial intelligence in materials science and engineering: Current landscape, key challenges, and future trajectories, Composite Structures 372, 119419 (2025). [5] P. Lyngby and K. S. Thygesen, Data-driven discovery of 2D materials by deep generative models, npj Computational Materials 8, 232 (2022). [6] Y. Zhao, E. M. D. Siriwardane, Z. Wu, N. Fu, M. AlFahdi, M. Hu, and J. Hu, Physics guided deep learning for generative design of crystal materials with symmetry constraints, npj Computational Materials 9, 38 (2023). [7] E. T. Chenebuah, M. Nganbe, and A. B. Tchagang, A deep generative modeling architecture for designing lattice-constrained perovskite materials, npj Computational Materials 10, 198 (2024). [8] C. Karpovich, E. Pan, and E. A. Olivetti, Deep reinforcement learning for inverse inorganic materials design, npj Computational Materials 10, 287 (2024). [9] N. Provatas and K. Elder, Phase-Field Methods in Materials Science and Engineering, 1st ed. (Wiley, 2010). [10] B. Li, J. Lowengrub, A. Ratz, and A. Voigt, Geometric Evolution Laws for Thin Crystalline Films: Modeling and Numerics, Communications in Computational Physics 6, 433 (2009). [11] D. Montes De Oca Zapiain, J. A. Stewart, and R. Dingreville, Accelerating phase-field-based microstructure evolution predictions via surrogate models trained by machine learning methods, npj Computational Materials 7, 3 (2021). [12] C. Hu, S. Martin, and R. Dingreville, Accelerating phasefield predictions via recurrent neural networks learning the microstructure evolution in latent space, Computer Methods in Applied Mechanics and Engineering 397, 115128 (2022). [13] V. Oommen, K. Shukla, S. Goswami, R. Dingreville, and G. E. Karniadakis, Learning two-phase microstructure evolution using neural operators and autoencoder architectures, npj Computational Materials 8, 190 (2022). [14] P. Tep and M. Bernacki, High-fidelity grain growth modeling: Leveraging deep learning for fast computations, Acta Materialia 301, 121486 (2025). [15] K. Alhada–Lahbabi, D. Deleruyelle, and B. Gautier, Machine learning surrogate for 3D phase-field modeling of ferroelectric tip-induced electrical switching, npj Computational Materials 10, 197 (2024). [16] C. Bonneville, N. Bieberdorf, P. Robbe, M. Asta, H. Najm, L. Capolungo, and C. Safta, Towards spatio-
temporal extrapolation of phase-field simulations with convolution-only neural networks (2026), 2601.04510. [17] I. Peivaste, A. Makradi, and S. Belouettar, Teaching artificial intelligence to perform rapid, resolution-invariant grain growth modeling via fourier neural operator, Computer Methods in Applied Mechanics and Engineering 440, 117945 (2025). [18] S. Fan, A. L. Hitt, M. Tang, B. Sadigh, and F. Zhou, Accelerate microstructure evolution simulation using graph neural networks with adaptive spatiotemporal resolution, Machine Learning: Science and Technology 5, 025027 (2024). [19] D. Lanzoni, F. Montalenti, and R. Bergamaschini, Deep learning for simulating the evolution of condensed matter systems at the continuum scale: Methods and applications, Journal of Physics: Condensed Matter 37, 403003 (2025). [20] X. Shi, Z. Chen, H. Wang, D.-Y. Yeung, W.-k. Wong, and W.-c. Woo, Convolutional LSTM Network: A Machine Learning Approach for Precipitation Nowcasting (2015), 1506.04214. [21] N. Ballas, L. Yao, C. Pal, and A. Courville, Delving Deeper into Convolutional Networks for Learning Video Representations (2016), arXiv:1511.06432. [22] K. Yang, Y. Cao, Y. Zhang, S. Fan, M. Tang, D. Aberg, B. Sadigh, and F. Zhou, Self-supervised learning and prediction of microstructure evolution with convolutional recurrent neural networks, Patterns 2, 100243 (2021). [23] P. Wu, A. S. Iquebal, and K. Ankit, Emulating microstructural evolution during spinodal decomposition using a tensor decomposed convolutional and recurrent neural network, Computational Materials Science 224, 112187 (2023). [24] D. Lanzoni, A. Fantasia, R. Bergamaschini, O. PierreLouis, and F. Montalenti, Extreme time extrapolation capabilities and thermodynamic consistency of physicsinspired neural networks for the 3D microstructure evolution of materials via Cahn–Hilliard flow, Machine Learning: Science and Technology 5, 045017 (2024). [25] A. Fantasia, D. Lanzoni, N. Di Eugenio, A. Monteleone, R. Bergamaschini, and F. Montalenti, A parametricallyConditioned Deep Learning Surrogate for Coherent Spinodal Decomposition, Advanced Theory and Simulations 9, e02144 (2026). [26] A. A. K. Farizhandi, O. Betancourt, and M. Mamivand, Deep learning approach for chemistry and processing history prediction from materials microstructure, Scientific Reports 12, 4552 (2022). [27] D. Lanzoni, M. Albani, R. Bergamaschini, and F. Montalenti, Morphological evolution via surface diffusion learned by convolutional, recurrent neural networks: Extrapolation and prediction uncertainty, Physical Review Materials 6, 103801 (2022). [28] K. Kutsukake, Review of machine learning applications for crystal growth research, Journal of Crystal Growth 630, 127598 (2024). [29] M. Petkovic, L. Vieira, and N. Dropka, Machine learning in crystal growth: A review of methods, data, and applications, Progress in Crystal Growth and Characterization of Materials 71, 100689 (2025).
13 [30] M. Lu, S. Rao, H. Yue, J. Han, and J. Wang, Recent advances in the application of machine learning to crystal behavior and crystallization process control, Crystal Growth & Design 24, 5374 (2024). [31] Y. Lu, X. Yang, B. Niu, T. Li, Z. Zheng, L. Jia, D. Chen, H. Qi, K. H. L. Zhang, M. Zhu, H. Zhang, and X. Lu, Achievement of high-quality gallium oxide epitaxial growth via machine learning, Advanced Functional Materials 36, e19854 (2026). [32] S. M. Allen and J. W. Cahn, A microscopic theory for antiphase boundary motion and its application to antiphase domain coarsening, Acta Metallurgica 27, 1085 (1979). [33] D. P. Kingma and J. Ba, Adam: A Method for Stochastic Optimization (2017), arXiv:1412.6980. [34] Y. Bengio, J. Louradour, R. Collobert, and J. Weston, Curriculum learning, in Proceedings of the 26th Annual International Conference on Machine Learning, ICML ’09 (2009) pp. 41–48.
SUPPLEMENTARY MATERIAL
Neural surrogates for crystal growth dynamics with variable supersaturation: explicit vs. implicit conditioning Matteo Rigoni,1 Daniele Lanzoni,1, 2 Francesco Montalenti,1 and Roberto Bergamaschini1 1
Dept. of Materials Science, University of Milano-Bicocca, 20125 Milano, Italy 2 Dept. of Physics, University of Genova, 16146 Genova, Italy
(a) NNseq
100
10-1
training validation
10-2
Loss
Loss
10-1
10-3
10-3 10-4
10-5
10-5
0
50
100
150
250
200
10-6
300
training validation
10-2
10-4
10-6
(b) NNpar
100
0
50
Training epoch
100
150
200
250
300
Training epoch
7.5k FIG. S1. Training and validation losses during training for the mini-sequence model NN5k (a) and the parameter 5 vs. NN5 5k 1.5k model NN∆µ vs. NN∆µ (b).
100 training validation
Loss
10-1 10-2 10-3 10-4 10-5 10-6
0
50
100
150
200
250
300
Training epoch FIG. S2. Comparison of the training and validation losses during training for the NNseq mini-sequence model of 3 (NN5k 3 ), 5 5k (NN5k 5 ) and 7 (NN7 ) input frames.
2
FIG. S3. Evolution sequence for the same case of Fig. 4 of the main manuscript (∆µ = 0.6) comparing the ground-truth PF solution (a) with the corresponding predictions from both NNseq (b,c,d) and NNpar (e) models. The pixel-by-pixel error δ is superimposed on the predicted images by colour bar. (f) Time evolution of the MAE between NN predictions and true PF solution.
3
FIG. S4. Evolution sequence for the same case of Fig. 5 of the main manuscript (∆µ = 0.6) comparing the ground-truth PF solution (a) with the corresponding predictions from both NNseq (b,c,d) and NNpar (e) models. The pixel-by-pixel error δ is superimposed on the predicted images by colour bar. (f) Time evolution of the MAE between NN predictions and the true PF solution.
4
max(MAE)
1
0.1
0.01 0.2
0.4
∆μ
0.6
0.8
0.2
0.4
∆μ
0.6
0.8
0.2
0.4
∆μ
0.6
0.8
0.6
0.8
max(MAE)
1
0.9 perc.
0.1 IQR
{ median
0.1 perc.
0.01 0.2
0.4
0.2
0.4
∆μ
0.6
0.8
0.2
0.4
0.6
0.8
0.2
0.4
∆μ
0.6
0.8
0.6
0.8
max(MAE)
1
0.1
0.01
∆μ
∆μ
0.2
0.4
∆μ
FIG. S5. Analysis of the distribution of MAE maxima on the test set of 2500 sequences used for Fig. 5 of the main manuscript, for all NNseq and NNpar trained models. For each of them, we report the overall distribution (first column) and those obtained by partitioning the test-set cases according to ∆µ (binning width of 0.05). Both violin plots and box-plots, reporting median, interquartile range (IQR) and the [0.1, 0.9] percentile range, are shown.
5
seq 1 120 τ
MAE=0.043
seq 2 0.008
MAE=0.026
seq 3 0.009
=0.25
60 τ
0.027
0.007
0.018
0.011
0.007
0.008
0.007
0.014
0.005
0.009
0.007
0.006
=0.50
40 τ
MAE=0.033
0.007
=0.75
PF FIG. S6. Comparison between true PF and NN-predicted profiles, for equivalent growth stage for simulations started from 3 different initial configurations (seq 1 is the same as Fig. 7(b-d) of the main manuscript) and different ∆µ. Reported profiles correspond to the φ = 0.5 contour lines. The MAE prediction error for each case is also indicated.
6
FIG. S7. Evolution sequences for crystal growth under the same conditions of the training set but on different domain sizes: (a) 256 × 256 (∆µ ≈ 0.58), (b) 512 × 512 (∆µ ≈ 0.23) and (c) 1024 × 1024 (∆µ ≈ 0.31) as obtained by models NN7.5k and 5 NN5k ∆µ . Pixel-by-pixel errors δ are superimposed on the evolution frames by color map. The corresponding evolutions of MAE over time are plotted.
7
(a) ∆μ=0.3
(b) ∆μ=0.6
0.9 perc.
0.5
IQR
{ median
0.5
max(MAE)
max(MAE)
0.1 perc. 0.1 0.05
0.01 0.005
0
0.1 0.05
0.01 0.05
0.1
0.15
initial coverage θ0
0.2
0.25
0.005
0
0.05
0.1
0.15
0.2
0.25
initial coverage θ0
FIG. S8. Dependence of the MAE maxima distribution as a function of the initial profile coverage θ0 for the same test set of 50 evolution sequences used in Fig. 9(e) of main manuscript but with (a) ∆µ = 0.3 and (b) ∆µ = 0.7 Box-plots show the median, interquartile range (IQR) and [0.1, 0.9] percentile range. The gray shaded area corresponds to the training set coverage range.