Surface Stability Modeling with Universal Machine Learning Interatomic Potentials: A Comprehensive Cleavage Energy Benchmarking Study

TL;DR

This study benchmarks universal machine learning interatomic potentials for predicting cleavage energies across various materials, revealing data composition's critical role over model complexity and highlighting efficient architectures trained on diverse datasets.

Contribution

It provides the first comprehensive assessment of uMLIPs for cleavage energy prediction, emphasizing the importance of training data diversity over architectural complexity.

Findings

Models trained on non-equilibrium data achieve <6% error.

High accuracy in predicting stable surface terminations.

Simpler models trained on diverse data match complex architectures.

Abstract

Machine learning interatomic potentials (MLIPs) have revolutionized computational materials science by bridging the gap between quantum mechanical accuracy and classical simulation efficiency, enabling unprecedented exploration of materials properties across the periodic table. Despite their remarkable success in predicting bulk properties, no systematic evaluation has assessed how well these universal MLIPs (uMLIPs) can predict cleavage energies, a critical property governing fracture, catalysis, surface stability, and interfacial phenomena. Here, we present a comprehensive benchmark of 19 state-of-the-art uMLIPs for cleavage energy prediction using our previously established density functional theory (DFT) database of 36,718 slab structures spanning elemental, binary, and ternary metallic compounds. We evaluate diverse architectural paradigms, analyzing their performance across chemical compositions, crystal systems, thickness, and surface orientations. Our results reveal that training data composition dominates architectural sophistication: models trained on the Open Materials 2024 (OMat24) dataset, which emphasizes non-equilibrium configurations, achieve mean absolute percentage errors below 6% and correctly identify the thermodynamically most stable surface terminations in 87% of cases, without any explicit surface energy training. In contrast, architecturally identical models trained on equilibrium-only datasets show five-fold higher errors, while models trained on surface-adsorbate data fail catastrophically with a 17-fold degradation. Remarkably, simpler architectures trained on appropriate data achieve comparable accuracy to complex transformers while offering 10-100x computational speedup. These findings show that the community should focus on strategic training data generation that captures the relevant physical phenomena.

Tables

Table 1. Table 1: The 19 state-of-the-art uMLIPs considered in this study, representing a diverse set of architectures. The number of model parameters and the training sets are summarized. Training datasets vary from equilibrium-focused collections (Materials Project) to non-equilibrium configurations (OMat24) and surface-adsorbate systems (OC20), enabling assessment of how training data composition affects out-of-distribution generalization. Parameter counts range from 4.55M to 1.4B, highlighting the diversity of model complexity.

Model Family	Model Name	#Params	Training Data Description
eSEN	eSEN-30M-MP	30.1M	MPtrj: 1.58M equilibrium structures
	eSEN-30M-OMAT24	30.0M	OMat24: 101M non-equilibrium structures
	eSEN-30M-OAM	30.2M	Combined: OMat24 + MPtrj + sAlex (10.4M)
EquiformerV2	Eqf.V2-153M-OMAT24	153M	OMat24: Non-equilibrium configurations
	Eqf.V2-153M-OAM	153M	Combined: OMat24 + MPtrj + sAlex
	Eqf.V2-153M-OC20	153M	OC20: Surface-adsorbate interactions
UMA	UMA-m-1p1-OMAT24	1.4B	OMat24: Full non-equilibrium dataset
GRACE	GRACE-2L-OMAT	12.6M	OMat24: Non-equilibrium structures
	GRACE-2L-OAM	12.6M	Combined: OMat24 + MPtrj + sAlex
	GRACE-2L-MP	12.6M	MPtrj: Equilibrium MD trajectories
MACE	MACE-MP-0b3	4.69M	Materials Project: Equilibrium only
	MACE-MPA	9.06M	MP + Alexandria (subset)
	MACE-OMAT	9.06M	OMat24: Non-equilibrium focus
SevenNet	SevenNet-MF-OMPA	25.7M	Combined: MPtrj + sAlex + OMat24
ORB	ORB-v3-OMAT	25.5M	OMat24: Non-equilibrium
	ORB-v3-MPA	25.5M	MP + Alexandria (subset)
	ORB-v2	25.2M	MPtrj + full Alexandria
MatterSim	MatterSim-v1	4.55M	Proprietary workflow-generated data
AlphaNet	AlphaNet-OMA	4.65M	Combined: OMat24 + MPtrj + sAlex

Table 2. Table 2: Comprehensive performance evaluation of 19 uMLIPs for cleavage energy prediction. Models are grouped by architectural family and sorted within each family by mean APE from lowest to highest error. Error metrics quantify prediction accuracy across 36,718 surface structures, with MAPE > > 100% indicating the count of structures with Mean Absolute Percentage Error exceeding 100% (outliers). Stability agreement indicates the percentage of 3,699 bulk materials where models correctly identify the thermodynamically most stable surface termination (Term.), Miller plane only (Miller), or neither (None). Computational time represents single-point energy evaluation on slab structures. Bold values indicate the best performance across all models for each metric (lowest values for error metrics, including MAPE > > 100%, highest for stability agreement).

	Error Metrics						Stability Agreement (%)
Model	APE Mean	APE Med.	MAE	RMSE	R2	MAPE	Term.	Miller	None	Params	Time
	(%)	(%)	(meV/Å2)	(meV/Å2)		$>$100%					(ms)
