Machine Learning Interatomic Potentials

• Machine Learning Interatomic Potentials (MLIAPs) are surrogate models that learn atomic energy contributions from quantum-mechanical data to represent material energy landscapes.
• They utilize diverse methodologies including descriptor-based regressors and equivariant graph neural networks to balance accuracy and computational efficiency.
• MLIAPs enable large-scale molecular dynamics and phase-space exploration while addressing challenges like long-range interactions and model transferability.

Machine Learning Interatomic Potentials (MLIPs, or MLIAPs) are surrogate models for the potential-energy surface of materials. They learn energies, forces, and, in some settings, stresses or charges from quantum-mechanical reference data, typically DFT, with the aim of combining nearly ab initio accuracy with the efficiency required for large-scale molecular dynamics, structure relaxation, and phase-space exploration (Mishin, 2021, Tang et al., 10 Jun 2025, Brunken et al., 21 May 2026). In contemporary usage, the term covers descriptor-based regressors, neural-network potentials, equivariant graph neural networks, and more explicitly physics-based constructions that introduce charge, density, or ensemble electronic degrees of freedom (Mishin, 2021, Atlas, 5 Mar 2026).

1. Physical formulation and scope

Most MLIAPs adopt a local energy decomposition in which the total energy is written as a sum of atomic contributions,

$E_{\text{tot}} = \sum_{i=1}^N E_i(\mathcal{E}_i),$

where $\mathcal{E}_i$ is the local environment of atom $i$ within a cutoff radius, and forces are obtained as

$$\mathbf{F}_i = -\frac{\partial E_{\text{tot}}}{\partial \mathbf{R}_i}.$$

This construction preserves translational invariance, and, when the local environment representation is chosen appropriately, it can also enforce invariance to global rotation and permutation of identical atoms (Tang et al., 10 Jun 2025, Mishin, 2021).

The standard motivation is the long-standing tradeoff between ab initio molecular dynamics and empirical force fields. Ab initio MD, with forces from Kohn–Sham DFT, is accurate but scales badly with system size and is limited to roughly $10^2$–$10^3$ atoms and picosecond times, whereas empirical force fields are cheap and scalable but often inaccurate or non-transferable, especially in multicomponent systems with mixed ionic, covalent, and metallic bonding (Tang et al., 10 Jun 2025). MLIAPs were introduced precisely to break this tradeoff by learning the mapping from atomic coordinates to energies and forces from a database of DFT calculations (Mishin, 2021).

A seminal conceptual template is the Behler–Parrinello neural network potential, first proposed in 2007, in which each element is associated with a neural network acting on atom-centered descriptors and the total energy is assembled as a sum of atomic energies (Tang et al., 10 Jun 2025). This local decomposition remains the dominant organizing principle even in later architectures such as DeePMD, MACE, NequIP, Allegro, and DPA3, although the descriptors and regression layers differ substantially (Tang et al., 10 Jun 2025, Brunken et al., 21 May 2026, Liu et al., 9 Mar 2026).

2. Representations and model classes

The main distinction among MLIAPs is how the local environment is represented and how the energy map is parameterized. A concise taxonomy used across recent work is summarized below.

Family	Local representation	Examples named in the literature
Descriptor-based local models	atom-centered symmetry functions, SOAP, bispectrum, moment tensor descriptors, ACE, scalar many-body density descriptor	Behler–Parrinello, GAP, SNAP, MTP, ACE, qSNAP, tabGAP
Graph-based neural models	message passing embeddings, spherical harmonics, irreducible representations, invariant or equivariant node features	SchNet, PaiNN, NequIP, Allegro, MACE, ViSNet, eSEN, DPA3
Explicit electronic or density-based models	atom-in-molecule densities, ensemble weights, QEq charges, global charge redistribution	ECT-EAM, NequIP-LR

Descriptor-based models remain important because low-dimensional descriptors can be unusually data-efficient in complex alloys. In Mo–Nb–Ta–V–W, a computationally fast Gaussian approximation potential with two-body, three-body, and a scalar many-body density descriptor based on the embedded atom method significantly outperformed SOAP in data efficiency, accuracy, and speed, and could be tabulated into a “tabGAP” form (Byggmästar et al., 2022). Related application-specific studies with qSNAP also show that the descriptor complexity $2J_{\max}$ directly controls the accuracy–cost tradeoff, with the number of descriptors increasing from 105 at $2J_{\max}=4$ to 4186 at $2J_{\max}=10$ and the speed decreasing from $4 \times 10^5$ to $\mathcal{E}_i$0 atoms timestep$\mathcal{E}_i$1 (Baghishov et al., 6 Jun 2025).

