Machine-Learned Interatomic Potentials

• Machine-learned interatomic potentials (ML-IAPs) are data-driven models that decompose total energy into local contributions using invariant descriptors.
• They employ regression methods such as neural networks, Gaussian processes, and kernel methods to achieve near-DFT accuracy at much lower computational cost.
• Research focuses on optimal database design, active learning, and uncertainty quantification to enhance transferability and mitigate extrapolation risks.

Machine-learned interatomic potentials (ML-IAPs), often denoted MLIPs in the recent literature, are data-driven models of the potential energy surface used to provide computationally efficient predictions of energy and Newtonian forces for atomistic simulation. Their standard construction decomposes total energy into local atomic contributions, $E=\sum_i E_i$, maps local atomic environments $\mathbf{R}_i$ to descriptors $\mathbf{G}_i$ that are invariant to translations, rotations, and permutations, and fits a regression model to a quantum-mechanical reference database, most often generated with density functional theory (DFT). Within the interpolation regime defined by that database, ML-IAPs can achieve near quantum accuracy at a computational cost orders of magnitude below DFT, while active-learning and physically informed variants are designed to reduce failures associated with extrapolation to previously unseen environments (Mishin, 2021, Podryabinkin et al., 2016).

1. Formal structure and conceptual interpretation

A canonical ML-IAP expresses the total energy as a sum of atomic energies and evaluates each atomic energy from a descriptor of the local environment. In the review literature, this is written schematically as

$$\mathbf{R}_i \rightarrow \mathbf{G}_i \overset{\mathcal{R}}{\rightarrow} E_i,$$

with $\mathcal{R}$ implemented by neural networks, Gaussian process regression, kernel ridge regression, moment tensor potentials, SNAP, or related regressors. The same literature emphasizes that the predictive ability and transferability of an ML-IAP are strictly limited to the atomic environments represented in the reference database, which is why database design, active learning, and uncertainty monitoring are central rather than auxiliary components of the field (Mishin, 2021).

One influential conceptual interpretation is that high-accuracy ML-IAPs are not merely black-box surrogates. For elemental Ti, linearized MLIPs were presented as a generalization of the embedded atom method (EAM) and modified embedded atom method (MEAM), with the embedding energy treated beyond the uniform density approximation. In that treatment, pairwise descriptors take the form

$b_{n0}^{(i)}=\sum_j f_n(r_{ij}),$

and angular-dependent descriptors are written as

$$b_{nl}^{(i)}=\sum_{j,k} f_n(r_{ij})f_n(r_{ik})\cos^l(\gamma_{jik}), \qquad l\ge 1.$$

The central claim is that the high accuracy and transferability of MLIPs stem mainly from employing a large, systematic set of physically meaningful descriptors, especially angular-dependent ones, rather than from using sophisticated black-box functions alone (Takahashi et al., 2017).

A parallel line of work introduces hybrid or physically informed ML-IAPs in which the regression model predicts parameters of a physics-based potential rather than energies directly: $$\mathbf{R}_i \rightarrow \mathbf{G}_i \overset{\mathcal{R}}{\rightarrow} \mathbf{p}_i \overset{\Phi}{\rightarrow} E_i.$$ In the review formulation, this class is intended to improve transferability to unknown atomic environments while retaining the accuracy of data-driven fitting inside the database domain (Mishin, 2021).

2. Descriptor and model families

The field contains several distinct but partially overlapping model classes. A large fraction of early methodological work concerns linearly parameterized or polynomial models. In the active-learning formulation for linearly parameterized potentials, the energy of a configuration $x$ is written as

$E(x)=\sum_{j=1}^{m}\theta_j b_j(x),$

which covers models such as Moment Tensor Potentials (MTPs) and, in the cited discussion, other linearized families including SNAP and GAP expansions (Podryabinkin et al., 2016). For Ti, linearized MLIPs used polynomial functions of pairwise and angular descriptors; for Si–C, the UF3 formulation used explicit two-body and three-body terms expanded in cubic B-splines,

$\mathbf{R}_i$0