Mean Baseline (Reference)
Mean Baseline	48.7	28.6	33.8	46	0.000	2816	—	—	—	1	—
UMA
UMA-m-1p1-OMAT24	5.8	2.8	5.2	16	0.882	145	87.5	3.8	8.7	1.4B	215
EquiformerV2
Eqf.V2-OMAT24-MP-sAlex	5.2	2.5	4.8	9	0.961	24	84.0	4.8	11.2	153M	267
Eqf.V2-OMAT24	6.3	3.7	6.4	15	0.895	113	85.3	4.6	10.0	153M	153
Eqf.V2-OC20	71.1	62.7	68.1	87	-2.562	6641	16.2	8.8	75.1	153M	130
eSEN
eSEN-30M-OMAT24	5.7	2.9	5.0	15	0.898	140	86.6	3.6	9.8	30M	125
eSEN-30M-OAM	6.3	4.1	5.6	9	0.959	18	85.7	4.2	10.1	30.2M	121
eSEN-30M-MP	15.6	14.7	14.7	18	0.839	27	77.7	4.3	18.0	30.1M	121
GRACE
GRACE-2L-OAM	8.3	6.4	7.6	11	0.944	17	82.4	4.1	13.5	12.6M	8
GRACE-2L-OMAT	8.4	6.9	7.7	10	0.948	15	82.5	3.6	13.9	12.6M	8
GRACE-2L-MP	15.3	13.6	14.5	19	0.835	29	70.9	4.2	24.9	12.6M	9
ORB
ORB-v3-MPA	7.4	5.1	6.8	11	0.947	17	77.8	8.9	13.2	25.5M	11
ORB-v3-OMAT	7.8	6.0	7.3	11	0.946	23	78.2	6.2	15.5	25.5M	17
ORB-v2	9.2	6.3	8.8	14	0.906	39	69.3	9.2	21.5	25.2M	6
MACE
MACE-OMAT	7.8	6.1	6.9	10	0.955	15	81.2	3.9	14.9	9.06M	20
MACE-MPA	9.7	8.0	9.1	13	0.922	22	79.0	3.8	17.2	9.06M	20
MACE-MP-0b3	25.2	24.6	24.7	30	0.584	17	69.9	4.5	25.6	4.69M	21
SevenNet
SevenNet-MF-OMPA	7.1	5.0	6.6	10	0.953	18	83.6	3.8	12.5	25.7M	84
MatterSim
MatterSim-v1	19.3	17.5	17.4	21	0.782	126	76.2	4.1	19.7	4.55M	53
AlphaNet
AlphaNet-OMA	11.1	9.3	10.3	14	0.908	21	77.6	4.5	18.0	4.65M	21

Equations

$$E_{\text{cleavage}}=\frac{E_{\text{slab}}-n_{\text{bulk}}\cdot E_{\text{bulk}}}{2\cdot A_{\text{slab}}}$$

Full text

Surface Stability Modeling with Universal Machine Learning Interatomic Potentials: A Comprehensive Cleavage Energy Benchmarking Study

Ardavan Mehdizadeh1 and Peter Schindler1,∗

Machine learning interatomic potentials (MLIPs) have revolutionized computational materials science by bridging the gap between quantum mechanical accuracy and classical simulation efficiency, enabling unprecedented exploration of materials properties across the periodic table. Despite their remarkable success in predicting bulk properties, no systematic evaluation has assessed how well these universal MLIPs (uMLIPs) can predict cleavage energies, a critical property governing fracture, catalysis, surface stability, and interfacial phenomena. Here, we present a comprehensive benchmark of 19 state-of-the-art uMLIPs for cleavage energy prediction using our previously established density functional theory (DFT) database of 36,718 slab structures spanning elemental, binary, and ternary metallic compounds. We evaluate diverse architectural paradigms, analyzing their performance across chemical compositions, crystal systems, thickness, and surface orientations. Our results reveal that training data composition dominates architectural sophistication: models trained on the Open Materials 2024 (OMat24) dataset, which emphasizes non-equilibrium configurations, achieve mean absolute percentage errors below 6% and correctly identify the thermodynamically most stable surface terminations in 87% of cases, without any explicit surface energy training. In contrast, architecturally identical models trained on equilibrium-only datasets show five-fold higher errors, while models trained on surface-adsorbate data fail catastrophically with a 17-fold degradation. Remarkably, simpler architectures trained on appropriate data achieve comparable accuracy to complex transformers while offering 10–100× computational speedup. These findings fundamentally reframe MLIP development priorities: rather than pursuing increasingly complex architectures, the community should focus on strategic training data generation that captures the relevant physical phenomena.

††footnotetext: * 1 Northeastern University, Department of Mechanical and Industrial Engineering, Boston, MA 02115, USA

∗ Corresponding author: p.schindler@northeastern.edu*

1 Introduction

Cleavage and surface energy govern fundamental interfacial processes, including crack tip propagation 1, 2, the brittle-to-ductile behavior of crystalline solids 3, 4, elemental segregation at grain boundaries 5, molecular adsorption at surfaces 6, heterogeneous catalysis 7, 8, 9, electron emission devices 10, 11, 12, photocathodes 13, semiconductor devices 14, and battery material interfaces 15, 16. The cleavage energy is equal to the surface energy for symmetric slabs, which determines the equilibrium shape of nanoparticles and determines the stability of surfaces 17, 18. First principles calculations can accurately predict this critical property; however, its computational expense, combined with the vast structural and chemical space of materials surfaces, limits high-throughput screening efforts.

The grand challenge of computational materials science lies in accurately predicting the properties of materials to enable the discovery and design of new materials that can meet the demanding requirements of modern technological applications 19. Over the past decades, computational methods have revolutionized this landscape 20, leading to a data-driven discovery approach, the “fourth paradigm of science”, that uses massive datasets generated by experiments and simulations 21. Machine Learning Interatomic Potentials (MLIPs), also termed “foundational potentials”, have emerged as a transformative solution, bridging the long-standing gap between accuracy and computational efficiency by emulating density functional theory (DFT) calculations at a fraction of the computational cost 22, 23, 24, 25, 26, 27. This powerful tool has rapidly evolved and enabled remarkable breakthroughs, from million-atom molecular dynamics simulations of crack propagation to the rapid exploration of complex phase diagrams and the accelerated discovery of novel alloys with enhanced mechanical properties 28, 29, 30, 31.

Early MLIPs typically targeted one or a few chemistries, using hand‐engineered descriptors whose complexity grew combinatorially as new elements were added 32. Architecturally, the field has progressed from fixed descriptor-based functions 32, 33 to dynamic graph neural networks like SchNet 34 and DimeNet 35 that learn atomic interactions directly from raw coordinates through iterative message-passing 36. The advent of equivariant 37 and attention-based 38 neural networks has further improved data efficiency and long-range interaction modeling, while active-learning workflows 39 automatically expand datasets into underrepresented configuration spaces. The recent development of universal MLIPs (uMLIPs) has been enabled by massive computational datasets spanning the periodic table, including Materials Project40, Open Catalyst41, 42, Alexandria43, and Open Materials 202444 databases, which collectively provide millions of DFT calculations across diverse chemical environments. These datasets each contribute unique training regimes, fundamentally shaping uMLIP capabilities (see Supplementary Section 1 for comprehensive dataset descriptions).

