Towards “on-demand” van der Waals epitaxy with adaptive ensemble sampling atomistic workflows

Abstract

Traditional approaches to achieving targeted epitaxial growth involve exploring a vast parameter space of thermodynamic and kinetic drivers (e.g., temperature, pressure, chemical potential, etc.). This tedious and time-consuming approach becomes particularly cumbersome to accelerate synthesis and characterization of novel materials with complex dependencies on the local chemical environment, temperature and lattice strains, specifically for nanoscale heterostructures of layered 2D materials. We combine the strengths of next generation supercomputers at the extreme scale, machine learning and classical molecular dynamics simulations within an adaptive real time closed-loop virtual environment steered by Bayesian algorithms to enable asynchronous ensemble sampling of the synthesis space and apply it to the recrystallization of amorphous transition-metal dichalcogenide (TMDC) bilayers to form stacked moiré superstructures under various growth parameters. We show that such batch parallel Bayesian optimization-based online ensemble sampling frameworks for materials simulations can be promising for achieving and accelerating the on-demand epitaxy of van der Waals stacked moiré devices, paving the way towards a robust autonomous materials synthesis pipeline to enable discovery of unprecedented functionalities.

---

Introduction

Accelerating conventional synthesis efforts in van der Waals (vdW) heterostructures is an important aspect of discovering novel pathways that can harness unprecedented functionalities in 2D materials.¹ The recent surge in efforts towards realizing autonomous synthesis^2,3 and characterization^4,5 frameworks is expediting such discoveries. However, narrowing down the optimal thermodynamic and kinetic parameter space is a nontrivial automation challenge, especially when gaining finer control over tuning these emerging functionalities on demand is the ultimate goal. For instance, vertically stacked van der Waals (vdW) heterostructures⁶ with optimal interlayer twists are extremely promising for manifesting tunable and novel functionalities at the nanoscale. These include emerging quantum phenomena, e.g., correlated electronic phases,⁷ unconventional superconductivity,⁸ interlayer excitons,⁹ etc., which can enable unprecedented optical and electronic properties in nanoscale devices. Traditionally, top-down methods such as mechanical manipulations of as-grown epitaxial layers (via release, transfer, stamping,¹⁰ etc.) have been a preferred mode to produce heterostructures with tunable geometries. However, achieving fine control over interlayer orientation as well as ensuring high quality lattice and desired functionality at the wafer scale can be challenging,^1,11 especially in the context of automated platforms.

To overcome the bottlenecks associated with conventional top-down approaches, alternative synthesis routes often exploit thermodynamically and kinetically driven pathways^4,12 of rotational self-organization^13,14 in vdW epitaxial systems.^4,15–17 It has been indeed demonstrated that bottom-up growth strategies can be leveraged to spontaneously evolve a range of rotationally aligned homo- and heterostructures by tuning the driving forces, e.g., substrate interactions,^13,14 epitaxial strain,^18,19 electron irradiation^5,20 and/or thermal annealing.^21,22 Harnessing such growth drivers either individually or combinatorially can therefore pave the way towards scalable synthesis recipes for vdW layered structures with targeted twists or moiré patterns.

However, the efficacy of such bottom-up techniques to design heterostructures with desired functionality often hinges upon the ability to navigate the complex and intertwined parameter space of synthesizing vdW layered structures.²³ In addition to nucleation and growth, the prominence of solid-state phase transitions²⁴ and out-of-equilibrium pathways^25,26 involving the effects of defects, precursors and metastable phases¹⁴ only adds to the complexity of the design space. To this end, the recent rise of high-throughput experimentation techniques^27,28 integrated in a closed-loop fashion with advanced artificial intelligence (AI)/machine learning (ML) platforms has been shown to be promising in rapid and autonomous exploration² of vast synthesis–structure–property spaces.³

While active learning (AL) driven robotic discovery platforms^3,29,30 are in general tailored to facilitate synthesis aided by autonomous characterization of target materials through optimal design of experiments, several challenges³¹ still remain in achieving accelerated discovery of “new” functional materials. Many of these campaigns can benefit from deeper and accessible physical insights through the robust integration of theory, modeling and simulation capabilities^32,33 in the autonomous experimentation (AE) loop. Moreover, it has recently been demonstrated how virtual environments³⁴ can serve as useful testbeds for rapid refinement and systematic benchmarking of algorithms and rigorous uncertainty quantification shielding the prohibitively expensive experimental data acquisition loops.

In the context of accelerating the epitaxial growth of desired twisted vdW layered structures, existing self-driving thin-film synthesis platforms³⁵ could naturally exploit predictive insights from atomistic simulation methods. In the past, purely simulation driven investigations^36–39 have been able to propose novel strategies to design twisted metastable states in graphene and transition metal dichalcogenide (TMDC) layers by means of lattice strain and thermally activated self-assembly pathways, ultimately predicting and motivating subsequent experimental efforts.^12,22,40

To bridge and guide autonomous nanoscale experimental (AE) platforms with insightful molecular dynamics simulations, Bayesian active learning workflows⁴¹ exploiting the asynchronous throughput accessible at the exascale can act as a catalyst towards further acceleration in autonomous materials discovery. Moreover, the recent advances⁴² in high performance computing (HPC) have pushed the boundaries of state-of-the-art extreme scale computational frameworks,⁴³ yielding unprecedented rates of floating point operations.

However, current ML and UQ-driven workflows for atomistics, especially dynamical trajectory datasets, have not yet fully and efficiently leveraged the extreme scale computational resources and are mostly limited in their current role to serve as initial motivations and/or offline knowledge bases to implicitly kick start an AE platform. A few key challenges associated with these efforts include (a) orchestration and scheduling bottlenecks to handle a large volume of resources (node-hours), (b) navigating asynchronous on-line inference frameworks along with (c) on-the-fly data reduction—all three of which are paramount to achieving theory-experiment integrated closed-loop hpc-enabled AI workflows incorporating predictive materials theory and simulation tools in automated experiments to make them truly autonomous.

Inspired by recent experimental efforts to recrystallize MoS2 crystals from amorphous precursors,⁴⁴ here we computationally sample optimal growth parameters to recrystallize (i.e. to evolve and reorganize long range order) a variety of target moiré heterostructures on-demand from a partially disordered amorphous system. We achieve this by demonstrating a resource-adaptive, on-the-fly decision making and parallelizable asynchronous Bayesian ensemble learning framework for dynamical atomistic simulations. We apply active learning strategies using both conventional acquisition function (e.g. expected improvement) and direct posterior sampling (Thompson sampling) based Bayesian optimization (BO) and study the acceleration for asynchronous batch parallel acquisition in our on-the-fly ensemble molecular dynamics workflow environment. Excellent advantages over parallel random sampling as well as efficacy in identifying optimal growth space for evolving “on-demand” target structural phases are then shown.

