Predicting emergence of crystals from amorphous precursors with deep learning potentials

Main

When amorphous solids crystallize, the first phases to emerge are often not the thermodynamic ground states, but are instead metastable crystals^1,2,3,4,5. The formation of metastable phases during such crystallization processes is ubiquitous in nature, a well-known example being the crystallization of amorphous calcium carbonate to vaterite, aragonite or calcite, which have roles in biomineralization and the carbon cycle^6,7. The crystallization sequences of amorphous ice and interstellar molecules have importance in geology, astrophysics and astrobiology^8,9. In materials science, the amorphous-to-crystalline transformation underpins the development of new technologies, including phase-change memories¹⁰, nanocrystallized metallic glass soft magnets¹¹ and ceramics with finely tuned electronic or optical properties¹². In synthesis in particular, the crystallization of non-crystalline precursors plays a central role^8,13,14,15. The complexity of nucleation has hindered the formulation of polymorph-selective synthesis models¹⁶, but has encouraged the development of numerous experimental techniques to explore the energy landscape, such as solid-state reactions, chimie douce, deposition, rapid cooling or high-pressure methods^17,18. Unlike polymorphs, an amorphous precursor can be accessed with many of these techniques and used as a precursor to metastable crystals^13,14,19. This path to metastable crystals is an instance of the rule of Ostwald, explained by classical nucleation theory, whereby crystals that share local structural motifs with the amorphous precursor will have a lower barrier to nucleation upon annealing (Fig. 1a)^1,20,21. The structural connection renders the amorphous state a uniquely useful starting point for the predictive synthesis of crystals. To capture the energy landscape itself, state-of-the-art deep learning potentials trained on first-principles calculations are well poised to provide a reliable description at the scale of atoms^22,23,24. Guided by the ubiquity of the amorphous precursor and its local structural connection to crystals, here we introduce a computational approach, known as a²c, that achieves highly accurate predictions of the (meta)stable crystallization products of amorphous matter with the transformation of local motifs in disordered states into ordered structures, using universal deep learning potentials, at a scale impractical for ab initio methods.

**Fig. 1: Predicting the crystallization products of amorphous precursors.**

Results

With a²c we start by creating an atomistic model of the amorphous precursor that captures the essential short-range order via melt-and-quench molecular dynamics (MQMD)²⁵ (Fig. 1b). We hypothesize that if a large number of subcells of the amorphous atom configuration are relaxed under periodic-boundary conditions following downhill energy gradients using an accurate interatomic potential (Fig. 1c), a fraction of those with the seeding motif should relax into the adjacent basin, which nucleates in practice (as exemplified in Supplementary Fig. 1). The lowest-energy configurations reached with this process should be a close match to the crystal structure of the experimentally observed initial crystallization product of the amorphous material, irrespective of the product being stable or metastable. Initial crystallization products predicted using a²c for 12 inorganic systems are presented in Fig. 1d. In each case, a²c sifts through tens of thousands of local motifs in the amorphous precursor using deep learning-based interatomic interactions, and identifies the correct crystal structure matching the experimental product (Supplementary Fig. 2 and Supplementary Notes 1 and 2).

In terms of utility, a²c is complementary to existing crystal structure prediction methods, particularly those that start with randomized initial configurations^26,27,28,29, as it optimizes not for finding a broad range of viable structures, but for selectively finding polymorphs that are realizable in experiments via crystallization of the amorphous precursors. In Fig. 2, this distinction in utility is quantified by computing the ‘acceleration factor’, which is defined as the ratio of the precision obtained by a²c to that obtained by a random structure search (RSS)²⁸ in finding the actual crystallization products of the amorphous precursors (Supplementary Fig. 3). We find that this factor ranges from ~1.2- to ~6 fold in systems where both methods were able to identify the product. This increase in polymorph selectivity arises from a²c confining the search to the part of the energy landscape near the amorphous basin (Fig. 1), whereas RSS, by design, performs a more global structure search.

**Fig. 2: Benchmarking a²c in predicting the crystallization products of amorphous precursors.**

Encouraged by the examples discussed, we focus next on how a²c can complement experimental observations by presenting several case studies.

Case study 1

The BiBO₃ system has two metastable polymorphs, (I) and (II)³⁰. Polymorph (I) can be accessed by crystallizing the melt, and its structure has been solved using diffraction³¹. Polymorph (II) can only be synthesized by crystallizing the glassy precursor, and its structure remained unsolved until Shinozaki and colleagues³² showed that it belongs to space group C2. The lowest-energy a²c product matches the experimental diffraction pattern, providing a solution for the crystal structure of BiBO₃ (II) by including atomic positions (Supplementary Fig. 4).

