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

The study addresses the bottleneck in discovering new multinary inorganic crystals: the combinatorial explosion of possible structures makes crystal structure prediction (CSP) NP-hard, and routine CSP with density functional theory (DFT) is too computationally expensive to scale. Despite large databases like ICSD, experimental coverage of ternary and higher-order compounds remains sparse and biased toward easily synthesized chemistries (notably oxides). Computational prescreening with DFT has aided discovery across many domains, but throughput is constrained when starting from unknown structures. Machine-learning potentials (MLPs), which can approach DFT accuracy at orders-of-magnitude lower cost, offer a path to accelerate CSP; however, training them without prior structural knowledge is challenging. Building on prior work that used melt–quench molecular dynamics trajectories to self-generate diverse local environments for training Behler–Parrinello neural network potentials (NNPs), the authors develop SPINNER: a CSP framework coupling NNPs with random and evolutionary searches. SPINNER aims to reliably find equilibrium (lowest free-energy) structures for multicomponent systems rapidly and without empirical structural priors, thereby enabling broader, faster exploration of uncharted materials space.

Prior computational discoveries leveraged DFT to propose synthesizable compositions and basic properties (e.g., battery cathodes, nitride semiconductors, metal nitrides, MAX phases, high-Tc superconductors). However, these typically rely on known prototypes or limited evolutionary steps due to DFT cost, risking metastable outcomes. MLPs/NNPs have been applied successfully to structure prediction of crystals, nanoclusters, and surfaces, dramatically reducing cost while retaining accuracy. A key challenge is building training sets absent prior knowledge of crystal structures. The authors’ earlier work showed that melt–quench–anneal DFT-MD trajectories can self-start from random configurations, rapidly equilibrate liquids, and sample diverse local orders relevant to crystalline phases; NNPs trained on such disordered data can then serve as accurate surrogates for CSP. Other CSP approaches include genetic algorithms (USPEX), particle swarm optimization (CALYPSO), and data-mining prototype substitution, each with strengths but often limited by DFT throughput or prototype coverage, especially for non-oxides where databases are sparse.

Overview: SPINNER (Structure Prediction of Inorganic crystals using Neural Network potentials with Evolutionary and Random searches) performs CSP for a given composition and formula units Z by (1) training NNPs using DFT-MD melt–quench–anneal trajectories, (2) iterative NNP refinement using low-energy crystalline structures from short CSP runs, and (3) a main CSP up to 5000 generations using NNP-accelerated relaxations with evolutionary/random operations. Final low-energy candidates are validated and ranked by full DFT relaxations. Training NNPs: Initial structures are random atomic distributions (~80 atoms) at the target composition. DFT-MD (VASP, PBE) performs: premelt at 4000 K for 4.5 ps; liquid MD for 16 ps at an ad hoc melting temperature determined per ref. 40; quench from Tm to 300 K at 200 K/ps; anneal at 500 K for 4 ps to sample amorphous states. DFT settings ensure convergence to within 20 meV/atom (energy), 0.3 eV/Å (forces), and 10 kbar (stress). NNPs are trained with SIMPLE-NN using Behler–Parrinello symmetry functions (per pair: 8 radial, 18 angular; cutoffs 6.0 Å radial, 4.5 Å angular). Training proceeds until validation errors fall below ~20 meV/atom, 0.3 eV/Å, and 20 kbar for energy, forces, and stress. Iterative NNP refinement (refining-stage CSP): With the initial NNP, CSP runs for 50 generations; from structures within the bottom 100 meV/atom, select 10 (5 lowest-energy and 5 with lowest antiseed weights) and relax them with DFT (AMP) to forces <0.1 eV/Å and stress <20 kbar; sample positions/energies during relaxation to augment training. Repeat refinement at 100, 200, and 400 generations (excluding previously learned structures), totaling 40 crystalline structures added. A k-point mesh is adaptively updated during DFT relaxations for cell shape changes. Main CSP with NNPs: Up to 5000 generations using LAMMPS-based relaxations under pairwise distance restraints derived from the melt–quench trajectories (MQ distance constraints: minimum observed pair distances), preventing unphysical configurations and guiding toward trained domains. Every 1000 generations, a single-shot DFT energy of the current lowest-energy structure is compared with its NNP energy; if discrepancy >50 meV/atom, the NNP is retrained on 10 additional structures as in the refining stage. Final candidates within the bottom 50 meV/atom are fully DFT-relaxed with tight criteria (≤2 meV/atom energy, ≤3 kbar pressure). Structure generation and evolutionary operations: Four modes: random generation, permutation, lattice mutation, and crossover (used mainly in refining stage). Random generation uses RandSpg to choose a space group, lattice vectors, and Wyckoff positions; first-generation density matches amorphous MD; thereafter, volumes sample 70–130% of the prior generation’s best structure. MQ distance constraints are enforced in generation and relaxation. Crossover leverages atomic energies from the NNP to place cutting planes and adjust lattice/registries to preserve low-energy chemical units (e.g., polyhedra). Structural diversity control uses partial radial distribution function (PRDF) similarity to remove near-duplicates and an antiseed algorithm that biases inheritance probabilities against overrepresented motifs. In refining stage, operation ratios: 30% random, 50% crossover, 10% permutation, 10% lattice mutation; in the main CSP, crossover is excluded and ratios are 70% random, 20% permutation, 10% lattice mutation. The number of generated structures per generation equals twice the number of atoms. Inheritance selects from bottom 200 meV/atom (refining) or 100 meV/atom (main), mixing random and antiseed-weighted choices, and directly carries over the lowest-energy subset (100 meV/atom in refining; 50 meV/atom in main). DFT evaluation: VASP with PBE is the default for ranking; in analyses, SCAN and spin–orbit coupling (SOC) are used to probe cases with suspected PBE mis-ordering. AMP automates DFT relaxations. Additional tests include PBEsol where prior work used it. Performance and scaling: SPINNER is highly parallel; a full workflow (melt–quench MD, NNP training/refinement, 5000-generation CSP, DFT candidate relaxations) typically completes in 3–5 days on a 36-core node. Example resource split (Mg2SiO4): ~32% DFT-MD, 16% NNP training, 49% SPINNER (NNP relaxations), 3% DFT candidate relaxations.