Results

A real-time hpc-driven asynchronous and autonomous workflow for accelerated computational synthesis

Our hpc resource-driven workflow to enable accelerated computational autonomous synthesis is demonstrated in Fig. 1. Choosing recrystallization of an amorphous molybdenum disulfide (MoS₂) layer on top of a monolayer crystalline substrate (Fig. 1a) as a case study (amorph + crystal), we aim to predict the parameter space (Fig. 1c), which could transform into crystalline bi-layers with target interlayer twist angles (Fig. 1e). We hypothesize that for a given stoichiometry, the parameter space for the recrystallization process involves annealing temperature (T), in-plane lattice mismatch strain (ε_lat, assuming it to be tensile for the thin film) and shear (ε_sh), and we could identify the appropriate growth regime leading to a specific target moiré structure (Fig. 2). The choice of these parameters is motivated by the experimentally reported investigations of a variety of transformation and reconstruction pathways involving van der Waals layered epitaxial phases.^12,24,26,45 While several other kinetic and thermodynamic driving forces, e.g., chemical potential, electron beam effects, effect of edge crystallography and pressure, etc., will play crucial roles in the emergence of long-range order during recrystallization and in dictating the final orientation of phases, for practicality, we restrict ourselves to implicitly incorporating such effects through annealing temperature, hydrostatic (e.g. uniform lattice mismatch strains) and symmetry distortions i.e. shear components in the mechanical strain tensor to demonstrate the potential of such computational investigations in autonomous discovery campaigns. We note that by considering this parameter space, we only demonstrate a representative workflow, where a diverse set of local environments resulting from temperature and strain effects are covered during the amorphous to crystalline phase transformation dynamics; further consideration of drivers like chemical potential and pressure could also be similarly used to explore the structural evolution pathways relevant for typical growth conditions within our current framework.

Fig. 1 Asynchronous active batch sampling strategy for optimal driving parameters for target twist angles. The goal is to transform a given initial structure with an amorphous MoS₂ layer on top of a crystalline layer (amorph + Xtal), into a crystallized twisted bi-layer with a predefined target (a). Asynchronous ensemble optimization (b) strategies using Bayesian acquisition of candidate recrystallization parameters (c). On-the-fly concurrent evaluations of atomistic simulation via a resource-adaptive dynamic orchestration under a parameter space consisting of annealing temperature, lattice and shear strains. Additionally this paramter space for recrystallization can be expanded (as indicated in smaller fonts) to consider other synthesis relevant parameters like chemical potential, pressure, ion/species flux etc. can naturally be (d) to identify the driving parameters ({p_i}) to achieve recrystallization with the target twist angle (e).

Fig. 2 Recrystallization under temperature and strain effects. As temperature is quenched near 0 K (a and b), overall order starts to emerge in the top layer as shown by both the mean-per-atom q₆ (purple) and pair-count based order parameters (grey) in the bottom panel (c and d) for different snapshots in the MD trajectory. While only temperature driven (i.e., heating at 2500 K) trajectories lead to almost no significant interlayer twist—AB stacked 2H phase (c(iii))—applying a lattice strain (ε_lat) with a supercell shear of (ε_sh) 1.5% each can induce a recrystallized bilayer with a twist angle (θ) of around 8.7° (d(iii)).

In principle, this involves exploration of materials dynamics over a nontrivial landscape sensitive to small perturbations in temperature and strain. Accelerating such a search will require adaptive sampling techniques. We take a multi-level iterative acceleration approach by adopting (i) Bayesian optimization (BO) to predict next simulation parameters along with (ii) an on-the-fly asynchronous high-throughput sampling orchestrated on a leadership class single batch computational allocation (cf. Methods) bypassing scheduling bottlenecks. Such workflows, hence, ensure an integrated gain in terms of optimal resource usage and computational throughput through uncertainty-aware dynamic task execution based sampling strategies.⁴⁶

Active learning of molecular dynamics (MD) trajectories performed with a candidate parameter combination set (p_i ≡ {T_i, ε_{lat i}, ε_{sh i}}, where i ∈ {1, 2, …, n_batch}) also requires automated parsing through the large volume simulation output as well as analytics related to the quantity of interest (q.o.i), which in this current case is the interlayer twist angle. Our solution to address this problem is to analyze simulation data at a regular time interval on-the-fly as soon as the data become available through an in-memory message passing and data-parallel processing approach (cf. Methods). Hyperparameters of this adaptive workflow, e.g., batch size of candidate parameters (n_batch) to be evaluated with MD and the available allocation, are correlated (for instance, with the allocation required for MD evaluation) and are managed by adaptive scheduling so as to maintain a uniform number of tasks running during most of the wall clock limit. The capability (Fig. 1) to combine on-the-fly analytics of individual simulations with a concurrent asynchronous evaluation framework makes it possible to expand the global search horizon while exploiting uncertainty of the data through BO.

Tunable twisted MoS₂ bilayers through recrystallization

We start from an initially disordered top layer (amorph-MoS₂) generated with melt-quench MD runs (cf. Methods) placed on top of a 2H crystalline layer (Fig. 2c and d). Upon rapid heating (e.g., T = 2500 K in Fig. 2) followed by a slow quench until ∼0 K, a recrystallized MoS₂ top layer forms an AB stacked 2H bi-layer phase in Fig. 2a with gradually emerging order (cf. details in the Methods on order quantification using a pair-distance count based simplistic order parameter (gray lines in Fig. 2) which compares well with a conventional hexatic symmetry identifier like q₆ shown in purple in Fig. 2). The final ordered structure is drastically altered as 1.5% strain in terms of both the lattice strain (ε_lat) and shear (ε_sh) is induced (Fig. 2b). A moiré bilayer with a final interlayer twist angle of 8.7° with a supercell wavelength of 20.8 Å appears as a result of this change in the recrystallization conditions. This indicates that there exists a promising approach to bias recrystallization pathways towards tunable twisted bi-layers as final phases with the help of controlling temperature and strain induced parameters. Based on such observations as well as guided by past computational investigations of melt-quench crystallization analysis⁴⁷ of MoS₂ involving reactive force fields (cf. modeling details), we pick a range of 1500–3000 K for annealing temperature and lattice strains up to 5% for our autonomous search space. We note that using such approaches to control solid–solid phase transitions via self-organization of disordered 2D materials can emerge as cleaner and tunable alternatives to traditional CVD driven growth of twisted bilayers,⁴⁸ as they involve complex interplay of nucleation and growth kinetics dominated by chemically sensitive driving factors.

