A convolutional neural network is trained on a large dataset of suitably randomized film profiles and corresponding elastic energy densities , computed by the finite element method. The trained model provides quantitative predictions of for arbitrary profiles, surrogating its explicit calculation, and is used for the time integration of partial differential equations describing the evolution of strained films. The close match found between the neural network predictions and the “ground-truth” evolutions obtained by the finite element method calculation of , even after tens-of-thousands of integration time-steps, validates the approach. A substantial computational speed up without significant loss of accuracy is demonstrated, allowing for million-steps simulations of islands growth and coarsening. The intriguing possibility of extending the domain size is also discussed.
I. INTRODUCTION
In the last decade, machine learning (ML) techniques1–3 have been introduced in almost any field of science, from purely theoretical to applied disciplines and engineering.4–9 Among the many uses, these methods proved to be an effective way to speed up otherwise computationally heavy calculations. In Condensed Matter Physics and Materials Science, considerable effort has been dedicated both to the acceleration of atomic-scale models10–12 and, more recently, to the approximation of continuum models based on partial differential equations.13–18 A few recent studies even proposed novel strategies to bridge the gap between the two paradigms by exploiting ML.19,20
Here we show that neural networks (NN)2 can be conveniently exploited to accelerate continuum simulations of strained-film growth (as encountered, e.g., in lattice-mismatched heteroepitaxy21), a key field in Materials Science.22,23 Indeed, a wide range of behaviors depending on growth parameters can be observed and exploited in various applications. In some cases, controlling the surface roughness that may arise during the growth is the main issue.24,25 Other applications are instead based on the formation of three-dimensional (3D) nanostructures, such as islands or low-dimensional objects,26–29 yielding novel properties of high interest. Understanding the interplay between growth parameters and material properties is thus key and the development of reliable and fast simulation tools to support and anticipate costly experiments would offer great benefits. As atomistic approaches are generally too expensive for simulating the film evolution on realistic spatial and, most critically, temporal scales, continuum modeling becomes fundamental.30
The evaluation of requires in general the solution of the mechanical equilibrium problem by numerical approaches. Among these, the Finite Element Method (FEM) is particularly well suited for generic free-boundary problems, although efficient alternatives are available for more specific cases.40 While single static calculations are not an issue for nowadays computational resources, accurate FEM solutions still pose severe limitations when time-integrating equations such as Eq. (1), as needs to be continuously updated, even for millions of iterations, thus multiplying the computational burden on simulation times. A few studies in the literature pursued this strategy by using advanced numerical approaches with optimized time-integration schemes41–44 or resorted to computationally-cheap semi-analytical methods,37,45 at the price of an approximated solution for the elastic problem.
In this work, we explore the possibility of by-passing the explicit solution of the elastic problem using a convolutional NN (CNN)2,46–48 model to directly predict the contribution. CNNs involve numerical operations that are much cheaper than those required for a FEM solver and can be efficiently implemented on multi-threaded systems and GPUs, allowing for the approximation of complex, non-linear functions at a fraction of the computational cost of conventional approaches. Here, we train a CNN capable of providing a mapping of arbitrary surface profiles to their corresponding with the same accuracy as FEM. After validating the model and inspecting its extrapolation performances, we exploit it in the time integration of the thin-film evolution equations. The substantial speed-up offered by this deep learning approach enables the analysis of evolution regimes, which are usually expensive to tackle.
We remark that, in principle, an even stronger computational-time compression might be obtained by predicting the full driving force μ or the profile sequence directly.16,17 The present choice of focusing only on the term that causes the major computational bottleneck, however, ensures full transferability, while retaining a better physical insight beyond the “blind” prediction of the required quantity. Our approach is also complementary to other ML methods devised to tackle differential equation problems. For example, the Physics Informed NN (PINN) method developed in Ref. 14 (and recently applied to linear elasticity18) could potentially substitute FEM solvers for the mechanical equilibrium problem in the whole domain. In this respect, the NN approach here proposed represents a data-driven alternative that does not require re-training for new profiles and completely bypasses the solution step by directly predicting at the film surface.
II. METHODS
A. Strain calculation and dataset definition
We describe the film morphology by a 1D profile function y = h(x), evolving according to Eq. (1). With no loss of generality for the workflow, we refer to the prototypical case of a Ge film subject to a purely diagonal 3.99% eigenstrain ɛ*49,50 as the one produced by the (compressive) lattice mismatch in Ge/Si(001) epitaxy.26,51
A domain size of L = 100 nm is considered since the characteristic length scale of morphology is expected to be of the order of a few tens of nanometers, as the fastest growing mode has wavelength λ ≈ 23 nm,32,33,35 given the aforementioned parameters.
The NN dataset is composed of ∼80 000 [h(x), ρɛ(x)] couples, which can be found in Ref. 54, and has been split at random into training (85%) and validation (15%) sets. Given the non-linear relationship between h and , two types of profiles are considered: (1) pure sine-waves with wavelengths λ ∈ {L/6, L/5, L/3, L/2, L} and amplitudes in [10−3, 10] nm, representative of the Fourier components of a generic profile; (2) arbitrary undulated morphologies set by the kind of gradient noise known as Perlin noise55,56 with amplitudes in the [5 × 10−3, 15.4] nm range. An example of the latter is reported in Fig. 1.
B. Neural network architecture
A fully convolutional NN1,2,57 implemented within the PyTorch framework58 has been used to map the surface profile h(x) to the corresponding elastic energy density ρɛ(x). The term “Fully Convolutional” refers to the fact that the NN comprises only convolutions and element-wise functions. As a result, the input of the network is not restricted in length, thus allowing for the generalization to computational domains of arbitrary size.
This architecture incorporates translational equivariance2 by construction, which increases NN accuracy and reduces the number of examples required in training. Other symmetries can also be straightforwardly added by ad-hoc manipulations of the NN parameters and structure. Mirror symmetry is enforced by making convolution kernels reflection-invariant by construction and through data augmentation. While in principle only one of these strategies should be sufficient, models trained using both obtained a lower generalization error in our tests. The mean value of h(x), , is also removed before passing a profile as input to the NN to enforce the vertical shift invariance proper of the considered system. The CNN adopted circular padding,59,60 enforcing PBCs exactly and consistently with the FEM calculation scheme.
In terms of structure, the NN is composed of stacked, one-dimensional convolutional layers followed by tanh activation functions. Each layer contains 20 different kernels with a receptive field, that is, the spatial extent of the convolutional filter, of 21 “pixels,” each corresponding to a collocation point used to discretize the [h(x), ρɛ(x)] pairs. A constant spacing of dx = 1 nm has been chosen to ensure accurate tracking of the surface dynamics while still allowing for large enough time-steps. A total of 5 convolution-tanh blocks are used, with a nominal total of ≈26 000 parameters. Thanks to reflection invariance, however, their effective number is reduced.
III. RESULTS AND DISCUSSION
A. Neural network training and validation
Figure 2 reports the training and validation losses. These values steadily decrease during the 6000 training epochs, without any sign of overfitting. A comparison between the NN prediction and the FEM for the same validation Perlin-noise profile for Fig. 1 is reported in Fig. 3(a).
An additional proof of the reliability of the trained NN model is given in Figs. 3(b)–3(d), where we compare the NN prediction against the FEM solution for three profiles, qualitatively different from those present in the dataset (additional examples can be found in Ref. 54). The accurate reproduction provided by the NN asserts its generalization capabilities, as the NN was never exposed to such functional forms during training. Noticeably, the main discrepancies are present at the singular points corresponding to the sharp angles of the profile in panel (d), clearly indicating that faceted profiles exceed the extrapolation limits of the trained model.
B. Validation of time evolutions based on the NN approximation of ρɛ
The NN-predicted elastic energy densities can be fruitfully exploited for the numerical integration of the surface evolution dynamics stated in Eq. (1) at a reduced computational cost. In particular, an orders of magnitude speed-up with respect to FEM is found. An explicit-Euler scheme is implemented, with a fixed time-step δt. Given the present purpose of demonstrating the stability of the NN approach during time integration, we report the simulation time t by the number of time-steps. The absolute timescale can be recovered by multiplying the reported numbers by Mδt, here set to 5 × 10−3 in all simulations. Since a typical evolution sequence occurs over hundreds of thousands of time-steps, one-to-one comparisons with the ground-truth provided by iterating the explicit FEM calculation of over time are possible only for small systems and limited time intervals. The NN approximation capabilities have been further confirmed in a quantitative study limited to small-slope profiles for which semi-analytical expressions of the strain field are available.62
As a first example, we consider the simplest case of annealing (f = 0) of a thick film with isotropic surface energy density γ, such that μ = κγ + ρɛ (κ being the profile curvature). In this case, any perturbation of the film exceeding a certain critical wavelength is expected to grow exponentially according to the well-known ATG-instability model.32–35
Figure 4(a) reports representative stages of the growth of a sinusoidal perturbation, setting an unstable wavelength of 25 nm that was on purpose excluded from the training set, so as to return a fair test on the NN predictive capabilities. A semi-log plot of the peak-to-valley profile height Δh as a function of time is also reported in panel (b). The linear behavior confirms the exponential growth regime predicted by the analytic ATG theory. In this specific case, a one-to-one comparison with the FEM evolution is provided. Notice, however, that the number of FEM evaluations required for tracking the dynamics is approximately equal to the size of the whole dataset used to train the NN, making it inconvenient to extensively repeat a similar analysis. A substantial overlap is obtained over the tens of thousands of integration steps of the simulation. Local discrepancies around the minima and maxima are recognizable only in the last evolution stages, when non-linear contributions induce the formation of cusp-points in between the peaks, leading to numerical divergence. To quantify the discrepancy between the NN predicted profiles and the true ones obtained by means of FEM, we resort to the normalized RMSE σh, defined as in Eq. (5) but now evaluated for the height profiles hNN(x, t) and hFEM(x, t). In panel (c), the σh error is reported as a function of time. As it can be observed, σh progressively increases, up to reaching the threshold value at the end of the reported simulation time, that is, close to the divergence [see also the curve at t = 8 × 104 in Fig. 4(a)]. This increment is not due to a failure in the NN prediction of as made evident by the constant, below-threshold σρ error, evaluated at representative evolution stages. It is instead the result of small-error accumulation in time integration, exacerbated by the exponential nature of the process considered.
A more appropriate description of the film growth dynamics, which also heals the divergent behavior of the ATG model allowing for the study of island formation and coarsening, can be achieved by introducing the substrate–film wetting interaction, as discussed in Sec. II A. Since the ML approach only applies to the elastic contribution, the same NN model can be used despite this modification. In Fig. 5, we report an example of a simulation of island growth, including the wetting energy term and a small deposition flux of f/M = 2 × 10−5, such to follow quasi-equilibrium conditions. Here and in the following figures, a 5× scaling is applied to the height axis for the sake of readability. As we want to provide a one-to-one comparison with the full FEM solution, at least for part of the dynamics, we directly set the initial profile thickness above the critical value hc, with an arbitrary modulation by random Fourier components with λ ∈ [15, 100] nm to trigger unstable modes. As expected, in the early stages of the reported evolution, the corrugation amplifies similarly to the ATG model while, in the latest stages, that is, after 1 × 106 steps, stable islands grow and coarsen, resembling Stranski–Krastanov30 behavior. The corresponding true evolution by FEM was carried on for the first 40 000 integration steps (still corresponding to half of the FEM evaluations required for generating the dataset). Again, close correspondence between the NN and FEM evolutions is achieved, as made evident by the overlap of profiles at both 20 000 and 40 000 steps in Fig. 5. Correspondingly, the σh error remains well below the threshold over the whole tested interval, certifying the accuracy of the NN approach. While the one-to-one comparison with FEM simulations cannot be extended to the last stage at 1 × 106 steps, it must be noted that, in contrast to the divergent behavior of the ATG case, the dynamics of island growth is quenched by the wetting interactions so that the error accumulation is expected to have a lesser impact on the long-time stages, preserving the reliability of the predicted profiles. As an indirect confirmation, in Fig. 5, we also compare the NN prediction of with its FEM counterpart for the three reported evolution stages. Evidently, the NN provides accurate evaluations of , with low σρ error over the whole time interval, even after 1 × 106 steps.
C. Growth simulations and generalization to large domains and long times
After having demonstrated the reliability of NN evolutions for islanding dynamics, we can use the trained model to produce more realistic simulations of island growth, starting from the bare substrate and depositing Ge for several nms. In particular, in Fig. 6, we report the evolution sequence for three different f/M ratios (for full animations, see supplementary material). In these cases, a random-Fourier-components-noise of amplitude 10−2 nm is reset at every step as long as the profile remains flat, so to trigger any unstable mode. As expected, for high f values [Fig. 6(a)], the film grows flat since the time scale of the instability is much slower than the one of deposition. By raising the deposition flux (b), the corrugation formation becomes concurrent with growth, yielding a final morphology with several island-like features. Finally, in the more thermodynamic, low-flux condition (c), the unstable modes develop faster than deposition and coarsen into a single island. Due to the constant supply of material by deposition, the morphology of the final island continuously adapts toward the minimum of energy, thus resulting in a sequence of quasi-equilibrium configurations. We remark that this last simulation required 30 × 106 integration steps, giving impressive proof of the stability of the proposed method. Also, it exceeds the training set size by two orders of magnitude, thus single-handedly justifying the computational effort of the dataset building.
Studying the equilibrium morphology of heteroepitaxial islands is relevant for assessing the long-time limits of simulations and provides further proof of the NN capability of returning meaningful outputs. This task can be conveniently treated as a free energy minimization problem constrained to the assigned amount of material. Since μ = δF/δh, a steepest-descent scheme can be exploited, profiting from the cheap estimation of the elastic term provided by the NN model. As initial conditions, we consider flat profiles with a Gaussian dimple to trigger the reshaping. Different film thicknesses h0 are then considered to inspect the change in the equilibrium morphology as a function of the amount of material. The minimization procedure is halted when the maximum deviation of μ(x) from its mean value is lower than 10−4 eV/nm3. As shown in Fig. 7(a), for small enough h0, due to the stabilization effect of the wetting contribution, the equilibrium configuration is the flat film. Vice versa, when considering thicker initial profiles, an island morphology yields the minimum energy. The largest island in the plot is also shown to overlap with the one obtained for the same amount of material from the surface diffusion evolution of Fig. 6(c), thus giving further proof of consistency of the NN predictions. To confirm that the computed profiles have a quantitative foundation, it is sufficient to verify the match of the NN evaluation of the elastic contribution with the corresponding FEM solution, as the other terms in μ are identical by construction. In Fig. 7(b), we show such a comparison for the profiles of panel (a). In all cases, one-to-one correspondence can be observed, and the σρ values are well below the 0.02 threshold. Once again, we remark that performing the same optimization via FEM would be a computationally demanding task, given the number of iterations required (≳10 000).
One of the key advantages offered by the fully convolutional architecture is that it is not restricted to the input size used in training. In Fig. 8, we take profit of this for inspecting the coarsening dynamics of islands formed on a 1-μm-long, that is, 10 times larger than the training examples, Ge film surface (see the video in the supplementary material). Even larger dimensions have been tested, up to 1 mm, but are not shown here for graphical reasons and can be found in Ref. 54. Given the large domain and the number of iterations, up to 7 × 106 steps, comparisons with explicit FEM calculations of are only possible at a few representative stages. The NN approach returns remarkable accuracy, with below-threshold σρ values at all stages. It must be however pointed out that these generalization capabilities are limited to configurations in which the profile is composed of features separated by distances comparable to the training set periodicity. Indeed, configurations exceeding such constraints are not expected to be reproduced using a dataset with a fixed size, as the trained model lacks any knowledge of the intrinsically long-range effects of the elastic field. A more detailed investigation of this effect is discussed in Ref. 62.
IV. CONCLUSIONS
The present study demonstrates how a Convolutional Neural Network can be fruitfully applied to learn one of the main driving forces of strained film growth, that is, the surface elastic energy density. A speed-up of ∼4 orders of magnitude is found with respect to using the same FEM code exploited for the training, while retaining the same level of accuracy, with typical discrepancies of ∼1% or less. Notably, such a gain is entirely ascribed to the calculation of irrespective of the adopted time-integration scheme. Here, the basic forward-Euler method has been implemented to ensure a fair comparison. However, the proposed approach could be profitably integrated into semi-implicit schemes, for example, leveraging automatic differentiation techniques, even when this would be non-trivial with conventional methods. This promises to push even further the spatial and temporal scales accessible to numeric investigation.
As the computational cost of a NN comes from its architecture and not from the number of examples in training, future efforts on building bigger FEM datasets (e.g., including faceted morphologies) promise to push the accuracy of the method even further while keeping evaluation computational costs low. Clearly, one major limitation of the current work is the simplifying assumption of 2D modeling. Besides the obviously larger computational effort, generalization to 3D description should be straightforward, and it is here left for future work.
While the present study focused on the prediction of the elastic energy density for simulating the growth of thin, strained Ge films, applications are not limited to it. The same approach could be applied to predict other quantities, for example, the displacement and stress/strain fields or other surface properties. Another interesting perspective concerns the possibility of tackling different heteroepitaxial systems. Due to the non-linear dependence of on elastic constants, in principle, a new dataset and model should be set for each new material. The only exception would be when considering a variation in the eigenstrain only, which is a reasonable approximation to describe SiGe/Si alloys. In this case, thanks to the quadratic dependence of on ɛ* a simple rescaling is sufficient. Nonetheless, the current model could still be used as a starting point in a transfer-learning procedure. This is expected to reduce the number of examples required for the training when considering systems affine to Ge/Si. A more challenging possibility would be to condition the NN output on Y and ν explicitly and build an ad-hoc dataset exploring different elastic constant combinations, so as to build a comprehensive NN model for the of generic heteroepitaxial films.
Finally, the same paradigm presented here could also be used for tackling more complex physical effects (e.g., including surface and elastic anisotropy or plastic relaxation by dislocations) or could even be transferred to other dynamical problems, thus broadening the class of simulations that could be accelerated.
SUPPLEMENTARY MATERIAL
Videos of the complete profile evolution for the growth simulations in Fig. 6 are provided in the supplementary material as Video_Growth_High_Rate.avi for the fastest deposition rate of panel (a), Video_Growth_Intermediate_Rate.avi for the intermediate one of panel (b), and Video_Growth_Low_Rate.avi for the slowest one of panel (c). The video of the complete annealing simulation in Fig. 8 is also made available in the supplementary material as Video_Annealing_L1000.avi.
ACKNOWLEDGMENTS
R.B. and F.M. acknowledge financial support from the ICSC–Centro Nazionale di Ricerca in High Performance Computing, Big Data and Quantum Computing, funded by the European Union – NextGenerationEU. L.M. acknowledges the support offered by MCIN/AEI/10.13039/501100011033 under Project No. PID2020-115118GB-I00.
AUTHOR DECLARATIONS
Conflict of Interest
The authors have no conflicts to disclose.
Author Contributions
Daniele Lanzoni: Conceptualization (supporting); Investigation (equal); Software (lead); Writing – original draft (lead). Fabrizio Rovaris: Conceptualization (supporting); Investigation (equal); Software (equal); Writing – review & editing (equal). Luis Martín-Encinar: Investigation (equal); Software (supporting); Writing – original draft (equal). Andrea Fantasia: Investigation (supporting); Software (supporting). Roberto Bergamaschini: Conceptualization (supporting); Investigation (supporting); Supervision (equal); Writing – original draft (lead). Francesco Montalenti: Conceptualization (lead); Funding acquisition (lead); Supervision (equal); Writing – review & editing (equal).
DATA AVAILABILITY
The datasets generated and/or analyzed during the current study are openly available in Materials Cloud at https://doi.org/10.24435/materialscloud:5r-9j, Ref. 54. The underlying code for this study is available in GitHub at https://github.com/mosegroup/nn4surf.