Machine learning for materials

The problem every simulation must solve

Every atomistic simulation rests on one function: the potential energy $E(\mathbf r_1, \dots, \mathbf r_N)$ of a set of atoms as a function of their positions. Its negative gradient gives the forces, the forces give the motion, and the motion, sampled over time, gives every property a simulation can predict, from the melting point of a metal to the binding free energy of a drug. The chapter on potential energy surfaces explains where this function comes from: in the Born–Oppenheimer approximation, it is the ground-state energy of the electrons for fixed nuclei, a solution of the Schrödinger equation.

The difficulty is cost. Solving the electronic problem accurately, even with the approximations of density-functional theory (DFT), takes minutes to hours of computer time for a few hundred atoms, and the cost grows roughly as the cube of the number of electrons. A molecular dynamics simulation needs the forces millions of times. Since the 1920s, physicists and chemists have therefore tried to replace the quantum-mechanical calculation by a cheap formula that reproduces its results: an interatomic potential, also called a force field. The history of the field is a history of increasingly flexible formulas, each built to capture physics that its predecessors missed, culminating in machine-learned potentials that are flexible enough to learn the quantum-mechanical energy surface directly from data. Knowing this history explains why modern machine-learned interatomic potentials (MLIPs) look the way they do, because almost every design decision in them answers a lesson learned by an earlier generation.

The figure below is a timeline of the models discussed in this chapter. Filter by family and click an entry to highlight it.

Pair potentials: Mie, Lennard-Jones and Morse

The earliest potentials described the energy as a sum over pairs of atoms, each pair contributing a function of its separation alone:

Gustav Mie proposed in 1903 that the force between molecules combines a short-range repulsion and a longer-range attraction, each an inverse power of distance. In 1924 John Lennard-Jones, then at Bristol, fitted such potentials to the measured equation of state and viscosity of gases, and the form with exponents 12 and 6 that carries his name became standard:

The attractive term was justified a few years later by Fritz London’s quantum-mechanical theory of dispersion forces, which showed that the correlated fluctuations of two atoms’ electron clouds produce an attraction falling as $r^{-6}$. The repulsive exponent of 12 has no such basis; the true repulsion from overlapping electron clouds is closer to exponential, as in the Buckingham potential, but $r^{-12}$ is the square of $r^{-6}$ and was cheap to compute. For noble gases the Lennard-Jones potential is remarkably good, and it remains the model system of simulation, the “hydrogen atom” of liquid-state theory on which much of statistical mechanics has been tested.

For chemical bonds, Philip Morse introduced in 1929 a potential, $V(r) = D[1 - e^{-a(r - r_0)}]^2$, whose quantum mechanics is exactly solvable and which, unlike a harmonic spring, allows a bond to dissociate. Pair potentials of these kinds, with harmonic terms for bonds and angles and Coulomb terms for partial charges, became the basis of the classical biomolecular force fields developed from the 1970s, such as AMBER, CHARMM, GROMOS and OPLS, which still dominate simulations of proteins and solutions.

Pair potentials, however, fail for most solids, and the failures are instructive. A pair potential predicts that the energy of a crystal varies linearly with the number of neighbours, so that removing an atom costs exactly as much energy as the bonds it had. For metals this is badly wrong: the measured energy to form a vacancy is only about a third of the cohesive energy per atom, whereas pair potentials make the two nearly equal. Pair potentials also force the Cauchy relation between elastic constants, $C_{12} = C_{44}$ for cubic crystals, which real metals violate (gold has $C_{12}$ about four times $C_{44}$). And they cannot explain why carbon, silicon and germanium prefer the open diamond structure with only four neighbours, when a pair potential always favours close packing. The lesson, repeated throughout this history, is that the energy of a bond depends on its environment. That requires many-body terms.

Many-body potentials for metals and semiconductors

For metals, the decisive step was the embedded atom method (EAM) of Murray Daw and Michael Baskes at Sandia, published in 1983–1984 and inspired by density-functional ideas. Each atom is thought of as embedded in the electron density created by its neighbours, and its energy is a nonlinear function of that density:

Because $F$ is concave, roughly $-\sqrt{\bar\rho}$ in the related Finnis–Sinclair potential of 1984 (derived from the second-moment approximation of tight-binding theory), each additional neighbour adds less energy than the previous one. That single feature fixes the vacancy formation energy, breaks the Cauchy relation and gives realistic surface energies. EAM potentials are as cheap as pair potentials, and for forty years they have been the workhorse of simulations of metals and alloys, of dislocations, grain boundaries and radiation damage, with millions of atoms. Notice the form of the equation: the total energy is a sum of atomic energies, each a function of a descriptor of the atom’s local environment (here the scalar $\bar\rho_i$). Twenty-four years later, the Behler–Parrinello network would generalise exactly this structure.

Covalent materials needed angles. Frank Stillinger and Thomas Weber added an explicit three-body term in 1985 that penalises bond angles deviating from the tetrahedral value, making diamond silicon stable. Jerry Tersoff took a different route in 1988 with a bond-order potential, inspired by Linus Pauling’s idea that a bond weakens as an atom forms more of them: the attractive part of each pair interaction is multiplied by a bond order that decreases with the number and angles of competing neighbours. Donald Brenner’s reactive empirical bond-order potential (REBO, 1990) extended this to hydrocarbons and became, with its successor AIREBO, the standard for carbon nanotubes and graphene. Michael Baskes’s modified EAM (MEAM, 1992) added angular terms to the embedding density for directional bonding in metals and semiconductors.

The most ambitious classical potential is ReaxFF, introduced by Adri van Duin, William Goddard and colleagues in 2001. It computes bond orders continuously from distances, so bonds can form and break, combines them with valence-angle, torsion, van der Waals and Coulomb terms, and equilibrates atomic charges at every step. ReaxFF can simulate combustion, catalysis and corrosion, but at a price that foreshadows machine learning: dozens of functional terms and hundreds of parameters per chemical system, fitted by hand to quantum-mechanical data, with limited accuracy and transferability.

Each of these potentials encodes physical insight about a class of bonding in a functional form, and each is accurate only for the situations its form anticipates. Fitting a potential for a new material, or making an existing one more accurate, meant inventing new functional terms, a task that took experts years and often ended in compromises between properties.

Ab initio molecular dynamics

A different solution to the problem was to stop approximating the energy at all. In 1985 Roberto Car and Michele Parrinello showed how to run molecular dynamics with forces computed from density-functional theory on the fly, by treating the electronic wavefunctions as extra dynamical variables that follow the nuclei adiabatically. Ab initio molecular dynamics (AIMD), in the Car–Parrinello form or the Born–Oppenheimer form in which the electronic problem is solved to convergence at every step, made it possible to simulate liquid water, chemical reactions and materials under extreme conditions without any empirical fitting, with the accuracy and transferability of DFT.

Its cost confines it to a few hundred atoms and tens of picoseconds even on large supercomputers. Many important phenomena, such as nucleation, diffusion in solids, the dynamics of interfaces and anything involving thousands of atoms or nanoseconds, are out of reach. By the 2000s the field faced a clear gap: classical potentials were fast but inaccurate and hard to build, ab initio dynamics was accurate but slow. Machine learning was proposed to close it, by learning a potential that reproduces DFT energies and forces at a cost close to that of a classical potential.

The first neural-network potentials

The idea of fitting a potential energy surface with a neural network appeared in the mid-1990s. In 1995 Thomas Blank, Steven Brown, August Calhoun and Douglas Doren fitted a network to the energy of CO on a nickel surface, and over the following decade several groups used networks for low-dimensional surfaces of small molecules and of molecules interacting with surfaces. These early networks took a fixed list of coordinates as input. That made them useful for a single small system, but they could not describe a different number of atoms, they had to learn rotational and permutational symmetry from data, and they could not be applied to condensed phases.

Behler–Parrinello: the atomic decomposition

The breakthrough came in 2007, when Jörg Behler and Michele Parrinello published a neural-network potential for bulk silicon that could describe systems of any size. Their construction rests on three ideas that are still at the heart of almost every MLIP:

The atomic decomposition makes the energy extensive and the model applicable to any number of atoms. Locality makes the cost grow linearly with system size. The symmetry functions, radial and angular functions summed over neighbours and refined in Behler’s 2011 paper, guarantee the physical symmetries exactly, so the network does not have to learn them. Forces follow by differentiating the total energy, so they are consistent with it and energy is conserved in dynamics. The structure is that of EAM, with the single scalar density replaced by a vector of dozens of environment descriptors and the hand-chosen embedding function replaced by a flexible neural network. The chapter on neural-network potentials develops it in detail, and the figure there shows the architecture in action.

Kernel methods and systematic descriptors

In 2010 Albert Bartók, Mike Payne, Risi Kondor and Gábor Csányi at Cambridge introduced Gaussian approximation potentials (GAP), which keep the atomic decomposition but replace the neural network by Gaussian-process regression (the chapter on Gaussian processes). The atomic energy is a weighted sum of kernel functions that measure the similarity between the atom’s environment and environments in the training set. The approach brings a principled uncertainty estimate, and with the smooth overlap of atomic positions (SOAP) descriptor of 2013 it gave a mathematically clean way of comparing environments: smear the neighbours into a density, expand it in radial functions and spherical harmonics, and form rotation-invariant combinations of the coefficients. GAP potentials for silicon, carbon, tungsten and many other materials, some built from thousands of carefully chosen structures, established that machine-learned potentials could match DFT across a wide range of configurations, including defects, surfaces and liquids.

Linear models followed: the spectral neighbour analysis potential (SNAP, Aidan Thompson and colleagues at Sandia, 2015), linear in bispectrum components and fast enough for very large LAMMPS simulations, and Alexander Shapeev’s moment tensor potentials (MTP, 2016), systematically improvable and paired with active learning. In 2019 Ralf Drautz unified these approaches with the atomic cluster expansion (ACE), a complete basis for atomic environments built from products of one-particle functions, of which symmetry functions, SOAP, SNAP and MTP are special or truncated cases. The development is told in the chapters on descriptors and kernel and linear potentials.

Learned representations: message passing

Hand-designed descriptors have a fixed resolution and a fixed cutoff. An alternative, inspired by deep learning’s success at learning features from raw images, is to let the network learn the representation from the atomic positions themselves. Message-passing neural networks, formalised for chemistry by Justin Gilmer and colleagues at Google in 2017, treat a structure as a graph whose nodes are atoms and whose edges connect atoms within a cutoff. Each atom carries a feature vector that is updated repeatedly by aggregating “messages” from its neighbours, which depend on the neighbours’ features and the interatomic distances. Kristof Schütt and colleagues’ SchNet (2017) made the messages continuous functions of distance through learned filters and trained the whole model end to end on energies and forces. Oliver Unke and Markus Meuwly’s PhysNet (2019) added explicit electrostatics, and Johannes Gasteiger, Janek Groß and Stephan Günnemann’s DimeNet (2020) passed messages along bonds using angles.

Message passing changes two things. The representation is learned for the task rather than fixed in advance. And information travels one cutoff per layer, so after $T$ layers an atom’s energy depends on atoms up to $T$ cutoffs away. That enlarges the effective range, but it complicates parallel simulation of large systems, a trade-off discussed in the chapter on message-passing potentials.

Equivariance and the data-efficiency revolution

Invariant models, whose internal features are all rotation-invariant scalars, throw away directional information early and must recover it through many layers or many-body descriptors. Equivariant networks instead carry vectors and higher-order tensors through the network, transforming them correctly under rotations, and combine them with the tensor-product rules of angular momentum. Tensor field networks (Nathaniel Thomas and colleagues, 2018) and Cormorant (2019) introduced the idea; PaiNN (2021) added vector features efficiently; and NequIP, from Simon Batzner, Boris Kozinsky and colleagues at Harvard (published 2022), showed that an E(3)-equivariant message-passing potential could reach state-of-the-art accuracy with strikingly little data, sometimes a few hundred structures where invariant networks needed tens of thousands. Allegro (2023) made the equivariant approach strictly local, so it scales to hundreds of millions of atoms on many GPUs. MACE, by Ilyes Batatia, Dávid Péter Kovács, Gregor Simm, Christoph Ortner and Gábor Csányi (2022), combined equivariant message passing with ACE-style many-body messages, so that two layers suffice, giving a fast and accurate architecture that underlies several of today’s foundation models.

Universal potentials

Until about 2021, every MLIP was trained for one chemical system: silicon, water, a particular alloy. The large databases of DFT calculations assembled by high-throughput projects, above all the Materials Project, made a new goal conceivable: a single potential covering most of the periodic table. Chi Chen and Shyue Ping Ong’s M3GNet (2022), trained on relaxation trajectories from the Materials Project, was among the first such universal potentials; it was followed by CHGNet (2023), by the graph networks behind Google DeepMind’s GNoME project, which proposed hundreds of thousands of new stable crystals (2023), by MACE-MP-0 (2023), MatterSim, SevenNet, ORB and the EquiformerV2 models trained on the 110-million-calculation OMat24 data set (2024), among others. These models can relax structures, screen materials and run stable molecular dynamics for systems they were never specifically trained on, and they can be fine-tuned with a small amount of system-specific data for higher accuracy. Benchmarks such as Matbench Discovery compare them on the task of predicting crystal stability. Their strengths, limits and use are the subject of the chapter on universal potentials.