with a custom knot-placement strategy that increased spline resolution near the observed minimum in the radial distribution function around the C–Si bond distance (MacIsaac et al., 2024).

Kernel models remain a major branch. Gaussian Approximation Potentials (GAPs) use descriptors such as SOAP, while sparse Gaussian process regression (SGPR) models predict local atomic energies from kernel similarities. For organic CNH chemistry, the total energy was written as

$\mathbf{R}_i$1

with modified SOAP descriptors and inducing points, and multiple SGPR experts were combined by a Bayesian Committee Machine (BCM) to address chemical diversity and kernel-scaling limits (Park et al., 2024). SOAP-based GAPs also appear in studies of carbon phase behavior, iron, and minimalist training-set transferability (Marchant et al., 2022, Byggmästar et al., 2022, Robredo-Magro et al., 21 Nov 2025).

Bispectrum-based models define another widely used family. In SNAP, the local atomic environment is encoded by bispectrum components $\mathbf{R}_i$2. A linear per-atom model was used for aluminum,

$\mathbf{R}_i$3

with a ZBL short-range reference term (Kumar et al., 2023), whereas a quadratic form was used for carbon at extreme conditions,

$\mathbf{R}_i$4

to fit energies, forces, and stresses jointly to QMD data (Willman et al., 2022).

Low-dimensional descriptor models continue to be competitive, especially for chemically complex alloys. For ferromagnetic and liquid iron, machine-learned EAM, explicit two-body-plus-three-body terms, and combined EAM-plus-three-body terms were tabulated and evaluated by cubic spline interpolation, spanning three orders of magnitude in cost relative to SOAP-based GAP (Byggmästar et al., 2022). For Mo–Nb–Ta–V–W alloys, a tabulated GAP using two-body, three-body, and a scalar many-body density descriptor

$\mathbf{R}_i$5

reached $\mathbf{R}_i$6–$\mathbf{R}_i$7 meV/atom accuracy for crystals and $\mathbf{R}_i$8 meV/atom for liquids while being $\mathbf{R}_i$9–$\mathbf{G}_i$0 orders of magnitude faster than SOAP-based models (Byggmästar et al., 2022).

Recent software-centered work has broadened the model landscape toward equivariant graph neural networks. A 2025 MLIP library release included MACE, NequIP, and ViSNet as Flax modules implemented in JAX, while a 2025 study of minimalist MLIPs used Allegro, an $\mathbf{G}_i$1-equivariant neural network without message passing, alongside SOAP/GAP models (Brunken et al., 28 May 2025, Robredo-Magro et al., 21 Nov 2025). At the universal-model end, MACE-Osaka26 combined inorganic, molecular, and heavy-element datasets to produce an open-source model spanning 97 elements, described as the broadest elemental coverage to date (Kuroda et al., 3 Mar 2026).

3. Reference data, active learning, and database construction

Because ML-IAP behavior is database-limited, the generation and curation of reference configurations is a primary methodological problem. A central active-learning approach for linearly parameterized potentials uses the D-optimality criterion. Given a candidate configuration $\mathbf{G}_i$2, one computes

$\mathbf{G}_i$3

where $\mathbf{G}_i$4 is the current $\mathbf{G}_i$5 matrix of selected configurations. If $\mathbf{G}_i$6, the configuration is added to the training set. In that formulation, $\mathbf{G}_i$7 is the extrapolation grade, and values with all $\mathbf{G}_i$8 are interpreted as interpolation rather than extrapolation (Podryabinkin et al., 2016).

Other data-generation strategies target specific regions of the potential energy surface. For bulk crystalline matter near stationary points, an automated non-diagonal-supercell protocol was proposed to sample the small-displacement regime efficiently. Starting from a geometry-optimized primitive cell, one enumerates NDSCs for a target $\mathbf{G}_i$9-point grid, applies Gaussian random displacements with $$\mathbf{R}_i \rightarrow \mathbf{G}_i \overset{\mathcal{R}}{\rightarrow} E_i,$$0 Å, isotropic volume changes, and random shear strains through