Case study 2

The initially crystallizing metastable phases often have compositions close to those of their amorphous precursors, because decomposition into multiple phases with diverging compositions would be penalized by solid-state transport. Under such conditions, we anticipate decomposition of the well-known metallic glass a-Fe₈₀B₂₀ to fan out gradually to nearby metastable crystals⁵, which is found to be consistent with the phase landscape obtained using a²c (Extended Data Fig. 1). Our results show that crystallization of a-Fe₈₀B₂₀ involves topotactic decomposition to a series of related phases towards Fe and Fe₃B, in accord with the spinodal-like phase separation suggested for such glasses^5,33,34. The a²c predictions reconcile the available experimental knowledge into a phase decomposition pathway for Fe₈₀B₂₀ (Supplementary Note 3).

Case study 3

Polymorph selection by changing the morphology of an amorphous precursor is a promising but challenging research direction. Synthesis of sp³-hybridized cubic boron nitride (c-BN) from the sp²-hybridized hexagonal polymorph (h-BN) requires excessive pressures and temperatures, whereas starting with amorphous BN (a-BN) requires less extreme conditions^35,36,37. We find that the fraction of nominally a-BN subcells that crystallize into c-BN from a given a-BN cell increases with pressure, which is explainable by the continuous shift from sp² to sp³ in B–N bonding in a-BN with pressure (Extended Data Fig. 2 and Supplementary Fig. 5).

Anatase and rutile are the two polymorphs a²c predicted to selectively form from a-TiO₂ (Supplementary Fig. 3). With further analysis, we find that the frequency of subcells that crystallize into anatase versus rutile in a²c is influenced by the local ordering of TiO₆ octahedra in amorphous TiO₂ (Supplementary Figs. 6–8 and Supplementary Note 4). These findings agree with the sensitivity to the preparation conditions of the crystallization products of amorphous oxide films of multivalent cations such as Ti^38,39,40. We further confirm that the findings are consistent with all-atom simulations sampling amorphous–crystalline TiO₂ interfaces (Supplementary Figs. 9 and 10 and Supplementary Note 5), where anatase has the lowest interfacial energy with a-TiO₂, and rutile remains competitive. Unlike a²c, however, such interface simulations would quickly become impractical to run at scale for unexplored systems (Supplementary Methods). Nevertheless, a²c is not a replacement for all such simulations. For example, as a limitation, sensitivity of the product to oxygen substoichiometry or thickness in thin a-TiO₂ films, and subsequent appearance of brookite^38,40, could not be captured by a²c and requires simulations of the deposition and nucleation to understand (Supplementary Note 6). In other words, although the melt–quench used here yields reasonable amorphous configurations to apply to a²c and it is possible to vary the morphology to some degree with quench rate or density, the actual experimental structure of the amorphous precursor can vary based on the synthesis process, and an explicit modeling of that process itself may be needed in certain cases⁴¹. Similarly, capturing transitory phases on the path to the product, as in a-MgF₂ (refs. ^42,43), or understanding the real-time crystallization pathway of a-GeTe⁴⁴ requires all-atom simulations of the transformation.

Discussion

On the atomic level, the size of the atom displacements involved in the transformations observed with the a²c approach, from the amorphous precursor to atom arrangements containing ordered structures, is reminiscent of the atomic shifts occurring in displacive and dilatational polymorphic phase transitions in Buerger’s classification^18,19, which take place between polymorphs with structural relationships under limited mobility (for example, crystallization at lower temperatures). Unlike the reconstructive cases, these transitions require minimal bond breaking or rearrangement beyond distortions or rotation, and they occur via a low-energy interface due to the structural similarity of the involved polymorphs. Accordingly, these transitions proceed rapidly¹⁸, and any nucleus can grow easily into its neighborhood in the amorphous matrix. This mechanistic nature has parallels with what is modeled atomistically with a²c, hence explaining why its predictions are consistent with experimental observations in many different chemical spaces.

Our work has demonstrated that prediction of the products of amorphous-to-crystalline transformations with atomic-level explainability is possible by combining a local motif search with the general-purpose deep learning interatomic potentials in many inorganic systems. Predictive agreement with experimental polymorphs across a wide range of chemical spaces validated the computational framework. With this validation, the presented approach can be employed in research tasks such as resolving the outcome of such transitions in complex material systems, or in the in silico discovery of material polymorphs or glass-crystal composite systems that inherently have synthetic pathways starting from the amorphous state.

Methods

Generation of amorphous structures

