Accelerated identification of equilibrium structures of multicomponent inorganic crystals using machine learning potentials

The discovery of multicomponent inorganic compounds can provide direct solutions to scientific and engineering challenges, yet the vast uncharted material space dwarfs synthesis throughput. While the crystal structure prediction (CSP) may mitigate this frustration, the exponential complexity of CSP and expensive density functional theory (DFT) calculations prohibit material exploration at scale. Herein, we introduce SPINNER, a structure-prediction framework based on random and evolutionary searches. Harnessing speed and accuracy of neural network potentials (NNPs), the program navigates configurational spaces 102–103 times faster than DFT-based methods. Furthermore, SPINNER incorporates algorithms tuned for NNPs, achieving performances exceeding conventional algorithms. In blind tests on 60 ternary compositions, SPINNER identifies experimental (or theoretically more stable) phases for ~80% of materials. When benchmarked against data-mining or DFT-based evolutionary predictions, SPINNER identifies more stable phases in many cases. By developing a reliable and fast structure-prediction framework, this work paves the way to large-scale, open exploration of undiscovered inorganic crystals.

stable) phases for ~80% of materials within 5000 generations, entailing up to half a million structure evaluations for each composition.When benchmarked against previous data mining or DFT-based evolutionary predictions, SPINNER identifies more stable phases in the majority of cases.By developing a reliable and fast structure-prediction framework, this work opens the door to large-scale, unbounded computational exploration of undiscovered inorganic crystals.
Human history has evolved together with material innovation: steel production from heating iron with carbon triggered a shift from the Bronze to Iron age, and the growth of high-purity Si crystal ingots led to burgeoning of the computer age.In modern times, synthetic multicomponent materials are constantly developed to meet the demands of diverse applications.The Inorganic Crystal Structure Database (ICSD), including most of the experimentally synthesized inorganic compounds, has approximately 200,000 materials registered to date, and the data entry steadily increases by ~5,000 every year 1 .The sheer size and active expansion of the database reflect that new materials continue to drive scientific and engineering advances in fields such as electronics, energy harvesting/storage, and high-Tc superconductors [2][3][4][5][6][7] .
Despite the vast material library available today, it is far from being complete in ternary or higherorder (simply multinary hereafter) phases.Based on a rough estimate, only approximately 16% and 1% of ternary and quaternary compounds, respectively, are at least partially revealed 8 .This inspires reasonable hope that valuable materials can be discovered in the largely unexplored multinary domain.Furthermore, the current material repositories are chemically and synthetically biased.For instance, elemental occurrences within ternary compounds peak at oxygen with 22,476 counts in the ICSD, which is more than three times that of the next most frequent elements (Fe, Si, and S).This is mainly because oxygen is the most earth-abundant element and forms stable compounds with most metals.Thus, oxides are easier to synthesize than other compounds and benefit from synthetic recipes established throughout the long history.This indicates that the present material database is biased toward those with facile synthesis, possibly missing promising compounds that are unfamiliar today 9 .
Considering the very large gap between current experimental throughput and the number of unknown materials, together with rising synthesis barriers, exploring the uncharted chemical space solely by experiment would be inefficient.Alternatively, experimental endeavors can take advantage of computational prescreening based on the density-functional theory (DFT) calculations.This has been demonstrated by a multitude of recent publications in which the discovery of new materials was accelerated by DFT: cathodes for Li-ion batteries 10 , nitride semiconductors 11 , metal nitrides 9 , 18-electron compounds 12 , and boron-based MAX phases 13 .Throughout these works, DFT results were able to suggest compositions that are presumably stable under ambient conditions and thus have high synthesizability.In addition, DFT predictions of basic properties could steer experimental resources to materials appropriate for specific applications.
Despite great promise, the current computational exploration of unknown materials inherently suffers from low throughput; in contrast to high-throughput screening of materials cataloged in the ICSD 14,15 , the investigation into as-yet-unreported materials should start from the crystal structure prediction (CSP), preferentially for equilibrium phases with the lowest Gibbs free energy (evaluated at the DFT level) under given chemical composition and thermodynamic conditions.However, CSP is a classic NP-hard problem in which the computational load in identifying the global minimum exponentially increases with material complexity 16 .Compounded further by the high costs of DFT calculations, this results in extremely low throughput of DFT-based CSP.
Consequently, computational investigations are often limited to known prototypes 9,10 or short evolutionary steps 11,12 , which risks resulting in metastable structures instead of the ground state.
Reflecting this, the compositional territory has been scanned rather 'conservatively' over subsets close to known chemical families.If free energies can be evaluated much faster than with DFT, then more aggressive and far-reaching exploration will be viable on a large scale, raising the chance of discovering materials with novel functionalities.
Recently, machine-learning potentials have attracted considerable attention because they can provide DFT-level energies at a fraction of the cost, which finds immediate applications to CSP 17- 20 .In applying machine-learning potentials to CSP, however, one is faced with difficulties in choosing training sets without information on crystal structures.We recently suggested a way to resolve this problem by using melt-quench molecular dynamics (MD) simulations 21 .The liquid simulation can self-start from a random distribution with rapid equilibration; thus, a priori information on the crystal structure is not required.In addition, the MD trajectory automatically samples diverse local orders that may appear in stable or metastable crystalline structures.We demonstrated that the Behler-Parrinello-type neural network potential (NNP) 22 trained with the melt-quench-annealing trajectories can serve as a high-fidelity surrogate model in CSP.
Going beyond the previous achievement, we herein develop a structure prediction framework combining NNPs with evolutionary or random searches, named SPINNER (Structure Prediction of Inorganic crystals using Neural Network potentials with Evolutionary and Random searches).
Free from any empirical knowledge on material structures, the program identifies the global minimum in a brute-force style by harnessing the accuracy and speed of NNPs.In blind tests on ternary compounds with significant complexity and diverse crystal symmetries, SPINNER successfully identifies the experimental (or theoretically more stable) phases for ~80% of the cases.
The program also outperforms other favored approaches in the majority of test materials.The average computational throughput is ~4 days per one composition on a 36-core node, which is estimated to be 10 2 -10 3 times faster than pure DFT-based approaches.The outstanding performance of SPINNER allows for large-scale, unbounded computational exploration of undiscovered inorganic crystals.
The basic workflow of SPINNER is schematically presented in Fig. 1(a).In brief, for an input chemical composition (elements and stoichiometry), SPINNER first carries out a melt-quenchannealing simulation and trains an NNP with MD trajectories.To enhance the accuracy for ordered phases, the NNP is iteratively retrained over low-energy structures in the refining-state CSP (the upper part of Fig. 1(a)).In the main CSP proceeding up to 5000 generations (the lower part of Fig. 1(a)), SPINNER collects low-energy candidate structures within 50 meV atom −1 , which are finally sorted after full relaxations at the DFT level.For DFT calculations, we use the Vienna Ab initio Simulation Package (VASP) 23 with the Perdew-Burke-Ernzerhof (PBE) functional for exchangecorrelation energies 24 .Figure 1b shows a schematic of the present evolutionary algorithms that are based on random generation, crossover, permutation, and lattice mutation.The random generation is heavily used over mutations, which we find to be instrumental in multi-component systems.
Detailed descriptions of DFT calculations and evolutionary algorithms are provided in the Methods section.We stress that no empirical information (for instance, prototype data 25 , bonding topology 26 , or Pauling's rules 27 ) is used in the present CSP scheme.

Results
Selection of test materials.In benchmarking search algorithms tackling NP-hard problems, the global minimum is usually unknown.The situation is slightly different in CSP because the equilibrium structures that are experimentally resolved by diffraction analysis mostly equate to the global minimum for the given chemical formula and thermodynamic conditions.Therefore, blind tests with the compositions reported in the ICSD can serve as an ultimate evaluation of the performance of a CSP algorithm.Here we focus on ternary compounds because the corresponding database spans a wide range of chemistries and structural prototypes.
Within the ICSD, we first select ternary compounds measured at room temperature or below and at atmospheric pressure, as we are mainly concerned with materials that are stable under ambient conditions.We also consider only high-quality (R < 0.1) ordered crystals with welldefined structures and rule out molecular crystals as well as compounds including 3d transition elements (V-Zn), lanthanides and actinides.In the latter materials with 3d and f electrons, the magnetic ordering influences the free energy, but the current implementation of NNPs has yet to resolve the fine energy scales of magnetic interactions.(We note that meaningful developments are underway 28,29 .)When distinct crystal structures with the same composition are available in the ICSD, the most stable structure within the PBE functional is regarded as the reference.Among the filtered structures, we randomly select 50 materials under the condition that at least one crystal is selected from the 32 crystallographic point groups (omitting 6/m due to zero occurrence), which secures diversity of the test structures.To exclude structures that are too simple, we choose compounds with the formula units (Z) in the primitive cell being at least 4. (So, the numbers of atoms are at least 12).We additionally handpick 10 structures to further diversify local motifs and chemistry.The full list of the selected 60 materials is provided in the Supplementary Information.
Regarding the band gap (Eg), there are 18 metals, 13 semiconductors (0 < Eg ≤ 2 eV), and 29 insulators (Eg > 2 eV).(The finite band gaps are calculated within one-shot hybrid functional calculations 30,31 .)Most structures have a hull energy of 0 in the Materials Project database 32 (see Supplementary Table 1).

Structure prediction by SPINNER.
In running SPINNER, we assume the same Z value as in the ICSD structure.Fig. 2(a) shows ΔEmin, the offset of the PBE energy of the most stable structure after 5000 generations of structure search from that of the experimental structure.Since the pool size is 24-80 (see the Methods section), up to half a million structure relaxations are performed for each composition.The color code of the data indicates the earliest generation at which the minimum energy structure is identified (Ng).For 75% of the compositions (45 of 60), SPINNER predicts the reference (ΔEmin = 0) or a lower-energy (ΔEmin < 0) structure, which is mostly found within 1000 generations (38 of 45).The largest error occurs for Sr2Pt3In4 with a ΔEmin of 36 meV atom −1 .Figures at the bottom of Fig. 2(a) display the unit cells of some materials with ΔEmin = 0.In Fig. 2(a), six materials have negative ΔEmin down to −16 meV atom −1 , meaning that the structure found by SPINNER is more stable than the experimental phase within the PBE functional.
We note that both structures share almost the same local order, differing in ranges beyond ~4 Å .
Since the ICSD does not register any other experimental structures for these compositions, the energy ordering obtained by the PBE functional is likely incorrect, calling for further investigation.
In ref. 33 , the PBE functional was found to incorrectly stabilize metastable phases for some binary materials, which was partly addressed by the SCAN functional 34,35 .The empty squares in Fig. 2(a) are the ΔEmin's recalculated by the SCAN functional.Except for PbOsO3 and LiYSn, the ΔEmin values becomes positive.This confirms that the SCAN functional is more accurate than PBE in ordering low-energy phases.Since heavy elements such as Tl, Pb, and Bi possess large spin-orbit coupling (SOC), we also recalculate energies for Tl3PbCl5 and PbOsO3 with SOC (see sun crosses) and find that the ΔEmin for Tl3PbCl5 and PbOsO3 increases to −0.2 and 16 meV atom −1 , respectively.Thus, the incorrect energy orderings are mostly rectified by introducing more sophisticated energy functionals.(We note that the PBE functional without SOC still correctly reproduce the equilibrium phases in many compositions including Tl, Pb, and Bi.) Except for Na3PS4 and KAlCl4, experimental structures are not found among the final candidate structures for materials with ΔEmin < 0.Even though the energies of the ICSD structures are close to the global minimum, they compete with other metastable structures populating at similar energies, lowering the chance of appearing during given generations.This indicates the importance of the appropriate exchangecorrelation functional in preparing training data to efficiently identify the ground state.However, the SCAN functional or SOC significantly increases the computational time for the DFT part.
We assess the inference accuracy of NNP based on two metrics: The first is the absolute energy difference (ΔE0) between DFT and NNP for the experimental structure (relaxed within each method), which relates to how well NNP predicts the structure and energy of the reference structure.
The second metric (ΔĒ) is similar to the first one but it is averaged over final candidates within bottom 50 meV atom −1 .ΔĒ indicates how accurately NNP ranks the energies of stable and metastable phases.ΔE0 and ΔĒ are provided for every material in Supplementary Table 1.When averaged over materials with ΔEmin ≤ 0, ΔE0 and ΔĒ are 12.9 and 11.8 meV atom −1 , respectively, meaning the potential energy surfaces of DFT and NNP agree well around the equilibrium and low-energy metastable structures.This confirms that the trained NNPs are good surrogate models of DFT.The remaining errors are partly attributed to a low resolution in describing medium-to long-order beyond cutoff radii of descriptors, which are 6 and 4.5 Å for radial and angular parts, respectively.This can be improved by increasing the cutoff radius, but it adversely affects the computational efficiency due to rising costs in evaluating descriptors.The analysis on the failed cases (ΔEmin > 0) is given in the Discussion section.
To understand how SPINNER effectively addresses the NP-hardness, we analyze the evolutionary steps leading to the experimental structure for materials with ΔEmin ≤ 0 based on two parameters, Nt and Nm.As schematically shown in Fig. 2(b), Nt counts the generations from random seeding to the appearance of the equilibrium structure at Ng. Nm indicates the number of mutations within Nt steps.The distribution of Nt in Fig. 2(b) shows that Nt is 0 for almost half of the cases, meaning that the ground-state structure is obtained directly from relaxing a random structure, without any mutations involved.Even in cases undergoing finite mutations, Nm is mostly within 5.
Intriguingly, the randomly generated structures are very close to the global minimum.There are two reasons contributing to this.First, pair-wise minimum distances obtained from melt-quench-annealing trajectories are imposed as constraints in the structure generation and relaxation (see the Methods section).This condition filters out ~99% of the randomly generated structures that are unlikely to relax into low-energy structures (see the decision tree in Fig. 1(b)).It also effectively protects atomic configurations from evolving into unphysical structures during relaxation.Note that the minimum distance constraints often extend to 2-3 Å , which are much stronger than simple conditions preventing too short bonds.Second, we find that when the random structures leading to the global minimum are relaxed by DFT instead of the NNP, more than half of them relax to high-energy metastable structures instead of the global minimum.This implies that the energy landscape around the global minimum is different between DFT and the NNP (see Fig. 2(c)).That is, DFT adaptively forms chemical bonds that can locally stabilize the starting configuration; thus, the search can be stuck at a high-energy metastable structure.Having not machine-learned such chemical bonds, the NNP unwittingly smooths out the corresponding energy region, directing the structure to the global minimum.Figuratively, the NNP "catalyzes" the relaxation process to the equilibrium, eliminating or lowering energy barriers around certain metastable structures.This also effectively increases the volume of the configurations that can funnel down to the equilibrium structure.the lowest energy is found by SPINNER.The structures are adopted from ref. 9 (Li2HfN2, Sc2C3N6, Hf2SeN2, Na3OsN2), ref. 36 (Ti3O3N2, Zr3O3N2), ref. 37 (Li2SnN2, Li2SnF6, Li2TiN2), ref. 38 (KScS2), ref. 12 (HfIrP, ScPtBi, AlPtSb), ref. 39 (Hf3PB4, Zr2InB2), ref. 40 (WMo5B6, W4Mo4B15), ref. 41 (Sn5S4Cl2, Cd4SF6), and ref. 42  Comparison with other approaches.To benchmark SPINNER against other CSP methods, we select ternary structures from the literature that were theoretically predicted by either data-mining known prototypes 9,12,[36][37][38] or using DFT-based evolutionary approaches such as a genetic algorithm [39][40][41] or particle swarm optimization 42 .The materials are listed in Fig. 3(a), and understandably, they are all non-oxides (nitrides, borides, etc).None of these compositions has ever been synthesized as far as we are aware except for KScS2, Sc2C3N6, and ScPtBi (see below).
For these materials, we perform CSP with SPINNER up to 1000 generations with Z ranging from 2 to 8. The plot at the top of Fig. 3(a) shows ΔEmin in reference to the structure in the literature.
(The structures from the references are fully relaxed within the present DFT method.)SPINNER identifies lower-energy structures for the majority of cases (13 of 21) and the same structures for the rest.In ref. 41 , CSP was performed with the PBEsol functional 43 instead of PBE for Sn5S4Cl2 and Cd4SF6.As a comparison, we evaluate ΔEmin with PBEsol and find that compounds identified in this study are still more stable than those predicted in the reference by 31 and 4 meV atom −1 for Sn5S4Cl2 and Cd4SF6, respectively.The lower part of Fig. 3(a) shows Z values of the primitive cells with the lowest energies (solid squares) in comparison with those in the references (solid circles).Lower energies are found mostly at Z values higher than those in the literature.In addition, SPINNER finds the lowest-energy structure often at a larger Z than that of the identified primitive cell.This indicates the importance of trials with diverse Z numbers.The NNP is far more advantageous in varying Z than DFT owing to the linear scaling with respect to the number of atoms, in contrast to the cubic scaling of DFT.
While structures found by SPINNER have usually local structures similar to reference structures (such as KScS2 and W4Mo4B15 in Fig. 3(b)), different short-range orders appear in materials such as Li2TiN2, Cd4SF6 and TaCN3 (see Fig. 3(b)).We note that Cd4SF6 and TaCN3 were discovered by DFT-based evolutionary searches 41,42 .The failure of the previous works in identifying the present lower-energy structures might be attributed to small Z numbers and/or few generations.In Fig. 3(a), Na3OsN2 shows the biggest energy drop.Both structures are similar in local order, but the reference structure is distorted from that found by SPINNER.We comment that metastable structures can have distinct physical properties compared to the ground state.In the example of HfO2, the high-temperature tetragonal phase that is metastable by 57 meV atom −1 at 0 K possesses a far higher dielectric constant (70) than that for the low-temperature monoclinic phase ( 16) 44 .This underscores the importance of identifying the true global minimum for accurate prediction of material properties.
For the 13 structures newly identified in Fig. 3(a) by SPINNER, we check the existence of corresponding prototypes in the ICSD using AFLOW-XtalFinder 45 .We find that 10 of 13 materials do not have any matching prototypes (see star marks in Fig. 3(a)).This suggests that the current prototypes for the ternary phase are not sufficiently cataloged, which is understandable because materials such as nitrides have been synthesized much less frequently than oxides.Notably, the prototype of the SPINNER-identified Zr3O3N2 and Ti3O3N2 is commonly Ti3O5.We think that ref. 36 missed this prototype because they did not fully consider partial ion-exchanges replacing O with N. In passing, the three compositions for which SPINNER and previous evolutionary searches agree have corresponding prototypes in the ICSD.
Sc2C3N6 was recently synthesized with the structure predicted by data mining 46 , which was confirmed by SPINNER.For KScS2, there was an experimental report 47 preceding the theoretical prediction, but it was not recognized by ref. 38 .The experimental crystal structure is consistent with the present work, which is more stable than that for ref. 38 by 2.5 meV atom −1 (see Fig. 3(b)).
The crystal structure of ScPtBi was predicted and confirmed by synthesis in ref. 12 .Although it was synthesized as a multiphase, SPINNER does not identify any phases within 50 meV atom −1 other than the one predicted in ref. 12 .(Possibly, SOC would be necessary due to Pt and Bi.)

Discussion
Regarding the computational cost, the whole computation from the melt-quench-annealing MD to 5000 CSP steps takes 3-5 days on a 36-core Amazon ® server (CPU of Intel ® Xeon Platinum 8000series).(For 1000 evolution steps, it takes 2-3 days.)This excludes human time in managing data flow and making decisions, which can be executed instantly by automation (under development).
On average, the workload can be split into the DFT-MD (25%), training of the NNP (10%), SPINNER (60%), and DFT relaxation of crystals (5%).The computational time for the DFT part varies widely depending on the number of electrons.The SPINNER part is highly scalable on massively-parallel computers.As a specific example of LiWCl6 tested above, it takes 84 hours for 5000 generations with the pool size of 64 and Z = 4.When tested with the same computational resources and conditions, a DFT evolutionary program 48 could proceed up to only 6 generations under the setting suggested in the manual.Therefore, the entire 5000 generations would take several years with DFT.(We add that actual generations required for identifying the equilibrium phase can be different between DFT and NNP.)In terms of monetary price, the structure prediction of one material by SPINNER costs approximately 200 US dollars based on the current pricing policy of the Amazon ® server.This estimation indicates that it would be viable to construct databases of ~1000 materials targeted for a specific application at a reasonable cost and in a reasonable time frame.Despite impressive performance, SPINNER failed with 25% of test ternary materials, which merits further discussion.First, ΔE0 and ΔĒ averaged over materials with ΔEmin > 0 (see Supplementary Table 1) are 41.0 and 43.3 meV atom −1 , which are substantially larger than 12.5 and 11.9 meV atom −1 for structures with ΔEmin ≤ 0. Most notably, the NNPs for SnGeS3 and Sr2Pt3In4 show the largest ΔE0 (ΔĒ) of 75.5 (227.3) and 83.5 (140.0)meV atom −1 , respectively.
For these materials, locating the equilibrium structures would be very difficult, as they are high in the energy scale.The origin of the poor NNP quality requires further investigation.(Increasing the cutoff radii of descriptors was not helpful.)Second, among the failed cases, qualities of the NNP are acceptable in several materials, with accuracy metrics on par with those with ΔEmin ≤ 0. In these materials, the experimental structures would be found in principle by extending the evolutionary process beyond 5000 generations.We note that BaGe2S5, Na3SbO3, and YPdGe have rather flat or pointed cell shapes.Such low-dimension-like geometries are entropically low with small chances to occur in the random generation.Since conventional unit cells for these materials have more 3-D-like structures than the primitive cells, we additionally carry out a structure search with Z equal to that of the conventional cells (see Supplementary Table 1).In all three cases, SPINNER successfully identifies experimental or theoretically more stable structures (see diamonds in Fig. 2(a)).This again highlights the importance of multiple trials at different Z numbers.
Even though the training set is constructed within a specific stoichiometry, dynamic fluctuations during the melt-quench-annealing process extend transferability of NNPs.For example, employing the NNP developed for Mg2SiO4, we try CSP for MgSiO3 and SPINNER successfully identifies the experimental crystal structure (ICSD-ID of 196432).(Note that both materials have the same valence states for every element.)Since the NNP training set has a single stoichiometry, atomicenergy offsets among chemical species become arbitrary, resulting in constant energy shifts between the DFT and NNP energies for crystalline structures of MgSiO3 49 .Nevertheless, energy ordering among the low-energy structures is well maintained, leading to the successful identification of the equilibrium phase.Such extended transferability would be advantageous in saving training sets for a variable-composition crystal search 40 , which would require multiple MD simulations covering a range of compositions.
To test SPINNER on compounds other than ternary phases, we carry out similar studies on B, TiO2, P3N5, NbPd3, Li10GeP2S12 2 , and InGaZnO4 3 whose equilibrium structures are experimentally known.In all cases, ground-state structures are found within or around 5000 generations.In detail, B is well known for a complicated bonding nature and rich polymorphism 17,50 .The ground-state structure (α-B12) is found at 5055 generations for Z = 12 (primitive cell), while it is found at 748 steps for Z = 36 (conventional cell).Similar to a previous study 50 , the training error is relatively large (~30 meV atom −1 ) but the energy ordering among low-energy structures is good enough to predict the equilibrium phase correctly.For TiO2, we try Z = 8 and obtain all the major polymorphs such as rutile, anatase, and brookite among the final candidate structures.Even though the test sets are small compared to the ternary compounds, these results support that the present CSP framework works regardless of material complexity.However, quaternary or higher-order compounds may generally require longer generations than in the present work, which needs further investigation.melt-quench-annealing trajectories, we carry out CSP for 50 generations and select 10 structures: 5 lowest-energy structures and 5 structures with the lowest antiseed weights (see below) within bottom 100 meV atom −1 .These structures are further relaxed within DFT by using AMP 2 30 .The relaxation terminates when atomic forces and stress tensors are less than 0.1 eV/Å and 20 kbar, respectively.The atomic positions and energies are sampled during the DFT relaxation, which are augmented to the training set for refining the NNPs.We find that the initial cell structure sometimes undergoes significant deformations, which necessitates adapting the k-point mesh consistently.To this end, the k-point mesh is checked and modified every 10 relaxation steps according to the converged k-point spacing.Then, we continue CSP and retrain the NNPs at 100, 200, and 400 generations in similar ways but excluding structures already learned in the previous iterations.In total, 40 crystal structures are sampled to refine the NNPs.

Workflow of crystal structure prediction.
With the refined NNPs, we perform up to 5000 generations of the main CSP (the lower part of Fig. 1(a)).The structures are generated by random seeding, permutation, lattice mutation, and crossover algorithms.Detailed descriptions of each mode are provided in the next subsection.The generated structures are relaxed within the NNPs using the LAMMPS code 54 .We use the restraint option during structure relaxation, which prevents atom pairs from being closer than a certain distance limit through repulsive harmonic forces.The pair-wise distance limits are set to the minimum values found for the corresponding pair during the melt-quench-annealing process.The restraint option effectively prohibits the structures from evolving into untrained domain.In comparison, the structural relaxation by DFT methods typically takes 10 2 -10 3 times that of the NNP.
To prevent unphysical structures from appearing as candidate structures, we compare the singleshot DFT energy and NNP energy of the lowest-energy structure every 1000 generations.If the discrepancy is greater than 50 meV atom −1 , then the NNP is retrained over 10 structures selected similarly to the refining stage (see above).Last, we perform DFT structure relaxation for final candidates that lie within bottom 50 meV atom −1 using the AMP 2 .For accurate evaluation of energies, we adopt tight convergence criteria of 2 meV atom −1 and 3 kbar for energy and pressure, respectively.
Evolutionary algorithms.Following ref. 48, SPINNER utilizes evolutionary algorithms to find the global minimum in the configuration-energy space under the given composition and Z number.
From extensive tests, we find that random generation, crossover, permutation, and lattice mutation algorithms are particularly effective in the structure prediction of multinary phases, and so the current version of SPINNER is based on these four generation modes.In generating random structures, the space group, lattice vectors, and Wyckoff positions are selected randomly using the RandSpg package 55 .The atomic density in the first generation is set to that of the amorphous phase generated by the DFT-MD simulations.Then, the volume of each structure in the ith generation (i > 1) is chosen randomly between 70% and 130% of the volume of the lowest-energy structure in the previous generation.Distance constraints are imposed on all pairs of atoms such that distances between pairs are longer than the minimum values found for the corresponding pair during the melt-quench-annealing process.In applying the crossover algorithm, we note that the conventional approach often unnecessarily disrupts stable chemical units (for instance SiO4 tetrahedra in Mg2SiO4).To avoid this, SPINNER utilizes atomic energies of the NNP in choosing cutting planes such that the average atomic energy within the slab is sufficiently low.The average atomic energies are also used in adjusting lattice vectors and atomic registries in combining the two slabs.This effectively conserves chemically-stable local units.
During the evolutionary steps, the structural similarity is measured by the partial radial distribution function (PRDF) 56 and if two structures are determined to be too similar, then one of them is discarded from the pool.We also adopt an antiseed algorithm to diversify the structure pool 57 .The antiseed weight of ith structure (Ai) is defined as follows: where the summation runs over structures in the inheritable pool (bottom 200 meV atom −1 in the refining stage and 100 meV atom −1 during the main CSP), dia is the distance between the ith and ath structures measured by PRDF, and σ is a Gaussian width.Unlike the original scheme 57 , we use the antiseed weight in defining the probability (Pi ∝ exp(−Ai/σA) where σA is a parameter) for a structure to be inherited in the next generation by mutation (see below).
The number of structures generated by the operations is set to twice the number of atoms in the simulation cell.The inherited structures are chosen within bottom 200 and 100 meV atom −1 for the refining-stage and main CSP, respectively, and we choose half of them randomly and the rest are chosen according to the probability distribution of {Pi} defined in the above.In the refining stage, ratios of operations of random generation, crossover, permutations, and lattice mutations are 30%, 50%, 10%, and 10%, respectively.The crossover tends to generate diverse local orders, which is useful in training robust NNPs.However, the algorithm is less effective in finding the ground state due to the low symmetries of the generated structures, so we exclude the crossover operation in the main CSP in which random (70%), permutation (20%), and lattice mutation (10%) are employed.In addition to the structural pool generated by these operations, we carry low-energy structures in the bottom 100 and 50 meV atom −1

Fig. 1 .
Fig. 1.Schematic workflow of SPINNER.(a) Schematic depiction of the whole workflow of CSP.Top left: An initial NNP is trained over disordered structures sampled from melt-quench-annealing trajectories of a given composition obtained by DFT calculations.Top right: The initial NNP is iteratively refined over low-energy crystal structures obtained during 400 generations of the structure search.Bottom left: With the refined NNP, the structure search is carried out up to 5000 generations.The quality of the NNP is monitored every 1000 generations.Bottom right: Free energies of final candidates are evaluated by DFT calculations.(b) Left: Representation of the evolutionary

Fig. 2 .
Fig. 2. Search results of SPINNER against structures in the ICSD.(a) Distribution of the energy difference between the most stable structure found by SPINNER and the ICSD structure (ΔEmin).The color codes indicate the generation at which the lowest-energy structure was identified (Ng).The circle, square, and sun cross indicate the exchange-correlation functional and diamonds are results obtained with the conventional cell.Equilibrium structures of some compositions with ΔEmin = 0 are displayed below.(b) Definition and distribution of Nt and Nm for materials with ΔEmin ≤ 0; Nt denotes the number of generations that elapsed from random seeding until the ground state is identified, and Nm is the number of mutations involved in this process.Nt = 0 means that the ground state is directly obtained by relaxation from a random structure.(c) Schematic illustration of potential energy surfaces of the DFT and NNP.

Fig. 3 .
Fig. 3. Benchmark results of SPINNER against other methods.(a) Upper part shows energy difference between the structure predicted by SPINNER and the reference structure (ΔEmin) predicted by data mining or evolutionary algorithms.The lower part compares the formula unit in the unit cell (Z) between the reference and the value at which (Hf2CN4, TaCN3).The chemical formulas with star marks do not have corresponding structural prototypes in the ICSD.(b) Examples illustrating the structural difference between the crystal structures in references (top) and ones identified by SPINNER (bottom).
to the next generation in the refining stage and main CSP, respectively.50.Huang, S.-D., Shang, C., Kang, P.-L.&Liu, Z.-P.Atomic structure of boron resolved using machine learning and global sampling.Chem.Sci.9, 8644-8655 (2018).For the symbols under SPINNER, we refer to the main text.Under ΔEmin we write the energy difference between the most stable structure found by SPINNER within ICSD structure primitive cell size and the ICSD structure calculated by PBE functional.When needed, we additionally write the energy difference calculated by SCAN functional (marked by †), or by PBE when the spinorbit coupling is considered ( ‡).Also, the numbers with the * mark indicate the energy difference between the most stable structure found within ICSD structure conventional cell size and the ICSD structure.The unit for Ehull, Ng, ΔE0, ΔĒ, and ΔEmin are meV atom −1 .