Despite training on these comprehensive datasets, the rapid proliferation of uMLIPs has raised important questions about reliability, transferability, and generalizability across different applications. With their varied architectures, descriptors, and training databases, determining which uMLIP framework delivers optimal performance for specific use cases has become a critical challenge for the materials science community. While uMLIPs like MACE 45, ORB 46, 47, GRACE 48, and eSEN 49 offer the promising vision of a universal force field that could streamline simulations and accelerate materials discovery, their performance must be systematically evaluated across diverse atomic environments beyond their training data. Particularly crucial is assessing how models trained primarily on equilibrium bulk crystal structures perform when predicting properties of materials containing defects, surfaces, or under dynamic conditions such as phonon behavior and thermal expansion that involve near-equilibrium crystalline configurations.

Recent benchmarking efforts have systematically evaluated uMLIPs across various materials properties and tasks. The Matbench Discovery framework 50 established a comprehensive evaluation paradigm addressing four fundamental challenges: prospective benchmarking that mimics real discovery campaigns, thermodynamically relevant targets (convex hull distance rather than formation energies), informative classification metrics beyond regression accuracy, and scalability to large chemical spaces. Complementing this, the MLIP Arena platform 51 shifts focus from error-based metrics to physics-aware evaluation, revealing that gradient-based force predictions maintain conservative fields essential for stable simulations, while direct force prediction models exhibit center-of-mass drifts exceeding $10^{2}\;\text{\AA }$, six orders of magnitude larger than gradient-based approaches. Their hydrogen combustion benchmarks exposed critical stability issues, with CHGNet52, M3GNet, and EquiformerV2 (OC20) failing to complete 1 ns molecular dynamics trajectories, highlighting that strong regression performance does not guarantee physical reliability.

Property-specific benchmarks have revealed systematic patterns in uMLIP capabilities across diverse applications. For phonon anharmonicity, graph neural networks and artificial neural networks trained on ab initio molecular dynamics achieve agreement up to fifth-order irreducible derivatives 53. Reactive hydrogen dynamics at metal surfaces show that REANN and MACE provide an optimal balance between accuracy and computational efficiency, achieving sub-millisecond force evaluations while maintaining chemical accuracy 54. Harmonic phonon properties across 10,000 ab initio calculations demonstrate that models with non-conservative force predictions excel at equilibrium geometry optimization yet catastrophically fail at phonon properties, while gradient-based forces show better phonon accuracy 55. Comprehensive assessments across high-entropy alloys and disordered oxides revealed that equivariant uMLIPs offer 1.5–-2× improvements over non-equivariant architectures, though this advantage requires an order of magnitude more training data to achieve comparable accuracy for simpler systems 56. Analyses of over 2,300 MLIP models have identified fundamental trade-offs, with 35 challenging properties, particularly defect formation energies, elastic constants, and energy rankings, where fewer than 15% of models achieve errors below DFT-derived or practically-motivated thresholds (e.g., 0.1 eV for defect formation energies).

57. Intriguingly, these benchmarking efforts hint at a pattern we explore systematically: model performance on out-of-distribution tasks appears more strongly influenced by training data composition than architectural sophistication.

The computational cost of generating comprehensive surface data, which requires numerous orientations, terminations, and convergence tests per material, has historically prevented its systematic inclusion in existing training sets. Hence, the training data landscape was heavily skewed toward bulk configurations, with surface structures appearing only in specialized catalysis datasets (OC20/OC22) that focus on adsorbate-surface interactions rather than surface creation itself 41, 42. While Focassio et al. recently assessed earlier uMLIPs (MACE-MP-0, CHGNet-0.3.0, and M3GNet) for surface energy predictions, their evaluation was limited to elemental systems (73 elements, 1,497 surfaces) and revealed systematic underestimation requiring extensive fine-tuning, with root mean squared error improvements from 25% to 6% 58, 59. The performance of modern uMLIPs on chemically complex surface energies remains unexplored.

Here, we present a comprehensive benchmark of 19 state-of-the-art uMLIPs for cleavage energy prediction using our previously established high-throughput DFT database of 36,718 surface structures spanning elemental, binary, and ternary metallic compounds. Our benchmarking study, comprising over 36,000 surface evaluations for each of the 19 models (totaling over 1.3 million individual energy predictions), represents the largest systematic evaluation of uMLIP transferability to surface properties to date. Our evaluation encompasses diverse architectural paradigms, including equivariant transformers (EquiformerV2, UMA), graph neural networks (GRACE, eSEN), orbital-based representations (ORB), and multi-scale approaches (MatterSim), all tested on an out-of-distribution task for which they were not explicitly trained. We systematically analyze prediction accuracy across chemical compositions, crystal systems, thickness, and surface orientations to identify the structural and electronic factors that govern uMLIP transferability. Beyond accuracy metrics, we evaluate each model’s ability to correctly identify the thermodynamically most stable surface terminations, a critical requirement for materials design applications. This work establishes cleavage energy as a rigorous benchmark for assessing uMLIP generalization to surface phenomena and provides actionable insights for developing next-generation foundational potentials capable of accurately modeling interfacial properties across the periodic table.

This paper is organized as follows: Section 2 details our high-throughput DFT benchmark dataset of 36,718 surface structures and describes the 19 evaluated uMLIPs spanning diverse architectural paradigms, as illustrated in Figure 1. Section 3 presents comprehensive performance metrics, analyzing prediction accuracy across chemical compositions, crystal systems, slab thicknesses, and surface orientations, and discusses the implications of our findings for uMLIP development and deployment. Section 4 concludes with recommendations for developing next-generation universal potentials capable of accurately modeling surface and interfacial properties.

2 Methods

We present a comprehensive benchmarking framework to evaluate how well uMLIPs, trained primarily on bulk structures, can predict cleavage energies. Figure 1 illustrates our benchmarking workflow: we consider DFT-computed cleavage energies of surface structures from 3,716 bulk materials and use 19 uMLIPs representing different architectural families and training strategies for cleavage energy predictions, then compare the predictions to the DFT ground truth, finally, we conduct extensive analyses and quantitative assessment. The following subsections detail each aspect of this methodology, from the computational challenges that have historically limited surface property datasets (Section 2.1) through our statistical analysis framework (Section 2.6).

