The core problem: you need the potential-energy surface
Almost everything atomistic simulation computes — molecular dynamics trajectories, free energies, phase diagrams, reaction rates, spectra — is downstream of one object: the potential-energy surface (PES), the total energy E({R}) as a function of every atomic position, together with the forces F_i = −∂E/∂R_i that follow from it. Get the PES right and the dynamics follow. Get it wrong and no amount of sampling helps.
There have historically been two ways to obtain it, and they fail in opposite directions.
Quantum-mechanical methods solve for the electronic structure at each configuration. Density-functional theory is the workhorse, implemented in mature plane-wave codes [1] [2], and it is accurate enough to be treated as reference data across chemistry and materials science. Its problem is cost: the scaling with system size makes long trajectories of large systems infeasible, which is exactly the statement that motivated the first neural-network PES work — "the accurate description of chemical processes often requires the use of computationally demanding methods like density-functional theory, making long simulations of large systems unfeasible" [3].
Classical force fields replace the electrons with a fixed functional form — bonded terms, Lennard-Jones, point charges — fitted to a limited set of data. They are cheap enough for million-atom, microsecond simulations in engines built for exactly that [4] [5] [6], and they underpin a great deal of applied work including pharmaceutical modelling [7]. Their problem is the fixed functional form: accuracy is bounded by what the form can express, transferability outside the fitting regime is poor, and bond breaking requires a specially constructed reactive form such as REBO or ReaxFF [8] [9].
Machine-learning interatomic potentials (MLIPs) target the gap. They fit the PES directly from quantum reference data using a flexible regressor with no fixed functional form, aiming for near-DFT accuracy at near-classical cost. The two founding papers state the bargain precisely: a neural-network representation of a DFT PES that "provides the energy and forces as a function of all atomic positions in systems of arbitrary size and is several orders of magnitude faster than DFT" [3]; and a kernel model whose "models do not have a fixed functional form and hence are capable of modeling complex potential energy landscapes… systematically improvable with more data," saving "orders of magnitude in computational cost" on the long trajectories needed for finite-temperature properties [10].
Two structural constraints shape every architecture that followed. First, the energy must be extensive and size-transferable, which is achieved by decomposing the total energy into a sum of atomic contributions, each a function of that atom's local environment within a cutoff [3] [11]. Second, physical symmetry must be respected: the energy is invariant under translation, rotation, and permutation of identical atoms, while forces rotate with the system. Whether a model imposes these by construction or has to learn them from data turns out to be the single largest determinant of data efficiency.
Method lineage
Family 1: descriptor-based neural networks
The Behler–Parrinello construction encodes each atom's neighbourhood in hand-designed symmetry functions — radial and angular features that are invariant by construction — and feeds them to a small per-element neural network whose outputs sum to the total energy [3]. The approach was demonstrated on bulk silicon against both empirical potentials and DFT, and the review a few years later documented the transition from low-dimensional to genuinely high-dimensional systems, with MD of thousands of atoms and energies and forces "obtained with very high accuracy," while being candid that development "requires substantial computational effort as typically thousands of reference calculations are required" [11].
The known limitation of this family is locality. Because each atomic energy depends only on a finite cutoff sphere, long-range electrostatics and global charge transfer are invisible. Fourth-generation high-dimensional NNPs address this with a charge-equilibration scheme built on environment-dependent atomic electronegativities, correctly describing global charge distributions in systems that purely local models get wrong [12].
Family 2: kernel methods and Gaussian processes
Gaussian Approximation Potentials take the non-parametric route: rather than a network, a Gaussian-process regressor over a smooth representation of the atomic environment, systematically improvable as data accumulates [10]. The statistical machinery is standard GP regression specialized to atomistic data [13], and the practical appeal is that a GP supplies a predictive variance alongside its prediction — an uncertainty estimate for free, which turns out to matter for active learning. GAP has been benchmarked head-to-head against permutationally invariant polynomials and neural networks on the same PES-fitting tasks [14] and used to build general-purpose potentials for hard cases such as silicon across its phase diagram [15].
Family 3: systematic linear expansions
A parallel line asks for a complete, systematically convergent basis for the local environment rather than a black-box regressor. Moment tensor potentials pose it as a mathematical approximation problem and construct a class of potentials that is systematically improvable by design [16]. The spectral neighbor analysis method (SNAP) expands the neighbour density in bispectrum components with linear fitting [17] [18]. The atomic cluster expansion (ACE) generalizes these into a complete descriptor whose evaluation cost scales linearly with the number of neighbours irrespective of expansion order, converging smoothly to meV accuracy and — once nonlinear functions of the expansion are introduced — reaching accuracy comparable to state-of-the-art ML potentials [19]. ACE now has performant production implementations [20] [21] [22]. The attraction of this family is interpretability and control: the approximation error can be reasoned about rather than only measured.
Family 4: message-passing and graph neural networks
Hand-designed descriptors cap what a model can see. Graph neural networks instead learn the representation, treating atoms as nodes and neighbours as edges and iterating learned messages. SchNet introduced continuous-filter convolutional layers appropriate to atoms at arbitrary positions rather than on a grid, predicting properties across chemical space and producing energy-conserving force fields for MD [23] [24] [25]. Because plain distance-based messages discard angular information, follow-up architectures reintroduce it — directional message passing over molecular graphs [26] [27], and Newtonian-mechanics-inspired message construction [28].
Family 5: E(3)-equivariance — the step change
The decisive idea is that symmetry should be a property of the architecture, not something learned. Most earlier symmetry-aware models are invariant: they act only on scalars, which is correct for energy but throws away the directional structure of the environment. NequIP instead uses E(3)-equivariant convolutions acting on geometric tensors, so internal features rotate with the system and only the final energy is invariant — "a more information-rich and faithful representation of atomic environments" [29]. The empirical consequence is the headline result of the field: state-of-the-art accuracy "with up to three orders of magnitude fewer training data, challenging the widely held belief that deep neural networks require massive training sets," data-efficient enough to fit potentials to high-order quantum-chemical reference [29]. Symmetry is, in effect, a very strong prior — and priors buy data.
Two directions extend it. MACE attacks the cost side: standard MPNNs pass only two-body messages, tying expressivity to depth and therefore to communication cost. Using higher body-order messages — four-body — cuts the required message-passing iterations to two, giving a fast, highly parallelizable model at or beyond state-of-the-art accuracy on rMD17, 3BPA and AcAc [30], subsequently evaluated across amorphous carbon, universal materials modelling, small-molecule chemistry, large molecules and liquid water, and shown to be data-efficient enough to reproduce experimental vibrational spectra when trained on as few as 50 reference configurations [31]. Allegro attacks the scaling side: atom-centred message passing limits accessible length scales, so it builds a strictly local equivariant model from iterated tensor products of learned equivariant representations with no atom-centred message passing, keeping accuracy while scaling to large simulations [32]. Related equivariant designs cover tensorial properties and spectra [33], attention/transformer formulations [34] [35] [36], and the design space itself has been mapped as a unified framework [37] [38].
Family 6: transferable potentials across chemical space
The families above are typically fitted per-system. ANI targeted breadth instead: modified Behler–Parrinello symmetry functions built into atomic environment vectors, trained on data spanning both configurational and conformational space, giving "an accurate and transferable potential for organic molecules" — DFT accuracy at force-field cost [39]. That programme is fundamentally a data programme: 20 million off-equilibrium conformations for ANI-1 [40], then coupled-cluster and DFT property sets [41]. Its most instructive result is methodological — transfer learning from a large DFT-trained model onto a smaller gold-standard CCSD(T)/CBS set yields ANI-1ccx, approaching coupled-cluster accuracy on reaction thermochemistry, isomerization and drug-like torsions while remaining "billions of times faster" [42]. Accuracy ceiling and data volume can be decoupled. The same recipe now runs on equivariant backbones for organic chemistry [43] [44], supported by purpose-built datasets such as SPICE, over 1.1 million conformations of drug-like molecules, dimers, dipeptides and solvated amino acids at ωB97M-D3(BJ)/def2-TZVPPD [45].
Family 7: universal / foundation potentials
The current frontier extends transferability from organic chemistry to the whole periodic table, trained on the accumulated output of high-throughput DFT [46] [47]. M3GNet provides a universal graph deep-learning interatomic potential for the periodic table [48]; a universal NNP covering arbitrary combinations of 45 elements followed [49]; and CHGNet is pretrained on the Materials Project Trajectory Dataset — energies, forces, stresses and magnetic moments from over a decade of DFT on more than 1.5 million inorganic structures — with magnetic moments included explicitly so the model can represent electronic state alongside ionic rearrangement [50]. The emerging usage pattern mirrors the rest of machine learning: pretrain broadly, then fine-tune on a small system-specific set for quantitative accuracy [51].
A minimal model: the diatomic potential-energy curve
Before the high-dimensional machinery, it is worth being concrete about what a "potential" is, because the one-dimensional case contains most of the intuition.
For a diatomic molecule the PES collapses to a single coordinate, the internuclear separation r, and the standard closed form is the Morse potential:
V(r) = D_e · (1 − e^(−a(r − r_e)))²
with three parameters: r_e, the equilibrium bond length where the curve is at its minimum; D_e, the well depth, i.e. the dissociation energy measured from the bottom of the well; and a, which sets the curvature and therefore the vibrational frequency. Its two essential features are that it is anharmonic — asymmetric about the minimum, steeply repulsive at short range from nuclear and Pauli repulsion, and flattening to a finite asymptote as r → ∞ because the bond simply breaks — and that it is smooth and differentiable, so forces come from −dV/dr.
The contrast that makes the point is the harmonic approximation, the second-order Taylor expansion about the minimum:
V_harm(r) = k/2 · (r − r_e)², with k = 2·D_e·a²
The harmonic form is exact in the limit of infinitesimal displacement and is what a classical force field's bond term typically uses. It is also qualitatively wrong everywhere else: it is symmetric, it rises without bound, and it therefore cannot dissociate. That failure is the whole reason a fixed functional form is a liability, and the reason MLIPs are built without one.
Grounding and honesty. The Morse functional form itself is textbook and predates this substrate's coverage window, so it is presented here as the standard closed-form diatomic potential rather than attributed to a paper in the collected corpus — an illustrative model, explicitly flagged. What is grounded is its continued use as the working approximation to a diatomic PES parameterized exactly by the quantities named above: substrate work modelling diatomics under external electric fields approximates the potential-energy surface by "a Morse function," predicting response from field-free observables including "the equilibrium bond length, the bond dissociation energy, the polarizability and dipole moment functions, and the vibrational frequency," validated against direct calculations on nine diatomics [52]. The parameter values used in the figure below are illustrative teaching values, not a fit to any particular molecule.
What generalizes and what does not. Three things carry over from this picture to a real MLIP: energy is a smooth function of geometry, forces are its negative gradient, and the physically important regions are the ones away from equilibrium. Two things do not. First, dimensionality — a diatomic PES is a curve; an N-atom PES lives in 3N−6 internal coordinates, which is why the atomic-energy decomposition and a local cutoff are structural necessities rather than conveniences [3] [11]. Second, the parameterization — a Morse curve has three interpretable parameters, while an MLIP has 10⁴–10⁷ uninterpretable ones, which is exactly the trade that buys flexibility and costs extrapolation guarantees.
Practical threads
Training data and active learning
An MLIP is only as good as the configurations it was shown, and generating them is the dominant cost [11]. Random sampling is wasteful, so the field converged on active learning: run dynamics with the current model, detect configurations where it is untrustworthy, compute quantum reference only for those, retrain, repeat.
The on-the-fly formulation makes forces "either predicted by Bayesian inference or, if necessary, computed by on-the-fly quantum-mechanical calculations and added to a growing ML database, whose completeness is thus never required" — with progressively fewer quantum calls as chemical processes recur [53]. DP-GEN structures the same loop into explicit exploration, labelling and training stages, producing uniformly accurate PES models "with a minimal number of reference data" [54]. The pattern has since been applied to amorphous and liquid oxides [55], to reactive catalysis with Bayesian force fields [56], and sharpened into hyperactive variants that bias exploration toward informative configurations [57] [58] [59].
Uncertainty quantification
Active learning presupposes a usable uncertainty estimate, and it is also what makes an MLIP safe in production — the model needs to say when it has left its training distribution. Gaussian processes give a variance natively [10] [13]; neural networks need it added, typically via committees of independently trained models whose disagreement proxies for error. The literature has developed quantitative uncertainty metrics that demonstrably control error in network-driven chemical discovery [60], uncertainty-aware message-passing architectures [27], and cheap single-model estimators that avoid paying for an ensemble at inference [61].
Extrapolation and transferability limits
The honest framing of the field's central weakness: MLIPs interpolate superbly and extrapolate unpredictably, and unlike a physically motivated functional form they have no correct asymptotic behaviour to fall back on. Locality is the sharpest instance — purely local models cannot see long-range charge transfer, which motivated fourth-generation charge-equilibration schemes [12] — and strictly local architectures make the accuracy/length-scale trade explicit rather than hiding it [32]. Transferability is therefore an empirical property to be measured, which is what the universal-potential and fine-tuning literature is really testing [48] [50] [51] [31].
Integration into MD engines
A potential that cannot be called from a production integrator is a benchmark, not a tool. The path to adoption has been through the established engines: LAMMPS, explicitly designed so that "adding code for new interatomic potentials" is a supported extension point [4] [18]; large-scale parallel and GPU-accelerated MD infrastructure [5] [6]; and universal force engines that decouple the dynamics driver from the energy evaluator [62]. Downstream property calculations — phonons and thermal transport, for instance — then run unchanged on top [63].
Benchmarks, and their discontents
The field's standard metric is force and energy RMSE on held-out configurations [64] [65]. The most useful critical result in the substrate argues this is measuring the wrong thing: the practical use case is producing realistic MD trajectories, and a benchmark suite spanning water, organic molecules, a peptide and materials shows "how the commonly benchmarked force accuracy is not well aligned with relevant simulation metrics," demonstrating when and how state-of-the-art methods fail [66]. Low force error does not guarantee a stable trajectory, a correct radial distribution function, or a physical free-energy landscape. Purpose-built training and evaluation datasets [45] [41] and architecture evaluations across genuinely different domains [31] are the partial answer.
Frontier and open problems
Foundation models for atomistic simulation. The pretrain-then-fine-tune pattern is now established for materials [48] [50] [49] [51] and organic chemistry [43], and the open questions are the familiar ones in a new domain: what the scaling behaviour is, how to evaluate a general model without a general benchmark, and whether a universal potential fine-tuned on a small set beats a bespoke potential trained on the same budget.
Long-timescale dynamics. Making each force evaluation cheap does not solve the timescale problem — MD is still limited by femtosecond timesteps against phenomena that take microseconds and longer. The complementary approaches are enhanced sampling of rare events [67] and learning coarse-grained variables so that the accelerated model and the resolution reduction are trained together [68] [69].
Reactive chemistry. Bond breaking was the historical province of specially constructed reactive force fields [8] [9]. MLIPs handle it naturally in principle, since nothing in the fit presumes a bond topology, but reactive regions are rare in unbiased sampling — which is precisely why active learning coupled to enhanced sampling has become the standard recipe for catalysis and surface chemistry [56] [70] [71] [72].
Coupling to experiment, and to electronic structure. The convincing validations are increasingly experimental observables rather than held-out DFT: infrared and vibrational spectra from ML-driven MD [73] [31], the thermodynamics of liquid and solid water [74], and spectroscopy of aqueous interfaces including quantum nuclear effects [75]. A parallel line learns the electronic structure itself rather than only the energy — transferable electron densities [76], deep-learning DFT Hamiltonians [77] [78], improved exchange-correlation functionals [79], and direct neural solution of the many-electron Schrödinger equation [80]. Taken together these suggest the eventual target is not a faster surrogate for DFT but a learned replacement for parts of the electronic-structure stack.
Dig deeper in lmmol
Methods siblings in this collection:
- Physics-informed simulation — the neighbouring question of embedding physical law into learned models, where MLIPs are the atomistic-scale instance of the same bargain.
- Simulation optimization — what to do once each simulation is cheap enough to run many times, which is the regime MLIPs create.
For molecular entities named here, use lmmol's graph index, all proteins, or all diseases rather than guessing an entity URL.
Implementation & visualization hooks
ml: a reference implementation of the two-stage recipe that defines current practice — an equivariant backbone pretrained on a broad dataset, then fine-tuned on a small system-specific set — with an active-learning loop and a committee-disagreement uncertainty estimate driving which configurations get labelled [29] [30] [54] [60] [51].
lmvideo / diagramkit: two diagrams. First, the architecture lineage as a graph — descriptors → kernels → message passing → equivariance → foundation models — with each node annotated by the constraint it removed [3] [10] [24] [29] [48]. Second, the active-learning loop as a cycle: run MD, flag high-uncertainty configurations, label with DFT, retrain, showing where the quantum cost is actually incurred [53] [54].