What the history teaches

Looking back, each generation added flexibility while keeping the lessons of the previous one. From pair potentials came the need for environment dependence; from EAM, the atomic energy as a function of a local descriptor; from bond-order potentials and Stillinger–Weber, the importance of angles; from ReaxFF, the cost of building flexibility by hand; from AIMD, the target accuracy; from Behler and Parrinello, the combination of locality, symmetry and learning that every modern MLIP inherits; from GAP and ACE, systematic and complete descriptions of environments and the value of uncertainty; from message passing and equivariance, learned representations and data efficiency; and from universal models, the power of very large data sets. The principles that recur throughout, and that the rest of this subject develops, are three: physical symmetries must be built in exactly, locality makes potentials scalable, and forces must be the exact gradient of an energy.

References and further reading

The quantity being learned

A machine-learned interatomic potential is a regression model. Its input is a set of atoms, each with a chemical element $Z_i$ and a position $\mathbf r_i$, plus the periodic cell for a crystal. Its output is the potential energy $E$, and, by differentiation, the forces on every atom and the stress on the cell. It is trained on examples whose energies and forces were computed by an electronic-structure method, and its purpose is to reproduce that method at a cost lower by a factor of a thousand to a million.

The reference method fixes the best accuracy the model can ever reach. For materials it is almost always density-functional theory (DFT), founded on the Hohenberg–Kohn theorems (1964) and the Kohn–Sham equations (1965), with an approximate exchange–correlation functional. The generalised-gradient functional PBE (Perdew, Burke and Ernzerhof, 1996) is the most common choice in materials databases; it describes metals and many oxides reasonably, but overbinds or underbinds certain bonds, misses dispersion unless a correction such as D3 is added, and underestimates band gaps and reaction barriers. Meta-GGA functionals such as r²SCAN and hybrid functionals are more accurate and more expensive; for small molecules, coupled-cluster theory, CCSD(T), serves as the “gold standard”. An MLIP trained on PBE reproduces PBE, including its errors: a perfect PBE potential for water predicts water that is overstructured and freezes well above 0 °C, because PBE water does. Choosing the reference is therefore the first modelling decision, and it should be made with the property of interest in mind.

The reference energies must also be consistent. Every structure in a training set must be computed with the same functional, pseudopotentials, plane-wave cutoff, smearing and convergence criteria, otherwise the model is asked to learn noise. Two large databases computed with different settings cannot simply be merged, a practical problem that the builders of universal potentials face constantly (the chapter on universal potentials).

Energies, forces and stresses

An electronic-structure calculation on a structure of $N$ atoms yields one total energy, but also $3N$ force components, by the Hellmann–Feynman theorem, at almost no extra cost, and for a periodic cell the six independent components of the stress tensor. Forces are by far the richest source of information: a single 100-atom structure provides 300 force labels and one energy. They are local, each describing the slope of the energy surface at one atom, which matches the local structure of MLIPs. Training on forces is what makes modern potentials data-efficient.

Almost all MLIPs obtain forces by differentiating their energy with automatic differentiation, exactly as the tutorial in the chapter on building an MLIP does. Some recent models instead predict forces directly as an extra output, which is faster, but the predicted forces are then not the gradient of any energy: molecular dynamics with them does not conserve energy and can heat up or drift, and geometry optimisation can fail to converge. Studies of such non-conservative models have found that they can be useful for fast pre-screening but must be used with care in dynamics. The safe default remains: learn an energy, differentiate it.

The training loss usually combines all three kinds of label, each weighted to balance their very different magnitudes and numbers:

Locality and nearsightedness

Every practical MLIP makes the same structural assumption as the Behler–Parrinello network and the embedded-atom method before it: the total energy is a sum of contributions from each atom, and each contribution depends only on the neighbours within a cutoff radius $r_c$, typically 5 to 6 Å.

The justification is a property of quantum mechanics that Walter Kohn called the nearsightedness of electronic matter (1996): for fixed chemical potential, the local electronic properties at one point depend only weakly on perturbations far away, with an influence that decays exponentially in insulators and, more slowly, in metals at finite temperature. A change of position of an atom 10 Å away hardly changes the forces here. The cutoff is a compromise between capturing as much of this influence as possible and keeping the number of neighbours, and hence the cost and the amount of data needed, manageable: a 6 Å sphere in a dense solid contains around 80 atoms.

Locality fails for some physics, and knowing where is essential. Electrostatics between charged or polar species decays as $1/r$ and is not local; neither is the response of charge to a distant perturbation, as in charge transfer between separated parts of a battery interface. Dispersion decays as $r^{-6}$, slowly enough that its sum beyond 6 Å contributes noticeably to cohesive and surface energies. And message-passing networks, whose effective range extends over several cutoffs, blur the strict notion of locality. These limits motivate MLIPs with explicit long-range terms, discussed in the chapter on neural-network potentials.

Symmetry

The energy of an isolated set of atoms does not change if the whole set is translated, rotated or reflected, or if two atoms of the same element exchange labels. These are exact symmetries of the Schrödinger equation, and a potential that violates them makes unphysical predictions: a molecule that spins up spontaneously, a crystal whose energy depends on how it is oriented relative to the coordinate axes. Together they form the Euclidean group E(3) of translations, rotations and reflections, plus the permutation group.

Translation invariance is the easy part: models use only relative positions $\mathbf r_j - \mathbf r_i$. Permutation invariance follows from summing over neighbours with a shared function. Rotation is the hard part, and there are three ways to handle it. Invariant descriptors, such as distances, angles, symmetry functions and SOAP, contain only rotation-invariant quantities, so anything built on them is invariant. Equivariant architectures carry vectors and higher tensors internally and combine them only in ways that respect rotations, producing invariant energies at the end. Data augmentation feeds the model many randomly rotated copies of each structure and hopes it learns the symmetry approximately. The first two give exact symmetry; the third gives approximate symmetry and is used by a few fast models.

The figure below shows what happens without symmetry. On the left, forces are computed from an energy that depends only on distances, so when the molecule is rotated, every force rotates with it. On the right, a model that takes raw Cartesian coordinates as input, fitted (to first order) on the molecule in its original orientation, is given the rotated molecule. Its forces do not follow, and its energy changes, although nothing physical has changed.

A model that has seen only one orientation has no way to know that a rotated copy is the same structure. Learning that fact from data would require training examples in every orientation, for every environment: a waste of data and, inevitably, an approximation. Building the symmetry in makes the model exact in this respect and lets all its capacity go into learning the chemistry.

Smoothness and the cutoff

A potential used for molecular dynamics must be a smooth function of the positions. If the energy jumps, energy is not conserved when an atom crosses the discontinuity; if the force jumps, the integrator’s error grows. Both happen at the cutoff if neighbours simply stop counting when they cross $r_c$. MLIPs therefore multiply every neighbour’s contribution by a cutoff function that goes smoothly to zero at $r_c$ together with at least its first derivative, such as Behler’s $f_c(r) = \tfrac12[\cos(\pi r/r_c) + 1]$ for $r < r_c$, or a polynomial envelope. Descriptors and network activations must also be differentiable: a network built from ReLU units, whose derivative jumps at zero, gives forces with kinks, which is why smooth activations such as tanh, SiLU (also called swish) or softplus are used.

Body order and completeness

The classical way to organise an energy function is the many-body expansion, a sum of one-body (isolated atom) energies, two-body corrections for every pair, three-body corrections for every triple, and so on:

A pair potential stops at $V_2$; Stillinger–Weber adds $V_3$. Each term is a function of the relative geometry of its atoms. The body order of a model is the highest number of atoms that appear together in one term. Many MLIP descriptors can be classified this way: radial symmetry functions are two-body, angular ones three-body, the SOAP power spectrum is three-body, and the bispectrum used by SNAP is four-body. The atomic cluster expansion makes body order an explicit, adjustable parameter, and MACE builds high body orders efficiently inside each message.

Body order matters because low-order descriptions are incomplete. Two different arrangements of neighbours can have identical lists of distances, and even identical lists of distances and angles; a model whose input is limited to these cannot tell them apart and must give them the same energy. Pozdnyakov and colleagues showed in 2020 that such degenerate pairs exist for three-body descriptors such as SOAP and that they cause real errors. Higher body orders, or equivariant message passing, remove these degeneracies.

The figure below gives a two-dimensional picture of how invariants are built. The angular arrangement of neighbours around a central atom is described by a density $n(\theta)$, which is expanded in circular harmonics $e^{im\theta}$ (the two-dimensional analogue of spherical harmonics). Rotating the environment by an angle $\alpha$ multiplies each coefficient $c_m$ by $e^{-im\alpha}$: its phase turns $m$ times as fast as the environment, but its magnitude $|c_m|$ does not change. The magnitudes form the power spectrum, a rotation-invariant fingerprint of the environment, and the analogous construction in three dimensions with spherical harmonics is the SOAP descriptor. Invariants of higher order, such as products $c_{m_1}c_{m_2}c^*_{m_1+m_2}$ (the bispectrum), retain information about the relative phases that the power spectrum discards.

A general recipe

With these requirements in place, every MLIP in the following chapters can be read as a variation on one recipe. Describe each atom’s neighbourhood within a cutoff in a way that is smooth and respects symmetry, either with a fixed descriptor or with learned, equivariant features. Map that description to an atomic energy with a flexible regression model: a linear model, a kernel, a neural network. Sum the atomic energies. Differentiate the sum to obtain forces and stresses. Fit the parameters to reference energies and forces with a loss like the one above, and test on structures the model has not seen. The families differ in how they make each choice, and the next chapters cover them in turn: descriptors, kernel and linear models, neural-network potentials, and message-passing and equivariant networks.

References and further reading

The input problem

A neural network takes a vector of numbers. A configuration of atoms is a list of positions $\mathbf r_i$ and element types $Z_i$. The obvious encoding, concatenating the Cartesian coordinates $(x_1, y_1, z_1, x_2, \dots)$, is a bad one, for a reason that has nothing to do with machine learning: the energy of a molecule does not change when you move it, turn it round or relabel its atoms, but the coordinate vector changes completely. A model fed raw coordinates would have to learn, from data, that millions of different input vectors describe the same physical object. With finite data it never quite does.

The symmetries of physics

For an isolated system of atoms, the potential energy $E(\mathbf r_1, \dots, \mathbf r_N)$ is unchanged by:

• translation, $\mathbf r_i \to \mathbf r_i + \mathbf t$;
• rotation and reflection, $\mathbf r_i \to \mathbf R\mathbf r_i$ with $\mathbf R$ orthogonal (the group O(3); together with translations, E(3));
• permutation of identical atoms, swapping the labels of two atoms of the same element.

A function with $f(g\cdot x) = f(x)$ for every operation $g$ of a group is invariant. The forces are not invariant: rotate the molecule and the force vectors rotate with it. A function with $f(g\cdot x) = \rho(g)\,f(x)$, where $\rho(g)$ is the corresponding transformation of the output, is equivariant. Forces, dipole moments and stress tensors are equivariant. The good news is that if the energy model is invariant, forces computed as $\mathbf F_i = -\nabla_{\mathbf r_i}E$ are automatically equivariant.

There are three ways to give a model these symmetries:

1. Data augmentation: show it many rotated copies and hope it learns. Approximate, and wasteful.
2. Invariant descriptors: transform the structure into features that are invariant by construction, then use any regressor (this lesson).
3. Equivariant architectures: build networks whose internal features transform correctly under rotation at every layer (the lesson on graph neural networks and the lesson on interatomic potentials).

Locality and extensivity

A second physical principle shapes almost every model. Electronic structure is nearsighted (Kohn, 1996): in insulators and in metals at finite temperature, the local electronic properties depend only weakly on atoms far away. This justifies writing the total energy as a sum of atomic contributions, each depending only on the atom’s neighbourhood within a cut-off radius $r_c$ (typically 4–6 Å):

The decomposition makes the energy extensive (twice the system, twice the energy), makes the cost linear in $N$, and lets a model trained on small cells be applied to systems of millions of atoms. It is also an approximation: long-range electrostatics and dispersion need separate treatment, which modern potentials add explicitly.

Descriptors, from simple to complete

Pair distances. All interatomic distances are invariant to rotation and translation; sorting them makes them invariant to permutation. But a sorted list does not say which distances belong to which atom, and its length changes with the number of atoms. It is a useful teaching example rather than a practical descriptor.