- Blind test on 60 ternary compositions (Z ≥ 4): SPINNER identifies the experimental ground state (ΔEmin = 0) or a more stable structure (ΔEmin < 0) for 45/60 cases (75%), typically within ≤1000 generations for 38/45 successes; headline success described as ~80%. - Six cases showed ΔEmin down to −16 meV/atom, indicating SPINNER-found structures lower than experimental in PBE; recalculations with SCAN generally corrected the order (except PbOsO3 and LiYSn), and SOC increased ΔEmin for Tl3PbCl5 and PbOsO3 to −0.2 and +16 meV/atom. SCAN gave correct ordering across a 10-compound zero-ΔEmin test, confirming SCAN’s improved accuracy over PBE for energy ordering. - NNP fidelity: For materials with ΔEmin ≤ 0, mean absolute energy errors were ~12.9 meV/atom for the experimental structure and ~11.8 meV/atom averaged over candidates within 50 meV/atom, demonstrating good surrogate accuracy near low-energy basins. Larger errors (up to ~140–227 meV/atom) correlated with failures (e.g., SnGeS4, Sr2Pt3In). - Evolutionary path analysis: In nearly half the successful cases, the ground state was obtained by direct relaxation from a random structure (Ne = 0); when mutations occurred, Nm was typically ≤5 (permutation used in 34% and lattice mutation in 66% of such cases). MQ distance constraints filtered ~99% of random seeds and stabilized relaxations. Notably, NNP relaxations often funneled structures to the global minimum more reliably than DFT, likely due to smoother learned landscapes reducing trapping in metastable basins. - Metastable phase recovery: For Na3PS4, TlSbO3, and TlGaSe2 (high-quality ICSD metastables within 50 meV/atom in PBE), SPINNER recovered the reported metastable structures. - Benchmarking against literature predictions (21 non-oxide compositions from data-mining or DFT-based evolutionary searches): SPINNER found lower-energy structures in 13/21 and matched the rest, often at higher Z than previously considered; even when evaluated with PBEsol where applicable, SPINNER’s structures remained lower by 31 meV/atom (Sn5S4Cl2) and 4 meV/atom (Cd4SF6). For many cases (10/13), SPINNER’s lowest-energy structures had no matching prototypes in ICSD (per AFLOW-XtalFinder), underscoring prototype incompleteness. - Computational efficiency: End-to-end per-composition runtime is ~3–5 days on a 36-core node. For Mg2SiO4, SPINNER completed 5000 generations in ~58 hours (pool size 56, Z=4). A comparable DFT-only USPEX run reached only 31 generations under recommended settings; extrapolated, 5000 DFT generations would take years. - Algorithmic comparison with USPEX using identical NNPs (USPEX-NNP): SPINNER achieved a ~34 meV/atom lower minimum energy for Mg2SiO4, with similar gaps for Cd12O6, Sr2P7Br, and Mg12IrB. - Transfer learning: NNPs trained on one stoichiometry could guide CSP for related compositions (e.g., Mg2SiO4 NNP successfully found MgSiO3 ground state), with consistent energy ordering despite constant energy offsets due to single-stoichiometry training. - Generality: Beyond ternaries, SPINNER identified ground states for TiO2 (including all major polymorphs within 50 meV/atom at Z=4,8), P3N5, NbPd3, Li10GeP2S12, and InGaZnO4, and uncovered several new low-energy TiO2 polymorphs, suggesting applicability across complexities.