Bayesian batch acquisitions to sample a target twist

Given that we have a tunable epitaxial phenomenon to generate a twisted homobilayer which can be captured through atomistic simulations, we now turn towards an autonomous approach to expedite the process of discovering optimal parameter combinations leading to a target twisted final state, given on demand. While sequential acquisition based Bayesian optimization loops have been commonly deployed in self-driving laboratories,^3,30,41,49 in the context of building virtual environments, throughput in terms of both sampling and decision making can be increased through extreme scale computational resources. Furthermore, many of the state-of-the-art materials simulation capabilities are highly optimized into applications exploiting leadership class supercomputers.⁴²

Our approach therefore is to accelerate the discovery of recrystallization parameters through asynchronous iterative evaluations of a batch of candidates (p_i, i = 1, 2, …, n_batch) recommended by BO. Several algorithms have explored the idea of parallel optimization of the acquisition function with various approximation^50,51 to account for the uncertainties of the pending evaluations. First, we start with q-expected improvement (EI) (cf. Methods) as our acquisition function to sample a target twist angle, θ_target = 8.5°. In Fig. 3a, we show the variation in predictions with different batch sizes (n_batch). The regrets from the target for each of the acquisition events are shown in Fig. 3a and 5a. These acquisitions are driven by Bayesian uncertainty quantification metrics as each MD trajectory is spawned with a set of annealing temperature, lattice and shear strains under the recommendation of algorithms as defined in the Methods section. To better follow the progress of this regret minimization process, we added cumulative minimum lines for each case (cf. Fig. 3a and 5a). From this analysis, we can see that the cumulative minimum gradually drops discretely as sample trajectories are evaluated. For qEI driven acquisition, with a batch size = 5, the best regret (0.1°) in Fig. 3a leads a final average recrystallized moiré twist of 8.4°, i.e. target (8.5°) − regret (0.1°) as achieved on the 38th recommendation (5-th cycle) at a temperature of 1386 K, 1.2% lattice strain and 5.53% shear (Fig. S5).

Fig. 3 Sampling a target twist with batch acquisition based on q-expected improvement. Regret (absolute error) between the sampled twist angle and the predefined target for q-EI based acquisitions with different batch sizes (n_batch) in (a). Solid lines correspond to the cumulative minimum of regret over the iterations. Mean regret evolution in (b) for the best sample in each batch over BO cycle iterations for varying batch sizes. The confidence interval shows standard deviation over 3 repetitions.

Meanwhile, starting with an identical set of initial observation data, with a larger batch size (n_batch = 15), the workflow predicts a twisted phase closer to the target, θ = 8.46° (Fig. 3a and Fig. S5c and d) for temperature = 2077 K, 2.3% lattice strain and 5.76% shear at the 4-th cycle (58-th sample).

While it is obvious that larger batches can be harnessed to cover a larger search space (Fig. S5b and d) in faster wall clock time with scalable resources, the indirect exploitation of the posterior, as it is the case for EI as an acquisition function,⁵² might not always lead to efficient searches (Fig. 4a) with a larger batch size. We could define a per-cycle sampling efficiency metric as an average rate (computed normalizing the change of mean best regret per cycle over 4 cycles) describing how fast BO iteration cycles generate an MD trajectory with a target twisted phase. As seen in Fig. 4b, this is driven by the fact that beyond a point larger batches tend to make the search more explorative rather than efficiently exploiting model uncertainties accessed through implicit acquisition criteria like EI.⁵³

Fig. 4 Effect of batch size and sampling strategies. Initialized with an identical set of observations, performance of two different acquisition strategies with batch evaluation–Thompson sampling and expected improvement in (a). Overall efficiency of various batch sampling strategies with respect to size and acquisitions in (b). Inherent stochastic nature of Thompson sampling outperforms EI based acquisition which can become explorative beyond an optimal batch size.

To probe the effect of acquisition strategies, we adopt batch Thompson sampling (cf. Methods), shown in Fig. 5a, for the same target angle (θ = 8.5°) and identical prior observation set. It is interesting to note that Thompson sampling with a comparable batch size (i.e., n_batch = 15) leads to more consistent sampling near the target, yielding three twist angles (sample index) e.g., 8.39°(18) and 8.57°(47) at different parameter combinations (T₁₈ = 1362 K, ε_{lat 18} = 0.4%, ε_{sh 18} = 5.69%), (T₄₇ = 3195 K, ε_{lat 47} = 0.2%, ε_{sh 47} = 5.5%) than EI. This ability to scan through a diverse range of recrystallization parameters could be due to the probabilistic nature of Thompson sampling from the posterior, which leads to a better exploration–exploitation trade-off in batch mode⁵² than primarily indirect exploitation (of posterior properties) driven acquisitions through EI. Furthermore, we assess the advantage of our active learning-driven sampling approach over vanilla random sampling (Fig. 5b). Starting with an identical initial observation set like all the previous sampling demonstrations, random sampling with n_batch = 15 is unable to sample twist angles close to the target (8.5°) even up to 120 samples, with the best value being 8.2°.

Fig. 5 Comparison of adaptive batch sampling via BO with random sampling. Performance of two different acquisition strategies with batch evaluation–Thompson sampling and expected improvement (a). Solid lines correspond to the cumulative minimum of regret over the iterations. Initialized with an identical set of observations, random sampling (b) underperforms within a comparable sampling duration.

Discussions

Designing twisted bilayer interfaces on-demand

We thus far have demonstrated a capability to identify target recrystallization conditions through our Bayesian optimization workflow coupled with asynchronous ensemble evaluation strategies driven by scalable computational resources. While constrained in the scope of parameter space (e.g., more realistic conditions might involve more drivers like pressure, chemical potential, rate of heating, etc.), these results show the potential of accelerating targeted solid state phase transitions in vdW 2D material interfaces. To further demonstrate the efficacy of our workflow, we provide quantitative evidence towards enabling on-demand design of twisted moiré interfaces of MoS₂ bilayers. Using batch Thompson sampling as a mode to evaluate the parameter space comprising temperature, lattice strain and shear, we focus our workflow towards sampling twist angles 15° and 25° in addition to 8.5°, thereby covering a wide range of moiré wavelengths and superlattices. As shown in Fig. 6, the current approach is able to achieve recrystallization into the moiré lattice with predefined target twist angles, especially for structures with nontrivial target twist angles (e.g. 25°) (cf. Fig. S6 for more details).