2.1 Computational Challenges in Cleavage Energy Calculations

The computational demands of comprehensive cleavage energy calculations have historically limited their application in high-throughput materials screening. A rigorous assessment requires evaluating all symmetrically distinct crystallographic planes, which scales drastically with decreasing crystal symmetry; even restricting to low Miller indices generates substantial computational demands 60, 61. Furthermore, surfaces with Miller indices up to 3 frequently contribute to Wulff constructions, directly influencing nanoparticle shapes and catalytic properties 62, making their evaluation essential for materials design applications.

Recent high-throughput studies have demonstrated that each orientation may exhibit multiple terminations with different stoichiometries and atomic arrangements, particularly for oxide materials, where polar compensation mechanisms can generate dozens of distinct surface configurations per facet 44. The computational workflow for each surface involves iterative optimization of slab thickness (typically $15$–$25\,\text{\AA }$ to achieve bulk-like interiors), vacuum spacing ($>15\,\text{\AA }$ to prevent inter-slab interactions), and surface reconstruction patterns that may involve complex (2×2), $(\sqrt{3}\times\sqrt{3})$, or even (7×7) superstructures 63. The Open Catalyst 2022 dataset exemplifies the computational challenge of surface-adsorbate systems specifically, requiring over $2\times 10^{8}\,\mathrm{CPU\,hours}$ to generate 62,331 calculations, though the majority (69%) were computationally intensive adsorbate+slab relaxations rather than bare surfaces 42. Machine-learning approximations can accelerate such calculations by factors of $10^{4}$–-$10^{6}$, yet training these models requires extensive DFT data 64. For comprehensive characterization, including multiple orientations, terminations, and adsorbate configurations, computational costs can be enormous even for a single material, creating a critical bottleneck in developing universal MLIPs for surface-property prediction.

2.2 DFT-Derived Cleavage Energy Dataset

Several databases of surface and cleavage energies have been developed to date, each with distinct scopes and limitations. Tran et al. 60 provided one of the earliest systematic collections, computing surface energies for over 100 polymorphs of approximately 70 elemental crystals with Miller indices up to 2 for non-cubic and 3 for cubic systems, totaling approximately 1,500 surfaces. While comprehensive for pure elements, it lacks the chemical complexity of multicomponent systems needed to evaluate uMLIP performance on diverse materials. Palizhati et al. 65 extended surface property calculations to intermetallic compounds, generating 3,033 cleavage energies across 36 elements and 47 space groups, with 72% bimetallic, 25% trimetallic, and 3% monometallic surfaces. While this introduced compositional complexity beyond pure elements, the dataset focused exclusively on intermetallic systems. Shi et al. 66 recently developed a comprehensive dataset of 3,526 surface energies and work functions for binary magnesium intermetallics, focusing on low-index surfaces (Miller indices up to 1) across 188 binary Mg compounds. This work specifically addressed the importance of surface properties in understanding corrosion behavior and mechanical performance of Mg alloys.

The Open Catalyst datasets represent significant contributions with OC20 41 containing $\sim 1.28$ million DFT relaxations (265 million single-point calculations) focused on adsorbate-metal interactions with Miller indices $\leq 2$, while OC22 42 specifically addresses oxide electrocatalysts with 62,331 systems (19,142 slabs and 43,189 adsorbate+slabs) and Miller indices up to 3. Though these datasets are extensive, they primarily focus on adsorbate-surface interactions rather than systematic clean surface energetics across diverse material classes.

For our benchmark, we selected the high-throughput DFT database previously developed by our group 67. This database contains 36,718 surface structures spanning elemental, binary, and ternary metallic compounds, stored as pymatgen structure dictionaries 68, providing substantial chemical diversity and scale for comprehensive uMLIP evaluation. The database focuses on Miller index 1 surfaces, including both symmetric and asymmetric surface terminations, and focuses on metallic systems (band gap $<0.1$ eV), providing a well-defined material class for benchmarking.

The cleavage energy, defined as the energy required to create new surfaces by separating a crystal along specific crystallographic planes, was calculated using:

[TABLE]

where $E_{\text{slab}}$ and $E_{\text{bulk}}$ represent the total energies of the slab and bulk configurations, respectively, $n_{\text{bulk}}$ is the number of bulk unit cells contained in the slab, and $A_{\text{slab}}$ is the surface area. All DFT calculations employed the Perdew-Burke-Ernzerhof (PBE) exchange-correlation functional 69 with plane-wave basis sets and a kinetic energy cutoff of 550 eV, as detailed in the original work 67.

2.3 Dataset Preparation and Filtering

To ensure computational tractability and maintain consistency in the DFT reference calculations, we applied a systematic filtering protocol to the original dataset. We excluded materials containing both oxide or fluoride anions (O, F) and specific transition metals (Co, Cr, Fe, Mn, Mo, Ni, V, W), as these combinations typically require Hubbard U corrections (DFT+U) to accurately describe their strongly correlated d-electrons and localized electronic states. Since the original DFT database employed standard PBE functionals without U corrections, these materials would have unreliable reference energies, potentially introducing systematic errors in our benchmarking analysis. This filtering removed only two bulk materials out of 3,718 total (0.05%): mp-19094 (V2O4) and mp-555934 (V2F4), both vanadium compounds that would require DFT+U corrections. These materials corresponded to 17 surface structures (6 surfaces for V2O4 and 11 for V2F4), resulting in the retention of 36,718 surface structures for benchmarking.

The resulting database provides cleavage energies ranging from 1.0 to 397.3 meV/Å2, with an average of 100.0 meV/Å2 for symmetric and asymmetric slabs. The filtered dataset was organized with comprehensive metadata including material identifiers (Materials Project IDs), Miller indices, surface terminations, slab thickness ratios ($n_{\text{slab}}/n_{\text{bulk}}$), surface areas, and both bulk and slab atomic structures.

While uMLIPs are computationally efficient, where completing each evaluation takes milliseconds rather than minutes/hours required by DFT, systematic organization was essential for managing the large amount of predictions across 19 models.

2.4 uMLIP Benchmarking Protocol