$$\mathbf{R}_i \rightarrow \mathbf{G}_i \overset{\mathcal{R}}{\rightarrow} E_i,$$1

with $$\mathbf{R}_i \rightarrow \mathbf{G}_i \overset{\mathcal{R}}{\rightarrow} E_i,$$2 drawn componentwise from $$\mathbf{R}_i \rightarrow \mathbf{G}_i \overset{\mathcal{R}}{\rightarrow} E_i,$$3. This procedure was shown for Al, W, Mg, and Si, yielding phonon RMSE below $$\mathbf{R}_i \rightarrow \mathbf{G}_i \overset{\mathcal{R}}{\rightarrow} E_i,$$4 THz and elastic properties often within $$\mathbf{R}_i \rightarrow \mathbf{G}_i \overset{\mathcal{R}}{\rightarrow} E_i,$$5–$$\mathbf{R}_i \rightarrow \mathbf{G}_i \overset{\mathcal{R}}{\rightarrow} E_i,$$6 of DFT (Allen et al., 2022).

For chemically diverse systems, global structure search can replace manually curated enumeration. The Si–C UF3 model was trained predominantly on structures discovered by the GASP genetic algorithm over the entire silicon–carbon composition space, including both relaxed and unrelaxed structures, producing 41,145 GA structures plus 80 AIMD structures (MacIsaac et al., 2024). A different route is indirect learning: an accurate but slower teacher ML-IAP, such as GAP, generates both configurations and labels for a much larger secondary dataset, which is then used to train a faster student model such as MTP, without new quantum-mechanical calculations at the secondary stage (Morrow et al., 2021).

Recent work also formalizes data selection as a cost–accuracy optimization problem. In an application-specific qSNAP study for Be, 20,000 configurations were generated by information entropy maximization, and leverage-score sampling based on the hat matrix

$$\mathbf{R}_i \rightarrow \mathbf{G}_i \overset{\mathcal{R}}{\rightarrow} E_i,$$7

was used to identify informative training subsets. That study found that leverage sampling can reduce the required training-set size by up to an order of magnitude compared with random sampling, and that the optimal DFT precision should be matched to the intrinsic precision limit of the target ML-IAP rather than maximized indiscriminately (Baghishov et al., 6 Jun 2025).

4. Training objectives, benchmarking practice, and uncertainty quantification

Most modern ML-IAPs are trained against energies, forces, and often stresses or virials simultaneously. The review literature gives a representative loss combining these terms with regularization, while application-specific studies often use weighted least squares. In the Be qSNAP analysis, the loss was written as

$$\mathbf{R}_i \rightarrow \mathbf{G}_i \overset{\mathcal{R}}{\rightarrow} E_i,$$8

and one of the main findings was that, when lower-precision DFT data are used, increasing the force weight can lower both force and energy errors on high-precision test data because forces converge faster than total energies (Baghishov et al., 6 Jun 2025).

Model-specific fitting procedures differ substantially. Linearized Ti MLIPs used ridge regression on energies and forces (Takahashi et al., 2017). Carbon SNAP at extreme conditions used weighted linear regression through FitSNAP, with group weights optimized by a genetic algorithm in DAKOTA and final hyperparameters $$\mathbf{R}_i \rightarrow \mathbf{G}_i \overset{\mathcal{R}}{\rightarrow} E_i,$$9 Å and $\mathcal{R}$0 (Willman et al., 2022). The aluminum SNAP model used FitSNAP and DAKOTA to optimize both cutoff radius and training weights, yielding an optimal $\mathcal{R}$1 Å and $\mathcal{R}$2 (Kumar et al., 2023). BCM training for CNH molecules used SGPR experts and combined committee predictions by uncertainty-weighted aggregation, while the covariance loss

$\mathcal{R}$3

served as an uncertainty indicator for active sampling (Park et al., 2024).

Validation practice has broadened beyond train/test RMSE. In bcc Fe, a three-step protocol evaluated ML-IAPs by Pareto fronts of computational speed versus RMSE, a benchmark-property quality factor