The MQMD process in a²c, shown in Fig. 1b, aims to efficiently obtain reasonably representative amorphous configurations in a three-stage procedure. In the first stage, the N atoms comprising the system are distributed randomly in a cubic box. Except when explicitly modeling variations in density, the density of this box is determined heuristically with a rule of thumb of 15% (adopted from the survey by Cui and colleagues⁴⁵) volumetric expansion of the volume per atom of the known stable phase or averaging the same in the case of the presence of a phase mixture of that composition on the convex hull in our standard workflow. This rule of thumb is similar to those previously reported in the literature²⁵ and is reminiscent of the Lindemann criterion⁴⁶, and it yields consistent amorphous configurations with no cavities. To maintain throughput and efficiency, N is restricted to the exact integral representation of the composition closest to 100, which has been shown to yield sufficiently representative structures in inorganic systems²⁵. The random atomic positions are then briefly relaxed with a soft-sphere interaction potential to minimize the atomic overlaps. We used the soft-sphere interaction in JAX-MD⁴⁷, which includes an empirical quadratic repulsion at distances shorter than a predefined distance, and is zero elsewhere. This preparation stage ensures that the subsequent amorphous structures do not inherit fingerprints from melting of any particular crystal; instead they rapidly develop the proper short-range structural elements inherent to liquid and amorphous configurations in the equilibration and quench stages²⁵. The number of steps in overlap minimization is limited to 30, or stopped early when the shortest interatomic distance reaches above the 90% of summed ionic radii of the atoms of each pair. Limiting the number of steps in this process ensures that we do not inject any arbitrary local arrangement using the soft-sphere potential. The next stage in the workflow is a three-step MD protocol in which a graph neural network (GNN) force field is used throughout: equilibrate a melt at T_high, cool/quench this phase to T_low, and equilibrate at this temperature. T_high and T_low were manually selected for each system (Extended Data Table 1), and we ran each of the three stages of the protocol for 1,500, 2,500 and 1,500 steps, respectively. The term ‘equilibrate’ refers to constant-temperature MD runs and does not imply that an equilibrium state is reached (which would be particularly impractical and irrelevant for simulations mimicking low-temperature deposition). In MD, we use the NVT ensemble with a Nose–Hoover thermostat and a timestep of 2 fs. The atomic positions can be relaxed with a short gradient descent (GD) run in the last stage; this was kept optional as we did not find any impact on crystallization. The MQMD workflow is implemented in JAX-MD⁴⁷ and can be run in high throughput⁴⁸. The representativeness of the amorphous structures generated by our MQMD workflow was validated by comparing with experimental data in several systems (Supplementary Fig. 11 and Supplementary Note 7).

Crystallization workflow

The crystallization phase in a²c (Fig. 1c) accepts the amorphous structure model generated in the MQMD as input, and exhaustively creates all orthorhombic parallelepipeds (subcells) that can be generated from this structure, keeping their edges parallel to the unit vectors of the parent cubic box (an early example of a subcell-finding approach is provided in ref. ²⁶). These subcells are allowed to be of any length and at any position on the n_grid × n_grid × n_grid point net, where the first point is located at position (0, 0, 0), encompassing the parent structure. Our default for n_grid is 10, except for Fe₈₀B₂₀, as noted in Extended Data Table 1. The full grid has ~1.7 × 10⁵ subcells for n_grid = 10 and ~1.0 × 10⁷ subcells for n_grid = 20. Because all the degrees of freedom of each subcell are allowed to relax (as described later), we found that the present settings were able to yield a comprehensive enough set of subcells to search for crystal seeds in the systems studied here, even for finding non-orthorhombic products. However, in any system, this algorithmic step can be enhanced by adding additional degrees of freedom to the subcell creation, for instance by allowing triclinic cells or shifts.