Fig. 6 Achieving recrystallization with on demand final twist angles. General flexibility of the workflow is shown using batch Thompson sampling. Starting from the same set of initial samples, the asynchronous Thompson sampling could search three diverse set of moiré sublattices with target interlayer twists-8.5° (a and d), 15° (b and e) and 25° (c and f).

It will be further intriguing to probe the effects of edge reactions, stacking orders, perturbed local environments and defect-pinning/doping that results from atomic rearrangement during the melt-quench recrystallization as well as precursor-driven growth,⁵⁴ on their electronic structure via first principles. The advantage of the proposed online workflow orchestration strategy is that such effects can now be more readily investigated leveraging HPC platforms.

Moreover, we successfully identify multiple combinations of parameter corresponding to each of these diversely specified target twists (Fig. 6). For instance, interlayer twists close to 15° can be recrystallized at (i) 2725 K, 2.7% lattice strain and 5.5% shear as well as at (ii) 2107 K, 1.3% lattice strain and 5.69% shear. Likewise, we identify more than one set of optimal parameters for θ ∼ 25°—(i) {2797 K, 0.9% lattice strain and 5.57% shear} and (ii) {1379 K, 3.2% lattice strain and 5.51% shear}. We note that very small changes in the recrystallization parameters, especially shearing the supercell, lead to drastically different final moiré superlattices. For this reason, we keep the search window for ε_sh to vary from 5.5% to 5.75% strain levels, corresponding to absolute distances of 1 Å to 1.5 Å tilts of in-plane supercell dimensions. Furthermore, the wide variations (Fig. S4) in optimal temperature and lattice strain conditions associated with a specific final twisted bi-layer phase lead us to speculate that the recrystallization landscape could be rough having multiple locally varying optimal cases. Traditional grid search or sequential active learning protocols hence might not be adequate to exhaustively identify favorable recrystallization conditions for a given target twist (Fig. S1). Furthermore, the current approach of high-throughput classical force-field driven simulations exhibits linear scaling O(N) in comparison to cubic scaling ab initio/DFT simulations O(N³), making it much more competitive towards accessing sufficient timespans for the in silico recrystallization simulation recipe described here, within a matter of a few hours of wall clock time with hpc-scalable batch parallel BO as long as enough compute is available. However, we note that bridging the timescales of state-of-the-art atomistic simulations with what is relevant for experiments is a long standing multiscale challenge, for which we remark that our end-to-end developments of asynchronous on-the-fly active learning atomistics are ideally suited to upscale the dynamics with coarse grained models, e.g., kinetic Monte Carlo, phase field simulations, etc. The presented online ensemble sampling workflow for identifying recrystallization regimes of MoS2 bi-layers could be seamlessly extended towards various other TMDC systems (e.g., W/Mo–S/Se/Te etc.), multilayers or even heterostructures, which are of increasing interest due to their highly tunable quantum properties. In practice, such adoptions are not limited by any of the workflow features which are rather general and mostly take atomistic systems and force fields as inputs. If appropriate force-field descriptions are available for a TMDC system of choice, the current framework should be fully transferable. With the rise of highly accurate foundation MLIPs, the bottlenecks of force descriptions are rapidly disappearing. Adapting these methods to more stochastic or continuum simulation approaches is also quite straightforward.

Conclusion

We have demonstrated that combining extreme-scale computations, Bayesian optimisation and asynchronous automated computational workflows, a powerful building block to enable on-demand synthesis can be achieved, and that it allows time and resource efficient discovery of synthesis pathways for complex materials. The targeted problem of achieving ‘on-demand’ van der Waals epitaxy is a major technological challenge, specifically when controlling both the phase and the interlayer stacking orientations. We demonstrate that within 20 iterations the autonomous computational synthesis platform achieves the targeted moiré quantum heterostructure. The global sampling opportunities boosted by extreme-scale resource-driven batch acquisition workflows and on-the-fly analysis and UQ capabilities are ideally suited to pin down more sophisticated interplay of thermodynamic and kinetic factors driving self-organized and strain-engineered interfacial evolution during vdW epitaxy. Furthermore, we note that integrating the autonomous computational synthesis workflows with real experimental synthesis efforts can quickly narrow down the number of experimental attempts to grow a targeted quantum heterostructure. The current ensemble simulation framework can also be used to generate digital-twins on experimental measurements, such as X-ray diffraction (XRD) and reflective high-energy electron diffraction (RHEED). As such, integration of our framework with experimental molecular beam epitaxy (MBE) or pulsed laser deposition (PLD) systems will allow us to perform on-the-fly comparisons with experimental XRD/RHEED measurements, so as to guide experiments in real-time, allowing a real autonomous thin-film synthesis platform with theory-in-the-loop.

Methods

: adaptive real-time ensemble task management environment for extreme scale computing

Within a single HPC batch allocation, ensemble evaluations/jobs, which require MPI parallel resources at the individual level, often suffer from scheduling challenges with typical SLURM-based workload managers.⁵⁵ This is primarily due to the fact that hard ceilings (e.g., 100 for Frontier, Perlmutter/NERSC) are often imposed toward chaining together multiple job submissions in batch mode (e.g., for SLURM). Such setups are highly unsuitable when (i) the size of the ensemble exceeds the ceiling of a scheduler and (ii) dynamic heterogeneous tasks are at the core of the workflow; e.g., based on outcomes of the ensemble of MD simulations, active learning or UQ algorithms (e.g., BO) have to operate and further spawn the next generation of ensembles.