$\mathcal{R}$4

and per-atom uncertainty during extended-defect simulations. That study found that $\mathcal{R}$5 correlates with RMSE, but reductions in RMSE below about $\mathcal{R}$6 meV/atom do not necessarily improve benchmark-property accuracy, and it explicitly argued for indicators beyond RMSE when assessing transferability (Zhang et al., 2023).

Uncertainty quantification is likewise model-family dependent. In D-optimal active learning, the extrapolation grade $\mathcal{R}$7 can be computed without QM data and was reported to be strongly correlated with model error (Podryabinkin et al., 2016). In Fe defect studies, safe regimes were defined by per-atom uncertainty below $\mathcal{R}$8 meV/atom for GAP and by per-atom extrapolation grade below $\mathcal{R}$9 for ACE-derived models (Zhang et al., 2023). This benchmarking practice has become particularly important when the intended applications involve dislocations, fracture, shock loading, or warm dense matter rather than equilibrium crystalline properties alone.

5. Accuracy, transferability, and representative domains of application

The most direct demonstrations of ML-IAP capability come from property-specific case studies. For elemental Ti, a pairwise MLIP trained on 2,700 DFT configurations and about 430,000 data points achieved an energy RMSE of $b_{n0}^{(i)}=\sum_j f_n(r_{ij}),$0 meV/atom with 2,925 coefficients, whereas the angular-dependent MLIP reached $b_{n0}^{(i)}=\sum_j f_n(r_{ij}),$1 meV/atom with 35,245 coefficients, and the latter accurately reproduced elastic constants, bulk moduli, and phonon dispersion curves, unlike EAM/MEAM and unlike the pairwise model on some elastic properties (Takahashi et al., 2017).

For extreme-condition carbon, a quadratic SNAP model was trained on a QMD database covering pressures from $b_{n0}^{(i)}=\sum_j f_n(r_{ij}),$2 to $b_{n0}^{(i)}=\sum_j f_n(r_{ij}),$3 TPa, temperatures from $b_{n0}^{(i)}=\sum_j f_n(r_{ij}),$4 to $b_{n0}^{(i)}=\sum_j f_n(r_{ij}),$5 K, densities from $b_{n0}^{(i)}=\sum_j f_n(r_{ij}),$6 to $b_{n0}^{(i)}=\sum_j f_n(r_{ij}),$7 g/cm$b_{n0}^{(i)}=\sum_j f_n(r_{ij}),$8, and 636 structures comprising 124,907 unique atomic environments. In two-phase MD, the resulting melting lines for diamond, BC8, and SC were within about $b_{n0}^{(i)}=\sum_j f_n(r_{ij}),$9 K, or $$b_{nl}^{(i)}=\sum_{j,k} f_n(r_{ij})f_n(r_{ik})\cos^l(\gamma_{jik}), \qquad l\ge 1.$$0, of QMD across $$b_{nl}^{(i)}=\sum_{j,k} f_n(r_{ij})f_n(r_{ik})\cos^l(\gamma_{jik}), \qquad l\ge 1.$$1–$$b_{nl}^{(i)}=\sum_{j,k} f_n(r_{ij})f_n(r_{ik})\cos^l(\gamma_{jik}), \qquad l\ge 1.$$2 GPa, and the model enabled simulations of 1 million-atom polycrystalline diamond (Willman et al., 2022). A different carbon study used GAP-20 with nested sampling up to $$b_{nl}^{(i)}=\sum_{j,k} f_n(r_{ij})f_n(r_{ik})\cos^l(\gamma_{jik}), \qquad l\ge 1.$$3 TPa and found remarkably good agreement with experiment up to $$b_{nl}^{(i)}=\sum_{j,k} f_n(r_{ij})f_n(r_{ik})\cos^l(\gamma_{jik}), \qquad l\ge 1.$$4 GPa, but also identified nonphysical graphite layer spacings and high-pressure phases not confirmed by DFT, making it a canonical example of both impressive transferability and clear extrapolation failure (Marchant et al., 2022).