We evaluated 19 state-of-the-art universal MLIPs spanning diverse architectural paradigms and training strategies, none of which were exposed to our cleavage energy dataset during training. Table 1 summarizes the evaluated models, the number of their parameters, and the training datasets they utilized. These models represent the current landscape of foundation potentials in computational materials science, thus providing a comprehensive generalization test. The evaluated models were categorized into nine families based on their underlying architectures. The uMLIP training datasets each contribute unique aspects: equilibrium properties (MP), surface structures (OC20/OC22), dynamic trajectories (MPtrj70, 52), dimensional diversity (Alexandria), and non-equilibrium configurations (OMat24). See Supplementary Section 1 for comprehensive descriptions of all training datasets.

For each uMLIP, we performed single-point energy calculations on both the slab and bulk structures from our database, maintaining identical atomic configurations (positions and cell parameters) as used in the original DFT calculations. The total energies of both slab and bulk configurations predicted by uMLIPs are used to calculate the cleavage energy with Equation 1, comparing against the ground truth DFT value contained in the dataset. This approach enabled direct comparison between uMLIP-predicted and DFT-reference cleavage energies. Additionally, we recorded energy prediction times for each slab structure using uMLIPs to assess the computational efficiency of different model architectures.

To handle the computational demands of evaluating 19 models on $\sim$36,000 structures, we utilized a high-performance computing node equipped with four NVIDIA L40S GPUs, each featuring 45 GB of GDDR6 memory (135 GB total) and 18,176 CUDA cores running on the Ada Lovelace architecture. The system, powered by an Intel Xeon Gold 6448Y processor with 502 GB of system memory and running Rocky Linux 9.3 with CUDA 12.3, provided the necessary computational resources for our large-scale evaluation. We implemented a parallel processing framework with checkpoint-based recovery, distributing the 19 models across the available GPUs. The GPUs’ high memory bandwidth (864 GB/s per GPU) and compute capability enabled efficient processing, while intermediate results were saved at regular intervals (128 structures initially, then every 4,000 structures) to enable fault-tolerant execution.

2.5 Computational Implementation

All uMLIP evaluations were conducted using the official implementations of each model architecture with rigorously controlled software environments. The FAIRChem family of models employed two distinct versions: fairchem-core v1.10.0 for EquiformerV2 and eSEN architectures, and fairchem-core v2.3.0 for the UMA model, reflecting the major architectural changes between FAIRChem generations 41, 42. The ORB models utilized orb-models v0.5.4 47, 46, while MACE implementations employed mace-torch v0.3.13 45. SevenNet calculations were performed using sevenn v0.11.2.post1 71, and GRACE models were evaluated through grace-tensorpotential v0.5.1 with TensorFlow 2.16.2 backend 48. MatterSim predictions utilized the mattersim package v1.2.0 72, while AlphaNet employed a custom implementation (v0.0.1) built on PyTorch 2.1.2 with PyTorch Geometric dependencies 73. Each uMLIP was installed in an isolated conda environment to prevent dependency conflicts and ensure computational reproducibility. All models were executed on NVIDIA L40S GPUs with CUDA 12.3, utilizing model-specific optimizations where available. The software versions and environment specifications used for each uMLIP are detailed above to ensure computational reproducibility. The filtered dataset and all uMLIP predictions are publicly available (see Data Availability Statement section).

We developed a unified Python framework to orchestrate the computational workflow, abstracting model-specific initialization and prediction interfaces across all 19 uMLIPs. The framework managed GPU memory allocation and cleanup between predictions, tracked detailed timing information for performance analysis, and saved results in hierarchical directory structures organized by model family.

2.6 Statistical Analysis

We evaluated model performance using multiple error metrics to capture different aspects of prediction accuracy. Mean Absolute Error (MAE) provided the average magnitude of errors in meV/Ų, while the Absolute Percentage Error (APE) and Mean APE (MAPE) provide normalized errors relative to the true cleavage energy values, enabling comparison across different energy scales. Root Mean Square Error (RMSE) emphasized larger deviations by squaring errors before averaging, helping identify models with significant outliers. We also calculated the coefficient of determination ($R^{2}$) to assess the linear correlation between uMLIP predictions and DFT reference values. For computational efficiency analysis, we measured the wall-clock time required for each model to complete single-point energy calculations on slab structures, providing a practical comparison of model speeds across different architectures. The error distributions were examined to identify systematic biases or limitations in specific chemical environments, with particular attention to outliers that might indicate fundamental limitations in model transferability.

3 Results and Discussion

First, we analyze the overall performance of uMLIPs for cleavage energy prediction across our benchmark dataset of 36,718 surface structures. We then examine how training data composition influences model accuracy, followed by a detailed analysis of the best-performing models (see Supplementary Information for all other models). Finally, we investigate chemical and structural dependencies in prediction errors and discuss the computational efficiency trade-offs between different architectures.

3.1 Overall Performance of uMLIPs

We evaluated 19 state-of-the-art uMLIPs on cleavage energy prediction across 36,718 surface structures (see Supplementary Figure 1 for a complete comparison of all 19 models). Figure 2 presents a comprehensive overview of select model performance, revealing substantial variation in prediction accuracy across different uMLIP architectures.

The ridge plots (Figure 2a) illustrate the distribution of APEs for each model, with numerical annotations indicating the mode (most frequent error value) of each distribution. The $x$-axis truncation at 50% emphasizes the primary error distribution while excluding extreme outliers. Lower peak values indicate that models typically produce smaller errors, with UMA-m-1p1-OMAT24 achieving the lowest mode at 1.3% compared to MatterSim-v1 at 14.5%. The median APE values range from 2.85% (UMA-m-1p1-OMAT24) to 17.50% (MatterSim-v1), representing a six-fold variation in performance.

The dual-axis representation (Figure 2b) reveals that transformer-based (UMA) and equivariant (eSEN, EquiformerV2) architectures trained on OMat24 achieve median MAEs below 3.5 meV/Å2, while the multi-scale MatterSim model exhibits approximately five-fold higher median errors. This stark difference, despite MatterSim’s sophisticated multi-scale architecture, underscores the dominant influence of training data composition on model performance for out-of-distribution tasks (see Supplementary Tables 2–5 for complete performance metrics, outlier statistics, and summary statistics across all models).