Coulomb matrix. Rupp et al. (2012) described a molecule by the matrix $M_{ij} = Z_iZ_j/\lvert\mathbf r_i - \mathbf r_j\rvert$ (with $M_{ii} = 0.5Z_i^{2.4}$), using its sorted eigenvalues as an invariant fingerprint. Simple and historically important, but global, so it does not decompose into atomic contributions.

Atom-centred symmetry functions. Behler and Parrinello (2007) described each atom’s environment by a fixed set of invariant functions. The radial functions

count neighbours in soft shells centred at $r_s$; angular functions $G_i^{4}$ add three-body information through $\cos\theta_{ijk}$. The cut-off function $f_c$ goes smoothly to zero so that the energy, and its derivative the force, stay continuous as atoms enter and leave the sphere. Each atom gets a vector of perhaps 50–100 such numbers, fed to an element-specific neural network that outputs $\varepsilon_i$. This was the first machine-learned potential that scaled to condensed phases. The third row of the lab is a two-dimensional version of $G^2$, summed over atoms.

SOAP. The smooth overlap of atomic positions (Bartók, Kondor and Csányi, 2013) places a Gaussian on every neighbour to form a smooth neighbour density $\rho_i(\mathbf r)$, expands it in radial functions and spherical harmonics, $\rho_i(\mathbf r) = \sum_{nlm} c^{i}_{nlm} g_n(r) Y_{lm}(\hat{\mathbf r})$, and forms the rotationally invariant power spectrum $p^i_{nn'l} = \sum_m c^{i}_{nlm}\,(c^{i}_{n'lm})^*$. The SOAP kernel, a normalised dot product of power spectra, underlies the GAP potentials and much of kernel-based materials ML.

Completeness and ACE. A descriptor built from distances and angles (two- and three-body information) is not complete: Pozdnyakov et al. (2020) constructed distinct environments that such descriptors cannot tell apart. The atomic cluster expansion (ACE; Drautz, 2019) systematically includes correlations of any body order, from a basis that can be made complete, and connects symmetry functions, SOAP and modern equivariant networks in a single framework.

In practice

The radial functions in NumPy, with a numerical check of all three symmetries:

import numpy as np

def cutoff(r, rc):
    """Smooth cutoff: 1 at r = 0, 0 (with zero slope) at r = rc."""
    return np.where(r < rc, 0.5 * (np.cos(np.pi * r / rc) + 1), 0.0)

def g2(pos, etas=(4.0,), centres=(1.0, 1.5, 2.0, 2.5, 3.0), rc=4.0):
    """Behler-Parrinello radial symmetry functions, one row per atom."""
    d = np.linalg.norm(pos[:, None] - pos[None], axis=-1)
    np.fill_diagonal(d, 2 * rc)                      # an atom is not its own neighbour
    fc = cutoff(d, rc)
    return np.array([[np.sum(np.exp(-eta * (d[i] - rs) ** 2) * fc[i])
                      for eta in etas for rs in centres] for i in range(len(pos))])

rng = np.random.default_rng(0)
pos = rng.normal(scale=1.2, size=(8, 3))             # a random 8-atom cluster
G = g2(pos)

# a random rotation (QR of a Gaussian matrix), a translation and a permutation
Q, _ = np.linalg.qr(rng.normal(size=(3, 3)))
perm = rng.permutation(8)
pos2 = (pos @ Q.T + np.array([3.0, -1.0, 2.0]))[perm]
G2 = g2(pos2)

print("raw coordinates changed by up to", np.abs(pos2 - pos).max().round(3))
print("per-atom descriptors equal after undoing the permutation:", np.allclose(G2, G[perm]))
print("sum over atoms unchanged:", np.allclose(G2.sum(0), G.sum(0)))

raw coordinates changed by up to 5.164
per-atom descriptors equal after undoing the permutation: True
sum over atoms unchanged: True

The per-atom descriptors are equivariant to permutation (renumbering the atoms renumbers the rows) and their sum is invariant: exactly the structure of $E = \sum_i\varepsilon_i$. For production work, libraries such as DScribe and librascal compute ACSF, SOAP and related descriptors for molecules and periodic crystals, and ACEpotentials.jl and pacemaker fit ACE potentials.

References and further reading

What is special about crystals

A molecule is a finite set of atoms. A crystal is infinite: a lattice of translations $\mathbf T = n_1\mathbf a_1 + n_2\mathbf a_2 + n_3\mathbf a_3$ (integers $n_k$, lattice vectors $\mathbf a_k$) repeats a basis of atoms at fractional coordinates $\mathbf f_j$, so that every atom sits at $\mathbf r = \sum_k (f_{jk} + n_k)\mathbf a_k$. A crystal structure file (CIF, POSCAR) stores just one unit cell: the three lattice vectors, the species and fractional coordinates of the basis, and, implicitly, the rule that it repeats forever.

Two consequences matter for machine learning.

• The unit cell is not unique. Any supercell, or any other choice of lattice vectors spanning the same lattice, describes the same crystal. A representation must give the same answer for all of them: an extra symmetry beyond those of the lesson on symmetry and descriptors.
• Neighbours cross the cell boundary. An atom near one face of the cell has neighbours that are periodic images of atoms near the opposite face. A model that only looks inside the cell sees the wrong environment.

Composition-only models

Sometimes the structure is not known, for example when screening hypothetical compositions before any structure has been predicted. Composition-based models describe a material by its formula alone: hand-crafted statistics of elemental properties weighted by stoichiometry (the Magpie feature set: mean, range and deviation of electronegativity, atomic radius, valence electrons and so on), or learned element embeddings combined by a network over the elements present (ElemNet, Roost). They are fast and useful for screening, but cannot distinguish polymorphs: diamond and graphite are both “C”.

Crystal graphs

Xie and Grossman’s crystal graph convolutional neural network (CGCNN, 2018) introduced the representation most structure-based models still use. Each atom in the unit cell is a node, and edges connect it to every atom, including periodic images, within a cut-off radius $r_c$ (or to its 12 nearest neighbours). Because an atom can be bonded to several images of the same atom, the result is a multigraph with several edges between the same pair of nodes, each with its own length. Edge features are the distances expanded in Gaussians, $e_{ij,k} = \exp\!\big[-(d_{ij} - \mu_k)^2/\gamma^2\big]$; node features are learned embeddings of the element.

Things to notice in the lab:

• In the triangular lattice there is only one atom per cell, so every edge goes to a periodic image of itself: the graph is a single node with six self-loops of different image shifts. Without periodic images there would be no graph at all.
• Raising $r_c$ past each coordination shell adds a shell of edges at once; the histogram is the discrete version of the radial distribution function from the lesson on trajectory analysis.
• Strain splits shells that were degenerate: under uniaxial strain the six nearest neighbours of the triangular lattice become two plus four, and a cut-off placed between the new distances changes the coordination number abruptly. Smooth cut-off functions and distance expansions are what keep the model’s output continuous under such changes.
• Supercells give the same per-atom environments, so the per-atom features and an extensive sum over atoms do not depend on the choice of cell.

Later models refine the idea. MEGNet adds a global state vector (temperature, pressure); ALIGNN adds a second “line graph” whose nodes are bonds, so that messages carry bond angles; M3GNet and CHGNet add three-body interactions and are trained on energies, forces and stresses from relaxation trajectories, so they serve as universal interatomic potentials (the lesson on interatomic potentials).

Where the data come from

Structure-based models are trained on high-throughput DFT databases:

Two cautions apply to all of them. First, they are DFT numbers, usually with the PBE functional: band gaps are systematically underestimated, and formation energies carry functional-dependent errors and empirical corrections. A model trained on them learns PBE, not experiment. Second, databases are biased towards known chemistry and towards structures that someone chose to compute, which is exactly the regime where a discovery model should be least trusted (the lesson on generalisation).

Stability: the convex hull

The most important property for materials discovery is whether a compound can exist. A negative formation energy is not enough: the compound must also be stable against decomposing into other compounds of the same elements. Plot the formation energy per atom of every known phase in a chemical system against composition and take the lower convex hull. Phases on the hull are thermodynamically stable (at 0 K, in DFT); a phase above it decomposes into the hull phases at its composition, and its energy above hull $E_{\text{hull}}$ measures by how much. Many metastable materials with $E_{\text{hull}}$ of a few tens of meV/atom can be synthesised, so $E_{\text{hull}}$ is a likelihood, not a verdict. Predicting $E_{\text{hull}}$ is harder than predicting formation energy, because small errors in the energies of competing phases shift the hull. Large-scale campaigns such as GNoME (2023) combined graph networks with DFT verification in active-learning loops to propose hundreds of thousands of new hull-stable crystals; how many of these can be made remains an open experimental question.

In practice

ASE’s neighbour list returns exactly the edges of a crystal graph: source atom, target atom, distance and the integer image shift:

from ase.build import bulk
from ase.neighborlist import neighbor_list
import numpy as np

atoms = bulk("NaCl", "rocksalt", a=5.64)        # primitive cell: 1 Na + 1 Cl
i, j, d, S = neighbor_list("ijdS", atoms, cutoff=3.0)

print("atoms in cell:", len(atoms), " directed edges:", len(i))
for a in range(len(atoms)):
    nb = j[i == a]
    print(f"{atoms[a].symbol}: {len(nb)} neighbours at {np.unique(d[i == a].round(3))} Å, "
          f"all of them {set(atoms[nb].get_chemical_symbols())}")
print("distinct periodic-image shifts used:", len({tuple(s) for s in S}))

atoms in cell: 2  directed edges: 12
Na: 6 neighbours at [2.82] Å, all of them {'Cl'}
Cl: 6 neighbours at [2.82] Å, all of them {'Na'}
distinct periodic-image shifts used: 11

The two-atom primitive cell of rock salt gives each ion its six octahedral neighbours of the other species, and almost every edge reaches a different periodic image. pymatgen (Structure.get_neighbor_list) and matminer (composition and structure featurisers) are the standard companions.

References and further reading

Regression on fixed descriptors

The models in this chapter share a structure. Each atom’s environment is first turned into a fixed vector of numbers, a descriptor $\mathbf q_i$, that respects the symmetries described in the chapter on what an MLIP must learn. The atomic energy is then a function of that descriptor, chosen from a family whose parameters can be fitted reliably: either a linear combination of basis functions, or a kernel expansion that compares the environment with stored reference environments. Both families have a decisive practical property: the prediction is linear in the parameters being fitted, so training is a linear least-squares problem with a unique solution, no random initialisation and no learning-rate schedule. Both also come with principled uncertainty estimates. Their limitation is that the descriptor is fixed in advance, so the model can be no better than the description of the environment it is given.

Linear models

In a linear MLIP, each atomic energy is a weighted sum of basis functions of the environment, $B_k(i)$, and the total energy is the sum over atoms:

Stacking every energy and force label of the training set into a vector $\mathbf y$ and the corresponding features into a design matrix $\mathbf X$, the coefficients minimise a regularised squared error, exactly as in the chapter on learning from data:

with a diagonal weight matrix $\mathbf W$ that balances energies, forces and stresses. The figure below shows the idea in its simplest form: the energy curve of a single bond fitted as a linear combination of smooth radial basis functions, each multiplied by a cutoff function. The grey curves are the basis functions scaled by their fitted coefficients; their sum is the orange fit.

With too few basis functions the fit cannot follow the curve; with many functions and weak regularisation it can oscillate between data points; with adequate data and regularisation it is excellent inside the training range. Outside that range, at short bond lengths where no data were given, every setting fails, in different ways. That is the central fact about all machine-learned potentials: they interpolate well and extrapolate unreliably. Real potentials protect against the most dangerous case, atoms pushed too close together, by adding a physically motivated repulsive term (such as the ZBL screened-Coulomb potential used in radiation-damage simulations) or by making sure the training set contains close contacts.

Gaussian approximation potentials