SPINNER directly addresses the CSP scalability challenge by replacing most costly DFT relaxations with NNP surrogates trained from self-generated disordered data and iteratively refined on low-energy crystalline structures. The combination of MQ distance constraints and evolutionary/random search tailored to NNPs (including atomic energy-guided crossover, antiseed diversity control, and PRDF-based pruning) enables efficient navigation of configuration space, often allowing direct relaxation from filtered random seeds to the ground state. The observed smoothing of certain metastable basins in the NNP potential energy surface effectively accelerates descent to global minima, overcoming traps that can hinder DFT-driven searches. The blind-test success rate (~75–80%) across chemically diverse ternaries, recovery of known metastables, and superior performance versus data-mining and DFT-based evolutionary predictions validate the framework. Systematic checks with SCAN and SOC show that apparent discrepancies with experiment in PBE ordering are largely resolved with more accurate functionals, indicating that SPINNER’s search effectiveness is robust and the remaining disagreements are tied to the electronic-structure level of theory rather than the search method. The ability to vary Z easily (linear scaling of NNP vs cubic scaling of DFT) expands the search space and proved important for finding lower-energy phases missed by prior studies constrained to small Z. Transfer learning further extends utility across compositions. Overall, the results demonstrate that an NNP-accelerated evolutionary/random CSP workflow can provide reliable ground-state identification and large throughput, enabling broader materials exploration.

The work introduces SPINNER, a fast and reliable CSP framework that couples melt–quench–trained and iteratively refined NNPs with evolutionary and random search algorithms adapted to machine-learning potentials. In extensive blind tests on 60 ternary compositions, SPINNER recovered experimental ground states or more stable structures in most cases and outperformed prototype-based data mining and DFT-driven evolutionary methods, all while operating 10^2–10^3 times faster than DFT-only CSP. Additional successes include identification of known metastables, discovery of lower-energy structures lacking ICSD prototypes, and generalization to non-ternary systems. These advances pave the way for large-scale, open-ended computational searches for undiscovered inorganic crystals and the construction of targeted materials databases at practical cost. Future directions include: fully automated end-to-end workflows; improved descriptors and larger cutoffs to better capture medium- and long-range order without sacrificing efficiency; extensions to handle magnetic materials and long-range interactions; enhanced diversity metrics beyond PRDF; strategies to optimize Z and cell shapes; and scaling to quaternary and higher-order systems with adaptive stopping criteria.

- Failure rate of ~25% in the ternary blind test: in some cases, NNPs exhibited large energy errors (e.g., SnGeS4, Sr2Pt3In), hindering correct ranking; others likely require more than 5000 generations or different cell choices. - Medium-to-long-range order is imperfectly captured due to finite descriptor cutoffs (6 Å radial, 4.5 Å angular), trading accuracy for speed; increasing cutoffs improves accuracy but reduces efficiency. - PRDF-based similarity metrics lose sensitivity for large supercells or compositions beyond ternaries; more detailed atom-centered descriptors are desirable. - Crossover method relying on atomic energies is incompatible with MLPs incorporating explicit long-range interactions; generalizing to such models remains to be addressed. - Current NNPs do not capture magnetic ordering and related energy scales; materials with 3d or f electrons and magnetic effects were excluded. - PBE can mis-order phases; while SCAN and SOC often rectify ordering, final stability assessments can be functional-dependent. - Search efficiency can depend on unit-cell choice; initial use of conventional cells improved outcomes for some cases (e.g., BaGe2S5, Na3SbO3, YPdGe). - Single-stoichiometry training induces constant energy offsets across compositions in transfer learning, requiring small retraining sets for absolute energies if needed. - Quaternary or higher-order systems may require longer searches and additional methodological tuning.