To mitigate the above challenges, we implement all our active learning MD simulations within a recently developed adaptive and dynamic orchestrator framework ,⁴⁶ which has an executor back-end of , an exascale friendly HPC resource manager with a hierarchical graph-based scheduler allowing for a generalized flexible custom operation within a single batch allocation. ⁴⁶ benefits from the native python executor interface of ,⁵⁶ and the concurrent asynchronous programming model of core python through objects.⁵⁷ A continuous throughput is maintained via dynamically spawning and monitoring tasks. Furthermore, to enable real-time streaming of post-processed data (e.g., interlayer twist angle evolution) from large-scale atomistic trajectories, an in-memory data analysis protocol is used by exploiting the heterogeneous (GPU + CPU) architecture of exascale systems (e.g., frontier) via a round-robin MPI-communicator splitting approach (cf. ref. 46). As explained in the following sections, such an online adaptive framework enables efficient coupling between available computing resource chunks for ensemble evaluations guided by adaptive sampling methods.

Parallel batch sampling with Bayesian optimization

Traditionally optimal design of experiments via active learning based on Bayesian optimization⁵⁸ proceeds through sequential acquisitions. With the adaptive asynchronous capability to produce scalable ensemble evaluations in our workflow infrastructure, batch mode acquisition algorithms are natural choices to effectively utilize available hpc resources and accelerate the search towards target twist angles in the spirit of efficient global optimization⁵⁹ with the following acquisition functions.

q-Expected improvement. With a goal to accelerate the solution of the minimization problem for an expensive black-box function f, we have used parallel/batch expected improvement^60–63 in this study. It is computed by Monte Carlo sampling and is given bywhere ξⁱ(x) ∼ f(x) are samples drawn randomly from the joint posterior distribution of the model evaluated at batch points x, before optimizing the expectation which dictates the acquisition strategy rather than inferences which are directly based on the posterior, and top performing n_batch candidates from the optimization are finally obtained. Relevant to the last point, we also explored more direct posterior sampling strategies, e.g. Thompson sampling, as described in the following.

Batch Thompson sampling. Unlike several acquisition functions, which are mostly driven by exploitative strategies, Thompson sampling, due to its stochastic nature, is able to handle the exploration–exploitation trade-off in a better fashion. In our simulation, we adopt the batch parallel version of Thompson sampling first proposed by Kandasamy et al.⁵¹ We use a quasi-Monte Carlo sampler which uses Sobol sequences⁶⁴ to sample a user-defined function from joint-posterior over the entire batch of size n_batch and optimizes over this chosen function to find the next batch of candidates. Following Kandasamy et al.⁵¹ we choose simple regret as the user-defined function,

Upon finishing the evaluation of all candidates in the batch, the next cycle of the BO is resumed.

We implement all our BO algorithms using the library⁶⁵ for a single task Gaussian process approximated with the Matérn kernel.⁶⁶

Molecular dynamics simulations of recrystallization in MoS₂ bilayers

We perform all our MD simulations by wrapping the shared library utilities of the open source package ,⁶⁷ through drivers. Starting with a 2H bilayer MoS₂ phase (a rhombohedral supercell with 864 atoms), we first induce disorder in the top crystalline layer by selectively heating up the layer at a high temperature of 5000 K for 50 ps with a timestep of 0.5 fs. We use a classical reactive force field ReaxFF, specifically re-parametrized to accurately capture order-disorder transitions in MoS₂ during multiple melt-quench cycles.⁴⁷ After the top layer undergoes melting, the whole bilayer system is equilibrated at 300 K at atmospheric pressure using a Berendsen NPT thermostat for 50 ps to achieve the disordered top layer placed on top of a crystalline bottom layer—this served as the initial sample (amorph + Xtal) for the recrystallization simulations.

Finally to simulate the recrystallization process under various candidate input parameter combination of p_i ≡ {T_i, ε_lat, ε_sh}, the system is rapidly heated for 50 ps at a temperature T_i along with in-plane lattice rescaling (x, y) of (1 + ε_lat) (to model epitaxial strain) and an in-plane shear strain (xy) of ε_sh measured on the substrate crystal layer (cf. Fig. 1a). This is then followed by a slow quench of the bilayer system to reduce the temperature level from T_i to 0 K over 0.5 ns i.e. 1 × 10⁶ MD steps. The final step of rapid melt-slow quench simulations is performed using a Langevin thermostat (NVT). Periodic boundaries are maintained along in-plane (i.e. x and y) directions. Lennard-Jones reflective walls are implemented near the bottom (0 Å) and top (12 Å) edges. We use Ovito for all the visualizations and analysis of MD trajectories.

It is important to note that the annealing approach does not guarantee the lowest energy configurations, which is intentional for the present study to let the workflow explore high-energy metastable basins and discover a diverse range of on-demand moire motifs which might not necessarily be equilibrium states. In the context of synthesis and materials discovery, where tracking down non-equilibrium phases are of particular interest, such a melt-quench dynamics approach could specifically be useful. While it is possible for atoms to undergo rare diffusion events which are path dependent mostly in the long-time slow quench limit, our 0.5 ns quench protocol could still be relatively fast to favor such path dependent quenching behavior in the first order. This could be explored in more detail by studying long-time dynamics limits, which is beyond the current scope.

Constructing an order parameter. To quantify the emerging order in the MoS₂ bilayers (Fig. 2), we define a simple pair correlation based custom order parameter which estimates the sharpness of the first nearest neighbor peak distribution. First, the pair correlation function is computed with a cutoff of 3 Å based on the first nearest neighbors. Then, we count the number of pairs within a pair distance window ([2.4 Å, 2.6 Å]), where the first nearest neighbor peak/mode of the distribution for pristine MoS₂ is around 2.47 Å.

Computing interlayer twist angles on-the-fly

Adopting a streaming approach to analyze MD trajectories as they become available through a distributed in-memory message passing framework, we accelerate and perform on-the-fly data reduction. As described in more detail in ref. 46, apart from being scalable and ideally suited for integrated active learning workflows, such frameworks serve as general testbeds for a variety of atomistic analysis and visualization algorithms. For this study, we wrap through a simple and computationally light-weight iterative algorithm to compute the interlayer twists for atomistic snapshots along the recrystallization trajectory.

Analysis of the polycrystallinity and validity of the twist computation algorithm

To test the validity of our approach, we performed a more rigorous per-atom bond-orientation based interlayer twist computation approach. As can be seen (Fig. S7) from the randomly considered final quenched phases for a series of recrystallization simulations under various temperature, lattice and shear strain conditions, the twist distributions mostly tend to be unimodal leading to no speculation of multiple grains/polycrystalline nature. Barely some small clusters of other angles are observed, which are likely due to disorder and vacancy debris after the quench. Furthermore, the modes of per-atom twist histograms show excellent agreement with our much simpler single twist results for most of the cases. Note that due to hexagonal symmetry, angles that are numerically close to 0 deg. and 60 deg. are physically identical and hence the bimodal distributions are still single crystal phases.