For aluminum, a single SNAP ML-IAP was trained on DFT-MD data between $$b_{nl}^{(i)}=\sum_{j,k} f_n(r_{ij})f_n(r_{ik})\cos^l(\gamma_{jik}), \qquad l\ge 1.$$5 and $$b_{nl}^{(i)}=\sum_{j,k} f_n(r_{ij})f_n(r_{ik})\cos^l(\gamma_{jik}), \qquad l\ge 1.$$6 K at ambient density $$b_{nl}^{(i)}=\sum_{j,k} f_n(r_{ij})f_n(r_{ik})\cos^l(\gamma_{jik}), \qquad l\ge 1.$$7 g/cm$$b_{nl}^{(i)}=\sum_{j,k} f_n(r_{ij})f_n(r_{ik})\cos^l(\gamma_{jik}), \qquad l\ge 1.$$8, using 108- and 256-atom cells and about 100 snapshots per temperature. The reported prediction errors were $$b_{nl}^{(i)}=\sum_{j,k} f_n(r_{ij})f_n(r_{ik})\cos^l(\gamma_{jik}), \qquad l\ge 1.$$9 meV/atom MAE and $$\mathbf{R}_i \rightarrow \mathbf{G}_i \overset{\mathcal{R}}{\rightarrow} \mathbf{p}_i \overset{\Phi}{\rightarrow} E_i.$$0 meV/atom RMSE for energies, and $$\mathbf{R}_i \rightarrow \mathbf{G}_i \overset{\mathcal{R}}{\rightarrow} \mathbf{p}_i \overset{\Phi}{\rightarrow} E_i.$$1 eV/Å MAE and $$\mathbf{R}_i \rightarrow \mathbf{G}_i \overset{\mathcal{R}}{\rightarrow} \mathbf{p}_i \overset{\Phi}{\rightarrow} E_i.$$2 eV/Å RMSE for forces. Despite the training range ending at $$\mathbf{R}_i \rightarrow \mathbf{G}_i \overset{\mathcal{R}}{\rightarrow} \mathbf{p}_i \overset{\Phi}{\rightarrow} E_i.$$3 K, the model was used to examine thermal conductivity, diffusion coefficient, viscosity, sound velocity, and ion-ion structure factors up to about $$\mathbf{R}_i \rightarrow \mathbf{G}_i \overset{\mathcal{R}}{\rightarrow} \mathbf{p}_i \overset{\Phi}{\rightarrow} E_i.$$4 K, and diffusion trends were checked up to $$\mathbf{R}_i \rightarrow \mathbf{G}_i \overset{\mathcal{R}}{\rightarrow} \mathbf{p}_i \overset{\Phi}{\rightarrow} E_i.$$5 K (Kumar et al., 2023).

Transferability across chemical families has likewise been demonstrated in molecular systems. The BCM potential for organic nitrogen compounds partitioned CNH data into chemical-family experts such as amines, azoles, alkaloids, cyanides, hydrazines, imidazoles, and nitriles; it then showed high consistency with DFT on long-chain amines such as $$\mathbf{R}_i \rightarrow \mathbf{G}_i \overset{\mathcal{R}}{\rightarrow} \mathbf{p}_i \overset{\Phi}{\rightarrow} E_i.$$6 even though training used only short-chain amines. Reported test-set energy MAEs ranged from $$\mathbf{R}_i \rightarrow \mathbf{G}_i \overset{\mathcal{R}}{\rightarrow} \mathbf{p}_i \overset{\Phi}{\rightarrow} E_i.$$7 meV for azoles to $$\mathbf{R}_i \rightarrow \mathbf{G}_i \overset{\mathcal{R}}{\rightarrow} \mathbf{p}_i \overset{\Phi}{\rightarrow} E_i.$$8 meV for alkaloids, and force MAEs ranged from $$\mathbf{R}_i \rightarrow \mathbf{G}_i \overset{\mathcal{R}}{\rightarrow} \mathbf{p}_i \overset{\Phi}{\rightarrow} E_i.$$9 to $x$0 eV/Å across the listed groups (Park et al., 2024).