Equivariant graph neural networks encode local geometry through radial basis functions, spherical harmonics, tensor products, or message passing on atomic graphs. In mlip v2, the supported families include MACE, NequIP, ViSNet, and eSEN; in these models, node features are updated through equivariant message passing and atomwise energies are summed to a total energy, while forces are obtained by automatic differentiation in JAX (Brunken et al., 21 May 2026). Mixture-of-Experts extensions have now also been developed for MLIAPs: in DPA3-based models, sparse activation, shared experts, and element-wise routing yielded state-of-the-art accuracy on OMol25, OMat24, and OC20M, and the learned routing patterns aligned with periodic-table trends (Liu et al., 9 Mar 2026).

A separate line of work questions whether exact rotational equivariance and conservative-force constraints must always be hard-wired. Large PET models trained on OMat24, MPtrj, Alexandria, and SPICE show that rotationally unconstrained models can be superior in accuracy and speed when trained on large datasets, although inference-time averaging or symmetry projection may be required to recover observables consistent with the relevant physical symmetries (Bigi et al., 22 Jan 2026). By contrast, explicitly physics-based “latent space” formulations build the hidden variables from DFT itself: in the latent-space design program, atom-in-molecule densities are expanded as ensembles of isolated-atom charge and excitation states, leading to the ensemble charge-transfer embedded atom method as a DFT-consistent interatomic potential (Atlas, 5 Mar 2026).

3. Training data, losses, and workflow design

For MLIAPs, the training set is often the decisive design choice. The standard supervised target consists of DFT energies, forces, and sometimes stresses for equilibrium, distorted, and high-temperature structures (Mishin, 2021). Recent case studies make this operational point explicit. In the La–Si–P system, a DeePMD-kit potential was trained on 71,800 distorted crystals and 130,000 liquid configurations generated with VASP, PAW pseudopotentials, and GGA-PBE; the resulting model reached an Energy RMS error of approximately 12 meV/atom and a Force RMS error of approximately $\mathcal{E}_i$2–$\mathcal{E}_i$3 eV/Å on 1320 randomly selected configurations from 66 systems (Tang et al., 10 Jun 2025).

That study also illustrates a now-common iterative workflow. The initial model was trained with all configurations having equal usage probability; the authors then increased the usage probability of problematic phases, initialized subsequent training rounds from the previous iteration, and actively enriched the dataset with new DFT configurations in regimes where the model failed, such as crystal–liquid transition states relevant to melting (Tang et al., 10 Jun 2025). This is closely aligned with the older review view that MLIPs are powerful numerical interpolators whose transferability depends on whether the database spans the configurations sampled in production simulations (Mishin, 2021).

Several recent works address training-set construction and supervision more systematically. An information-entropy–driven generation strategy maximizes $\mathcal{E}_i$4, where $\mathcal{E}_i$5 is the covariance of bispectrum features, to produce a material-agnostic dataset with 30,349 configurations and 179,367 atoms; MLIAPs trained on that dataset were reported to be extremely robust over very broad swaths of configuration space, even without dataset fine-tuning or hyper-parameter optimization (Subramanyam et al., 2024). For application-specific potentials, leverage-score sub-sampling in qSNAP made it possible to approach the limiting error with a few hundred configurations and to reduce DFT cost by roughly an order of magnitude relative to random sampling, while the same study showed that the usefulness of higher-precision DFT depends on the intrinsic error floor of the chosen model class (Baghishov et al., 6 Jun 2025).

Loss design is also evolving beyond plain energy–force regression. A physics-informed weakly supervised framework introduced PITC and PISC losses based on a first-order Taylor expansion of the energy and on conservative-force path consistency, and reported reduced energy and force errors, often lower by a factor of two, especially in sparse-data regimes and in fine-tuning settings where only energies are available (Takamoto et al., 2024). At the workflow level, AMLP combines large-language-model agents for DFT setup recommendation, automated input generation, AIMD output parsing, MACE-ready HDF5 construction, and ASE-based validation; on acridine polymorphs, straightforward fine-tuning of a foundation model achieved mean absolute errors of approximately 1.7 meV/atom in energies and approximately 7.0 meV/Å in forces (Lahouari et al., 25 Sep 2025).

4. Long-range interactions, charge transfer, and interpretability

A persistent limitation of standard local MLIAPs is their difficulty with long-range electrostatics, charge transfer, and compositional heterogeneity. This is especially pronounced in systems where multiple charge states share the same composition, or where explicit global redistribution of electronic charge controls the energetics (Maruf et al., 23 Mar 2025, Brunken et al., 21 May 2026).