These subcells can be filtered further using maximum atom counts per cell or stoichiometry constraints in cases where compositions distant from the parent amorphous starting material are not of interest or not plausible. For example, in Fe₈₀B₂₀, all subcells are kept up to an atom count of 16, whereas for MgF₂, subcells having this target stoichiometry and up to 21 atoms are kept. With the stoichiometry and atom count constraints, the subcell set size per system often falls in the range of 10³–10⁴ for most systems with n_grid = 10, and is one to two orders of magnitude larger for n_grid = 20. Details on each system studied are provided in Extended Data Table 1. Next, all geometric degrees of freedom of these subcells, including atomic positions and the cell size and shape, are optimized using an NequIP-based equivariant GNN potential (see ‘Deep learning methods’) to describe the interatomic interactions, as well as the fire algorithm⁴⁹ for at least 150 steps. The lowest-energy structures found by this procedure at each composition are the selectively predicted stable or metastable crystallization products, which should match the experiments. Matching is confirmed in each case, first by a combination of visual inspection and pymatgen’s structure matching and/or space group analyzer methods⁵⁰, at various tolerance settings. Our lowest-energy bound is strict; that is, the experimental crystallization products were either the absolute lowest-energy a²c predictions or practically degenerate within a few milli-electronvolts per atom, which is on par or lower than the variability in density functional theory (DFT) convergence and/or structure optimization. For the studied systems, we provided the generated amorphous structure, the subcells extracted using a²c, the index of a particular initial subcell from the amorphous configuration that subsequently relaxed into the target phase, the final force-field structure, and, for further validation, a DFT-obtained final structure of the exact relaxation path starting from the given initial subcell. DFT relaxations were performed with MP2020 settings⁵¹ (see ‘Data availability’).

Deep learning methods

General-purpose machine learning interatomic potentials for various deep learning architectures are emerging to model condensed phases^{22,23,24,52,53,54}, with certain architectures being chemically universal (that is, a single potential is applicable to the entire periodic table and under a variety of modeling conditions). The particular instance of the general-purpose (commonly referred to as ‘universal’) deep learning potential used in this work was developed as part of the GNoME project²⁴. The potential adopts the equivariant GNN architecture of NequIP⁵⁵, and uses three layers of message passing, even irreps up to l_max of 2 with multiplicities of 64, 32 and 16 for l = 0, 1 and 2, and a two-hidden-layer radial multilayer perceptron with 64 neurons acting on a radial basis of eight Bessel functions. We used a local interatomic cutoff radius of 5 Å in graph generation. This particular model has ~2.4 million parameters trained on energies and forces, with a loss function composed of equally weighted Huber losses where the energy term is normalized per atom and the force term is normalized by the total number of force components, 3N, each with a Huber-loss delta of 0.01. The potential was trained using the Adam optimizer with a batch size of 64 structures and a learning rate of 5 × 10⁻³ using a linear learning rate decay. We trained our potential on a dataset of ~7.4 million DFT calculations. The dataset is made up of structures across a wide array of chemical spaces, sampled from ab initio random-search-generated structures²⁸ and substitutionally generated structures as described in the paper outlining the GNoME project²⁴. From the substitution data, the first, second and twentieth steps of each relaxation trajectory are sampled. The particular model reaches a test error of 0.04 eV per atom on the MatBench Discovery benchmark⁵⁶. We further showed that this potential is robust in modeling amorphous systems, achieving a testing error of 38 meV per atom with respect to DFT for thousands of amorphous structures encompassing nearly all binary chemical systems⁴⁸. We found this potential to be accurate in predicting the relative stabilities of polymorphs in the chemical systems relevant to this work, with an estimated error of 9 meV per atom and a stability ranking correlation coefficient of 0.8 or above (Supplementary Fig. 12 and Supplementary Note 8). This same potential was used in all the results reported in this work, including both MQMD and crystallization workflows. All MQMD and crystallization simulations were run on a single P100 or A100 graphics processing unit.

Crystal structure prediction with random structure search

Global exploration methods that utilize random configurations as the starting points of their complex search procedures have a long history in crystal structure prediction^{26,28,29,57,58}. Here we used a recent and widely adopted version—random structure search (RSS)—and the respective AIRSS package developed by Pickard and colleagues²⁸ to generate around 10,000 structures with bond-distance constraints at target compositions, and relaxed them with our force field. We applied a two-step process to categorize the resulting structures. First we checked if they could be matched to one of the known polymorphs that may nucleate from the amorphous atom configurations using pymatgen’s StructureMatcher⁵⁰, and, if not, we categorized them with respect to their space group. Precision was calculated as the true positive count (in this case, this value is 1 if the known crystallization product is found by the method, or 0 otherwise) divided by the number of all unique polymorphs found within the energy threshold. We used an energy threshold of up to 0.05 eV per atom of the lowest-energy crystal in each system when counting polymorphs. For TiO₂, the number of unique polymorphs found with RSS is roughly similar to the earlier machine learning-based RSS work by Reinhardt and colleagues⁵⁹.

Structural analysis and post-processing

Pymatgen was used extensively to process and analyze amorphous and crystal structures⁵⁰. VESTA was used to generate the structure visuals⁶⁰. For smoothness, radial distribution functions (RDFs) and the bond angle distributions for the amorphous structures were averaged over the last 200 steps of the low-temperature equilibration run, in ten-step increments.