Matbench Discovery -- A framework to evaluate machine learning crystal stability predictions

1 Introduction

• Materials discovery spans roughly 10^10 possible quaternary compounds permitted by electronegativity and charge-balancing rules, with still less exploration of higher-order compositions.The passage contrasts this search space with approximately 10^5 experimentally tested and 10^7 simulated combinations.Source: §1 · p. 1
• The potential materials search space includes applications such as more efficient solar cells, lighter and longer-lived batteries, and smaller, more efficient transistor gates.The passage frames materials discovery as an optimization problem of mixing and arranging atoms against a merit function capturing material properties.Source: §1 · pp. 1–2
• ML can extract trends from large datasets and handle high dimensionality, multiple objectives, uncertainty, and noisy or sparse data.These capabilities complement traditional computational materials-science methods.Source: §1 · p. 2
• Density-functional theory offers a fidelity–cost compromise but consumes substantial computing resources, motivating efforts to reduce its computational burden.The passage reports up to 45% of core hours at Archer2 and over 70% of materials-science allocation time at NERSC devoted to DFT.Source: §1 · p. 2
• ML models are less accurate and reliable but orders of magnitude faster than ab-initio simulation, making them suitable for triaging high-throughput searches before more expensive DFT.Universal models trained on diverse datasets cover 90 or more application-relevant elements.Source: §1 · p. 2
• No large-scale benchmark yet measures ML’s ability to accelerate materials discovery, leaving the best-suited methodologies unclear.Prior prospective studies used one-shot predictors, iterative Bayesian optimizers, or universal interatomic potentials on specific systems and custom datasets, without systematic standardized comparison.Source: §1 · p. 2
• Universal interatomic potentials outperform the other tested methodologies in accuracy and robustness under a framework designed to simulate real-world ML-guided discovery campaigns.The framework aims to rank models by their ability to accelerate materials discovery.Source: §1 · p. 2
• The authors hope standardized benchmarks will let researchers with limited materials-science backgrounds contribute to model architecture and methodology development.They present this as a possible pathway to aid progress in materials science.Source: §1 · p. 2

2 Evaluation Framework for Materials Discovery

• The framework addresses prospective testing, thermodynamic-stability targets, decision-relevant metrics, and scalability to broad chemical spaces and test sets larger than training data.It argues that realistic covariate shift and classification performance better reflect discovery use than retrospective splits or global regression errors alone.Source: §2 · pp. 2–3

2.1 Matbench Discovery

• Matbench Discovery permits practitioner-available training data but tests prediction of relaxed-structure convex-hull distance from unrelaxed structures, avoiding a discovery workflow that depends on prior DFT relaxation.Convex-hull distance is used instead of formation energy because it better informs stability decisions; composition-only models may predict formation energies accurately while performing poorly on decomposition enthalpy.Source: §2.1 · pp. 3–4

2.1.1 Materials Project Training Set

• The canonical training set is the Materials Project v2022.10.28 snapshot of energies, forces, stresses, and magnetic moments across ionic steps, with flexibility to use subsets and auxiliary targets.The benchmark treats use of relaxation byproducts as a fair comparison because those data share the workflow’s overall acquisition cost, and it allows system-level data curation and optimizer choices to count toward performance.Source: §2.1.1 · p. 5

2.1.2 WBM Test Set

• After protostructure cleaning and deduplication, the WBM test set contains 215,488 structures, including 32,942 labeled stable at or below the MP convex hull.Its five substitution batches enable OOD analysis as substitutions move structures away from MP coverage; the benchmark’s primary metrics aggregate batches, and its labels use the MP hull despite fewer than half of the roughly 33,000 below-hull candidates remaining on the joint MP+WBM hull.Source: §2.1.2 · p. 6

2.1.3 Limitations of this Framework

• The framework does not capture synthesis-relevant disorder or dynamic stability, and its energy-focused metrics may not reflect UIP performance on other tasks.Using the MP hull as the positive-label reference avoids path-dependent cumulative precision metrics, but can overestimate performance; the overestimation decreases as the hull becomes better sampled.Source: §2.1.3 · pp. 6–7

2.2 Models