Defect mechanics has been a particularly stringent proving ground. In bcc Fe, GAP, PACE, MTP, and neural-network models were compared for dislocations and cracks. All optimized models reached energy RMSE $x$1 meV/atom; mean relative benchmark-property errors $x$2 were reported; the screw-dislocation Peierls barrier was predicted in the range $x$3–$x$4 meV per Burgers vector $x$5, compared with DFT values of $x$6–$x$7 meV/$x$8; and cleavage on the pre-existing crack plane was found to be the zero-temperature atomistic fracture mechanism of pure bcc Fe under mode-I loading, independent of ML package and training database (Zhang et al., 2023).

Large-scale materials applications extend beyond equilibrium thermodynamics. A SNAP ML-IAP for Ni–Mo alloys was used in LAMMPS to simulate polycrystalline Ni and Ni-10%Mo with average grain diameters from $x$9 to $E(x)=\sum_{j=1}^{m}\theta_j b_j(x),$0 nm and system sizes up to about 454,600 atoms. It reproduced Hall–Petch and inverse Hall–Petch behavior, linked strengthening at very small grain size to reductions in atomic von Mises strain by solute doping and annealing, and showed that annealing can weaken large-grain polycrystals because of large atomic movements in grain-boundary regions (Li et al., 2023). For Si–C, the UF3-based model was reported to capture peritectic decomposition of SiC and carbon dimer formation during SiC surface reconstruction, phenomena that the paper states cannot be captured with prior classical interatomic potentials for that material (MacIsaac et al., 2024).

6. Software ecosystems, scalability, and universal deployment

Methodological progress has been accompanied by increasingly integrated software stacks. The active-learning paper that introduced D-optimal selection published software, test cases, and usage examples at the mlip repository associated with MTP fitting and on-the-fly training (Podryabinkin et al., 2016). SNAP development for carbon and aluminum relied on FitSNAP, with DAKOTA used for hyperparameter or weight optimization (Willman et al., 2022, Kumar et al., 2023). The BCM workflow for CNH compounds integrated SGPR models with molecular dynamics through the AutoForce Python package (Park et al., 2024). The NDSC phonon-data workflow was accompanied by scripts and datasets intended for reproducibility and insertion into broader MLIP pipelines (Allen et al., 2022).

A more explicit unification effort appears in the 2025 MLIP library white paper. That library is written entirely in JAX, exposes three architectures—MACE, NequIP, and ViSNet—supports ASE and JAX-MD wrappers, and provides pre-trained organics models. On SPICE2-derived data, models were trained for 220 epochs with reported training times of 245 h for NequIP, 158 h for ViSNet, and 266 h for MACE on NVIDIA H100 GPUs. Runtime benchmarks reported, for example, $E(x)=\sum_{j=1}^{m}\theta_j b_j(x),$1 ms/step for MACE-large on 1UAO in JAX+JAX-MD versus $E(x)=\sum_{j=1}^{m}\theta_j b_j(x),$2 ms in Torch+ASE, and the paper summarized this as an approximately $E(x)=\sum_{j=1}^{m}\theta_j b_j(x),$3–$E(x)=\sum_{j=1}^{m}\theta_j b_j(x),$4 speed-up for JAX+JAX-MD over Torch+ASE (Brunken et al., 28 May 2025).

AI2Pot extends this unification toward high-throughput training and large-scale MD deployment. It supports MTP and NEP within a PyTorch/PyTorch-Lightning frontend, but re-engineers descriptor construction, energy, force, and virial evaluation as hand-crafted C++/CUDA operators shared between training and inference. Models can be deployed in ASE and LAMMPS, and benchmarked throughput on a single NVIDIA RTX 4090 reached 923,500 atom-steps/s for NEP and 355,300 atom-steps/s for MTP at 3.4 million atoms. Reported full training times on a Ge–Sb–Te dataset were 1.3 h for NEP and 0.75 h for MTP with batch size 32 (Liu et al., 8 Jul 2026).