Gaussian approximation potentials (GAP), introduced by Albert Bartók, Mike Payne, Risi Kondor and Gábor Csányi in 2010, replace the fixed basis by a kernel: a function $k(\mathbf q, \mathbf q')$ that measures how similar two atomic environments are. The atomic energy is a weighted sum of kernels between the environment and a set of $M$ representative environments $\mathbf q_m$ chosen from the training data:

The connection to Gaussian processes (the chapter on uncertainty) is that the kernel defines a prior over functions: it says how strongly the energies of two environments are correlated before any data are seen. Conditioning on the training data then gives a posterior mean, the prediction above, and a posterior variance, an estimate of the model’s own uncertainty. The training labels are total energies, sums of atomic energies, and forces, derivatives of those sums, but because both are linear in the $\alpha_m$ they can be used directly in the regression.

A full Gaussian process with $N$ atomic environments needs to invert an $N\times N$ matrix, which becomes impossible for the hundreds of thousands of environments in a typical training set. GAP uses sparse Gaussian-process regression, in the projected-process form described by Snelson and Ghahramani: the function is represented through $M \ll N$ representative environments, chosen by clustering or by CUR matrix decomposition, and the weights solve

where $\mathbf K_{MN}$ contains kernels (and their derivatives, for force labels) between representative environments and training data, and the diagonal matrix $\boldsymbol\Lambda$ holds the expected error of each label, which plays the role of the regularisation. Setting these noise levels per type of structure, smaller for crystals than for liquids for example, is how GAP users tell the fit which data to trust most.

Switch the figure to kernel mode to see a one-dimensional Gaussian process. Between training points the predictive band is narrow; beyond them it widens, honestly reporting that the model does not know. This uncertainty is the main attraction of kernel potentials for building training sets, because new structures can be chosen where the predicted variance is largest (active learning).

GAP models have been built for many materials, often with combined descriptors: a two-body term, a three-body term and SOAP, each with its own kernel. The GAP-18 potential for silicon, trained on a carefully designed database of crystal phases, surfaces, defects, liquid and amorphous structures, reproduces DFT across this whole range and was used to study the structure of amorphous silicon. GAP-17 for carbon describes amorphous carbon, graphite, diamond and their transformations, and GAP models of tungsten were developed for the dislocation and defect physics of fusion-reactor materials. The cost of prediction scales with $M$ and with the cost of computing SOAP, which makes GAP accurate but relatively slow, typically one to two orders of magnitude more expensive than an EAM potential.

The spectral neighbour analysis potential

The spectral neighbour analysis potential (SNAP), developed by Aidan Thompson and colleagues at Sandia in 2015, takes the linear route with a powerful descriptor. The density of neighbours around each atom is mapped onto the surface of a four-dimensional sphere and expanded in hyperspherical harmonics; the bispectrum components, triple products of the expansion coefficients contracted with Clebsch–Gordan coefficients, are invariant to rotation and capture four-body correlations. The atomic energy is a linear combination of bispectrum components (or, in quadratic SNAP, of their products). Because the model is linear, fitting is fast and evaluation is cheap: SNAP was designed for and runs efficiently in the LAMMPS code, on GPUs and on the largest supercomputers, and has been used for tungsten, tantalum, indium phosphide and carbon at extreme conditions, including billion-atom simulations.

Moment tensor potentials

Alexander Shapeev’s moment tensor potentials (MTP, 2016) build the basis from moment tensors of the neighbourhood, tensors of increasing rank formed from the neighbour vectors weighted by radial functions:

Increasing the level systematically enlarges the basis, so accuracy can be improved in a controlled way, and the MTP family has a mathematical completeness result: any smooth, symmetric function of the environment can be approximated to arbitrary accuracy as the level grows. MTPs are usually trained with an active-learning scheme based on the D-optimality criterion: during a simulation, the “extrapolation grade” of each new configuration measures how far outside the span of the training set it lies, and configurations beyond a threshold are sent for DFT calculation and added to the training set. This scheme, implemented in the MLIP package by Podryabinkin and Shapeev (2017), made MTPs a popular choice for alloys and for crystal-structure prediction.

The atomic cluster expansion

In 2019 Ralf Drautz introduced the atomic cluster expansion (ACE), which put all these descriptors into a single framework. ACE starts from a one-particle basis, products of radial functions and spherical harmonics, and projects the neighbour density onto it:

The key computational idea is the density trick. Writing an explicit sum over all triples or quadruples of neighbours would cost a power of the number of neighbours. ACE instead sums once over neighbours to form the atomic base, then multiplies base functions, so the cost of a high body-order basis function is the cost of a product, not of a sum over clusters. Exactly the same trick appears in the angular-expansion figure of the chapter on what an MLIP must learn: the power spectrum $|c_m|^2$ is a two-neighbour (three-body) invariant computed from single sums.

ACE is complete: as the radial and angular resolution and the body order grow, it can represent any smooth symmetric function of the environment. Symmetry functions, SOAP (body order 3), SNAP’s bispectrum (body order 4) and MTP’s contractions all correspond to particular choices or truncations within it. Linear ACE potentials are fitted like SNAP; nonlinear ACE, implemented in the pacemaker code, passes several linear ACE sums through an embedding function reminiscent of EAM, which improves accuracy for metals. Finally, ACE provided the theory for MACE, which builds ACE-like many-body features inside equivariant message passing (the chapter on message-passing potentials).

Choosing among them

Systematic comparisons, such as the 2020 study of Zuo, Ong and colleagues that trained GAP, MTP, neural-network, SNAP and quadratic SNAP potentials on the same data for six elements, found that the more flexible models (MTP, GAP, neural networks) reached lower errors, while SNAP was cheapest to evaluate and qSNAP sat between them. The broad trade-offs are these. Linear models (SNAP, MTP, linear ACE) are fast to fit and fast to run, have a convex training problem and well-defined extrapolation measures, and scale to enormous simulations; their accuracy is limited by the size of the basis, which grows quickly with the number of elements. Kernel models (GAP) are among the most accurate for a given amount of data and give calibrated uncertainty, at a higher evaluation cost that grows with the number of representative environments. Both are excellent for single- or few-element materials with carefully constructed training sets. For many-element chemistry and very large data sets, the learned representations of neural networks and message-passing models, discussed in the next chapters, have become the more common choice.

References and further reading

The high-dimensional neural network potential

The neural-network potential of Jörg Behler and Michele Parrinello (2007) is the template from which most of the field grew, and it is worth understanding completely. Its three ingredients, introduced in the history and justified in the chapter on what an MLIP must learn, are the atomic decomposition of the energy, a symmetric description of each atom’s neighbourhood, and one feed-forward neural network per chemical element that maps that description to an atomic energy. Behler called it a high-dimensional neural network potential (HDNNP) because, unlike the networks of the 1990s, it describes systems with any number of atoms and therefore any number of dimensions.

The figure follows a configuration through the model. Each atom’s neighbours within the cutoff are summarised in a descriptor vector; the same network turns every descriptor into an atomic energy; the atomic energies add up to the total.

Atom-centred symmetry functions

The descriptors Behler used, systematised in his 2011 paper, are atom-centred symmetry functions (ACSFs). Every function is a sum over neighbours, so it does not depend on their order, and it uses only distances and angles, so it does not depend on rotation or translation. A smooth cutoff function $f_c(r)$, equal to 1 at short range and falling smoothly to zero at $r_c$, multiplies every term.

A typical potential uses 30 to 100 symmetry functions per element. Their parameters are chosen so that the radial functions resolve the neighbour shells of the materials involved and the angular functions distinguish the relevant bond angles, and the inputs are standardised (shifted and scaled to zero mean and unit variance over the training set) before entering the network. Choosing them is partly an art, and a poor choice limits the accuracy of the whole potential; the tutorial in the chapter on building an MLIP shows a set that is too smooth to resolve the neighbour shells of copper failing to fit forces, and a sharper set succeeding. The chapter on descriptors compares ACSFs with SOAP and ACE.

The atomic network and the chain rule for forces

The atomic network is a small multilayer perceptron, typically two or three hidden layers of 20 to 100 neurons with smooth activation functions such as tanh, softplus or SiLU (swish). With 50 inputs and two hidden layers of 30 neurons it has about 2500 weights per element, a tiny network by deep-learning standards, evaluated once per atom per time step.

Forces require the derivative of the total energy with respect to every atomic position. Each atomic energy depends on positions only through descriptors, so the chain rule splits the derivative into a network part and a descriptor part:

The first factor comes from backpropagating through the network, exactly as in the chapter on backpropagation; the second is the derivative of the symmetry functions, which is analytic. Because forces are derived from a single energy function, they are conservative.

Training

The network weights are fitted by minimising a loss combining energy and force errors, as in the chapter on what an MLIP must learn, with Adam or, in Behler’s original work, a Kalman-filter optimiser. Including forces in the loss requires differentiating the force expression with respect to the weights, a second-order derivative that automatic-differentiation frameworks handle by keeping the computational graph of the first derivative (create_graph=True in PyTorch). The data are split into training, validation and test sets; the validation error decides when to stop training and which hyperparameters to keep, and the test error is reported.

The figure below trains a small atomic network in your browser. The system is a two-dimensional cluster of five atoms interacting through a smoothly truncated Lennard-Jones potential; each atom is described by eight radial symmetry functions; the network has two hidden layers. To keep the browser fast, it is trained on energies only. Watch the three panels: the learning curve, the parity plot of predicted against true energies for test structures, and the parity plot of the forces, computed exactly as the negative gradient of the trained network.

Two lessons appear quickly. With 30 training structures the network fits the training energies well but the test error stays high, a clear case of overfitting; with 480 structures training and test errors track each other. And even when the energies are learned well, the forces are noticeably worse, because a network trained only on energies is never told the slope of the surface. Real potentials are therefore trained on forces, which also supply $3N$ labels per structure instead of one.

Many elements

For a system with several elements, each element has its own network, and the descriptors must distinguish neighbours by element: a radial function counting oxygen neighbours of a hydrogen atom is different from one counting hydrogen neighbours. With $S$ elements, radial functions grow with $S$ and angular functions with the number of element pairs, roughly $S^2/2$, so descriptor sets for five or more elements become unwieldy, and every network must be trained with data covering all the element combinations it will meet. Remedies include element-weighted symmetry functions that weight each neighbour by a property of its element instead of creating separate functions per element (Gastegger and colleagues, 2018), and learned element embeddings, the approach taken by message-passing networks and universal potentials.

ANI: transferable potentials for organic molecules

Justin Smith, Olexandr Isayev and Adrian Roitberg’s ANI-1 (2017) applied the Behler–Parrinello architecture, with modified angular symmetry functions, to organic molecules containing H, C, N and O, trained on about 20 million DFT calculations of conformations of small molecules. It was among the first network potentials designed to transfer to molecules outside its training set. ANI-1x (2018) built its data set by active learning, adding only conformations on which an ensemble of networks disagreed, and reached similar accuracy with a fraction of the data. ANI-1ccx (2019) then used transfer learning: a network pre-trained on a large DFT data set was retrained on a much smaller set of coupled-cluster CCSD(T) energies, producing a potential that approaches coupled-cluster accuracy for reaction thermochemistry and torsion profiles at a cost billions of times smaller. ANI-2x (2020) extended the element set to S, F and Cl. The ANI models showed three ideas that are now standard: large chemically diverse data sets, active learning to build them, and transfer from cheap to expensive levels of theory.

Deep Potential

The Deep Potential of Linfeng Zhang, Jiequn Han, Han Wang, Roberto Car and Weinan E (2018) keeps the atomic decomposition but replaces hand-chosen symmetry functions by a learned description. In the widely used DeepPot-SE (“smooth edition”) form, each neighbour’s relative position is encoded in a row of a matrix $\mathcal R_i$ (containing a smooth function of distance $s(r_{ij})$ and the scaled relative vector), and an embedding network maps $s(r_{ij})$ to a learned feature vector for each neighbour, forming a matrix $\mathcal G_i$. The descriptor

is invariant to rotation (the products $\mathcal R_i\mathcal R_i^{\mathsf T}$ contain only dot products of neighbour vectors), to translation and to permutation (sums over neighbours), and its features are learned rather than chosen. A fitting network then maps $\mathcal D_i$ to the atomic energy.

Deep Potential’s influence comes as much from its software as from its architecture. DeePMD-kit made training and running potentials straightforward and integrated them with LAMMPS and GPUs, and DP-GEN automated a concurrent-learning loop: train an ensemble of models, run molecular dynamics with them, select configurations where the models disagree, compute them with DFT, and repeat. With these tools Deep Potential models were built for water and ice, metals and alloys, and battery materials; a Deep Potential water model was used to compute the phase diagram of water over a wide range of temperatures and pressures, and a 2020 simulation of about 100 million atoms with ab initio accuracy on the Summit supercomputer won the Gordon Bell Prize.

Four generations: adding long-range physics

Behler has classified network potentials into generations. The first generation were the low-dimensional networks of the 1990s. The second generation, the HDNNP described above, is local: each atomic energy depends only on the cutoff sphere. This is sufficient for many materials, but it misses electrostatics beyond the cutoff, which matters for ionic materials, interfaces and molecules in solution.

The third generation adds a second set of networks that predict an environment-dependent charge on each atom, also from symmetry functions; the electrostatic energy of these charges is computed exactly, with Ewald summation in periodic systems, and added to the short-range energy. Charges are learned from a reference such as Hirshfeld charges, or adjusted to reproduce DFT dipole moments. The charges are still local, however, and cannot respond to distant changes in the system. The fourth generation (4G-HDNNP, Ko, Finkler, Goedecker and Behler, 2021) determines the charges by a global charge equilibration: networks predict an electronegativity for each atom, and the charges are found by minimising a charge-dependent energy for the whole system subject to total charge. This captures non-local charge transfer, for example between a cluster and a distant dopant, which no local model can describe. Other approaches to long-range physics include explicit dispersion corrections, learned long-range descriptors, and the latent Ewald summation used in some message-passing models.

References and further reading

From fixed descriptors to learned representations

The potentials of the previous two chapters describe each atom with a descriptor chosen before training. Choosing it well requires expertise, its resolution is fixed, and for many elements it becomes very long. Deep learning’s central lesson, learned on images and text, is that representations learned from data for the task at hand usually beat hand-designed ones. Message-passing neural networks bring this lesson to atoms. A structure is treated as a graph whose nodes are atoms and whose edges connect atoms closer than a cutoff; every atom carries a feature vector that starts as a learned embedding of its element and is refined, layer by layer, by exchanging messages with its neighbours. The final features are mapped to atomic energies and summed, so the atomic decomposition, locality within a finite range, and forces by differentiation all survive.

Justin Gilmer and colleagues described the general framework in 2017 as a sequence of message and update steps:

SchNet: continuous-filter convolutions

SchNet, from Kristof Schütt, Klaus-Robert Müller, Alexandre Tkatchenko and colleagues (2017), was the first widely used message-passing potential. Graph networks for molecules had treated bonds as discrete edges, but in a potential energy surface interactions change continuously with distance. SchNet makes the message a continuous-filter convolution: each neighbour’s features are multiplied, element by element, by a filter that is itself a neural-network function of the interatomic distance.

The radial basis expansion plays the same role as the radial symmetry functions of Behler–Parrinello networks, but what is done with it is learned. SchNet uses only distances, so all its features are rotation-invariant. Trained on energies and forces of molecular dynamics trajectories of small organic molecules (the MD17 benchmark of Chmiela and colleagues), it reached chemical accuracy and set a template that PhysNet (2019), with explicit charges and dispersion, and many later models followed.

Receptive field and locality

Each layer lets an atom collect information from neighbours within the cutoff $r_c$. After $T$ layers, information can have travelled $T$ hops, so an atom’s features, and hence its energy, depend on atoms up to about $T r_c$ away: its receptive field. With $r_c = 5$ Å and five layers, that is 25 Å, which can include thousands of atoms.

A large receptive field is a mixed blessing. It lets the model capture interactions somewhat beyond one cutoff and gives each layer a chance to refine the picture. But it makes the model less local, which hurts in two ways. Large simulations are parallelised by dividing space among processors, each needing copies of atoms within the interaction range of its region (ghost atoms); for a message-passing model that range is $T r_c$, so the ghost regions, and the communication, become very large. And a model whose energy depends on a 25 Å sphere needs training data that represent such large environments. These costs motivated the strictly local Allegro architecture described below.

Beyond distances

An invariant network that uses only distances must reconstruct angles and other geometric information indirectly through several layers, and some arrangements are hard or impossible for it to distinguish. DimeNet (Johannes Gasteiger, Janek Groß and Stephan Günnemann, 2020) addressed this with directional message passing: messages are attached to bonds rather than atoms, and a message from bond $kj$ to bond $ji$ depends on the angle $\angle kji$, expanded in a spherical basis. Angles enter explicitly, at the cost of operating on all triplets of atoms. A more general and efficient solution is to let the features themselves carry direction, which requires understanding how quantities transform under rotation.

How quantities transform under rotation

Physical quantities fall into types according to how they change when the system is rotated. A scalar, such as energy or charge, does not change. A vector, such as a force, a dipole or a bond direction, rotates with the system: $\mathbf v \mapsto \mathbf R\mathbf v$. A second-rank tensor, such as the stress or a quadrupole, transforms as $\mathbf T \mapsto \mathbf R\mathbf T\mathbf R^{\mathsf T}$. Group theory organises these types into irreducible representations of the rotation group SO(3), labelled by an integer $l = 0, 1, 2, \ldots$ and of dimension $2l + 1$. Type $l = 0$ is a scalar, $l = 1$ a vector, $l = 2$ the five independent components of a symmetric traceless tensor, and so on. Under a rotation $\mathbf R$, a type-$l$ quantity is multiplied by a $(2l+1)\times(2l+1)$ Wigner D-matrix $\mathbf D^{(l)}(\mathbf R)$.

The natural way to produce type-$l$ quantities from geometry is with spherical harmonics $Y_{lm}(\hat{\mathbf r})$, $m = -l, \dots, l$, evaluated on the unit vector between two atoms. They are the three-dimensional generalisation of the circular harmonics $e^{im\theta}$ in the figure of the chapter on what an MLIP must learn, and they transform exactly as type-$l$ objects. An equivariant network keeps, for each atom, features of several types, $\mathbf h_i = \bigoplus_l \mathbf h_i^{(l)}$, and uses only operations that commute with rotation, so that rotating the input rotates every feature by its own $\mathbf D^{(l)}$.

Two equivariant quantities are combined with the tensor product, the operation familiar from the addition of angular momenta in quantum mechanics:

The example under the equation shows that this is not exotic: the dot and cross products of vectors are tensor products. Nonlinearities need care, because applying, say, a SiLU to each component of a vector would break its rotation behaviour; equivariant networks apply nonlinearities only to scalars, or use gates, multiplying each higher-$l$ feature by a nonlinear function of a scalar. Libraries such as e3nn implement these operations generally.

PaiNN and NequIP

PaiNN (Schütt, Unke and Gastegger, 2021) is the simplest equivariant design: each atom carries scalar features and vector features ($l \le 1$), and messages update both, using the neighbour direction to create vectors from scalars and dot products to create scalars from vectors. Directional information flows through the network without needing angles explicitly.

NequIP (Simon Batzner, Albert Musaelian, Boris Kozinsky and colleagues; preprint 2021, published 2022) uses features of general $l$ and builds each message as a tensor product of the neighbour’s features with the spherical harmonics of the bond direction, weighted by a learned radial function:

The most striking result was data efficiency. On molecular benchmarks, NequIP trained on a few hundred structures matched or exceeded invariant models trained on tens of thousands, and its learning curves, the test error against training-set size, were steeper. Physically, built-in directional information means the network does not have to discover geometry from distances, so every training example teaches more. The price is computational: tensor products with higher $l$ are expensive, and equivariant models are slower per atom than invariant ones of similar size.

MACE: many-body messages

In a two-body message, each neighbour contributes separately and the messages are summed: the message is a sum of pair terms. Many-body correlations, such as the joint arrangement of three neighbours, must then be built up over several layers. MACE (Ilyes Batatia, Dávid Péter Kovács, Gregor Simm, Christoph Ortner and Gábor Csányi, 2022) instead makes each message many-body, using the construction of the atomic cluster expansion on learned features. It first forms an equivariant atomic base by a single sum over neighbours, then takes tensor products of this base with itself up to a correlation order $\nu$, and finally forms a message from these products:

Because of the density trick, the products cost far less than explicit sums over triples or quadruples of neighbours, and with correlation order 3 each message already describes four-body correlations. As a result, MACE needs only two message-passing layers to reach high accuracy, which keeps its receptive field small (two cutoffs) and makes it fast for an equivariant model. MACE achieved state-of-the-art accuracy on standard benchmarks with modest cost, and it became the architecture of the MACE-MP-0 foundation model and of many system-specific potentials. Batatia and colleagues also showed that many earlier architectures, from Behler–Parrinello networks and ACE to NequIP, fit into a single “design space” of equivariant, atom-centred models differing in body order, symmetry type and number of layers.

Allegro: strictly local and scalable

Allegro (Albert Musaelian, Simon Batzner and colleagues, 2023) takes the opposite approach to the receptive-field problem. It keeps equivariant features but attaches them to pairs of atoms $ij$ rather than to atoms, and updates them using only the neighbours of the central atom $i$, never passing messages between atoms. Every pair energy therefore depends only on the cutoff sphere of atom $i$: the model is strictly local, like a Behler–Parrinello network, but with learned equivariant features. Strict locality makes it straightforward to distribute a simulation over many GPUs with small ghost regions, and Allegro demonstrated simulations of protein–water systems with tens to hundreds of millions of atoms at high accuracy.

Equivariant transformers and the trade-offs

Attention, the mechanism of transformers, has also been made equivariant: in Equiformer and EquiformerV2 (Yi-Lun Liao, Tess Smidt and colleagues, 2023–2024), attention weights between neighbours are computed from invariant combinations of equivariant features, and the values being mixed are equivariant. EquiformerV2 used computational shortcuts (rotating each edge into a frame aligned with the bond so that tensor products become cheaper) to reach higher $l$, and it underlies several of the largest models trained on the Open Catalyst and OMat24 data sets. At the other extreme, some fast models, such as ORB, drop exact equivariance and learn approximate symmetry from rotated training data, trading guarantees for speed.

The practical choice depends on the task. For a single system with limited DFT data, the data efficiency of equivariant models (NequIP, MACE, Allegro) is decisive. For very large simulations, strictly local models (Allegro, Behler–Parrinello, ACE) parallelise best. For many-element chemistry, learned element embeddings make message-passing models scale gracefully with the number of elements, which is why universal potentials are almost all of this kind (the chapter on universal potentials). Mature open-source implementations exist for nearly all of them, including SchNetPack, NequIP and Allegro, MACE, DeePMD-kit, TorchANI and fairchem, most with LAMMPS or ASE interfaces.

References and further reading

The training set is the model

An MLIP knows only what its training data show it. Within the region of configuration space covered by the data it interpolates, and its errors are typically a few meV per atom; outside that region it extrapolates, and its predictions can be wrong by electronvolts. The architecture decides how efficiently the model learns from data, but the data decide what it can describe. Most of the effort in building a reliable potential, and most of the expertise, goes into the training set.

A good training set covers every kind of configuration the intended simulations will visit, and some margin beyond. For a material, this usually means the relevant crystal phases at their equilibrium volumes and under compression and tension; distorted cells with randomly displaced (“rattled”) atoms; snapshots of molecular dynamics at and above the target temperatures, including the liquid; point defects, surfaces and, where relevant, grain boundaries and interfaces; and configurations with atoms closer together than they ever are at equilibrium, which teach the model the repulsive wall. For chemistry in solution or at interfaces, it includes the molecules, the solvent, and configurations along the reaction paths of interest. The GAP-18 silicon potential, for example, was trained on a database assembled by hand from 18 kinds of structures, from crystal phases and vacancies to surfaces, liquid and amorphous silicon.

Molecular dynamics snapshots are the commonest source of configurations, but consecutive frames are almost identical and add little information, as the chapter on generalisation explains. Data are therefore selected for diversity: frames are compared in descriptor space and chosen by farthest-point sampling or by CUR decomposition, so that each new structure differs as much as possible from those already chosen. Random structure searching, normal-mode sampling for molecules, and high-temperature or biased sampling extend coverage to configurations that equilibrium dynamics rarely visits.

Consistent reference calculations

Every label in the training set must come from the same electronic-structure method with the same settings: the same functional and dispersion correction, pseudopotentials or projector-augmented-wave data sets, plane-wave cutoff, density of k-points, smearing and spin treatment. Differences that seem minor, such as a k-point mesh that is converged for one cell size but not another, introduce noise of several meV per atom that no model can fit. Total energies from DFT also contain an arbitrary constant that depends on the pseudopotentials; MLIPs handle it by subtracting a reference energy for each element, either the energy of an isolated atom computed with the same settings or a per-element average fitted by linear regression, and learning only the remaining interaction energy.

Data are stored in formats that keep positions, cell, energy, forces and stress together. The extended XYZ format, read and written by ASE, is the common choice: one block per structure, with the energy and stress in the comment line and the forces as extra columns.

The loss and how it is weighted

As the chapter on what an MLIP must learn set out, the loss combines squared errors in energies per atom, forces and stresses, each with a weight. The weights decide the balance between them, and the right balance depends on the application: molecular dynamics depends on forces, phase stability on energy differences, and constant-pressure simulations on stresses. Typical schemes weight forces heavily at first, so that the overall shape of the surface is learned quickly from the abundant force labels, and increase the energy weight towards the end of training, so that energy differences between structures are accurate. Training uses mini-batches of structures, the Adam optimiser with a decaying learning rate, and early stopping on the validation error; many codes also keep an exponential moving average of the weights, which gives smoother and slightly more accurate final models.

The figure below shows why forces matter so much, on the smallest possible example: an ensemble of five small networks fitted to the energy of a diatomic molecule as a function of bond length.

With four training points and energies only, the networks agree at the four points and disagree everywhere else. Adding the forces at the same four points, the slopes of the curve, pins the shape down between them. Remove the short distances from the training range and the repulsive wall vanishes from all five models: a real simulation with such a potential would let atoms fall onto each other, the commonest catastrophic failure of a new MLIP.

Learning curves and data efficiency

The most informative single diagnostic of a model and a data set is the learning curve: the test error plotted against the number of training structures on log–log axes. For a well-behaved model it is a straight line, the error falling as a power of the data, $\epsilon \propto N^{-\beta}$, until it levels off at the noise floor of the reference data or the limit of the model’s expressiveness. Comparing learning curves rather than single numbers separates two questions: how good a model can eventually become, and how much data it needs to get there. Equivariant message-passing models, as the chapter on message passing describes, have steeper learning curves than invariant ones, which is the precise meaning of their data efficiency. Learning curves also tell you whether more data would help: a curve still falling steeply means the model is data-limited.

Uncertainty

A potential that knew when it was wrong would be far safer to use, and several methods provide an estimate of uncertainty. The simplest and most widely used is an ensemble (or committee): train several models that differ in their random initialisation or in the subset of data they see, and use the spread of their predictions as an uncertainty. Where the models agree they are usually right, and where they disagree they are usually wrong, as the figure shows outside the training range. Gaussian-process potentials such as GAP provide a predictive variance directly (kernel and linear potentials), linear models provide extrapolation grades, and cheaper approximations exist for neural networks, such as estimating the variance from only the last layer.

Uncertainty estimates are rarely perfectly calibrated: an ensemble spread of 50 meV/Å does not mean that the error is 50 meV/Å, only that it is likely to be larger than where the spread is 10. They are best used as a ranking: to flag configurations in a simulation that need checking, and to choose which new structures to compute in active learning. Most modern training pipelines, such as DP-GEN, FLARE and the active-learning modes of MTP and ACE codes, are built around exactly this loop: simulate, detect uncertain configurations, compute them with DFT, add them, retrain.

Validation beyond the test error

A low test error, a few meV per atom in energy and a few tens of meV/Å in force, is necessary but not sufficient. A random test set is drawn from the same distribution as the training set, so it says little about configurations that a long simulation will eventually reach and that the training set missed. A potential can have excellent average errors and still produce nonsense in dynamics. Fu and colleagues (2023) showed, across many models and systems, that force accuracy is only loosely correlated with simulation stability and with the accuracy of derived properties. Validation must therefore test physics, as Morrow, Gardner and Deringer (2023) set out:

Stability and conservation. Run long molecular dynamics at the highest temperature of interest, in the NVE ensemble to check energy conservation and in the thermostatted ensemble to check that structures do not collapse, explode or drift into unphysical states. A “hole” in the potential, a region where the predicted energy is spuriously low, is often found only this way.

Static properties. Compare equilibrium lattice constants, the energy–volume curve and bulk modulus, elastic constants, phonon dispersions, vacancy and interstitial formation energies, surface energies and stacking-fault energies with DFT. Phonons are a sensitive test of the second derivatives of the surface, which the training never fits directly.

Dynamic and thermodynamic properties. Compare radial distribution functions, densities and diffusion coefficients with ab initio molecular dynamics, and, where possible, with experiment for the target quantity: a melting point, a thermal expansion coefficient, a phase transition pressure.

Processes of interest. Compute energies along the processes the potential will be used for: migration barriers of vacancies and ions, reaction paths, the energy of a dislocation core. These are often outside the bulk of the training data and are where extrapolation errors hide.

The tutorial in the chapter on building an MLIP applies the first three kinds of test to a small copper potential.

Common failure modes

A handful of problems account for most failures. Missing short-range data produces holes in the repulsive wall and atoms collapsing onto each other; it is fixed by adding compressed structures, dimers and close contacts, or a physical repulsive baseline. Extrapolation to unseen environments, such as a surface when only bulk was trained, or a new temperature range, gives silent errors; it is detected by uncertainty estimates and fixed by active learning. An inadequate reference method, such as PBE for dispersion-bound layered materials, gives a potential that is faithful to the wrong answer. Missing long-range physics affects ionic materials, polar interfaces and charged defects, and requires the long-range extensions of neural-network potentials. Electronic degrees of freedom that the model does not see, such as magnetic order, spin states or the charge state of a defect, make the same atomic configuration correspond to several different energies, which no function of positions alone can fit. And data leakage, testing on frames from the same trajectory as the training frames, makes errors look smaller than they are; test sets should come from independent simulations.

References and further reading

What we will build

This chapter builds a complete machine-learned interatomic potential for copper, from raw data to molecular dynamics, using the ideas of the preceding chapters: the atomic decomposition and symmetry functions of neural-network potentials, the energy-and-force loss of what an MLIP must learn, and the validation strategy of training and validation. Every line of code below was run while writing this chapter, on a laptop CPU, and every number quoted is the actual output. The whole exercise takes about an hour of computer time.

The steps are: generate reference data; turn each atom’s environment into symmetry-function descriptors; define the atomic neural network; train it on energies and forces; use it as an ASE calculator and test it beyond its test error; and finally train an equivariant MACE model on the same data, to see how a modern package compares with a model written by hand.

The scripts need Python 3.10 or later with NumPy, PyTorch and ASE (pip install numpy torch ase), and the last step needs pip install mace-torch. To keep the exercise self-contained and fast, the reference energies and forces come from ASE’s EMT potential, a simple many-body potential for metals, instead of from DFT. EMT plays the role of the “quantum-mechanical” reference here; in real work you replace it by a DFT calculator (GPAW, VASP, Quantum ESPRESSO or CP2K through their ASE interfaces) and nothing else in the pipeline changes.

Step 1: reference data

A potential can only describe configurations like those it was trained on, so the data must cover what the simulations will visit. Here we want copper crystals at temperatures up to 1200 K and under ±3 % strain. The script starts from 32-atom cubic cells of face-centred cubic copper at three lattice strains, rattles the atoms slightly, and runs short Langevin molecular dynamics with EMT at four temperatures, saving a snapshot every 20 steps. That gives 480 configurations with 32 atoms each, 15 360 atomic environments and 46 080 force components. Shuffling and splitting gives 400 structures for training and 80 for testing.

# Step 1: a reference data set of copper configurations, labelled with the EMT potential
import numpy as np
from ase.build import bulk
from ase.calculators.emt import EMT
from ase.md.langevin import Langevin
from ase.md.velocitydistribution import MaxwellBoltzmannDistribution
from ase.io import write
from ase import units

rng = np.random.default_rng(0)
frames = []
for T in (300, 600, 900, 1200):
    for strain in (-0.03, 0.0, 0.03):
        atoms = bulk("Cu", "fcc", a=3.61 * (1 + strain), cubic=True).repeat((2, 2, 2))   # 32 atoms
        atoms.rattle(stdev=0.05, seed=int(rng.integers(1_000_000)))
        atoms.calc = EMT()
        MaxwellBoltzmannDistribution(atoms, temperature_K=T, rng=rng)
        dyn = Langevin(atoms, 2 * units.fs, temperature_K=T, friction=0.02, rng=rng)
        for _ in range(40):                       # 40 snapshots, 20 steps apart
            dyn.run(20)
            a = atoms.copy()
            a.info["energy"] = atoms.get_potential_energy()
            a.arrays["forces"] = atoms.get_forces()
            frames.append(a)
rng.shuffle(frames)
write("cu_train.extxyz", frames[:400])
write("cu_test.extxyz", frames[400:])
print(len(frames), "configurations;", len(frames[0]), "atoms each")

The data are written in extended XYZ format, which stores positions, cell, energy and forces together and is read by every MLIP code.

Step 2: describing atomic environments

Each atom is described by 36 numbers: 24 radial and 12 angular symmetry functions over its neighbours within a 5 Å cutoff, which in fcc copper includes the first three neighbour shells (at 2.55, 3.61 and 4.42 Å). The radial functions are Gaussians in the neighbour distance at eight centres and three widths, summed over neighbours and multiplied by the smooth cutoff function. The angular functions are of Behler’s G5 type, $2^{1-\zeta}\sum_{j,k}(1 + \lambda\cos\theta_{ijk})^\zeta e^{-\eta(r_{ij}^2 + r_{ik}^2)}f_c(r_{ij})f_c(r_{ik})$, at three values of $\zeta$, both signs of $\lambda$ and two radial decays.

The descriptors must be differentiable with respect to positions, because forces will be obtained by differentiating through them. The code therefore builds the neighbour list once with ASE (which handles periodic images), and computes all distances and angles from the position tensor inside PyTorch, so that automatic differentiation can follow every operation. The function triplets lists every pair of neighbours of the same central atom, which is what the angular terms sum over.

The choice of widths matters more than any other setting in this tutorial. The first version of this script used radial Gaussians with widths between 1.25 and 5 Å, a common-looking choice. Trained identically, it reached a force error of 296 meV/Å after 25 epochs; the version below, with widths between 0.25 and 0.7 Å, reached 34 meV/Å at the same point. The figure, which runs the plotting code in the page, shows why: the broad functions overlap so much that they cannot tell the first, second and third neighbour shells apart, so the descriptor barely changes when an atom moves between them.

Step 3: the atomic network

The model is the Behler–Parrinello architecture in its simplest form, since there is only one element: one network, shared by every atom, maps a standardised descriptor vector to an atomic energy, and the energy of a structure is the sum over its atoms. The network has two hidden layers of 32 neurons with SiLU activations, a smooth function so that forces are smooth too; that is 2 273 weights in total. The descriptors are standardised with the mean and standard deviation over the training set, stored with the model, and the output bias starts at the mean energy per atom, so that training begins from a sensible baseline.

Step 4: training on energies and forces

Forces are computed exactly as the negative gradient of the predicted energy with torch.autograd.grad. During training this gradient must itself be differentiated with respect to the weights, which is why create_graph=True is passed: PyTorch keeps the graph of the force computation so that the force error can be backpropagated. To process many structures at once, each mini-batch of ten structures is joined into one disconnected graph, with atom indices offset so that the neighbour lists of different structures do not mix, and structure energies are collected with index_add.

The loss weights the squared energy error per atom by 10 and the squared force error per component by 1, in units of eV and eV/Å. Forces dominate the information, 96 numbers per structure against one energy, so a moderate energy weight keeps energies accurate without starving the force fit. Adam with a cosine learning-rate schedule runs for 120 epochs. Here is the complete script for steps 2 to 4:

"""A Behler-Parrinello neural-network potential for copper, from scratch in PyTorch.

Steps: (2) symmetry-function descriptors, (3) an atomic network, (4) training on energies
and forces. Run make_data.py first to create cu_train.extxyz and cu_test.extxyz.
"""
import time
import numpy as np
import torch
from ase.io import read
from ase.neighborlist import neighbor_list

torch.set_default_dtype(torch.float64)
torch.manual_seed(0)

R_C = 5.0                                                   # cutoff radius (Å)
ETA = torch.tensor([1.0, 3.0, 8.0])                         # radial sharpness (1/Ų): widths 0.7 to 0.25 Å
R_S = torch.linspace(2.0, 4.8, 8)                           # radial centres spanning the neighbour shells (Å)
ZETA = torch.tensor([1.0, 2.0, 4.0])                        # angular sharpness
LAMBDA = torch.tensor([-1.0, 1.0])                          # angular sign
ETA_A = torch.tensor([0.05, 0.25])                          # angular radial decay (1/Ų)
N_DESC = len(ETA) * len(R_S) + len(ETA_A) * len(ZETA) * len(LAMBDA)   # 24 radial + 12 angular = 36


def cutoff(r):
    """Smooth cutoff: equal to 1 at r = 0, and to 0 with zero slope at r = R_C."""
    return 0.5 * (torch.cos(torch.pi * r / R_C) + 1.0) * (r < R_C)


def triplets(i):
    """All pairs (a, b) of neighbour entries that belong to the same centre atom."""
    order = np.argsort(i, kind="stable")
    counts = np.bincount(i)
    starts = np.cumsum(counts) - counts
    A, B = [], []
    for s, c in zip(starts, counts):
        idx = order[s:s + c]
        a, b = np.triu_indices(c, 1)
        A.append(idx[a]); B.append(idx[b])
    return np.concatenate(A), np.concatenate(B)


def featurise(atoms, offset=0):
    """Neighbour lists for one structure; positions stay free so forces can be taken later."""
    i, j, S = neighbor_list("ijS", atoms, R_C)
    A, B = triplets(i)
    return dict(i=i + offset, j=j + offset, shift=S @ atoms.cell.array, A=A, B=B)


def descriptors(pos, i, j, shift, A, B):
    """Radial (G2) and angular (G5-type) symmetry functions for every atom, differentiable in pos."""
    n = pos.shape[0]
    d = pos[j] - pos[i] + shift                             # neighbour vectors, periodic images included
    r = d.norm(dim=1)
    fc = cutoff(r)
    g2 = torch.exp(-ETA[:, None, None] * (r - R_S[None, :, None]) ** 2) * fc     # (eta, Rs, pairs)
    G2 = torch.zeros(n, len(ETA) * len(R_S)).index_add(0, i, g2.reshape(-1, len(r)).T)
    cos = (d[A] * d[B]).sum(1) / (r[A] * r[B])               # cosine of each angle j-i-k
    rad = torch.exp(-ETA_A[:, None] * (r[A] ** 2 + r[B] ** 2)) * fc[A] * fc[B]                 # (eta_a, pairs)
    ang = 2 ** (1 - ZETA)[:, None, None] * (1 + LAMBDA[None, :, None] * cos) ** ZETA[:, None, None]   # (zeta, lambda, pairs)
    g4 = ang[None] * rad[:, None, None, :]                                                     # (eta_a, zeta, lambda, pairs)
    G4 = torch.zeros(n, N_DESC - len(ETA) * len(R_S)).index_add(0, i[A], g4.reshape(-1, len(A)).T)
    return torch.cat([G2, G4], dim=1)


class BPNN(torch.nn.Module):
    """Total energy = sum over atoms of NN(G_i); one network for the single element Cu."""
    def __init__(self, g_mean, g_std, e0, hidden=32):
        super().__init__()
        self.register_buffer("g_mean", g_mean)
        self.register_buffer("g_std", g_std)
        self.net = torch.nn.Sequential(
            torch.nn.Linear(N_DESC, hidden), torch.nn.SiLU(),
            torch.nn.Linear(hidden, hidden), torch.nn.SiLU(),
            torch.nn.Linear(hidden, 1))
        with torch.no_grad():
            self.net[-1].bias.fill_(e0)                     # start from the mean energy per atom

    def forward(self, pos, graph, struct, n_struct):
        G = descriptors(pos, *(graph[k] for k in ("i", "j", "shift", "A", "B")))
        e_atom = self.net((G - self.g_mean) / self.g_std).squeeze(-1)
        return torch.zeros(n_struct).index_add(0, struct, e_atom)   # energy of each structure


def batch(frames):
    """Join several structures into one disconnected graph."""
    pos, E, F, struct, parts, off = [], [], [], [], [], 0
    for k, a in enumerate(frames):
        g = featurise(a, off)
        g["A"] = g["A"] + sum(len(p["i"]) for p in parts)    # triplet indices point into pair arrays
        g["B"] = g["B"] + sum(len(p["i"]) for p in parts)
        parts.append(g)
        pos.append(a.positions); E.append(a.get_potential_energy()); F.append(a.get_forces())
        struct += [k] * len(a); off += len(a)
    graph = {k: torch.tensor(np.concatenate([p[k] for p in parts])) for k in ("i", "j", "shift", "A", "B")}
    n_atoms = torch.tensor([len(a) for a in frames], dtype=torch.float64)
    return dict(pos=torch.tensor(np.concatenate(pos)), graph=graph, struct=torch.tensor(struct),
                n=len(frames), n_atoms=n_atoms, E=torch.tensor(E), F=torch.tensor(np.concatenate(F)))


def predict(model, b, create_graph=False):
    pos = b["pos"].clone().requires_grad_(True)
    E = model(pos, b["graph"], b["struct"], b["n"])
    F = -torch.autograd.grad(E.sum(), pos, create_graph=create_graph)[0]   # F = -dE/dr, exactly
    return E, F


def rmse(model, b):
    E, F = predict(model, b)
    e = ((E - b["E"]) / b["n_atoms"]).detach()
    return 1000 * e.pow(2).mean().sqrt().item(), 1000 * (F - b["F"]).detach().pow(2).mean().sqrt().item()


if __name__ == "__main__":
    train_frames, test_frames = read("cu_train.extxyz", ":"), read("cu_test.extxyz", ":")
    test = batch(test_frames)
    batches = [batch(train_frames[k:k + 10]) for k in range(0, len(train_frames), 10)]   # fixed mini-batches
    G = torch.cat([descriptors(b["pos"], *(b["graph"][k] for k in ("i", "j", "shift", "A", "B"))) for b in batches])
    e0 = torch.cat([b["E"] / b["n_atoms"] for b in batches]).mean()
    model = BPNN(G.mean(0), G.std(0) + 1e-8, e0.item())
    opt = torch.optim.Adam(model.parameters(), lr=2e-3)
    sched = torch.optim.lr_scheduler.CosineAnnealingLR(opt, T_max=120)
    w_E, w_F = 10.0, 1.0             # weights for (eV/atom)² and (eV/Å)² errors
    t0 = time.time()
    for epoch in range(1, 121):
        for k in torch.randperm(len(batches)).tolist():
            b = batches[k]
            E, F = predict(model, b, create_graph=True)
            loss = w_E * (((E - b["E"]) / b["n_atoms"]) ** 2).mean() + w_F * ((F - b["F"]) ** 2).mean()
            opt.zero_grad(); loss.backward(); opt.step()
        sched.step()
        if epoch in (1, 10, 25, 50, 80, 120):
            eE, eF = rmse(model, test)
            print(f"epoch {epoch:3d}   test RMSE  energy {eE:6.2f} meV/atom   forces {eF:6.1f} meV/Å   ({time.time() - t0:4.0f} s)")
    torch.save(model.state_dict(), "cu_bpnn.pt")

Running it prints the test error at a few points during training:

The final model predicts the energies of unseen structures to 0.3 meV per atom and their forces to 17 meV/Å, against force magnitudes of about 700 meV/Å in the data: an error of a few per cent. Training took about 33 minutes of wall time on a laptop CPU.

Step 5: using the potential, and testing it properly

A test error is not a validation. The last script wraps the network as an ASE calculator, so that it can drive any ASE workflow, and then asks three questions that the test set cannot answer. Does it reproduce the equation of state of the crystal, including volumes slightly outside the training range? Does it predict the equilibrium lattice constant? And does molecular dynamics with it conserve energy in a cell larger than any it was trained on?

"""Step 5: use the trained network as an ASE calculator, then test it beyond the test RMSE."""
import numpy as np
import torch
from ase.build import bulk
from ase.calculators.calculator import Calculator, all_changes
from ase.calculators.emt import EMT
from ase.md.verlet import VelocityVerlet
from ase.md.velocitydistribution import MaxwellBoltzmannDistribution
from ase import units
from bpnn import BPNN, featurise, descriptors

class BPNNCalculator(Calculator):
    implemented_properties = ["energy", "forces"]
    def __init__(self, model, **kw):
        super().__init__(**kw)
        self.model = model.eval()
    def calculate(self, atoms=None, properties=("energy",), system_changes=all_changes):
        super().calculate(atoms, properties, system_changes)
        g = {k: torch.tensor(v) for k, v in featurise(atoms).items()}        # neighbour list rebuilt every call
        pos = torch.tensor(atoms.positions, requires_grad=True)
        e_atom = self.model.net((descriptors(pos, g["i"], g["j"], g["shift"], g["A"], g["B"]) - self.model.g_mean) / self.model.g_std)
        E = e_atom.sum()
        F = -torch.autograd.grad(E, pos)[0]
        self.results = {"energy": E.item(), "forces": F.numpy()}

state = torch.load("cu_bpnn.pt")
model = BPNN(state["g_mean"], state["g_std"], 0.0)
model.load_state_dict(state)
calc = BPNNCalculator(model)

# (a) energy-volume curve and lattice constant, compared with the reference potential
print(" a (Å)   E_EMT (eV/atom)   E_NN (eV/atom)")
a_list = np.linspace(3.45, 3.80, 8)
for a in a_list:
    at = bulk("Cu", "fcc", a=a, cubic=True).repeat((2, 2, 2))
    at.calc = EMT(); e1 = at.get_potential_energy() / len(at)
    at.calc = calc;  e2 = at.get_potential_energy() / len(at)
    print(f" {a:.3f}   {e1:10.4f}        {e2:10.4f}")
for name, c in (("EMT", EMT()), ("NN", calc)):
    es = []
    for a in a_list:
        at = bulk("Cu", "fcc", a=a, cubic=True).repeat((2, 2, 2)); at.calc = c; es.append(at.get_potential_energy())
    p = np.polyfit(a_list, es, 3); r = np.roots(np.polyder(p)); a0 = r[np.isreal(r)].real; a0 = a0[(a0 > 3.4) & (a0 < 3.9)][0]
    print(f"{name}: equilibrium lattice constant {a0:.3f} Å")

# (b) NVE molecular dynamics: does the total energy stay constant? are forces close to EMT along the way?
atoms = bulk("Cu", "fcc", a=3.61, cubic=True).repeat((3, 3, 3))          # 108 atoms: larger than any training cell
MaxwellBoltzmannDistribution(atoms, temperature_K=600, rng=np.random.default_rng(1))
atoms.calc = calc
dyn = VelocityVerlet(atoms, 2 * units.fs)
etot, ferr = [], []
for step in range(250):                                                    # 500 steps of 2 fs = 1 ps
    dyn.run(2)
    etot.append(atoms.get_potential_energy() + atoms.get_kinetic_energy())
    if step % 25 == 0:
        ref = atoms.copy(); ref.calc = EMT()
        ferr.append(np.sqrt(np.mean((atoms.get_forces() - ref.get_forces()) ** 2)))
etot = np.array(etot) / len(atoms)
print(f"NVE, 108 atoms, 1 ps: total-energy drift {1000 * (etot[-1] - etot[0]):+.2f} meV/atom, "
      f"fluctuation {1000 * etot.std():.2f} meV/atom; T at end {atoms.get_temperature():.0f} K")
print(f"force RMSE against EMT along the trajectory: {1000 * np.mean(ferr):.0f} meV/Å")

The output:

 a (Å)   E_EMT (eV/atom)   E_NN (eV/atom)
 3.450       0.0677            0.0724
 3.500       0.0226            0.0227
 3.550      -0.0015           -0.0012
 3.600      -0.0067           -0.0066
 3.650       0.0046            0.0038
 3.700       0.0303            0.0298
 3.750       0.0686            0.0677
 3.800       0.1181            0.1137
EMT: equilibrium lattice constant 3.590 Å
NN: equilibrium lattice constant 3.590 Å
NVE, 108 atoms, 1 ps: total-energy drift +0.03 meV/atom, fluctuation 0.01 meV/atom; T at end 279 K
force RMSE against EMT along the trajectory: 7 meV/Å

The network reproduces the equilibrium lattice constant exactly and the energy–volume curve to within a meV per atom across the strains it was trained on (lattice constants between about 3.50 and 3.72 Å). At 3.45 and 3.80 Å, outside the training range, the errors grow to 5 meV per atom, the first sign of extrapolation that the chapter on training and validation warns about: the model was never shown such strongly compressed or expanded crystals. The molecular dynamics run, in a 108-atom cell three times larger than the training cells, conserves the total energy to 0.03 meV per atom over a picosecond, which shows that the forces are smooth and exactly conservative, and the forces along the trajectory agree with EMT to 7 meV/Å, better than the test error because room-temperature configurations are near the centre of the training distribution.

Step 6: the same data with MACE

Writing a potential by hand shows every moving part. For production work, mature packages implement far more capable architectures with careful defaults, GPU support and interfaces to LAMMPS. The same data set can be used to train an equivariant MACE model (message-passing potentials) from the command line. MACE expects the reference values under explicit keys, so the data are first copied into new files with the energy and forces stored as REF_energy and REF_forces:

from ase.io import read, write
for name in ("train", "test"):
    frames = read(f"cu_{name}.extxyz", ":")
    for a in frames:
        a.info["REF_energy"] = a.get_potential_energy()
        a.arrays["REF_forces"] = a.get_forces()
        a.calc = None
    write(f"cu_{name}_mace.xyz", frames)

and then trained with

mace_run_train --name="cu_mace" --train_file="cu_train_mace.xyz" --valid_fraction=0.1 \
    --test_file="cu_test_mace.xyz" --energy_key="REF_energy" --forces_key="REF_forces" \
    --E0s="average" --model="MACE" --num_channels=32 --max_L=1 --correlation=3 \
    --r_max=5.0 --batch_size=10 --max_num_epochs=30 --swa --start_swa=20 \
    --default_dtype="float64" --device=cpu --seed=1

This is a deliberately small MACE (32 channels, vector features up to $L = 1$, correlation order 3) trained for only 30 epochs on a CPU. The --swa flag turns on what MACE calls stage two: from epoch 20 the energy weight in the loss rises from 1 to 1000 while the force weight stays at 100. The first stage is dominated by forces, which carry many numbers per configuration, but forces say nothing about how the energy changes under uniform strain, where every atom stays at a symmetric site with zero force; without stage two the energy–volume curve is learned poorly, even when the forces are excellent.

Run on a laptop CPU, the 30 epochs took about an hour and three quarters. The log shows the two stages clearly. At epoch 18, before stage two, the validation errors were 9 meV/Å in the forces but 78 meV/atom in the energy. Ten epochs into stage two the energy error had fallen to 3 meV/atom, at some cost in the forces. On the held-out test set the final model reaches 3.3 meV/atom and 30.6 meV/Å, against 0.29 meV/atom and 16.6 meV/Å for the hand-written network, which trained four times as many epochs.

The physical test then tells a different story. Fitting the energy–volume curve of the perfect crystal over 3.45–3.75 Å gives a lattice constant of 3.727 Å from MACE, against 3.597 Å from EMT over the same range, and at $a = 3.65$ Å, inside the range of the training data, the perfect-crystal energy is wrong by 46 meV/atom. The reason is the data set, not the architecture. Step 1 generated configurations at only three volumes, strains of −3 %, 0 and +3 %, and the test set was drawn from the same three, so a small test error certifies accuracy at those three volumes and nothing in between. The Behler–Parrinello network, whose fixed radial functions change smoothly with distance, happened to interpolate correctly between them; the more flexible message-passing model did not. The remedy is to sample volume continuously, or at many strains, and to include perfect and slightly strained crystals, exactly as the chapter on training data recommends. It is the most important lesson of this tutorial: a test error measured on data drawn like the training data is not a validation of the physics.

The trained model is saved as cu_mace.model and is used from ASE with from mace.calculators import MACECalculator; calc = MACECalculator(model_paths="cu_mace.model", device="cpu"), or from LAMMPS for large simulations. Fine-tuning a foundation model instead of training from scratch, as described in the chapter on universal potentials, uses the same command with --foundation_model pointing to a pre-trained model.

From here to a research potential

Turning this exercise into a potential for research changes four things. The reference becomes DFT, computed with consistent settings, which makes each configuration seconds to hours of computer time instead of milliseconds, so every configuration must earn its place. The data set grows by active learning: train an ensemble, run molecular dynamics, compute the configurations on which the models disagree, retrain (active learning). The system gains elements, which for a Behler–Parrinello potential means element-resolved descriptors and one network per element, and for MACE needs no change at all. And validation extends to the properties the potential will be used for: defects, surfaces, phonons, diffusion, melting, whatever the research question requires. The structure of the code, data in, descriptors or learned features, an atomic model, a loss on energies and forces, an ASE calculator out, stays exactly as here.

References and further reading

From one material to the periodic table

Until about 2021, every machine-learned potential was a specialist. A potential for silicon knew nothing about oxygen; a potential for water could not describe a salt solution unless salt had been in its training set. Building one took months of DFT calculations, training and validation, and the result applied to one chemical system. This made MLIPs powerful tools for groups studying a particular material in depth, but useless for the broad questions of materials discovery: which of ten thousand candidate compounds are stable, which conduct lithium ions quickly, which might catalyse a reaction.

Two developments changed this. High-throughput DFT projects had accumulated enormous numbers of calculations across inorganic chemistry, most importantly the Materials Project, which since 2011 has computed the structures and energies of well over a hundred thousand compounds with consistent settings. And message-passing networks with learned element embeddings (the previous chapters) scale gracefully to many elements, because adding an element adds one embedding vector rather than a new set of descriptors and a new network. Training such a network on a large fraction of the Materials Project’s data produced the first universal interatomic potentials: models that give reasonable energies and forces for almost any inorganic material out of the box. Their role is often compared to that of foundation models in language processing: a large model pre-trained on broad data, used directly for many tasks or adapted to a specific one with little additional data.

The data sets

A universal potential is only as broad as its training data, and the history of these models is in large part a history of data sets.

The MPtrj data set, assembled by the authors of CHGNet in 2023, collected the intermediate structures of the Materials Project’s geometry relaxations: about 1.6 million structures from about 146 000 materials, with energies, forces, stresses and magnetic moments, mostly computed with the PBE and PBE+U functionals. Because relaxation trajectories start from reasonable guesses and move towards equilibrium, MPtrj is dominated by near-equilibrium structures. Most universal potentials up to 2024 were trained on it.

Larger and more diverse data followed. The Alexandria database (Schmidt, Marques and colleagues) provides millions of additional DFT calculations of hypothetical crystals. Google DeepMind’s GNoME project generated DFT data for hundreds of thousands of newly predicted stable crystals and their relaxations. Meta’s OMat24 data set (2024) contains about 110 million DFT calculations, deliberately including structures far from equilibrium, rattled, strained and sampled from high-temperature molecular dynamics, precisely to cover the configurations that relaxation trajectories miss. The Open Catalyst data sets OC20 and OC22 cover adsorbates on catalytic surfaces. For organic and biomolecular chemistry, the SPICE data set provides about a million calculations at a higher level of DFT on drug-like molecules, peptides and their complexes.

Mixing these data sets is not trivial. They use different functionals (PBE, PBE+U, r²SCAN, ωB97M-D3), different pseudopotentials and different settings, so their energies are not directly comparable. Training on several at once requires either restricting to compatible subsets, learning separate output “heads” for each level of theory, or careful correction schemes.

The models

M3GNet (Chi Chen and Shyue Ping Ong, 2022) was one of the first universal potentials: a graph network with explicit three-body interactions inside its message passing, trained on Materials Project relaxation data across 89 elements. It showed that a single model could relax structures and run dynamics across inorganic chemistry. CHGNet (Bowen Deng, Gerbrand Ceder and colleagues, 2023) was trained on MPtrj and also predicts magnetic moments, which act as a proxy for the charge states of transition-metal ions, important in battery materials where the oxidation state of a metal changes during cycling.

GNoME (Amil Merchant, Ekin Dogus Cubuk and colleagues at Google DeepMind, 2023) used graph networks in an active-learning loop: candidate structures were generated by substitution and random search, filtered by the network, computed with DFT, and added to the training data, over several rounds. The process proposed about 2.2 million structures predicted to lie on or near the convex hull of stability, of which roughly 380 000 were stable relative to all known materials, and the networks trained along the way became universal potentials themselves. The claim of novelty and stability for so many structures prompted discussion about how many are genuinely new and synthesisable, a useful reminder that computational stability is a necessary but not a sufficient condition for a real material.

MACE-MP-0 (Ilyes Batatia, Gábor Csányi and a large collaboration, 2023) is a MACE model (message-passing potentials) trained on MPtrj. Its paper demonstrated stable molecular dynamics, out of the box, across a remarkable range of systems, from water and ice to catalysis, batteries, metal–organic frameworks and high-entropy alloys. It is distributed in several sizes as an ASE calculator, and has become a common starting point for fine-tuning.

Several further models appeared in 2024. MatterSim (Han Yang and colleagues at Microsoft Research) was trained on structures sampled across elements, temperatures up to 5000 K and pressures up to 1000 GPa, aiming at finite-temperature and high-pressure properties. SevenNet is a NequIP-based universal model with a parallel algorithm for large-scale dynamics. ORB (Orbital Materials) is a fast model that is not built to be exactly equivariant but learns rotational symmetry approximately from rotated training data. EquiformerV2 models trained on OMat24 (the eqV2 models of Meta’s fairchem project) reached the top of the stability-prediction benchmarks at the time of their release. For organic molecules, MACE-OFF (trained on SPICE), AIMNet2 and the ANI family play the same role. The field moves quickly, and newer models appear every few months; the principles below matter more than any ranking.

How universal potentials are benchmarked

A natural test of a universal potential is whether it can do the job materials discovery needs: predict which hypothetical crystals are thermodynamically stable. Matbench Discovery (Janosh Riebesell, Rhys Goodall, Alpha Lee and colleagues) frames exactly this. Models are given about 250 000 unrelaxed hypothetical crystal structures from the WBM data set, generated by element substitution; they must relax each structure and predict its energy relative to the convex hull of known phases, and are scored on how well they classify structures as stable or unstable, with metrics such as the F1 score and the discovery acceleration factor (how much more often a model’s predicted-stable structures are truly stable than random picks would be). Later versions added a test of thermal conductivity, which depends on third derivatives of the energy (anharmonic phonons) and probes the shape of the potential energy surface far more severely than energies do.

These benchmarks were valuable because they compared models on the same task with the same data, and they showed universal potentials outperforming one-shot property predictors at screening. But a good leaderboard score does not guarantee accurate dynamics for a particular system: the benchmark rewards energies near equilibrium for typical crystals, while a simulation of a liquid, an interface or a defect may depend on regions of the energy surface that the benchmark never tests. The validation principles of the chapter on training and validation still apply.

Systematic errors

Universal potentials have characteristic errors that every user should know. The most widespread is softening of the potential energy surface. Models trained largely on near-equilibrium data, such as relaxation trajectories, systematically underestimate the curvature of the energy surface: predicted forces are too small, phonon frequencies too low, elastic constants and bulk moduli too small, and surface and defect energies often too low. Deng and colleagues showed in 2024 that this softening appears consistently across universal models and can be largely corrected by fine-tuning with a modest number of structures sampled away from equilibrium. Data sets such as OMat24, which deliberately include rattled and high-temperature structures, were built partly to address it.

Universal potentials also inherit the errors of their reference. A model trained on PBE overestimates lattice constants slightly and misses most dispersion, unless a dispersion correction is added (MACE-MP-0 can be combined with the D3 correction); PBE+U data mix two Hamiltonians, so energies of compounds of the same elements computed with and without $U$ must be reconciled. And coverage is uneven: chemistries, structures and conditions that are rare in the training data, such as molecular crystals, liquids, surfaces with adsorbates, very high pressures or radiation damage, may be described much less well than common oxides and intermetallics.

Using and fine-tuning a universal potential

Most universal potentials are distributed as Python packages with an ASE calculator, so they slot into any ASE-based workflow: relaxation, molecular dynamics, nudged elastic band calculations, phonons. The following example relaxes face-centred cubic copper with MACE-MP-0 (it was run while writing this chapter; the first call downloads the model):

from ase.build import bulk
from ase.filters import FrechetCellFilter
from ase.optimize import BFGS
from mace.calculators import mace_mp           # pip install mace-torch

calc = mace_mp(model="medium", default_dtype="float64", device="cpu")
atoms = bulk("Cu", "fcc", a=3.7, cubic=True)
atoms.calc = calc
BFGS(FrechetCellFilter(atoms), logfile=None).run(fmax=1e-3)   # relax positions and cell
print(atoms.cell[0, 0])                                       # 3.626 Å; experiment 3.615 Å

The predicted lattice constant of 3.626 Å is close to the PBE value for copper, about 0.3 % above experiment, exactly as expected for a model that reproduces PBE. The same calc object can drive molecular dynamics with ASE’s integrators, or be used from LAMMPS through the MACE interface for larger systems.

For production work on a specific system, the recommended route is usually fine-tuning: start from the pre-trained universal model and continue training on a few hundred to a few thousand DFT calculations of the system of interest, computed with settings compatible with the foundation model’s training data. The pre-trained model supplies a good description of chemistry in general, and the fine-tuning data correct its systematic errors for the configurations that matter. To avoid “catastrophic forgetting”, in which fine-tuning degrades the model elsewhere, MACE offers a multi-head scheme that keeps training on a sample of the original data alongside the new data. Fine-tuned foundation models often match or exceed potentials trained from scratch on the same system with far less data.

A sensible workflow for a new project is therefore: screen and explore with a universal potential; validate it on a handful of DFT calculations of configurations that matter; if the errors are too large, fine-tune on system-specific data, possibly selected by active learning; and validate the fine-tuned model with the property-based tests of the previous chapter before trusting long simulations.

References and further reading

The question every lab faces

You can afford twenty more syntheses, or two hundred more DFT calculations. Which ones? Choosing by intuition, on a grid or at random wastes most of the budget. Sequential experimental design uses everything measured so far to choose the next experiment, so that each one is as informative or as promising as possible. Two closely related goals use the same machinery.

• Active learning chooses data to make a model good everywhere it will be used, for example DFT calculations to train an interatomic potential (the lesson on interatomic potentials).
• Bayesian optimisation (BO) chooses experiments to find the best material, the composition or process conditions that maximise a property, in as few trials as possible.

The loop

Both are a loop of four steps:

1. Fit a surrogate model to all data so far, one that predicts both a mean $\mu(\mathbf x)$ and an uncertainty $\sigma(\mathbf x)$. Gaussian processes (the lesson on Gaussian processes) are the classic choice for small data; ensembles of networks or of trees for larger data.
2. Compute an acquisition function $a(\mathbf x)$ that scores how useful it would be to measure at $\mathbf x$.
3. Run the experiment (or calculation) at $\mathbf x_{\text{next}} = \arg\max_{\mathbf x} a(\mathbf x)$.
4. Add the result to the data and repeat.

Acquisition functions: exploration against exploitation

Every acquisition function balances two motives: exploit regions where the model already predicts high values, and explore regions where it is uncertain and something better may hide.

Upper confidence bound (UCB). Optimism in the face of uncertainty:

With $\kappa = 0$ it is pure exploitation and tends to get stuck on the first good peak it finds (try it in the lab: the broad peak near $x = 0.28$ captures it). Large $\kappa$ explores widely. Srinivas et al. (2010) proved convergence guarantees for suitably growing $\kappa$.