Beyond energy prediction accuracy, we evaluated each model’s ability to identify the thermodynamically most stable surface terminations (Figure 2c), a critical requirement for practical materials design applications. For each of the 3,699 bulk materials in our dataset, we compared whether the uMLIPs correctly identified the lowest-energy surface among all possible crystallographic orientations (up to a Miller index of 1) and atomic terminations, as predicted by DFT. We classified the following three outcomes: (1) exact termination match, where the model correctly identified both the Miller index (crystallographic plane) and the specific atomic termination (i.e., which atoms form the outermost layer) that yields the minimum cleavage energy; (2) Miller match only, where the model identified the correct crystallographic plane but selected the wrong atomic termination among the possible terminations for that plane; and (3) no match, where the model predicted an entirely different crystallographic orientation as most stable one. The analysis across 3,699 materials reveals a clear performance hierarchy: UMA (87.5%) $>$ eSEN (86.6%) $>$ EquiformerV2 (85.3%) for exact termination matches. Remarkably, all OMat24-trained models achieve greater than 76% exact termination match accuracy without explicit surface energy training, validating that exposure to non-equilibrium configurations during training enables models to learn bond-breaking energetics relevant to surface stability.

To provide context for model performance, we establish a trivial baseline using the mean cleavage energy across the entire database (100.0 meV/Ų) as a constant prediction for all structures. This mean baseline has a MAE of 33.8 meV/Ų and MAPE of 48.7%, with 2,816 structures (7.7%) showing errors exceeding 100%. The top-performing models like UMA-m-1p1-OMAT24 achieve up to nine-fold reduction in error compared to this mean baseline (see Supplementary Table 10 for complete baseline performance metrics and improvement factors).

Quantitative analysis reveals an observable correlation between prediction accuracy and termination identification capability (as shown in Supplementary Figure 5). Models with lower MAE consistently achieve higher termination match rates ($R^{2}=0.542$), with the relationship strengthening for low-energy surfaces (0–50 meV/Ų) where small energy differences make correct ranking most challenging. This correlation confirms that accurate energy predictions translate directly to reliable identification of stable surface configurations.

Note that surface properties obtained from DFT, such as the cleavage energy, can depend on slab thickness and may further exhibit oscillations depending on whether the slab contains an even or odd number of atomic layers. A well-established strategy to mitigate this issue is to perform the bulk total energy calculation using a conventional unit cell that is reoriented along the relevant Miller index, thereby ensuring consistency in the $k$-point sampling between bulk and slab calculations.63 This protocol was also adopted in constructing the DFT database employed for benchmarking in this work 67. By contrast, it is less clear whether uMLIP predictions capture such effects, since uMLIPs do not involve a $k$-point mesh, and the corresponding training datasets typically rely on primitive bulk unit cells rather than conventional, reoriented cells aligned with the surface orientation.

3.2 Outlier Analysis

While median errors provide the primary performance metric, examining outliers reveals important failure modes across architectures (as shown in Supplementary Table 5). We define outliers as predictions with MAPE exceeding 100%, representing cases where models predict cleavage energies at least twice the DFT reference value. Most models exhibit remarkably low outlier rates ($<0.5\%$ of structures), with EquiformerV2-OC20 as a notable exception (18.1% outlier rate), confirming its unsuitability for general surface predictions. Among well-performing models, an interesting pattern emerges: models trained on OMat24 alone (eSEN-30M-OMAT24, EquiformerV2-153M-OMAT24, UMA-m-1p1-OMAT24) show slightly higher outlier counts (140, 113, and 145 structures, respectively) compared to their counterparts trained on combined datasets, yet maintain lower mean outlier APE (181–252% vs 695–878%). This suggests that pure OMat24 training produces more frequent but less severe failures, while combined dataset training yields fewer but more catastrophic errors when models fail. The median error of outliers provides additional insight—negative values around $-125$ meV/Ų for OMat24-trained models indicate systematic underestimation for these challenging cases, likely corresponding to high-energy surfaces or unusual coordination environments underrepresented in training data. In contrast, models trained on combined datasets show positive median outlier errors ($\sim 30$ meV/Ų), suggesting overestimation when encountering unfamiliar configurations. These outlier characteristics, while affecting less than 0.5% of predictions for most models, could be critical for applications requiring reliable uncertainty quantification or when screening for extreme surface properties.

3.3 Impact of Training Data versus Architecture

Interestingly, we observe that the choice of training dataset exerts a more profound influence on cleavage energy prediction accuracy than architectural design. This observation highlights the critical importance of training data curation for out-of-distribution tasks.

Consider the eSEN architecture family: eSEN-30M-OMAT24 (trained on OMat24) achieves a median APE of 2.9%, while the architecturally identical eSEN-30M-MP (trained on MPtrj) shows 14.7% median APE, a five-fold degradation despite using the same model architecture. Similarly, within the EquiformerV2 family, the model trained on OMat24 (3.7% median APE) dramatically outperforms its OC20-trained counterpart (62.7% median APE), representing a 17-fold improvement solely through training data selection.

This pattern persists across all architectural families examined. The GRACE models show median APE ranging from 6.9% (OMat24-trained) to 13.6% (MPtrj-trained), while maintaining identical 2-layer GNN architectures. Even the sophisticated MACE architecture cannot overcome poor training data: MACE-OMAT achieves 6.1% median APE compared to MACE-MP’s 24.6%, despite MACE’s advanced equivariant message-passing design.

The superiority of OMat24-trained models likely stems from their emphasis on non-equilibrium configurations. By sampling structures with Gaussian perturbations ($\sigma=0.5\,\text{\AA }$) and at elevated temperatures (up to 3000K), OMat24 exposes models to atomic environments with forces and stresses spanning orders of magnitude beyond equilibrium values, precisely the regime relevant for bond-breaking during cleavage. In contrast, models trained exclusively on equilibrium structures (MP) or surface-adsorbate interactions (OC20) lack exposure to the intermediate states critical for cleavage energy prediction.

This finding has profound implications for uMLIP development strategies. While architectural innovations, from simple graph networks to sophisticated equivariant transformers, contribute to performance improvements, our results demonstrate that a well-designed training dataset can provide an order-of-magnitude greater impact on out-of-distribution generalization. Hence, computational resources may yield the greatest returns if directed toward systematic dataset design that captures relevant physical phenomena.