At the universal-model scale, MACE-Osaka26 merged HE26 heavy-element data with MPtrj and OFF23 using total-energy alignment and produced a 97-element model trained on 32 NVIDIA A100 GPUs across 4 nodes. Reported test errors included an MPtrj energy MAE improvement from 33.7 to 29.3 meV/atom and an OFF23 RMSE reduction from 25.2 to 12.3 meV/atom relative to the cited Osaka24 baseline, while heavy-element compounds in HE26 achieved $E(x)=\sum_{j=1}^{m}\theta_j b_j(x),$5 meV/atom energy MAE and $E(x)=\sum_{j=1}^{m}\theta_j b_j(x),$6 meV/Å force MAE (Kuroda et al., 3 Mar 2026). This development places universal ML-IAPs alongside application-specific models as a distinct software and data-management paradigm.

7. Limitations, controversies, and evolving scope

The central limitation of ML-IAPs remains extrapolation. The 2021 review states the problem directly: outside the database domain, extrapolation to atomic environments not present or underrepresented in the training set leads to unpredictable or even unphysical results, because no underlying physical model constrains the potential in the way traditional analytic forms do (Mishin, 2021). Carbon phase-diagram studies illustrate this concretely: GAP-20 reproduced graphite, diamond, and melting behavior well up to about $E(x)=\sum_{j=1}^{m}\theta_j b_j(x),$7 GPa despite being trained mainly at zero pressure, yet at higher pressure it produced multiple minima in graphite interlayer spacing and stable strained phases not confirmed by DFT (Marchant et al., 2022).

A second controversy concerns what counts as sufficient validation. Several studies warn that low train/test errors are not enough. The Fe defect paper showed that significant quantitative predictions of dislocation glide paths and $E(x)=\sum_{j=1}^{m}\theta_j b_j(x),$8 can be sensitive to database, ML package, and cutoff radius, and argued that RMSE must be complemented by benchmark properties and uncertainty indicators (Zhang et al., 2023). The Be qSNAP study similarly argued that MLIP quality depends jointly on training-set selection, reference-data precision, energy-versus-force weighting, and model complexity, and used Pareto analysis rather than a single accuracy measure to identify efficient operating points (Baghishov et al., 6 Jun 2025).

A third controversy concerns whether ML-IAPs are restricted to interpolation-only use. Typical usage guidelines emphasize exhaustive training sets and caution against applying models to situations not represented in training. A 2025 study explicitly challenged that position by constructing SOAP/GAP and Allegro models from minimalist on-the-fly datasets—247 structures for BaTiO$E(x)=\sum_{j=1}^{m}\theta_j b_j(x),$9, 488 for BiFeO$\mathbf{R}_i$00, 1300 for PbZrO$\mathbf{R}_i$01, and 463 for HfO$\mathbf{R}_i$02, all with default or near-default hyperparameters—and showing qualitatively and quasi-quantitatively correct predictions for vortex–antivortex electric textures, ferroelectric switching paths, competing polymorphs, and the oIII ferroelectric HfO$\mathbf{R}_i$03 polymorph. The same study also noted that accuracy drops for states very far from the training domain (Robredo-Magro et al., 21 Nov 2025). This suggests that the operative question is not simply whether extrapolation occurs, but how aggressively it is monitored, constrained, and validated.

The likely long-term direction is therefore plural rather than singular. Active learning based on extrapolation grades, physically informed hybrids that couple ML with physics-based potentials, optimal data-generation schemes such as NDSC or leverage sampling, distilled teacher–student workflows, universal many-element foundation-style models, and high-performance deployment frameworks all address different points on the same design surface: accuracy, transferability, interpretability, and computational cost (Podryabinkin et al., 2016, Mishin, 2021, Morrow et al., 2021, Kuroda et al., 3 Mar 2026). The modern ML-IAP field is best understood not as a single method, but as a family of atomistic surrogate models whose success depends on matching representation, database, uncertainty control, and simulation target with unusual precision.