• The initial benchmark release evaluates 13 models using varied energy, force, stress, and magnetic-moment training targets, including gradient-based or direct force and stress prediction.Source: §2.2 · p. 7
• Universal interatomic potential models trained on energies, forces, and stresses occupy the leaderboard’s top positions, with a pronounced regression-metric gap over seven energy-only models.The table ranks models by F1 score.Source: Table 1 · p. 8
• Voronoi RF, CGCNN, and MEGNet perform worse than the dummy on regression metrics but better on some classification metrics; Voronoi RF combines Magpie and relaxation-robust Voronoi features.This illustrates why regression metrics alone can be misleading.Source: Table 1 · p. 8, §2.2 · p. 9
• CGCNN+P augments relaxed training structures with five perturbed structures mapped to the same DFT final energy, but its mitigation strategy did not improve F1 score or DAF.Unlike the other structure-based GNNs considered, it learns a step-function mapping each potential-energy valley to its local minimum.Source: Table 1 · p. 8, §2.2 · p. 8
• ALIGNN processes interatomic bond and line graphs for three-body angles, but its global pooling makes it ill-suited to force-field applications; ALIGNN-FF removes that pooling.ALIGNN was trained on the MP-crystals-2022.10.28 dataset for this benchmark.Source: §2.2 · p. 8
• Wrenformer replaces Wren’s message-passing layers with QKV self-attention, reducing memory use and scaling to structures with more than 16 Wyckoff positions.Its coordinate-anonymized Wyckoff inputs allow enumeration of spacegroup and Wyckoff-position combinations for a composition and maximum unit-cell size.Source: §2.2 · p. 9
• BOWSR combines symmetry-constrained Bayesian optimization with a surrogate energy model for iterative potential-energy-landscape search, and its high sample count makes it the most expensive tested model.The benchmark uses the pretrained Eform_MP_2019 MEGNet model as its energy model.Source: §2.2 · p. 9

3 Results

• DAF 6.33 for EquiformerV2 + DeNS on each model’s top 10k predictions approaches the task’s stated optimum, and it leads all nine reported metrics.Across the full dataset, its DAF is 5.04; the maximum possible DAF is approximately 6.5.Source: §3 · p. 9
• UIPs outperform energy-only models that predict from unrelaxed structures, learning force and stress labels to emulate DFT relaxation before predicting final energies.The benchmark reports that energy-only models can still perform well on F1 and DAF despite substantially worse regression metrics.Source: §3 · p. 9
• Regression accuracy does not guarantee better stability classification: CGCNN+P has better regression metrics than CGCNN but performs worse as a classifier.The authors attribute this to CGCNN+P’s sharper concentration of predictions near the 0 eV/atom threshold, where even small errors can flip classifications.Source: §3 · p. 9
• UIPs offer higher precision on campaigns validating about 20k or fewer materials, while CHGNet initially exceeds several models that perform better across the full test set.Energy-only models’ precision falls to 0.6 or less in the first 5k screened materials.Source: Figure 2 · p. 10, §3 · p. 10
• The “triangle of peril” marks regions where rolling MAE exceeds the distance to the stability threshold, making classifications vulnerable to error-driven flips.All models have lower rolling error on the plot’s left than its right; the authors associate this imbalance with greater susceptibility to false positives than false negatives.Source: Figure 3 · p. 11, §3 · p. 11

4 Discussion

• Several models achieve a discovery acceleration factor above 2.5 across the full dataset, reaching up to 6 on each model’s 10k most stable predictions; UIPs show superior accuracy and extrapolation performance.The authors connect UIP performance to modeling forces, which helps chart a path toward DFT-relaxed structures before final energy prediction.Source: §4 · p. 11
• The benchmark addresses only a limited subset of UIP applications, so the authors call for additional task-based evaluations focused on prospective settings, relevant targets, informative metrics, and scalability.As the convex hull becomes more comprehensively sampled, the fraction of unknown stable structures and the enrichment of test sets are expected to decline.Source: §4 · p. 11
• Larger training datasets have improved reported UIP performance, but model quality remains bounded by the quality and level of theory of training data.The authors recommend generating large quantities of training data at higher-than-PBE fidelity.Source: §4 · pp. 11–12
• Future discovery models should address properties across temperature and pressure and help predict synthesis pathways, including reaction barriers and non-crystalline, out-of-equilibrium materials.The authors identify temperature-dependent dynamical properties as an area suited to interatomic potentials.Source: §4 · p. 12

6 Author Contributions

• Riebesell led methodology, software, data curation, validation, analysis, and drafting, and trained CGCNN, CGCNN+P, Wrenformer, and Voronoi RF.Source: §6 · p. 12
• Goodall contributed conceptualization, software, validation, analysis, and writing; Benner contributed software and trained ALIGNN and MACE.Source: §6 · p. 12
• Chiang investigated MACE and contributed analysis and editing; Deng curated MPtrj, trained CHGNet, and contributed editing.Source: §6 · p. 12
• Ceder, Asta, Lee, Jain, and Persson provided supervision; Ceder and Asta acquired funding, and Persson also contributed editing and funding.Source: §6 · p. 12