The $R^{2}$ values across models trained on OMat24 consistently exceed 0.88, as shown in Table 2 (see also Supplementary Figures 2b and 4 for a complete $R^{2}$ comparison and pairwise metric correlations, respectively), confirming strong linear correlation with DFT predictions. In stark contrast, the negative $R^{2}$ for EquiformerV2-OC20 ($-2.562$) indicates predictions worse than using the mean value, mathematically confirming complete failure for out-of-distribution generalization from catalyst-specific to general surface data. We hypothesize that a fundamental issue is that OC20-trained models cannot accurately predict bulk energies, as the OC20 dataset contains only surface-adsorbate configurations without any bulk structures, rendering cleavage energy calculations, which require both bulk and surface energy predictions, inaccurate. This statistical validation underscores that training data relevance outweighs architectural sophistication for successful transfer to new property domains.

3.4 Impact of Architecture on Computational Efficiency

However, the choice of model architecture critically affects the computational cost during inference, with a moderate speed-accuracy trade-off. As shown in Table 2 (and Supplementary Figure 3), prediction times vary up to 45-fold across architectures: ORB-v2 and GRACE models require only 6–8 ms per structure, while transformer-based models like UMA-m-1p1-OMAT24 (215 ms) and EquiformerV2-MP-sAlex (267 ms) are 36× and 45× slower, respectively. This speedup translates to computational durations of days instead of months for high-throughput screening of the order of tens of millions of surfaces. Remarkably, simpler architectures like GRACE-2L-OAM achieve reasonable accuracy (8.3% MAPE) at 8 ms per structure, providing a 33-fold speed advantage over the most accurate EquiformerV2-MP-sAlex while maintaining sub-10% errors. SevenNet-MF-OMPA offers a middle ground with competitive accuracy (7.1% MAPE) at 84 ms per structure. The fastest models (ORB at 6–17 ms, GRACE at 8–9 ms, MACE at 20–21 ms) still maintain sub-10% MAPE when trained on OMat24, suggesting that for large-scale screening, these simpler architectures provide an optimal balance. This 10–45× speed advantage, combined with comparable accuracy when trained on appropriate data, questions whether architectural complexity is justified for routine materials screening applications.

3.5 Detailed Analysis of Optimal Model Performance

Among all evaluated models, UMA-m-1p1-OMAT24 demonstrated superior performance with a mean absolute error (MAE) of 5.2 meV/Ų and MAPE of 5.8%, and most importantly, achieved the highest termination match rate of 87.5%, correctly identifying the thermodynamically most stable surface configuration for nearly 9 out of 10 materials. Figure 3 provides a detailed analysis of this top-performing model across 36,718 surface structures.

The parity plot (Figure 3a) reveals excellent agreement between uMLIP-predicted and DFT-calculated cleavage energies, with $R^{2}=0.882$. The hexagonal binning visualization highlights the high density of accurate predictions along the ideal prediction line, with most structures concentrated in the 0–150 meV/Ų range. Notably, the model maintains good correlation even for high-energy surfaces above 200 meV/Ų, though with increased scatter, suggesting robust extrapolation to extreme cases.

The error distribution analysis (Figure 3b) reveals a slight systematic underestimation bias, with the distribution centered at a mean error of $-3.8$ meV/Ų and a median of $-2.1$ meV/Ų. The leptokurtic profile with high kurtosis indicates a sharp central peak with heavier tails than a Gaussian distribution, suggesting the model makes highly accurate predictions for most structures while producing larger errors for a small subset of challenging cases. Despite these outliers, the narrow central peak with most errors falling within $\pm 25$ meV/Ų demonstrates consistent model behavior across the majority of surface chemistries and crystallographic orientations.

As shown in Figure 3c (see Supplementary Table 1 for detailed statistics across all models and energy ranges), the elevated MAPE ($\sim$15%) for low-energy surfaces (0–50 meV/Å2) is primarily a mathematical artifact of relative error calculation—small absolute errors become large percentage errors when divided by small cleavage energy values. This affects only 2,813 structures (7.7% of the dataset). Conversely, high-energy surfaces (350–400 meV/Å2) maintain MAPE below 15%, demonstrating robust performance across the full energy spectrum, though this bin contains only 48 structures (0.1% of the dataset), limiting statistical significance.

While we present detailed analysis for the top-performing UMA-m-1p1-OMAT24 model here, similar comprehensive performance analyses, including parity plots, error distributions, and energy-dependent MAPE for all other 18 evaluated models, are provided in Supplementary Figures 6–23.

3.6 Chemical and Structural Dependencies

To understand the sources of prediction errors, we analyzed model performance of UMA-m-1p1-OMAT24 across different chemical and structural categories (Figure 4; similar decomposition analyses for all other models are provided in Supplementary Figures 24–41).

Analysis of element-specific performance (Figure 4a) reveals striking patterns in the model’s predictive accuracy. Alkaline earth metals demonstrate exceptional accuracy ($\mathrm{Be}$: 1.1%, $\mathrm{Mg}$: 3.8%, $\mathrm{Ca}$: 3.6%, $\mathrm{Sr}$: 3.4%), along with select transition metals ($\mathrm{Hf}$: 1.4%, $\mathrm{Re}$: 1.2%, $\mathrm{Ta}$: 1.9%, $\mathrm{Ni}$: 3.2%). In stark contrast, halogens exhibit the poorest performance ($\mathrm{Cl}$: 39.3%, $\mathrm{Br}$: 80.2%, $\mathrm{I}$: 74.6%), while alkali metals show variable accuracy ($\mathrm{Na}$: 6.2%, $\mathrm{K}$: 9.7%, $\mathrm{Cs}$: 34.0%). Lanthanides display intermediate but elevated errors, particularly $\mathrm{Eu}$ at 26.7%.

These patterns provide insights into the limitations of current training datasets. The excellent performance for transition metals likely reflects their prevalence in the OMat24 dataset, while the systematically higher errors for halogens and heavy alkali metals suggest these chemistries are underrepresented and/or exhibit complex bonding characteristics difficult to capture by the uMLIPs. The anomalously high error for Eu among lanthanides may stem from its variable oxidation states ($+2/+3$), indicating that elements with complex f-orbital occupancy remain challenging for uMLIPs (see Supplementary Table 8 for mean errors across the ten most common elements in the dataset).