One response is to augment short-range equivariant networks with explicit charge equilibration. NequIP-LR combines a NequIP-style E(3)-equivariant message-passing network for short-range interactions with a global charge equilibration scheme based on predicted atomic electronegativities. Its total energy is decomposed into a short-range term from a second equivariant network and a Coulomb term computed from Gaussian charges, with Ewald summation used in periodic systems (Maruf et al., 23 Mar 2025). Across datasets including $\mathcal{E}_i$6 on MgO(001), benzotriazole on Cu(111), $\mathcal{E}_i$7, and a protonated carbon chain, it outperformed both short-range equivariant and long-range invariant baselines in energy and force prediction, and it handled multiple charge states that local composition-based models cannot distinguish (Maruf et al., 23 Mar 2025).

A related but more general software-level strategy appears in mlip v2, where all architectures support partial charge prediction, explicit charge conservation, global charge conditioning, and a modular long-range soft-core Coulomb term following PhysNet (Brunken et al., 21 May 2026). The same framework also integrates Hessian-augmented training, NPT via a Monte Carlo barostat, and NEB, illustrating how long-range physics and advanced simulation algorithms are being folded into unified MLIP ecosystems rather than added as ad hoc post-processing (Brunken et al., 21 May 2026).

At the opposite end of the design spectrum, the latent-space approach fixes the latent variables by quantum mechanics rather than learning them numerically. There, the system density is decomposed into atom-in-molecule densities, each atom-in-molecule density is written as an ensemble over isolated atomic charge and excitation states, and the same ensemble weights appear in both the density and energy functionals (Atlas, 5 Mar 2026). This makes the representation interpretable by construction: the hidden variables correspond to explicit physical states, and charge transfer becomes a reweighting of ensemble components rather than an implicit feature in a black-box embedding (Atlas, 5 Mar 2026). A plausible implication is that long-range electrostatics and interpretability need not be separate research programs in MLIAPs; they can arise from the same explicit electronic degrees of freedom.

5. Validation regimes and application domains

Recent work increasingly treats held-out energy and force errors as necessary but insufficient. Transferability is now commonly assessed through observables that are not directly targeted in the loss: energy–volume curves, pair-correlation functions, melting points, nucleation pathways, phonons, defect energetics, NEB barriers, and long-time MD stability (Tang et al., 10 Jun 2025, Brunken et al., 21 May 2026, Robredo-Magro et al., 21 Nov 2025).

The La–Si–P DeePMD model is a clear illustration. It reproduced energy–volume curves for elemental, binary, and ternary crystalline phases over $\mathcal{E}_i$8–$\mathcal{E}_i$9 of the equilibrium volume and matched liquid pair-correlation functions at 2500 K across metallic, covalent, and mixed compositions (Tang et al., 10 Jun 2025). It was then used to predict melting by solid–liquid coexistence, giving $i$0 K for LaSiP$i$1 against an experimental 1308 K and $i$2 K for La$i$3SiP$i$4 against a DSC onset near 1330 K, while also enabling 4096-atom nucleation simulations that showed LaP crystallization from La–Si–P melts (Tang et al., 10 Jun 2025). These are not mere train/test metrics; they are thermodynamic and kinetic observables emerging from nanosecond-scale MD.

For molecular crystals, AMLP and MACE were validated on acridine polymorphs by comparing lattice energies, structure relaxations, NVE energy conservation, NVT radial distribution functions, and orientational order. Three committees achieved approximately 2 meV/atom energy MAE and approximately 7 meV/Å force MAE, reproduced DFT geometries with mean RMSD around $i$5–$i$6 Å, and showed NVE energy conservation on the order of $i$7 over 40 ps in supercells of 1104–2944 atoms (Lahouari et al., 25 Sep 2025). The same study also demonstrated that committees with similar static errors can differ sharply in dynamical robustness, as one of the three MACE models destabilized several polymorphs above room temperature (Lahouari et al., 25 Sep 2025).

A more aggressive claim about extrapolation comes from minimalist MLIPs trained on only a few hundred configurations in ferroelectric oxides. With default SOAP-based GAP and Allegro models trained on 247 to 1300 DFT structures, the authors reported correct unstable phonon modes for cubic reference phases, a vortex–antivortex electric texture in BaTiO$i$8, a BiFeO$i$9 switching barrier reproduced to within 1 meV/f.u., and the emergence of AFE$$\mathbf{F}_i = -\frac{\partial E_{\text{tot}}}{\partial \mathbf{R}_i}.$$0, AFE$$\mathbf{F}_i = -\frac{\partial E_{\text{tot}}}{\partial \mathbf{R}_i}.$$1, ferrielectric, and ferroelectric polymorphs in PbZrO$$\mathbf{F}_i = -\frac{\partial E_{\text{tot}}}{\partial \mathbf{R}_i}.$$2 and HfO$$\mathbf{F}_i = -\frac{\partial E_{\text{tot}}}{\partial \mathbf{R}_i}.$$3 that were absent from training (Robredo-Magro et al., 21 Nov 2025). This does not remove the standard caveat that MLIAPs are primarily interpolative, but it does suggest that carefully chosen local information can sometimes support limited serendipitous discovery.