• MACE achieves a DAF of 3.5 on the full test set versus 3.85 on the unique-protostructure subset.The passage attributes the roughly 10% difference largely to stable-structure prevalence: 16.7% in the full set versus 15.3% in the subset.Source: Metrics on full test set and for 10k… · p. 18
• Table 1 uses unique WBM protostructures to better reflect extrapolation to out-of-domain materials, while acknowledging possible deduplication bias.Protostructures are matched using prototype labels and chemical system; structures with symmetry-preserving lattice changes can count as duplicates.Source: Metrics on full test set and for 10k… · p. 18
• Table S2 evaluates each model on the 10k test materials it predicts to be most stable, with each model selecting a different slice.The passage presents this as a simulation of campaigns limited to 10k validations.Source: Metrics on full test set and for 10k… · p. 18

• A material is classified as stable when its predicted energy above hull is below the stability threshold.Because all models predict formation energy, they are insensitive to changes in threshold t.Source: ROC Curves · p. 18
• The ROC curves plot the fraction of unstable structures classified as stable against the fraction of stable structures classified as stable.These are the false positive rate and true positive rate, respectively.Source: ROC Curves · p. 18
• The diagonal no-skill line represents a dummy model that randomly ranks material stability.Source: ROC Curves · p. 18

• All models perform well far below the convex hull, but struggle most around x = 0, where performance varies more above the hull.Occasional underpredictions occur for materials far above the hull.Source: Parity Plots · p. 18
• Wrenformer shows many severe energy overpredictions at x = 0, dominated by poorly represented spacegroup 71 A2BC compounds trapped in local minima.These compounds are unsuitable training examples for its global-minimum energy prediction conditioned on protostructure, and their inclusion leads to errors on similar correctly relaxed test structures.Source: Parity Plots · p. 18
• MACE has severe underpredictions at x = 0 with no identified common composition or crystal-symmetry pattern.Source: Parity Plots · p. 18
• MACE also exhibits rare, reproducible post-relaxation energy failures; the largest reported derailed prediction is −1022 eV/atom for wbm-3-31970.The received M3GNet-trained checkpoint produced several hundred such cases, while the authors’ MPtrj-trained checkpoint produced 44 among 250k test crystals.Source: Parity Plots · p. 18
• Table S1 covers all 257k WBM test structures, with duplicate protostructures included, finding 16.7% stable structures (42,825 of 256,963) and a maximum possible DAF ceiling of 6.Table S2 reports stability metrics for the 10k materials each model predicts most stable; each model selects a different test-set slice.Source: Table S1 · p. 19, Table S2 · p. 19

Hull Distance Box plot

• BOWSR has the largest median hull-distance error, while Voronoi RF has the largest interquartile range; models systematically underpredict hull distance.The proposed explanation is that MP training data favors stable, low-formation-energy structures, unlike the distributionally distinct WBM test set.Source: Hull Distance Box plot · p. 20, Figure S4 · p. 23

Classification Histograms using Model-Predicted Energies

• CGCNN+P’s predictions more closely match the DFT hull-distance distribution than CGCNN’s, yet can classify worse by placing predictions nearer the 0 eV/atom stability threshold.Near that threshold, even small errors can flip a stability classification.Source: Classification Histograms using… · p. 20

Measuring extrapolation performance from WBM batch robustness

• UIPs show less degradation across later WBM substitution batches than Wrenformer and Voronoi RF, while MEGNet and CGCNN retain higher rolling MAE with little batch degradation.The contrast suggests composition or coarse structural features alone may be insufficient for an extrapolatable PES map; Wrenformer may improve with sufficiently diverse training data, but its lack of predictive forces and stresses may limit some applications.Source: Measuring extrapolation performance… · p. 23

Largest Errors vs. DFT Hull Distance

• Among the 200 largest-error structures, model predictions form two groupings when colored by disagreement, and high-error predictions tend toward hull distances that are too small.The figure suggests model committees could screen large databases for possible DFT errors.Source: Largest Errors vs. DFT Hull Distance · p. 23, Figure S7 · p. 26

Exploratory Data Analysis

• WBM’s convex-hull-distance distribution is much narrower than its formation-energy distribution, making stability prediction a more discriminative and challenging task than formation-energy prediction.The reported WBM formation-energy spread approaches 10 eV/atom, while hull-distance spread is less than 1 eV/atom; the authors describe stability prediction as ranking subtle energy differences among chemically similar compounds.Source: Exploratory Data Analysis · pp. 23, 26, Figure S11 · p. 30
• WBM differs chemically from MP: only 11% of its structures contain oxygen versus 53% in MP, and WBM is dominated by ternaries at 75%.WBM also more frequently contains noble-transition metals, post-transition metals, lanthanides, actinides, and metalloids relative to MP.Source: Figure S8 · p. 27, Figure S9 · p. 28, Figure S10 · p. 29