Author contributions

SB designed the core real-time adaptive learning workflow and on-the-fly MD simulations of recrystallization, implemented batch TS sampling and wrote the manuscript draft. AG co-wrote and designed the BO approach. AB co-wrote the manuscript, implemented BO and analyzed the results. PVB discussed overall results and co-wrote the manuscript. PG wrote, supervised and conceived the project.

Conflicts of interest

There are no conflicts to declare.

Data availability

Data and code for this article, including example workflows and analysis notebooks, are available at Zenodo at https://doi.org/10.5281/zenodo.18444064.

Supplementary information (SI) is available. See DOI: https://doi.org/10.1039/d6dd00049e.

Acknowledgements

This research is sponsored by the INTERSECT Initiative as part of the Laboratory Directed Research and Development Program of Oak Ridge National Laboratory, managed by UT-Battelle, LLC, for the U.S. Department of Energy under contract DE-AC05-00OR22725 via the QCAD project. This research used resources of the Oak Ridge Leadership Computing Facility at the Oak Ridge National Laboratory, which is supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC05-00OR22725, via the Innovative and Novel Computational Impact on Theory and Experiment (INCITE) program. Resources of the Compute and Data Environment for Science (CADES) was used via a user project at the Center for Nanophase Materials Sciences (CNMS), which is a US Department of Energy, Office of Science User Facility at the Oak Ridge National Laboratory.

References