Crystal system analysis (Figure 4b, detailed in Supplementary Table 7) demonstrates that prediction accuracy correlates with structural symmetry. Hexagonal and orthorhombic structures are predicted most accurately ($\mathrm{MAPE}\approx 5\%$), while triclinic and trigonal systems present greater challenges ($\mathrm{MAPE}$ between $10$–-$15\%$). This pattern persists across all models, suggesting that low-symmetry structures represent a fundamental challenge rather than an architecture-specific limitation.

Slab thickness analysis reveals robust model performance across typical computational thickness ranges (Figure 4c). UMA-m-1p1-OMAT24 maintains consistent accuracy for slabs up to 25 Å (MAPE: $\sim 6\%$), demonstrating reliable predictions across typical DFT slab thicknesses. Performance degrades for very thick slabs ($\geq 25$ Å), showing elevated errors (10.3% for UMA), likely reflecting receptive field limitations of the graph neural network architectures, where messages from the top side of the slab no longer reach the bottom side. Notably, all OMat24-trained models show similar trends (Supplementary Table 9), with MAPE increasing by approximately two-fold for slabs $\geq 25$ Å, suggesting this represents a systematic prediction challenge.

Additional structural analysis reveals that compositional complexity and surface symmetry have a limited impact on prediction accuracy for UMA-m-1p1-OMAT24. Unary compounds show slightly higher MAPE (9.8%) compared to binary (6.6%) and ternary (5.4%) systems (Supplementary Table 6), potentially reflecting the smaller sample size of elemental structures (236 vs. 10,185 and 26,297, respectively), making errors more sensitive to outliers. Surface symmetry shows a negligible effect, with symmetric and asymmetric slabs exhibiting comparable error rates, indicating models handle both slab types equally well. This pattern persists across all evaluated models except EquiformerV2-OC20, which exhibits catastrophic failure for unary compounds (273% MAPE) while maintaining its already poor performance for binary and ternary systems (see Supplementary Figures 43–60 for composition and symmetry analyses of all 19 models).

4 Conclusions

Our comprehensive benchmark of 19 uMLIPs on 36,718 surface structures reveals fundamental insights that reframe the challenge of uMLIP development. Rather than pursuing increasingly complex architectures, our results demonstrate that strategic training data generation, capturing the full spectrum of atomic environments, is paramount for accurate property prediction.

We show that training data diversity trumps architectural sophistication. The five-fold to seventeen-fold improvement from dataset selection (OMat24 vs. MPtrj) overshadows differences between architectural families. Moreover, exposure to non-equilibrium configurations proves essential for learning bond-breaking physics, with OMat24-trained models achieving sub-6% MAPE without explicit surface training, suggesting that sufficiently diverse bulk training data incorporating non-equilibrium states can enable accurate surface property prediction. Hence, models can identify stable configurations without expensive surface-specific training data. Finally, simpler architectures (ORB, MACE) offer 10–100× speed advantages over transformer-based models while maintaining comparable accuracy when trained on appropriate data, questioning whether architectural complexity is justified for routine materials screening.

While our study provides the most comprehensive surface property benchmark to date, the following limitations are worth noting. Our evaluation uses fixed DFT geometries, thus not testing the surface relaxation capabilities of uMLIPs. Further, the dataset focuses on materials with a near-zero bandgap and surface with a Miller index up to 1; extending this to higher indices and semiconducting or insulating materials would provide additional generalization insights. Additionally, all DFT references use the PBE functional, and exploring functional-dependent performance remains for future work (e.g., $r^{2}$SCAN). We have recently also developed a new equivariant model with symmetry breaking and uMLIP-informed features to directly train on the same database (predicting both the cleavage energy and work function). The best model in that work achieves a test-MAPE of 3.2%, which, surprisingly, is not even two times better than UMA-m-1p1-OMAT24, which was not directly trained on our dataset 74.

The impact of this work extends beyond surface energy prediction. By establishing cleavage energy as a rigorous out-of-distribution benchmark, we provide a framework for evaluating uMLIP transferability to other challenging properties. Our finding that training data quality dominates over architectural complexity suggests reallocating computational resources from architecture development to strategic dataset generation. Recent advances in simulation infrastructure, such as the TorchSim package 75 enabling batched GPU-accelerated simulations, could further accelerate uMLIP evaluation, potentially reducing benchmarking time by orders of magnitude. Additional future approaches may rely on a fine-tuning approach using MatterTune 76, or the recently established fully automated framework to capture the relevant potential energy surface 77.

Future directions should prioritize expanding training datasets to include underrepresented chemistries (particularly halogens and f-block elements), low-symmetry structures, and slab structures. In addition, deliberate sampling of transition states, metastable configurations, and other high-energy structures relevant to the intended applications will be essential for improving model transferability and robustness. The success of OMat24-trained models across diverse architectures demonstrates that current MLIP architectures are sufficient; the bottleneck lies in training data that captures relevant physics. As the field moves toward foundational potentials for materials design, our work establishes that a systematic evaluation on out-of-distribution tasks is essential for developing truly generalizable uMLIPs capable of predicting the full spectrum of materials properties.

Data Availability Statement

The filtered dataset and predictions for all slab cleavage energies for every uMLIP model will be available, open-source, on Zenodo at doi.org/10.5281/zenodo.16970767 after peer-review.

Conflict of Interest

The authors declare that they have no conflicts of interest relevant to the content of this manuscript.

Acknowledgments

P.S. acknowledges the start-up funding from the Department of Mechanical and Industrial Engineering at Northeastern University and funding through the TIER1 Seed Grant. This work was completed in part using the Explorer cluster, supported by Northeastern University’s Research Computing team.

Author Contributions

A. M. contributed to Writing – Original Draft (Lead), Writing – Review & Editing (Lead), Conceptualization (Equal), Software (Lead), Investigation (Lead), Data Curation (Lead), and Visualization (Lead). P. S. contributed to Writing – Original Draft (Supporting), Writing – Review & Editing (Supporting), Conceptualization (Equal), Visualization (Supporting), Resources (Lead), Supervision (Lead), and Funding acquisition (Lead). Both authors have read and approved this manuscript.