Practical deployment also increasingly includes acceleration strategies. ML-MIX introduces spatial mixing of an expensive and a cheap potential inside LAMMPS, with core, blending, and buffer regions and constrained linear fitting of the cheap model to the expensive reference in relevant regions of configuration space (Birks et al., 26 Feb 2025). On 8,000-atom systems, the method achieved up to an 11x speedup without sacrificing accuracy on minimum energy paths or defect diffusion, and for larger domains the achievable speedup approached the ratio between cheap and expensive potential speeds (Birks et al., 26 Feb 2025). This suggests that the “MLIP bottleneck” is no longer only a question of model architecture; it is also a question of spatial allocation of fidelity.

6. Limitations, controversies, and active directions

The core limitation emphasized across the literature is still extrapolation. The 2021 review framed MLIPs as high-dimensional interpolators whose accuracy inside the training domain can reach $$\mathbf{F}_i = -\frac{\partial E_{\text{tot}}}{\partial \mathbf{R}_i}.$$4 meV/atom relative to DFT, but whose extrapolation is often poor because the learned model is mathematical rather than explicitly physical (Mishin, 2021). Recent work complicates, rather than overturns, that view.

One controversy concerns whether stronger physical constraints are always beneficial. PET results indicate that rotationally unconstrained models can match or surpass fully constrained models in the large-data regime, especially when inference-time rotational averaging or symmetry projection is used to restore physically consistent observables (Bigi et al., 22 Jan 2026). A different controversy concerns how much training coverage is really required: the minimalist ferroelectric study shows that somewhat trivial and easy-to-compute models can predict complex structural behaviors qualitatively and quasi-quantitatively correctly, whereas more traditional guidance emphasizes exhaustive coverage and warns against out-of-domain use (Robredo-Magro et al., 21 Nov 2025, Mishin, 2021). A cautious synthesis is that extrapolation capacity is not binary; it depends on descriptor dimensionality, data geometry, and the extent to which the target phenomenon reuses local motifs already present in training.

A more specific and currently sharp limitation is systematic energy underprediction. A large-scale benchmark over more than 12 million calculations and more than 150,000 inorganic crystals found that nine frontier MLIAPs systematically underpredict total energy, formation energy, and especially energy above hull, even though more than 90% of the test structures are in the training data (Nong et al., 21 Jul 2025). The mean absolute errors for $$\mathbf{F}_i = -\frac{\partial E_{\text{tot}}}{\partial \mathbf{R}_i}.$$5 exceed approximately 30 meV/atom even for the best model, and the errors become especially severe for structures with high symmetry degrees of freedom, where subtle atomic displacements strongly affect the energy landscape (Nong et al., 21 Jul 2025). The authors attribute this to insufficient handling of lattice symmetry and Wyckoff-site symmetry, and argue for explicit DOF encoding or symmetry-regularized loss functions (Nong et al., 21 Jul 2025). This directly challenges the practice of equating low single-shot energy MAE on DFT-relaxed structures with reliable thermodynamic stability prediction.

Another persistent limitation is bias inherited from the reference method. In La–Si–P, melting temperatures were systematically underpredicted by approximately 5–20%, and the authors explicitly linked part of this error to the GGA-PBE reference itself (Tang et al., 10 Jun 2025). More generally, any MLIAP trained to DFT reproduces DFT’s systematic errors unless the model or training protocol adds corrective structure.

Active development is therefore moving along multiple fronts at once: mixture-of-experts scaling for broader chemical coverage and parameter efficiency (Liu et al., 9 Mar 2026); charge-aware and long-range models for systems with explicit electronic response (Maruf et al., 23 Mar 2025); weakly supervised and physics-informed losses for sparse-data regimes (Takamoto et al., 2024); unified open-source infrastructures for training, inference, NPT, NEB, Hessians, and benchmarking (Brunken et al., 21 May 2026); and more explicitly physics-based latent spaces that connect atomic and electronic length scales through ensemble densities and DFT constraints (Atlas, 5 Mar 2026). The cumulative trend is not toward a single universal architecture, but toward a layered ecosystem in which representation, data generation, physical constraints, simulation algorithmics, and software design are all treated as first-class components of the potential itself.