1. J. Li, X. Yang, Z. Zhang, W. Yang, X. Duan and X. Duan, Towards the scalable synthesis of two-dimensional heterostructures and superlattices beyond exfoliation and restacking, Nat. Mater., 2024, 23, 1326–1338.
2. N. J. Szymanski, B. Rendy, Y. Fei, R. E. Kumar, T. He, D. Milsted, M. J. McDermott, M. Gallant, E. D. Cubuk, A. Merchant and others, An autonomous laboratory for the accelerated synthesis of novel materials, Nature, 2023, 624, 86–91.
3. S. Ament, M. Amsler, D. R. Sutherland, M.-C. Chang, D. Guevarra, A. B. Connolly, J. M. Gregoire, M. O. Thompson, C. P. Gomes and R. B. Van Dover, Autonomous materials synthesis via hierarchical active learning of nonequilibrium phase diagrams, Sci. Adv., 2021, 7, eabg4930.
4. C. Liu, Y.-C. Lin, M. Yoon, Y. Yu, A. A. Puretzky, C. M. Rouleau, M. F. Chisholm, K. Xiao, G. Eres, G. Duscher and others, Understanding substrate-guided assembly in van der Waals epitaxy by in situ laser crystallization within a transmission electron microscope, ACS Nano, 2021, 15, 8638–8652.
5. K. M. Roccapriore, M. G. Boebinger, O. Dyck, A. Ghosh, R. R. Unocic, S. V. Kalinin and M. Ziatdinov, Probing electron beam induced transformations on a single-defect level via automated scanning transmission electron microscopy, ACS Nano, 2022, 16, 17116–17127.
6. A. K. Geim and I. V. Grigorieva, Van der Waals heterostructures, Nature, 2013, 499, 419–425.
7. Y. Cao, V. Fatemi, A. Demir, S. Fang, S. L. Tomarken, J. Y. Luo, J. D. Sanchez-Yamagishi, K. Watanabe, T. Taniguchi, E. Kaxiras and others, Correlated insulator behaviour at half-filling in magic-angle graphene superlattices, Nature, 2018, 556, 80–84.
8. Y. Cao, V. Fatemi, S. Fang, K. Watanabe, T. Taniguchi, E. Kaxiras and P. Jarillo-Herrero, Unconventional superconductivity in magic-angle graphene superlattices, Nature, 2018, 556, 43–50.
9. P. Rivera, H. Yu, K. L. Seyler, N. P. Wilson, W. Yao and X. Xu, Interlayer valley excitons in heterobilayers of transition metal dichalcogenides, Nat. Nanotechnol., 2018, 13, 1004–1015.
10. K. Kim, M. Yankowitz, B. Fallahazad, S. Kang, H. C. Movva, S. Huang, S. Larentis, C. M. Corbet, T. Taniguchi, K. Watanabe and others, van der Waals heterostructures with high accuracy rotational alignment, Nano Lett., 2016, 16, 1989–1995.
11. X. Xu, Y. Pan, S. Liu, B. Han, P. Gu, S. Li, W. Xu, Y. Peng, Z. Han, J. Chen and others, Seeded 2D epitaxy of large-area single-crystal films of the van der Waals semiconductor 2H MoTe2, Science, 2021, 372, 195–200.
12. J.-H. Baek, H. G. Kim, S. Y. Lim, S. C. Hong, Y. Chang, H. Ryu, Y. Jung, H. Jang, J. Kim, Y. Zhang and others, Thermally induced atomic reconstruction into fully commensurate structures of transition metal dichalcogenide layers, Nat. Mater., 2023, 22, 1463–1469.
13. P. Sutter, S. Wimer and E. Sutter, Chiral twisted van der Waals nanowires, Nature, 2019, 570, 354–357.
14. P. Sutter, R. Ibragimova, H.-P. Komsa, B. A. Parkinson and E. Sutter, Self-organized twist-heterostructures via aligned van der Waals epitaxy and solid-state transformations, Nat. Commun., 2019, 10, 5528.
15. P. Gehring, B. F. Gao, M. Burghard and K. Kern, Growth of high-mobility Bi2Te2Se nanoplatelets on hBN sheets by van der Waals epitaxy, Nano Lett., 2012, 12, 5137–5142.
16. J. Kim, C. Bayram, H. Park, C.-W. Cheng, C. Dimitrakopoulos, J. A. Ott, K. B. Reuter, S. W. Bedell and D. K. Sadana, Principle of direct van der Waals epitaxy of single-crystalline films on epitaxial graphene, Nat. Commun., 2014, 5, 4836.
17. Y.-C. Lin, N. Lu, N. Perea-Lopez, J. Li, Z. Lin, X. Peng, C. H. Lee, C. Sun, L. Calderin, P. N. Browning and others, Direct synthesis of van der Waals solids, ACS Nano, 2014, 8, 3715–3723.
18. Y. Wan, J.-K. Huang, C.-P. Chuu, W.-T. Hsu, C.-J. Lee, A. Aljarb, C.-W. Huang, M.-H. Chiu, H.-L. Tang, C. Lin and others, Strain-directed layer-by-layer epitaxy toward van der Waals homo-and heterostructures, ACS Mater. Lett., 2021, 3, 442–453.
19. Y. Lu, J. Chen, M. J. Coupin, S. Sinha and J. H. Warner, Lattice-Mismatch-Driven Small-Angle Moiré Twists in Epitaxially Grown 2D Vertical Layered Heterostructures, Adv. Mater., 2022, 34, 2205403.
20. A. Maksov, O. Dyck, K. Wang, K. Xiao, D. B. Geohegan, B. G. Sumpter, R. K. Vasudevan, S. Jesse, S. V. Kalinin and M. Ziatdinov, Deep learning analysis of defect and phase evolution during electron beam-induced transformations in WS2, npj Comput. Mater., 2019, 5, 12.
21. K. M. Omambac, H. Hattab, C. Brand, G. Jnawali, A. T. N'Diaye, J. Coraux, R. van Gastel, B. Poelsema, T. Michely, F.-J. Meyer zu Heringdorf and others, Temperature-controlled rotational epitaxy of graphene, Nano Lett., 2019, 19, 4594–4600.
22. T. A. de Jong, T. Benschop, X. Chen, E. E. Krasovskii, M. J. de Dood, R. M. Tromp, M. P. Allan and S. J. Van der Molen, Imaging moiré deformation and dynamics in twisted bilayer graphene, Nat. Commun., 2022, 13, 70.
23. C. Jenks, H. N. Lee, J. Lewis, C. R. Kagan, P. Nealey, P. Braun, J. E. Holladay, Y. Gao, D. S. Sholl and B. A. Helms, others, Basic Research Needs for Transformative Manufacturing, 2020.
24. W. Li, X. Qian and J. Li, Phase transitions in 2D materials, Nat. Rev. Mater., 2021, 6, 829–846.
25. J. M. Woods, Y. Jung, Y. Xie, W. Liu, Y. Liu, H. Wang and J. J. Cha, One-step synthesis of MoS2/WS2 layered heterostructures and catalytic activity of defective transition metal dichalcogenide films, ACS Nano, 2016, 10, 2004–2009.
26. M. Bianchini, J. Wang, R. J. Clément, B. Ouyang, P. Xiao, D. Kitchaev, T. Shi, Y. Zhang, Y. Wang, H. Kim and others, The interplay between thermodynamics and kinetics in the solid-state synthesis of layered oxides, Nat. Mater., 2020, 19, 1088–1095.
27. A. G. Kusne, T. Gao, A. Mehta, L. Ke, M. C. Nguyen, K.-M. Ho, V. Antropov, C.-Z. Wang, M. J. Kramer, C. Long and others, On-the-fly machine-learning for high-throughput experiments: search for rare-earth-free permanent magnets, Sci. Rep., 2014, 4, 6367.
28. J. R. Hattrick-Simpers, J. M. Gregoire and A. G. Kusne, Perspective: composition–structure–property mapping in high-throughput experiments: turning data into knowledge, APL Mater., 2016, 4(5), 053211.
29. J. Chen, S. R. Cross, L. J. Miara, J.-J. Cho, Y. Wang and W. Sun, Navigating phase diagram complexity to guide robotic inorganic materials synthesis, Nat. Synth., 2024, 1–9.
30. H. Sheng, J. Sun, O. Rodríguez, B. B. Hoar, W. Zhang, D. Xiang, T. Tang, A. Hazra, D. S. Min, A. G. Doyle and others, Autonomous closed-loop mechanistic investigation of molecular electrochemistry via automation, Nat. Commun., 2024, 15, 2781.
31. J. Leeman, Y. Liu, J. Stiles, S. B. Lee, P. Bhatt, L. M. Schoop and R. G. Palgrave, Challenges in High-Throughput Inorganic Materials Prediction and Autonomous Synthesis, PRX Energy, 2024, 3, 011002.
32. B. G. Sumpter, R. K. Vasudevan, T. Potok and S. V. Kalinin, A bridge for accelerating materials by design, npj Comput. Mater., 2015, 1, 1–11.
33. A. Ghosh, M. Ziatdinov, O. Dyck, B. G. Sumpter and S. V. Kalinin, Bridging microscopy with molecular dynamics and quantum simulations: An atomAI based pipeline, npj Comput. Mater., 2022, 8, 74.
34. K. J. Kanarik, W. T. Osowiecki, Y. Lu, D. Talukder, N. Roschewsky, S. N. Park, M. Kamon, D. M. Fried and R. A. Gottscho, Human–machine collaboration for improving semiconductor process development, Nature, 2023, 616, 707–711.
35. B. P. MacLeod, F. G. Parlane, T. D. Morrissey, F. Häse, L. M. Roch, K. E. Dettelbach, R. Moreira, L. P. Yunker, M. B. Rooney, J. R. Deeth and others, Self-driving laboratory for accelerated discovery of thin-film materials, Sci. Adv., 2020, 6, eaaz8867.
36. S. Zhu, P. Pochet and H. T. Johnson, Controlling rotation of two-dimensional material flakes, ACS Nano, 2019, 13, 6925–6931.
37. S. Bagchi, H. Johnson and H. B. Chew, Rotational stability of twisted bilayer graphene, Phys. Rev. B, 2020, 101, 054109.
38. S. Dai, Y. Xiang and D. J. Srolovitz, Twisted bilayer graphene: Moiré with a twist, Nano Lett., 2016, 16, 5923–5927.
39. S. Bagchi, H. T. Johnson and H. B. Chew, Strain-controlled dynamic rotation of twisted 2D atomic layers for tunable nanomechanical systems, ACS Appl. Nano Mater., 2020, 3, 10878–10884.
40. N. P. Kazmierczak, M. Van Winkle, C. Ophus, K. C. Bustillo, S. Carr, H. G. Brown, J. Ciston, T. Taniguchi, K. Watanabe and D. K. Bediako, Strain fields in twisted bilayer graphene, Nat. Mater., 2021, 20, 956–963.
41. A. G. Kusne, H. Yu, C. Wu, H. Zhang, J. Hattrick-Simpers, B. DeCost, S. Sarker, C. Oses, C. Toher, S. Curtarolo and others, On-the-fly closed-loop materials discovery via Bayesian active learning, Nat. Commun., 2020, 11, 5966.
42. F. Alexander, A. Almgren, J. Bell, A. Bhattacharjee, J. Chen, P. Colella, D. Daniel, J. DeSlippe, L. Diachin, E. Draeger and others, Exascale applications: skin in the game, Philos. Trans. R. Soc. A Math. Phys. Eng. Sci., 2020, 378, 20190056.
43. S. M. Mniszewski, J. Belak, J.-L. Fattebert, C. F. Negre, S. R. Slattery, A. A. Adedoyin, R. F. Bird, C. Chang, G. Chen, S. Ethier and others, Enabling particle applications for exascale computing platforms, Int. J. High Perform. Comput. Appl., 2021, 35, 572–597.
44. R. Wang, N. Koch, J. Martin and S. Sadofev, Recrystallization of MBE-Grown MoS2 Monolayers Induced by Annealing in a Chemical Vapor Deposition Furnace, Phys. Status Solidi RRL, 2023, 17, 2200476.
45. A. Weston, Y. Zou, V. Enaldiev, A. Summerfield, N. Clark, V. Zólyomi, A. Graham, C. Yelgel, S. Magorrian, M. Zhou and others, Atomic reconstruction in twisted bilayers of transition metal dichalcogenides, Nat. Nanotechnol., 2020, 15, 592–597.
46. S. Bagchi, others, MatEnsemble: An exascale compatible asynchronous integrated on-the-fly active learning frameowrk for materials ensemble simulations, 2025, https://github.com/Q-CAD/MatEnsemble.
47. I. Ponomarev, T. Polcar and P. Nicolini, New reactive force field for simulations of MoS2 crystallization, J. Phys. Chem. C, 2022, 126, 9475–9481.
48. M. Xu, H. Ji, L. Zheng, W. Li, J. Wang, H. Wang, L. Luo, Q. Lu, X. Gan, Z. Liu and others, Reconfiguring nucleation for CVD growth of twisted bilayer MoS2 with a wide range of twist angles, Nat. Commun., 2024, 15, 562.
49. B. J. Shields, J. Stevens, J. Li, M. Parasram, F. Damani, J. I. M. Alvarado, J. M. Janey, R. P. Adams and A. G. Doyle, Bayesian reaction optimization as a tool for chemical synthesis, Nature, 2021, 590, 89–96.
50. J. M. Hernández-Lobato, J. Requeima, E. O. Pyzer-Knapp and A. Aspuru-Guzik, Parallel and distributed Thompson sampling for large-scale accelerated exploration of chemical space, International conference on machine learning, 2017, pp. 1470–1479.
51. K. Kandasamy, A. Krishnamurthy, J. Schneider and B. Póczos, Parallelised Bayesian optimisation via Thompson sampling, International conference on artificial intelligence and statistics, 2018, pp. 133–142.
52. J. González, Z. Dai, P. Hennig and N. Lawrence, Batch Bayesian optimization via local penalization, Artificial intelligence and statistics, 2016, pp. 648–657.
53. A. De Palma, C. Mendler-Dünner, T. Parnell, A. Anghel and H. Pozidis, Sampling acquisition functions for batch Bayesian optimization, arXiv, 2019, preprint, arXiv:1903.09434,  DOI:10.48550/arXiv.1903.09434.
54. M. A. M. Filho, W. Farmer, C.-L. Hsiao, R. B. Dos Santos, L. Hultman, J. Birch, K. Ankit and G. K. Gueorguiev, Density functional theory-fed phase field model for semiconductor nanostructures: The case of self-induced core–shell InAlN nanorods, Cryst. Growth Des., 2024, 24, 4717–4727.
55. A. B. Yoo, M. A. Jette and M. Grondona, Slurm: Simple linux utility for resource management, Workshop on job scheduling strategies for parallel processing, 2003, pp. 44–60.
56. D. H. Ahn, N. Bass, A. Chu, J. Garlick, M. Grondona, S. Herbein, H. I. Ingólfsson, J. Koning, T. Patki, T. R. Scogland and others, Flux: Overcoming scheduling challenges for exascale workflows, Future Gener. Comput. Syst., 2020, 110, 202–213.
57. B. Quinlan, Futures – Execute computations asynchronously, PEP 3148, 2009, available: https://www.python.org/dev/peps/pep-3148/.
58. J. Močkus, On Bayesian methods for seeking the extremum, Optimization techniques IFIP technical conference: Novosibirsk, July 1–7, 1974, 1975, pp. 400–404.
59. D. R. Jones, M. Schonlau and W. J. Welch, Efficient global optimization of expensive black-box functions, J. Global Optim., 1998, 13, 455–492.
60. S. Daulton, S. Ament, D. Eriksson, M. Balandat and E. Bakshy, Unexpected improvements to expected improvement for Bayesian optimization, Proceedings of the 37th International Conference on Neural Information Processing Systems, Red Hook, NY, USA, 2023.
61. D. Ginsbourger, R. Le Riche and L. Carraro, Computational intelligence in expensive optimization problems, Springer, 2010, pp. 131–162.
62. H. J. Kushner, A New Method of Locating the Maximum Point of an Arbitrary Multipeak Curve in the Presence of Noise, J. Basic Eng., 1964, 86, 97–106.
63. D. Ginsbourger, R. L. Riche and L. Carraro, A Multi-points Criterion for Deterministic Parallel Global Optimization based on Gaussian Processes, 2008.
64. J. Snoek, H. Larochelle and R. P. Adams, Practical bayesian optimization of machine learning algorithms, Adv. Neural Inf. Process. Syst., 2012, 25.
65. M. Balandat, B. Karrer, D. Jiang, S. Daulton, B. Letham, A. G. Wilson and E. Bakshy, BoTorch: A framework for efficient Monte-Carlo Bayesian optimization, Adv. Neural Inf. Process. Syst., 2020, 33, 21524–21538.
66. C. K. Williams and C. E. Rasmussen, Gaussian processes for machine learning, MIT press Cambridge, MA, 2006, vol. 2.
67. A. P. Thompson, H. M. Aktulga, R. Berger, D. S. Bolintineanu, W. M. Brown, P. S. Crozier, P. J. In’t Veld, A. Kohlmeyer, S. G. Moore, T. D. Nguyen and others, LAMMPS-a flexible simulation tool for particle-based materials modeling at the atomic, meso, and continuum scales, Comput. Phys. Commun., 2022, 271, 108171.

---