Machine Learning Force Fields

Glossary

Introduction

In 1964, physicist Richard Feynman famously remarked “that all things are made of atoms and that everything that living things do can be understood in terms of the jigglings and wigglings of atoms”. 1 As such, an atomically resolved picture can provide invaluable insights on biological (and other) processes. The first molecular dynamics (MD) study of a protein in 1977 by McCammon et al. 2 did not consider explicit solvent molecules and was limited to less than 10 ps of simulation. Still, it challenged the (at that time) common belief that proteins are essentially rigid structures3 and, instead, suggested that the interior of proteins behaves more fluid-like. Since then, systems consisting of more than a million atoms have been studied4, simulation times extended to the millisecond regime5, and the study of entire viruses in atomic detail made possible6, 7. Recently, a distributed computing effort even allowed to study the viral proteome of SARS-CoV-2 for a total of 0.1 seconds of simulation time 8.

To perform MD simulations, typically, the Newtonian equations of motion are integrated numerically, which requires knowledge of the forces acting on individual atoms at each time step of the simulation9. In principle, the most accurate way to obtain these forces is by solving the Schrödinger equation (SE), which describes the physical laws underlying most chemical phenomena and processes10. Unfortunately, an analytic solution is only possible for two-body systems such as the hydrogen atom. For larger chemical structures, the SE can only be solved approximately. However, even with approximations, an accurate numerical solution is a computationally demanding task. For example, the CCSD(T) method (coupled cluster with singles, doubles and perturbative triples), which is widely regarded as the “gold standard” of chemistry (as its predictions compare well with experimental results)11, scales $\propto N^{7}$ with the number of atoms $N$. Strictly speaking, the true scaling of the CCSD(T) method is $\mathcal{O}(n^{7})$, where $n$ is the number of basis functions used for the wave function ansatz. Depending on the desired accuracy and which atoms are present (more precisely, how many electrons are in their shells), $n$ can vary greatly. However, the number of atoms is usually a good proxy. Because of this, it is unfeasible to calculate the forces for many different configurations of large chemical systems, which is required for running MD simulations, with accurate methods. Instead, simple empirical functions are commonly used to model the relevant interactions. From these so-called force fields (FFs), atomic forces can be readily derived analytically.

Most conventional FFs model chemical interactions as a sum over bonded and non-bonded terms12, 13. The former can be described with simple functions of the distances between directly bonded atoms, or angles and dihedrals between atoms sharing some of their bonding partners. The non-bonded terms consider pair-wise combinations of atoms, typically by modeling electrostatics with Coulomb’s law (assuming a point charge at each atom’s position) and dispersion with a Lennard-Jones potential14. Due to the computational efficiency of these terms, such classical FFs allow to study systems consisting of many thousands of atoms. However, while offering a qualitatively reasonable description of chemical interactions, the quality of MD simulations, and hence the insights that can be obtained from them, are ultimately limited by the accuracy of the underlying FF.15 This is particularly problematic when polarization, or many-body interactions in general, are of significant importance, as these effects are not adequately modeled by the terms described above. While it is possible to construct polarizable FFs16, 17, 18, 19 or account for other important effects, e.g. anisotropic charge distributions20, 21, to improve accuracy (in exchange for computational efficiency), it is not always clear a priori when such modifications are necessary. Beyond that, conventional FFs require a preconceived notion of bonding patterns and thus cannot describe bond breaking or bond formation. While there exist reactive FFs offering an approximated description of reactions22, 23, 24, they are often not sufficiently accurate for quantitative studies or restricted to specific types of reactions. Mixed quantum mechanics/molecular mechanics (QM/MM) treatments25 pose an alternative solution: Here, the SE is only solved for a small part of the system where high accuracy is required or reactions may take place, e.g. the active site of an enzyme. Meanwhile, all remaining atoms are treated at the FF level of accuracy. Although this is more efficient than a pure quantum-mechanical approach, it is often necessary to model a large number of atoms at the QM level for converged results26, which is still highly computationally demanding.

In a “dream scenario” for computational chemists and biologists, it would be possible to treat even large systems at the highest levels of theory, which would require prohibitively large computational resources in the real world. Machine learning (ML) methods could help to achieve this dream in a principled manner, thus closing the gap between the accuracy of ab initio methods and the efficiency of classical FFs (Fig. 1). ML methods aim to learn the functional relationship between inputs (chemical descriptors) and outputs (properties) from patterns or structure in the data. Ideally, a trained learning machine would then reflect the underlying effective “rules” of quantum mechanics.27 Practically, ML models can take a shortcut by not having to solve any equations that follow from the physical laws governing the structure-property relation. Because of this unique ability, ML methods have been enjoying growing popularity in the chemical sciences in recent years. They allow to explore chemical space and predict the properties of compounds with both unprecedented efficiency and high accuracy.28, 29, 30, 31, 32, 27 ML has also been used to augment and accelerate traditional methods used in molecular simulations, e.g. for sampling equilibrium states33, 34 and rare events35, computing reaction rates36, exploring protein folding dynamics37 and other biophysical processes38, 39, 40, 41, 42, Markov state modeling43, 44, 45, 46, 47, 48, 49, and coarse-graining50, 51, 52, 53 (for a recent review on applications of ML in molecular simulatons, see Ref. 54). Recent advances made it even possible to predict molecular wave functions, which can act as an interface between ML and quantum chemistry55, 56, as knowledge of the wave function allows to deduce many different quantum mechanical observables at once. ML can also be combined with more traditional semi-empirical methods, for example by predicting accurate repulsive potentials for density functional tight-binding approaches.57 Instead of circumventing equations, ML methods can also help when solving them: They have been used to find novel density functionals58, 59, 60 and solutions of the Schrödinger equation61, 62. Other promising applications include the generation of molecular structures to tackle inverse design problems,63, 64, 65, 66, 67 or planning chemical syntheses68.

For constructing ML-FFs, suitable reference data to learn the relevant structure-property relation include energy, forces, or a combination of both, obtained from ab initio calculations. Contrary to conventional FFs, no preconceived notion of bonding patterns needs to be assumed. Instead, all chemical behavior is learned from the reference data. This allows to reconstruct the important interactions purely from atomic positions without imposing a restricted analytical form on the interatomic potential and enables a natural description of chemical reactions. For example, it is now possible to construct ML-FFs for small molecules from CCSD(T) reference data close to spectroscopic accuracy and with computational efficiency similar to conventional FFs69, 70. This has enabled studies that would be prohibitively expensive with conventional methods of computational chemistry and allowed to obtain novel chemical insights (see Fig. 2).

Other properties than energies and forces can be predicted as well: For example, dipole moments, which are a measure for how polar molecules are, can be used to calculate infrared spectra from MD simulations71, 72, 73 and allow a comparison to experimental measurements. Other prediction targets could be used to screen large compound databases for molecules with desirable properties several orders of magnitude faster than with ab initio methods. The HOMO/LUMO gap, which is important for materials used in electronic devices such as solar cells74, is only one prominent example of many potentially interesting prediction targets.

This review will focus on the construction of ML-FFs for the usage in MD simulations and other applications (for details on how to set up such simulations or how to extract physical insights from them, refer to Refs. 75, 76, 77). The text is structured as follows: Section 2 reviews fundamental concepts of chemistry (2.1) and machine learning (2.2) relevant to the construction of ML-FFs and discusses special considerations when the two are combined (2.3). As this article is intended for both chemists and machine learning practitioners, these sections provide all readers with the necessary background to understand the remainder of the review. Experts in either field may want to skip over the respective sections, as they discuss fundamentals. The next part (Section 3) serves as a step-by-step guide and reference for readers that want to apply ML-FFs in their own research. Here, the best practices for constructing ML-FFs are outlined, possible problems that may be encountered along the way (and how to avoid them) are discussed and an overview of several software packages, which may be used to accelerate the construction of ML-FFs, is provided. Section 4 lists several applications of ML-FFs described in the literature and highlights physical and chemical insights made possible through the use of ML. The review is concluded in Section 5 by a discussion of obstacles that still need to be overcome to extend the applicability of ML-FFs to an even broader context.

Mathematical and conceptual framework

Section 2.1 reviews important chemical concepts such as the potential energy surface and invariance properties of physical systems, which are essential for constructing physically meaningful models. It is meant as a short summary of the most important physical principles and fundamental chemical knowledge for readers with a primarily ML-focused background who are interested in constructing ML-FFs. On the other hand, to offer readers with a chemical background a first orientation, an overview of two important methodologies in Machine Learning, namely kernel-based learning approaches and artificial neural networks, is given in Section 2.2. Finally, Section 2.3 lists constraints related to the physical invariances mentioned earlier and gives examples of models for constructing ML-FFs and how they implement these constraints in practice.

The Schrödinger equation (SE)78, which describes the interaction of atomic nuclei and electrons, is sufficient for understanding most chemical phenomena and processes10. Unfortunately, it can only be solved analytically for very simple systems, such as the hydrogen atom. For more complex systems like molecules, exact numerical solutions are often impractical due to a steep increase of computational costs as a function of system size. For this reason, numerous approximation schemes have been devised to enable insights into more complicated chemical systems. Virtually all of these are based on the Born-Oppenheimer (BO) approximation79, which decouples electronic and nuclear motion so that the latter can be neglected. It is assumed that electrons adjust instantaneously to changes in the nuclear positions, which is motivated by the observation that atomic nuclei are heavier than electrons by several orders of magnitude, thus moving on a vastly different timescale. Hence, the nuclear positions appear almost stationary to the electrons and therefore enter the resulting “electronic SE” only parametrically: The energy of the electrons depends on the external potential caused by the nuclei, which in turn is fully determined by their positions and nuclear charges. By summing electronic energy and Coulomb repulsion between nuclei, the total potential energy of the system is obtained, which is one of the most important properties of molecules. Alongside entropic contributions, it determines the relative stability of different compounds, whether reactions are exothermic or endothermic, and can even serve as proxy for more complex properties. For example, the potency of some drugs can be estimated from their binding energy to biomolecules80.

By introducing a parametric dependency between energy and nuclei, the BO approximation implies the existence of a functional relation $f:\{Z_{i},\mathbf{r}_{i}\}_{i=1}^{N}\mapsto E$, which maps the nuclear charges $Z_{i}$ and positions $\mathbf{r}_{i}$ of $N$ atoms directly to their potential energy $E$. This function, called the potential energy surface (PES), governs the dynamics of a chemical system, similar to a ball rolling on a hilly landscape. Minima (“valleys”) on the PES correspond to stable molecules and significant structural changes (or even chemical reactions) occur when a system crosses over a transition state (“ridge”) from one minimum into another (Fig. 3).

Knowledge of the PES therefore also allows to predict how a system evolves over time. For example, by studying a thermal ensemble of molecules starting from the same minimum on the PES, it is possible to determine which fraction of them will reach different minima and in what time frame, allowing to assess their reactivity and which products are formed. It is also possible to deduce the macroscopic thermodynamic properties of a system by studying how it behaves at an atomic level. In such molecular dynamics (MD) simulations, a classical treatment of nuclear dynamics is sometimes sufficiently accurate. In case of significant nuclear delocalization, which may occur in systems with light atoms, strong bonds, or for shallow potential energy landscapes81, nuclear quantum effects (NQEs) must be included as well. Even then, methods like path-integral MD establish a one-to-one correspondence between the properties of a quantum object and a classical system with an extended phase space, eliminating the need to solve the nuclear SE82, 83, 84.

At each time step of a dynamics simulation, the forces $\mathbf{F}_{i}$ acting on each atom $i$ must be known so that the equations of motion can be integrated numerically (e.g. using the Verlet algorithm85). They can be derived from the PES by using the relation $\mathbf{F}_{i}=-\nabla_{\mathbf{r}_{i}}E$, i.e. the forces are the negative gradient of the potential energy $E$ with respect to the atomic positions $\mathbf{r}_{i}$ (see also Section 2.1.2). Forces can also be used to perform geometry optimizations, e.g. to find special configurations of atoms which correspond to critical points on the PES. For example, the height of a reaction barrier can be computed from the energy difference between the saddle point (transition state) and either of the two minima (equilibrium structures) which are connected by it.

Although the BO approximation simplifies the SE, even approximate solutions can be computationally demanding, in particular for large systems containing many degrees of freedom. Thus, it is often unfeasible to derive ab initio energies and forces for each time step of an MD simulation. For this reason, analytical functions, i.e. force fields (FFs), are typically used to represent the PES, circumventing the problem of solving an equation altogether. The difficulty is then shifted to finding an appropriate functional form and parametrization of the FF. ML methods automate this demanding and time-consuming process by learning an appropriate function from reference data.

1.2 Invariances of physical systems

Closed physical systems are governed by various conservation laws that describe invariant properties. They are fundamental principles of nature that characterize symmetries that must not be violated. As such, conservation laws provide strong constraints that can be used as guiding principles in search of physically plausible ML models. The basic invariances of molecular systems are directly derived from Noether’s theorem86, which states that each conserved quantity is associated with a differentiable symmetry of the action of a physical system. Typical conserved quantities include the total energy (following from temporal invariance), as well as angular and linear momentum (roto-translational invariance). Energy conservation imposes a particular structure on vector fields in order for them to be valid force fields with corresponding potentials. Namely, forces must be the negative gradient of the potential energy with respect to atomic positions. This relation ensures that when atoms move, they always acquire the same amount of kinetic energy as they lose in potential energy (and vice versa), i.e. the total energy is constant (the work done along closed paths is zero). The conservation of linear and angular momentum implies that the potential energy of a molecule only depends on the relative position of its atoms to each other, i.e. it does not change with rigid rotations and/or translations. Another invariance (not derived from Noether’s theorem) follows from the fact that, from the perspective of the electrons, atoms with the same nuclear charge appear identical to each other. They can thus be exchanged without affecting the energy and forces, which makes the PES symmetric with respect to permutations of some of its arguments. To ensure physically meaningful predictions, ML-FFs must be made invariant under the same transformations as the true PES by introducing appropriate constraints.

2 Machine learning foundations

A question that frequently arises for researchers new to the field of ML concerns the difference of ML modeling to plain interpolation in the noise free regression case. After all, the Shannon sampling theorem gives bounds for the number of “training samples” needed to reconstruct a band-limited signal exactly87. Since the regression tasks considered in this review use ab initio data as reference, they can be considered practically noise-free. Furthermore, PESs are usually smooth, i.e. there is a well-defined frequency cutoff in the spectrum of this “signal”. Thus, both requirements for Shannon interpolation are satisfied and it should in principle be possible to reconstruct FFs via interpolation of the training samples without error, provided there are enough of them.

This is where ML diverges from signal interpolation theory. In practice, there is often not enough data available to fully capture all the necessary information for a perfect reconstruction. In that case, the goal of ML methods is not to recover the training data, but rather to estimate the true process with its underlying regularities that also describes all new and unseen data – this is often denoted as generalization. The key to generalization is selecting a model based on the well known principle of Occam’s razor, i.e. the notion that simpler hypotheses are more likely to be correct88. The capacity of the model can be controlled using the bias–variance trade-off89 (a compromise between expressiveness and complexity) and is practically done by exercising model selection techniques (see Section 2.2.3) such as cross-validation that leave out part of the data from the ML training process and use it later to obtain a valid estimate of the generalization error.90, 30 The reason why regularization is often needed is that ML algorithms are universal approximators, i.e. they can approximate any continuous function on a closed interval arbitrarily close. Since for a finite amount of reference data infinitely many such functions are thinkable, a regularization mechanism is often needed to select a preferably simple function from the vast space of possibilities.

ML methods typically rely on the fact that nonlinear problems, such as predicting energy from nuclear positions, can be “linearized” by mapping the input to a (often higher-dimensional) “feature space” (see Fig. 4).91, 92, 93, 90 Note that such feature spaces are explicitly constructed for kernel-based learning methods (see Section 2.2.1) or learned respectively for deep learning models 94 (see Section 2.2.2). Kernel-based methods achieve this by taking advantage of the so-called kernel trick95, 96, 92, 97, 98 which allows implicitly operating in a high-dimensional feature space without explicitly performing any computation there. In contrast, artificial neural networks (NNs) decompose a complex non-linear function into a composition of linear transformations with learnable parameters connected by nonlinear activation functions. With increasingly many of such nonlinear transformations organized in “layers” (deep NNs), it is possible to efficiently learn highly complex feature spaces.

While NNs tend to require more training data to reach the same accuracy as kernel methods (see Fig. 5),99 they typically scale better to larger data sets. In general, neither method is strictly superior over the other100 and both have advantages and disadvantages that must be weighed against each other for a specific application. Recently, it has even been discovered that in the limit of infinitely wide layers, deep NNs are equivalent to kernel methods, which shifts the main differentiating factor between both methodologies to how they are constructed and trained101, 102 and makes deep NNs accessible to kernel-based analysis methods 103, 104.

In the following, kernel methods and neural networks are described in more detail to highlight the most important properties that differentiate both methodologies.

The representer theorem states that the functional relation

where $\epsilon$ denotes measurement noise, can be optimally approximated as a linear combination

where hyperparameter $d$ is the degree of the polynomial and $\langle\cdot,\cdot\rangle$ is the dot product, or the Gaussian kernel given by

with hyperparameter $\gamma$ controlling its width/scale and $\lVert\cdot\rVert$ denoting the $L^{2}$-norm (see Refs. 90, 30, 125, 97 for more examples of kernel functions).

The structure and number of dimensions of the associated Hilbert space $\mathcal{H}$ depends on the choice of $K(\mathbf{x},\mathbf{x^{\prime}})$ and dimension of the inputs $\mathbf{x}$ and $\mathbf{x^{\prime}}$. As an example, consider the polynomial kernel (Eq. 3) with degree $d=2$ and two-dimensional inputs. The corresponding homogeneous $(c=0)$ polynomial mapping is given by $\phi:(x_{1},x_{2})\mapsto(x_{1}^{2},\sqrt{2}x_{1}x_{2},x_{2}^{2})$, so the associated $\mathcal{H}$ is three-dimensional. While in this case, it is still possible to compute $\phi$ and evaluate the inner product of two points $\phi(\mathbf{x})$ and $\phi(\mathbf{x^{\prime}})$ explicitly, the advantage of using kernels becomes apparent when the Gaussian kernel (Eq. 4) is considered. Rewriting Eq. 4 as

and expanding the third factor in a Taylor series $e^{2\gamma\langle\mathbf{x},\mathbf{x^{\prime}}\rangle}=\sum_{d=0}^{\infty}\frac{1}{d!}\left(2\gamma\langle\mathbf{x},\mathbf{x^{\prime}}\rangle\right)^{d}$ reveals that the Gaussian kernel is equivalent to an infinite sum over (scaled) polynomial kernels (see Eq. 3) and the associated $\mathcal{H}$ is $\infty$-dimensional. Fortunately, by using the kernel function $K(\mathbf{x},\mathbf{x^{\prime}})$, it is possible to operate in $\mathcal{H}$ implicitly and evaluate $\hat{f}(\mathbf{x})$ (Eq. 2) without computing the mapping $\phi$. This is often referred to as the kernel trick95, 96, 92, 90, 125, 97.

It remains the question how the coefficients $\alpha_{i}$ in Eq. 2 are determined. One way to do so is by adopting a Bayesian, or probabilistic, point of view.126, 127 Here, it is assumed that the reference data $\{(y_{i};\mathbf{x}_{i})\}_{i=1}^{M}$ are generated by a Gaussian process (GP), i.e. drawn from a multivariate Gaussian distribution. For simplicity, it can be assumed that this distribution has a mean of zero, as other values can be generated by simply adding a constant term. Further, the possibility that the reference data might be contaminated by noise (for example due to uncertainties in measuring $y_{i}$) is accounted for explicitly. Typically, Gaussian noise is assumed, i.e.

where $\lambda$ is the variance of the normally distributed noise $\mathcal{N}$. In the GP picture, the choice of $K(\mathbf{x},\mathbf{x^{\prime}})$ expresses an assumption about the underlying function class. For example, choosing the Gaussian kernel implies that $f(\mathbf{x})$ does not change drastically over a length scale controlled by $\gamma$ (see Eq. 4). As such, a particular kernel function $K$ corresponds to an implicit regularization, i.e. an assumption about the underlying smoothness properties of the function to be estimated.128 The challenge lies in finding a kernel that represents the structure in the data that is being modeled as good as possible.128, 103 Many kernels are able to approximate continuous functions on a compact subset arbitrarily well,128, 129 but a strong prior has the advantage of restricting the hypothesis space, which drastically improves the convergence of the learning task with respect to the available training data.130

Under these conditions, it is now possible to rigorously answer the question “given the data $\mathbf{y}=\left[y_{1}\cdots y_{M}\right]^{\top}$, how likely is it to observe the value $y_{*}$ for input $\mathbf{x}_{*}$?” As $y_{*}$ is generated by the same GP as the reference data, the conditional probability $p(y_{*}|\mathbf{y})$ can be expressed as

where $\mathbf{I}_{M}$ is the identity matrix of size $M$, $\mathbf{K}$ is the $M\times M$ kernel matrix131, 90 with entries $K_{ij}=K(\mathbf{x}_{i},\mathbf{x}_{j})$ and $\mathbf{K}_{*}=\left[K(\mathbf{x}_{*},\mathbf{x}_{1})\cdots K(\mathbf{x}_{*},\mathbf{x}_{M})\right]$. In other words, Eq. 7 expresses a probability distribution over possible predictions, where its mean value

is the most likely estimate for $y_{*}$ (given the reference data) and its variance

provides information about how strongly other likely predictions vary from the mean. Note that while Eq. 9 can be used as uncertainty estimate for a particular prediction, it should not be confused with error bars. The optimal coefficients $\boldsymbol{\alpha}=\left[\alpha_{1}\cdots\alpha_{M}\right]^{\top}$ in Eq. 2 are thus given by

or simply $\boldsymbol{\alpha}=\mathbf{K}^{-1}\mathbf{y}$ in the noise-free case ($\lambda=0$). However, even in the absence of noise, it can be beneficial to choose a non-zero $\lambda$ to obtain a regularized solution. The addition of $\lambda>0$ to the diagonal of $\mathbf{K}$ increases numerical stability and has the effect of damping the magnitude of the coefficients, thereby increasing the smoothness of $\hat{f}(\mathbf{x})$. The downside is that the known reference values $y_{i}$ are only approximately reproduced. This, however, also decreases the chance of overfitting and can lead to better generalization, i.e. increased accuracy when predicting unknown values.

Matrix factorization methods like Cholesky decomposition132 are typically used to efficiently solve the linear problem in Eq. 10 in closed form. However, this type of approach scales as $\mathcal{O}(M^{3})$ with the number of reference data and may become problematic for extremely large data sets. Iterative, e.g. gradient-based, solvers can reduce the complexity to $\mathcal{O}(M^{2})$.133 Once the coefficients have been determined, the value $y_{*}$ for an arbitrary input $\mathbf{x}_{*}$ can be estimated according to Eq. 2 with $\mathcal{O}(M)$ complexity (a sum over all $M$ reference data points is required).

However, the approximation above gives rise to an over-determined system with fewer parameters than training points and therefore reduced model capacity. Strictly speaking, the involved matrix does not satisfy the properties of a kernel matrix anymore, as it is neither symmetric nor positive semi-definite. To obtain a kernel matrix that still maintains these properties, the Nyström 134 approximation

It should be mentioned that kernel regression methods are known under different names in the literature of different communities. Due to their relation to GPs, some authors prefer the name Gaussian process regression (GPR). Others favor the term kernel ridge regression (KRR), since determining the coefficients with Eq. 10 corresponds to solving a least-squares objective with $L^{2}$-regularization in the kernel feature space $\phi$ and is similar to ordinary ridge regression142. Sometimes, the method is also referred to as reproducing kernel Hilbert space (RKHS) interpolation, since Eq. 2 “interpolates” between known reference values (when coefficients are determined with $\lambda=0$, all known reference values are reproduced exactly). All these methods are formally equivalent and essentially differ only in the manner the relevant equations are derived. There are small philosophical differences, however: For example, in the KRR and RKHS pictures, $\lambda$ in Eq. 10 is a regularization hyperparameter that has to be introduced ad hoc, whereas in the GPR picture, $\lambda$ is directly related to the Gaussian noise in Eq. 6. The expansion coefficients obtained from Eq. 10 can change drastically depending on the choice of $\lambda$, so this is an important detail. Further, while Eq. 9 can be used to compute uncertainty estimates for all kernel regression methods, the GPR picture allows to relate it directly to the variance of a Gaussian process.

The most important concepts discussed in this section are summarized visually in Fig 6.

2.2 Artificial neural networks

Provided that the number of dimensions of the “hidden layer” $\mathbf{h}$ is large enough, this arrangement can approximate any mapping between inputs $\mathbf{x}$ and outputs $\mathbf{y}$ to arbitrary precision, i.e. it is a general function approximator150, 151.

In theory, shallow NNs as shown above are sufficient to approximate any functional relationship.151 However, deep NNs with multiple hidden layers are often superior and were shown to be more parameter-efficient.152, 153, 154, 155 To construct a deep NN, $L$ hidden layers are combined sequentially

mapping the input $\mathbf{x}$ to several intermediate feature representations $\mathbf{h}_{l}$, until the output $\mathbf{y}$ is obtained by a linear regression on the features $\mathbf{h}_{L}$ in the final layer. For PES construction, typically, the NN maps a representation of chemical structure $\mathbf{x}$ to a one-dimensional output representing the energy. Contrary to the coefficients $\boldsymbol{\alpha}$ in kernel methods (see Eq. 10), the parameters $\{\mathbf{W}_{l},\mathbf{b}_{l}\}_{l=1}^{L+1}$ of an NN cannot be fitted in closed form. Instead, they are initialized randomly and optimized (usually using a variant of stochastic gradient descent) to minimize a loss function that measures the discrepancy between the output of the NN and the reference data.156 A common choice is the mean squared error (MSE), which is also used in kernel methods. During training, the loss and its gradient are estimated from randomly drawn batches of training data, making each step independent of the number of training data $M$. On the other hand, finding the coefficients for kernel methods scales as $\mathcal{O}(M^{3})$ due to the need of inverting the $M\times M$ kernel matrix. Evaluating an NN according to Eq. 16 for a single input $\mathbf{x}$ scales linearly with respect to the number of model parameters. The same is true for kernel methods, but here the number of model parameters is tied to the number of reference data $M$ used for training the model (see Eq. 2), which means that evaluating kernel methods scales $\mathcal{O}(M)$. As the evaluation cost of NNs is independent of $M$ and only depends on the chosen architecture, they are typically the method of choice for learning large datasets. A schematic overview of the mathematical concepts behind NNs is given in Fig. 8.

2.3 Model selection: How to choose hyperparameters

In addition to the parameters that are determined when learning an ML model for a given dataset, e.g. the weights $\mathbf{W}$ and biases $\mathbf{b}$ in NNs or the regression coefficients $\boldsymbol{\alpha}$ in kernel methods, many models contain hyperparameters that need to be chosen before training. They allow to tune a given model to the prior beliefs about the dataset/underlying physics and thus play a significant role in how a model generalizes to different data patterns. Two types of hyperparameters can be distinguished: The first kind influences the composition of the model itself, such as the type of kernel or the NN architecture, whereas the second kind affects the training procedure and thus the final parameters of the trained model. Examples for hyperparameters are the width (number of neurons per layer) and depth (number of hidden layers) of an NN, the kernel width $\gamma$ (see Eq. 4), or the strength of regularization terms (e.g. $\lambda$ in Eq. 10).

The range of valid values is strongly dependent on the hyperparameter in question. For example, certain hyperparameters might need to be selected from the positive real numbers (e.g. $\gamma$ and $\lambda$, see above), while others are restricted to positive integers or have interdependencies (such as depth and width of an NN). This is why hyperparameters are often optimized with primitive exhaustive search schemes like grid search or random search in combination with educated guesses for suitable search ranges, or more sophisticated Bayesian approaches157. Common gradient-based optimization methods can typically not be applied effectively. Fortunately, for many hyperparameters, model performance is fairly robust to small changes and good default values can be determined which work across many different datasets.

Before any hyperparameters may be optimized, a so-called test set must be split off from the available reference data and kept strictly separate. The remainder of the data is further divided into a training and a validation set. This is done because the performance of ML models is not judged by how well they predict the data they were trained on, as it is often possible to achieve arbitrarily small errors in this setting. Instead, the generalization error, i.e. how well the model is able to predict unseen data, is taken as indicator for the quality of a model. For this reason, for every trial combination of hyperparameters, a model is trained on the training data and its performance measured on the validation set to estimate the generalization error. Finally, the best performing model is selected. To get better statistics for estimates of the generalization error, instead of splitting the remaining data (reference data excluding test set) into just two parts, it is also possible to divide it into $k$ parts (or folds). Then, a total of $k$ models is trained, each using $k-1$ folds as the training set and the last fold as validation set. This method is known as $k$-fold cross validation158, 30.

As the validation data influences model selection (even though it is not used directly in the training process), the validation error may give too optimistic estimates and is no reliable way to judge the true generalization error of the final model. A more realistic value can be obtained by evaluating the model on the held-out test set, which has neither direct nor indirect influence on model selection. To not invalidate this estimate, it is crucial not to further tweak any parameters or hyperparameters in response to test set performance. More details on how to construct ML models (including the selection of hyperparameters and the importance of keeping an independent test set) can be found in Section 3. The model selection process is summarized in Fig. 9.

3 Combining machine learning and chemistry

The need for ML methods often arises from the lack of theory to describe a desired mapping between input and output. A classical example for this is image classification: It is not clear how to distinguish between pictures of different objects, as it is unfeasible to formulate millions of rules by hand to solve this task. Instead, the best results are currently achieved by learning statistical image characteristics from hundreds of thousands of examples that were extracted from a large dataset representing a particular object class. From that, the classifier learns to estimate the distribution inherent to the data in terms of feature extractors with learned parameters like convolution filters that reflect different scales of the image statistics.101, 98, 94 This working principle represents the best approach known to date to tackle this particular challenge in the field of computer vision.

On the other hand, the benchmark for solving molecular problems is set by rigorous physical theory that provides essentially exact descriptions of the relationships of interest. While the introduction of approximations to exact theories is common practice and essential to reduce their complexity to a workable level, those simplifications are always physical or mathematical in nature. This way, the generality of the theory is only minimally compromised, albeit with the inevitable consequence of a reduction in predictive power. In contrast, statistical methods can be essentially exact, but only in a potentially very narrow regime of applicability. Thus, a main role of ML algorithms in the chemical sciences has been to shortcut some of the computational complexity of exact methods by means of empirical inference, as opposed to providing some mapping between input and output at all (as is the case for image classification). Notably, recent developments could show that machine learning can provide novel insight beyond providing efficient shortcuts of complex physical computations.159, 62, 59, 105, 55, 33, 70, 160

Force field construction poses unique challenges that are absent from traditional ML application domains, as much more stringent demands on accuracy are placed on ML approaches that attempt to offer practical alternatives to established methods. Additionally, considerable computational cost is associated with the generation of high-level ab initio training data, with the consequence that practically obtainable datasets with sufficiently high quality are typically not very large. This is in stark contrast with the abundance of data in traditional ML application domains, such as computer vision, natural language processing etc. The challenge in chemistry, however, is to retain the generality, generalization ability and versatility of ML methods, while making them accurate, data-efficient, transferable, and scalable.

To increase data efficiency and accuracy, ML-FFs can (and should) exploit the invariances of physical systems (see Section 2.1.2), which provide additional information in ways that are not directly available for other ML problems. Those invariances can be used to reduce the function space from which the model is selected, in this manner effectively reducing the degrees of freedom for learning,69, 161 i.e. making the learning problem easier and thus also solvable with a fraction of data. As ML algorithms are universal approximators with virtually no inherent flexibility restrictions, it is important that physically meaningful solutions are obtained. In the following, important physical constraints of such solutions and possible ways of their realization are discussed in detail. Furthermore, existing kernel-based methods and neural network architectures tailored for the construction of FFs and how they implement these physical constraints in practice are described.

A necessary requirement for ML-FFs is that, in the absence of external forces, the total energy of a chemical system is conserved (see Section 2.1.2). When the potential energy is predicted by any differentiable method and forces derived from its gradient, they will be conservative by construction. However, when forces are predicted directly, this is generally not true, which makes deriving energies from force samples slightly more complicated. The main challenge to overcome is that not every vector field is necessarily a valid gradient field. Therefore, the learning problem cannot simply be cast in terms of a standard multiple output regression task, where the output variables are modeled without enforcing explicit correlations.

A big advantage of predicting forces directly is that they are true quantum-mechanical observables within the BO approximation by virtue of the Hellmann-Feynman theorem162, 163, i.e. they can be calculated analytically and therefore at a relatively low additional cost when generating ab initio reference data. As a rough guideline, the computational overhead for analytic forces scales with a factor of only around $\sim$1–7 on top of the energy calculation.164. In contrast, at least $3N+1$ energy evaluations would be necessary for a numerical approximation of the forces by using finite differences. For example, at the PBE0/DFT (density functional theory with the Perdew-Burke-Ernzerhof hybrid functional) level of theory165, calculating energy and analytical forces for an ethanol molecule takes only $\sim$1.5 times as long as calculating just the energy (the exact value is implementation-dependent), whereas for numerical gradients, a factor of at least $\sim$10 would be expected.

As forces provide additional information about how the energy changes when an atom is moved, they offer an efficient way to sample the PES, which is why it is desirable to formulate ML models that can make direct use of them in the training process. Another benefit of a direct reconstruction of the forces is that it avoids the amplification of estimation errors due to the derivative operator that would otherwise be applied to the PES reconstruction (see Fig. 10).105, 70, 58

A crucial requirement for ML-FFs is the rotational and translational invariance of the potential energy, i.e. $E\left(\mathbf{R}\right)=E\left(\mathcal{R}\mathbf{R}+\mathcal{T}\right)$, where $\mathcal{R}$ and $\mathcal{T}$ are rigid rotations and translations and $\mathbf{R}$ are the Cartesian coordinates of the atoms. As long as the representation $\mathbf{x}(\mathbf{R})$ of chemical structure chosen as input for the ML model itself is roto-translationally invariant, ML-FFs inherit its desired properties and even the gradients will automatically behave in the correct equivariant way due to the outer derivative $\partial(\mathbf{x}(\mathcal{R}\mathbf{R}+\mathcal{T})/\partial\mathbf{R}=\mathcal{R}\partial\mathbf{x}(\mathbf{R})/\partial\mathbf{R}$. One example of appropriate features to construct a representation $\mathbf{x}$ with the desired properties are pairwise distances. For a system with $N$ atoms, there are $\binom{N}{2}$ different pairwise distances, which results in reasonably sized feature sets for systems with a few dozen atoms. Apart from very few pathological cases, this representation is complete, in the sense that any possible configuration of the system can be described exactly and uniquely117. However, while pairwise distances serve as an efficient parametrization of some geometry distortions like bond stretching, they are relatively inefficient in describing others, e.g. rotations of functional groups. In the latter case, many distances are affected even for slight angular changes, which can pose a challenge when trying to learn the geometry-energy mapping. Complex transition paths or reaction coordinates are often better described in terms of bond and torsion angles in addition to pairwise distances. The problem is that the number of these features grows rather quickly, with $\binom{N}{3}$ and $\binom{N}{4}$, respectively. At that rate, the size of the feature set quickly becomes a bottleneck, resulting in models that are slow to train and evaluate. While an expert choice of relevant angles would circumvent this issue, it reduces some of the “data-driven” flexibility that ML models are typically appreciated for. Note that models without roto-translational invariance are practically unusable, as they will start to generate spurious linear and/or angular momentum during dynamics simulations.

In the BO approximation, the potential energy of a chemical system only depends on the charges and positions of the nuclei. As a consequence, the PES is symmetric under permutation of atoms with the same nuclear charge. However, symmetric regions are not necessarily sampled in an unbiased way during MD simulations (see Fig. 11). Consequently, ML-FFs that are not constrained to treat all symmetries equivalently may (due to the uneven sampling) predict different results when permuting atoms.

While it is in principle possible to arrive at a ML-FF that is symmetric with respect to permutations of same-species atoms indirectly via data augmentation168, 29 or by simply using datasets that naturally include all relevant symmetric configurations in an unbiased way, there are obvious scaling issues with this approach. It is much more efficient to impose the right constraints onto the functional form of the ML-FF, such that all relevant symmetric variants of a certain atomic configuration appear equivalent automatically. Such symmetric functions can be constructed in various ways, each of which has advantages and disadvantages.

Assignment-based approaches do not symmetrize the ML-FF per se, but instead aim to normalize its input, such that all symmetric variants of a configuration are mapped to the same internal representation. In its most basic realization, this assignment is done heuristically, i.e. by using inexact, but computationally cheap criteria. Examples for this approach are the Coulomb matrix28 or the Bag-of-Bonds31 descriptors, that use simple sorting schemes for that purpose. Histograms169, 107 and some density-based117, 170, 171 approaches follow that same principle, although not explicitly. All of those schemes have in common that they compare the features in aggregate as opposed to individually. A disadvantage is that dissimilar features are likely to be compared to each other or treated as the same, which limits the accuracy of the prediction. Such weak assignments are better suited for datasets with diverse conformations rather than gathered from MD trajectories that contain many similar geometries. In the latter case, the assignment of features might change as the geometry evolves, which would lead to discontinuities in the prediction and would effectively be treated by the ML model as noise (see $\epsilon$ in Eq. 1).

Arguably the most universal way of imposing symmetry, especially if the functional form of the model is already given, is via invariant integration over the relevant symmetry group $f_{\text{sym}}(\mathbf{x})=\int_{\pi\in\mathcal{S}}f(\mathbf{P}_{\pi}\mathbf{x})$. Typically, $\mathcal{S}$ would be the permutation group and $\mathbf{P}_{\pi}$ the corresponding permutation matrix that transforms each vector of atom positions $\mathbf{x}$. Some approaches 117, 175, 176 avoid this implicit ordering of atoms in $\mathbf{x}$ by adopting a three-dimensional density representation of the molecular geometry defined by the atom positions, albeit at the cost of losing rotational invariance, which then must be recovered by integration. Invariant integration gives rise to functional forms that are truly symmetric and do not require any pre- or post-processing of the in- and output data. A significant disadvantage is however, that the cardinality of even basic symmetry groups is exceedingly high, which affects both training and prediction times.

This combinatorial challenge can be solved by limiting the invariant integral to the physical point group and fluxional symmetries that actually occur in the training dataset. Such a sub-group of meaningful symmetries can be automatically recovered and is often rather small164. For example, each of the molecules benzene, toluene and azobenzene have only 12 physically relevant symmetries, whereas their full symmetric groups have orders $6!6!$, $7!8!$ and $12!10!2!$ symmetries respectively.

3.2 (Symmetric) Gradient Domain Machine Learning ((s)GDML)

Gradient domain machine learning (GDML) is a kernel-based method introduced as a data efficient way to obtain accurate reconstructions of flexible molecular force fields from small reference datasets of high-level ab initio calculations105. Contrary to most other ML-FFs, instead of predicting the energy and obtaining forces by derivation with respect to nuclear coordinates, GDML predicts the forces directly. As mentioned in Section 2.3.1, forces obtained in this way may violate energy conservation. To ensure conservative forces, the key idea is to use a kernel $\mathbf{K}\left(\mathbf{x},\mathbf{x}^{\prime}\right)=\nabla_{\mathbf{x}}K_{E}\left(\mathbf{x},\mathbf{x}^{\prime}\right)\nabla_{\mathbf{x}^{\prime}}^{\top}$ that models the forces $\mathbf{F}$ as a transformation of an unknown potential energy surface $E$ such that

A limitation of the GDML method is that the structural descriptor $\mathbf{x}$ is not permutationally invariant, because the values of its entries (inverse pairwise distances) change when atoms are re-ordered. An extension of the original approach, sGDML69, 164 (symmetric GDML), additionally incorporates all relevant rigid space group symmetries, as well as dynamic non-rigid symmetries of the system at hand into the kernel, to further improve its efficiency and ensure permutational invariance. Usually, the identification of symmetries requires chemical and physical intuition about the system at hand, which is impractical in an ML setting. Here, however, a data-driven multi-partite matching approach is employed to automatically recover permutations of atoms that appear within the training set164. A matching process finds permutation matrices $\mathbf{P}$ that realize the assignment between adjacency matrices $(\mathbf{A})_{ij}=\|\vec{r}_{i}-\vec{r}_{j}\|$ of molecular graph pairs $G$ and $H$ in different energy states

and thus between symmetric transformations177. The resulting approximate local pairwise assignments are subsequently globally synchronized using transitivity as the consistency criterion174 to eliminate impossible assignments. By limiting this search to the training set, combinatorially feasible, but physically irrelevant permutations $\tau$ are automatically excluded (ones that are inaccessible without crossing impassable energy barriers). Such hard symmetry constraints (derived from the training set) greatly reduce the intrinsic complexity of the learning problem without biasing the estimator, since no additional approximations are introduced.

3.3 Gaussian approximation potentials (GAPs)

Gaussian approximation potentials (GAPs)178 were originally developed for materials such as bulk crystals, but were later also applied to molecules179. They scale linearly with the number of atoms of a system and can accommodate for periodic boundary conditions. Similar to high-dimensional neural network potentials114 (see Section 2.3.4), GAPs decompose each system into atom-centered environments $i$ such that its energy can be written as the sum of atomic contributions

with $\mathbf{r}_{ij}=\mathbf{r}_{j}-\mathbf{r}_{i}$ and $\mathbf{r}_{i}$ being the position of atom $i$. A smooth cutoff function is applied to the pairwise distances $\lVert\mathbf{r}_{ij}\rVert$ to ensure that the contributions $E_{i}$ are local and no discontinuities are introduced when atoms enter or leave the cutoff radius. Even though such a decomposition is inherently non-unique and no labels for atom-wise energies are available in the reference data, they can still be approximated by a Gaussian process: The sum over atomic environments can be moved into the kernel function, yielding a kernel for systems $\mathbf{x}$ and $\mathbf{x}^{\prime}$ with $N$ and $N^{\prime}$ atoms, respectively:

Thus, reference energies for the whole system are sufficient for the model to learn a suitable energy decomposition into atomic environments.

Several descriptors and kernels for GAPs have been developed based on a local “atomic density” $\rho(\mathbf{r})=\sum_{j}\delta(\mathbf{r}-\mathbf{r_{j}})$. Initially, Bartók et al. 178 proposed to employ local atomic coordinates projected onto a 4D hyper sphere. Since this projection can represent the volume of a 3D sphere, the introduction of an additional radial basis can be avoided. To achieve rotational invariance, the bispectrum of 4D spherical harmonics of these coordinates was used as a descriptor. Alternatively, the SOAP (smooth overlap of atomic positions) kernel117 is defined as the integral over rotations $\mathcal{R}$ of atomic densities

Given smoothed local densities $\rho(\mathbf{r})=\sum_{j}\exp(-\gamma\|\mathbf{r}-\mathbf{r_{j}}\|^{2})$, it has been shown that the SOAP kernel is equivalent to the linear kernel over the SO(3) power spectrum and bispectrum for $n=2$ and $n=3$, respectively117. Both approaches are invariant to permutation of neighboring atoms as well as the rotation of the local environment. Further representations include best matches of the atomic densities over rotations176 and kernels for symmetry-adapted prediction of tensorial properties180, 181.

3.4 Neural Network Potentials

The first neural network potentials (NNPs) used a set of internal coordinates, e.g. distances and angles, as structural representation to model the PES182, 183, 184, 185, 186. While being roto-translationally invariant, internal coordinates impose an arbitrary order on the atoms and are thus not reflecting the equivalence of permuted inputs. As a result, the NNP might assign different energies to symmetrically equivalent structures. Beyond that, the number of atoms influences the dimensionality of the input $\mathbf{x}$, limiting the applicability of the PES to chemical systems of the same size. Decomposing the energy prediction in the spirit of a many-body expansion circumvents these issues187, 188, 189, however, it scales unfavorably with system size and number of chemical species, because each term in the many-body expansion has to be modeled by a separate NN.

Behler and Parrinello 114 were the first to propose so-called high-dimensional neural networks potentials (HDNNPs), where the total energy of a chemical system is expressed as a sum over atomic contributions $E=\sum_{i}E_{i}$, predicted by the same NN (or one for each element). The underlying assumption is that the energetic contribution $E_{i}$ of each atom depends mainly on its local chemical environment. As all atoms of the same type are treated identically and summation is commutative, the output does not change when the input is permuted. Due to the decomposition into atomic contributions, systems with varying numbers of atoms can be predicted by the same NNP. In principle, this framework also enables transferability between system sizes, e.g. a model can be trained on small systems, but applied to predict energies and forces for larger systems. However, this requires sufficient sampling of the local environments to remove spurious correlations caused by the training data distribution, as well as corrections for long-range effects.

The introduction of HDNNPs inspired many NN architectures that can be broadly categorized into two types. Descriptor-based NNPs116, 190, 191, 192 rely on fixed rules to encode the environment of an atom in a vector $\mathbf{x}$, which is then used as input for an ordinary feed-forward NN (see Eq. 16). These architectures include many variants of the original Behler-Parrinello network, such as ANI193 and TensorMol194. On the other hand, end-to-end NNPs195, 196, 159, 197, take nuclear charges and Cartesian coordinates as input and learn a suitable representation from the data.

Many end-to-end NNPs have been inspired by the graph neural network by Scarselli et al. 198 and were later collectively cast as message-passing neural networks (MPNNs)197. In this type of model, molecules are regarded as undirected graphs, where atoms are represented by nodes and interactions between them as edges. By exchanging information between nodes along edges (message-passing), complex chemical interactions can be modeled. A prominent example is the Deep Tensor Neural Network (DTNN)159. Since its introduction, this approach has been refined to create new architectures, such as SchNet109, 199, HIP-NN115 or PhysNet108. End-to-end NNPs that do not directly fall into the category of MPNNs are covariant compositional networks that are able to employ features of higher angular momentum200, 201, 202 as well as models using a pseudo-density as input171.

Because no fixed rule is used to construct descriptors, end-to-end NNPs are able to automatically adapt the environment representations $\mathbf{x}$ to the reference data (in contrast to the descriptor-based variant). However, as long as $\mathbf{x}$ is invariant with respect to translation, rotation, and permutation of symmetry equivalent atoms, both types of NNPs adhere to all physical constraints outlined in Section 2.3.1. NNPs are commonly used to predict energies, while conservative forces are obtained by derivation. Despite being energy-based, it is still possible to incorporate information from ab initio forces by including them in the loss term optimized during training. At this point, it should be noted that the requirement for continuously differentiable models excludes the use of certain activation functions, for example the popular ReLU activation203, when constructing ML-FFs based on neural networks. To avoid discontinuities in the forces, activation functions used for NNPs must always be smooth.

The first descriptor-based NNP introduced by Behler and Parrinello 114 uses atom-centered symmetry functions (ACSFs)116 consisting of two-body terms

to encode information about the chemical environment of each atom $i$. Here, $r_{ij}$ is the distance between atoms $i$ and $j$, $\theta_{ijk}$ the angle spanned by atoms $i$, $j$ and $k$ centered around $i$, and the summations run over all atoms within a cutoff distance $r_{\rm cut}$. As the atom order is irrelevant for the values of $G^{2}_{i}$ and $G^{3}_{i}$ and only internal coordinates are used to calculate them, all physical invariants are satisfied. A cutoff function such as

ensures that $G^{2}_{i}$ and $G^{3}_{i}$ vary smoothly when atoms enter or leave the cutoff sphere and the parameters $\eta$, $r_{s}$, $\zeta$, and $\lambda(=\pm 1)$ determine to which distances, or combinations of angles and distances, the ACSFs are most sensitive. When sufficiently many $G^{2}_{i}$ and $G^{3}_{i}$ with different parameters are combined and stored in a vector $\mathbf{x}_{i}$, they form a “fingerprint” of the local environment of atom $i$. This environment descriptor is then used as input for a neural network for predicting the energy contributions $E_{i}$ of atoms $i$ and the total energy $E=\sum_{i}E_{i}$ is obtained by summation.

Since the ACSFs only use geometric information, they work best for systems containing only atoms of one element, for example crystalline silicon114. To describe multi-component systems, typically, the symmetry functions are duplicated for each combination of elements and separate NNs are used to predict the energy contributions for atoms of the same type204. Since the combinatorial explosion can lead to a large number of ACSFs for systems containing many different elements, an alternative is to modify the ACSFs with element-dependent weighting functions205. Most descriptor-based NNPs, such as ANI193 and TensorMol194, use variations of Eqs. 22 and 23 (sometimes allowing parameters of ACSFs to be optimized during training) to construct the environment descriptors $\mathbf{x}_{i}$. Different ways to encode the structural information are possible, for example using three-dimensional Zernike functions190, or the coefficients of a spherical harmonics expansion192, but the general principle remains the same. Also, while most descriptor-based NNPs use separate parametrizations for different elements, it is also possible to use a single NN to predict all atomic energy contributions192. The common feature for all variations of this approach is that the functional form of the environment descriptor is predetermined and manually designed.

A potential drawback of the previously introduced ACSFs is that they must be chosen by an expert before training the neural network. If the choice of symmetry functions is poor, for example when the resulting descriptor is (nearly) identical for two very different structures, the expressive power of the neural network and the achievable accuracy are limited a priori. Additionally, a growing number of input dimensions can quickly become computationally expensive, both for calculating the descriptors and for evaluating the NN. This is especially the case when modeling multi-component systems, where commonly orthogonality is assumed between different elements (which increases the number of input dimensions) or the descriptors are simply weighted by an element-dependent factor (which may limit the structural resolution of the descriptor).

In contrast, end-to-end NNPs directly take atomic types and positions as inputs to learn suitable representations from the reference data. Similar to descriptor-based NNPs, many end-to-end NNPs obtain the total energy $E$ as a sum of atomic contributions $E_{i}$. However, those are predicted from learned features $\mathbf{x}_{i}$ encoding information about the local chemical environment of each atom $i$. This allows them to adapt the features based on the size and distribution of the training set as well as the chemical property of interest during the training process. The idea is to learn a mapping to a high-dimensional feature space, so that structurally (and energetically) similar atomic environments lie close together and dissimilar ones far apart.

Within the deep tensor neural network framework159, this is achieved by iteratively refining the atomic features $\mathbf{x}_{i}$ based on neighboring atoms. The features are initialized to $\mathbf{x}_{i}^{0}=\mathbf{e}_{Z_{i}}$, where $\mathbf{e}_{Z}$ are learnable element-dependent representations that are updated for $T\in$ steps. This procedure is inspired by diffusion graph kernels206 as well as the graph neural network model by Scarselli et al. 198. Many end-to-end networks have adapted this approach which can be written in general as

where the summation runs over all atoms within a distance $r_{\rm cut}$ and a cutoff function $f_{\rm cut}$ ensures smooth behavior when atoms cross the cutoff. Here, the “atom-wise” function $\mathcal{G}^{t}$ is used to refine the atomic features after they have been updated with information from neighboring atoms through the interaction-function $\mathcal{F}^{t}$. Usually, the interatomic distance $r_{ij}$ is not used directly as input to $\mathcal{F}^{t}$, but expanded in a set of uniformly spaced radial basis functions159, 109, 108 to form a vectorial input $\mathbf{g}(r_{ij})$. Both $\mathcal{F}^{t}$ and $\mathcal{G}^{t}$ functions are NNs with the specific implementations varying between different end-to-end NNP architectures. As only pair-wise distances are used and the order of atoms is irrelevant due to the commutative property of summation, the features $\mathbf{x}_{i}$ obtained by Eq. 25 are automatically roto-translationally and permutationally invariant (and thus also the energy predictions).

Gilmer et al. 197 have cast graph networks of this structure as message-passing neural networks and proposed a variant that uses a set2set decoder207 instead of a sum over energy contributions to achieve permutational invariance of the energy. SchNet109 takes an alternative view of the problem and models interactions between atoms with convolutions. The convolution filters need to be continuous (to have smooth predictions) but are evaluated at finite points, i.e. the positions of neighboring atoms. To ensure rotational invariance, only radial convolution filters are used, leading again to an interaction function that is a special case of Eq. 25.

While the previously introduced approaches aim to learn as much as possible from the reference data, several models have been proposed to better exploit chemical domain knowledge. The hierarchical interacting particle neural network (HIP-NN)115 obtains the prediction as a sum over atom-wise contributions $E_{i}^{t}$ that are predicted after every update step $t$. A regularizer penalizes larger energy contributions in deeper layers, i.e. enforcing a declining, hierarchical prediction of the energy. PhysNet108 modified the energy function to include explicit terms for electrostatic and dispersion interactions,

Best practices and Pitfalls

A number of careful modeling steps are necessary to construct an ML-FF for a particular problem of interest (Fig. 13). Even before starting this process, some forethought is appropriate due to certain limitations of ab initio methods themselves. This section gives an overview about all steps necessary to construct an ML-FF from scratch and highlights possible “pitfalls”, i.e. issues that may occur along the way, in particular when the recommended practices are not followed. First, some preliminary considerations, which should be taken before starting with the construction of an ML-FF, are discussed (3.1). Next, basic principles for choosing an appropriate ML method for a specific task are given (3.2). Then, the importance of high quality reference data, different strategies to collect it (3.3), and how the data has to be prepared (3.4), are outlined. This is followed by an overview of how to train an ML model on the collected data (3.5) and guidelines for using the trained ML-FF in a production setting, e.g. for running MD simulations (3.6). Finally, popular software packages for constructing ML-FFs are briefly described and code examples are given (3.7).

Before running any ab initio calculations to collect data for training ML models, it is advisable to think about the limitations of the chosen level of theory itself. The issues discussed here are problem-specific and often not unique to ML-FFs, but PES reconstruction in general. As such, a comprehensive list is not possible, but a few examples are given below.

On the spectrum of quantum chemistry methods, ML-FFs fit into the niche between highly efficient conventional FFs210 and accurate, but computationally expensive ab initio methods.211 Efficiency-wise, they are still inferior to classical FFs, because their functional forms are considerably more complex and thus more expensive to evaluate. Even the fastest ML-FFs are still one to three orders of magnitude slower164, 212, 213. On the other end, ML-FFs are lower bounded by the accuracy of the reference data used for training, which means that the underlying ab initio method will always be at least equally accurate. In practical terms, this means that in order to be useful, ML-FFs need to offer time savings over directly running ab initio calculations and an improved accuracy compared to conventional FFs. For this purpose, the full procedure of data generation, training and inference must be taken into account, as opposed to just regarding inference speed, which will be much quicker than ab initio methods. While this consideration sounds trivial at first, it is still advisable to think about whether constructing an ML-FF really is economical. For example, if the goal is to run just a single short MD trajectory, the question is how much data is necessary for the model to reach the required accuracy. Some models may require several thousands of training points to produce accurate enough predictions, even for fairly small molecules. Then, when factoring in the overall time required for going through the process of creating the ML model, testing it, and running the MD simulation, it might be more efficient to simply run an ab initio MD simulation in the first place. Further, not every ML method is equally applicable or appropriate for all systems due to methodical and/or conceptual constraints. Such limitations are discussed in greater detail in Section 3.2.

Many ab initio methods use a single Slater determinant to express the wave function of a system. The problem with this approach is that different determinants may be dominant in different regions of the PES, leading to a poor description of the wave function if the wrong determinant is chosen. Especially when many calculations are performed for various strongly distorted geometries, for example when a reaction is studied and bonds need to be broken, it may happen that the solution “jumps” discontinuously from one electronic state to another, leading to inconsistent reference data. When an ML model is trained on such a data set, it will try to find a compromise between the inconsistencies and its performance typically be unsatisfying. It is therefore advisable to check for possible multireference effects prior to generating data and, if necessary, switch to a multireference method (for a comprehensive review on multireference methods, see Ref. 214).

The models discussed in this review all assume that energy contributions are local to some degree. This assumption is either introduced explicitly by a cutoff radius, or it enters the model through the use of a specific structural descriptor. For example, by using inverse distances to encode chemical structures for kernel methods (as is done e.g. in GDML, see Section 2.3.2), relative changes between close atoms are weighed more strongly when comparing two conformations. While assuming locality is valid in many practical applications, there exist many cases where this assumption breaks down. An example are extensive conjugated $\pi$-systems, where a rotation around certain bonds might break the favorable interaction between electrons, leading to a “non-local” energy contribution. If such effects exist, an appropriate model should be chosen, for example the cutoff radius may need to be larger than usual or a different structural descriptor must be picked.

2 Choosing an appropriate ML method

Several different variants of ML-FFs have been discussed in Section 2.3 and many more are described in the literature. Although all these methods can be applied to construct ML-FFs for any chemical system, some methods might be more promising than others for certain tasks. For researchers who want to apply ML methods to a specific problem for the first time, the abundance of different models to choose from may be overwhelming and it might be difficult to find an appropriate choice.

In the following, possible applications of ML-FFs are broadly categorized based on simple questions about the task at hand. For each case, advantages and disadvantages of individual models are discussed to provide help and guidance to the reader for identifying an appropriate model for their use case.

When in doubt which method to use, a rule of thumb could be to prefer kernel methods when there are less than $\sim 10^{3}$–$10^{4}$ training points and NN-based approaches otherwise (but this may also be a matter of preference).

Depending on the desired accuracy, the amount of ab initio reference data which can be collected within a reasonable time frame may be vastly different. For example, if reference calculations are performed at the DFT level of theory, it is often feasible to collect several thousands of data points, even for relatively large molecules. On the other hand, if CCSD(T) accuracy and a large basis set is required, already a few hundred reference calculations for small molecules can require a considerable amount of computing time. Although it is of course always desirable to perform as few reference calculations as possible, for some tasks, collecting a large data set is unavoidable. For example, if a model should be able to predict a variety of different molecules containing many different elements, the relevant chemical space must be sampled sufficiently.

In general, kernel-based models tend to achieve good prediction accuracies even with few training points, whereas NNs often need more data to reach their full potential (although there may be exceptions for both model variants, see also Fig. 5). Further, the optimal model parameters for kernel models can be determined analytically (see Eq. 10), which, at least for small datasets, is typically faster than training a NN via (a variant of) stochastic gradient descent. However, when the data set size $M$ is very large, solving Eq. 10 analytically can become prohibitively expensive as it scales $\mathcal{O}(M^{3})$ (and requires $\mathcal{O}(M^{2})$ memory to store the kernel matrix). Further, evaluating kernel models scales with $\mathcal{O}(M)$ (see Eq. 2), whereas the cost of evaluating NN-based methods has (as long as the number of parameters does not have to be increased for larger datasets) constant complexity. For this reason, NNs tend to be more suitable for large datasets. Note that there are approximations which improve the scaling of kernel methods (so they can be applied even to very large datasets) at the cost of accuracy (see Eqs. 11-13).

To be applicable to multiple systems, a model must decompose its prediction into atomic contributions. Models that use no such decomposition must either use a fixed size descriptor or several separate models need to be trained.

Some ML-FFs only need to be able to predict systems with a fixed composition and number of atoms, for example to study the dynamics of a single molecule, whereas other applications require the ability to predict different systems with varying size, e.g. when clusters consisting of a different number and kind of molecules are studied with the same model.

While all ML-FFs can be applied in the first case, the latter requires either that the length of chemical descriptors is independent of the number of atoms, or that model predictions are decomposed into local contributions based on fixed-size fingerprints of atomic environments (which naturally makes them extensive). Most NNPs (see Section 2.2.2) and many kernel methods, e.g. GAPs (see Section 2.3.3) or FCHL106, 107, use such a decomposition and can be applied to differently sized chemical systems without issues. An exception are e.g. (s)GDML models (see Section 2.3.2), which encode chemical structures as vectors of inverse distances between atomic pairs. Consequently, the length of the descriptor changes with the number of atoms and the model can only be applied to a single type of system. In some special cases, it may be possible to choose a maximum descriptor length and pad descriptors of smaller molecules with zeros, but this may introduce other problems and/or reduce the accuracy.

If only a handful of elements is relevant, all models are equally suitable. When a large number of elements needs to be considered, the model should be able to encode and use information about atom types efficiently.

As long as an ML-FF is only applied to single-component systems (consisting of a single element), for example elemental carbon or silicon, all relevant information is contained in the relative arrangement of atoms and nuclear charges need not be encoded explicitly. However, as soon as there are multiple atom types (as is common for most applications of ML-FFs), the model must have some way to distinguish between them. A notable exception are some models such as (s)GDML, which use inverse pairwise distances as structural descriptor. Here, information about atom types is implicitly contained, because specific entries always correspond to the same combination of atom types.

Many local descriptors of atomic environments only use geometric information in the form of distances and angles between pairs and triplets of atoms (see Eqs. 22 and 23). To include information about atomic types, geometric features have to be included separately for every possible combination of elements, leading to a drastic increase of descriptor size (descriptors for kernel machines based e.g. on SOAP117 or FCHL106, 107 also grow in size when the number of atom types is increased). Many descriptor-based NNPs further use separate NNs to predict atomic contributions of different elements (see Fig. 12). A disadvantage of these approaches is that the number of terms in the descriptor increases combinatorially with the number of elements covered by the model (in particular if three-body or even four-body terms are used), which impacts the computational cost of training and evaluating the model. Also, larger amounts of training data may become necessary for good results. As long as only a few elements need to be considered, these downsides are not an issue, but if a model for a significant fraction of the periodic table is required, a more efficient representation is desirable. Most end-to-end NNPs employ so-called element embeddings (see Fig. 12), which do not become more complex when the number of elements is increased. This has the additional benefit of potentially increasing the data efficiency of the model by utilizing alchemical information. Another alternative is to introduce element-dependent weighting functions (instead of duplicating terms in ACSF descriptors)205.

Models for very large target systems should be able to exploit chemical locality, so that reference calculations for fragments can be used as training data. Additionally, this allows trivial parallelization of predictions over multiple machines.

Often, ML-FFs are used to study small or medium-sized molecules. In such cases, all models are equally applicable. For very large systems containing many atoms however, some methods have particularly advantageous properties. For example, it might be infeasible to run ab initio calculations for the full target system. In this case, being able to fragment the system into smaller parts, for which reference calculations are affordable, is very useful.

To be trainable on such fragments, ML-FFs must introduce an explicit assumption about chemical locality by introducing a cutoff radius. Every method that decomposes predictions into a sum of local atomic contributions can thus be trained in this way. ML-FFs without cutoffs on the other hand need reference data for the complete system (see above). Another advantage of local models is that their predictions are embarrassingly parallel: The contributions of individual atoms can be calculated on separate machines (storing a copy of the model), each requiring only information about neighboring atoms within the cutoff radius. Apart from possible efficiency benefits, this may even become necessary if the computations to handle all atoms do not fit into the memory of a single machine (for example when the system of interest consists of millions of atoms215). Note that while not all ML methods to construct FFs can be parallelized in this way, most models contain mostly linear operations, which are amenable to other parallelization methods, e.g. by utilizing GPUs (graphics processing units).

At this point, a subtle difference between cutoffs used in NNPs of the message-passing type (see Section 2.3.4) and descriptor-based NNPs (as well as kernel machines based on local atomic environments) should be pointed out. In message-passing schemes, information between all atoms within the cutoff radius is exchanged over $T$ iterations, thus the effective cutoff radius increases by a factor of $T$. This means that in order to distribute the computation over multiple machines, it is either necessary to communicate updates to other machines after each iteration, or a sufficiently large subdomain needs to be stored on all machines.

If strong long-range contributions to the energy are present, it is advisable to either use a model without cutoffs, or augment the pure ML approach by explicitly including physical interaction terms.

As described earlier, many ML-FFs introduce cutoffs to exploit chemical locality. An obvious downside of this approach is that all interactions beyond the cutoff cannot be represented. For uncharged molecules without strong dipole moments, relevant interactions are usually sufficiently short-ranged that this is not problematic. However, when strong long-ranged (e.g. charge-dipole) interactions are important, cutoffs may introduce significant errors. Models such as (s)GDML, which consider the whole chemical structure without introducing cutoffs, do not suffer from this issue in principle.

While it is possible to simply increase the cutoff distance until more long-ranged contributions can be neglected, this decreases the computational and data efficiency of models which were designed with cutoffs in mind. A better alternative could be to include the relevant physical interaction terms explicitly in the model. For example, TensorMol194 and PhysNet108 include such correction terms by default, but other models can be augmented in a similar fashion. Although not strictly necessary, even models without cutoffs may profit from such terms by an increased data efficiency.

3 Data collection

A fundamental component of any ML model is the reference data. While its architecture and other technical details are responsible for the potential accuracy of a model, the choice of reference data and its quality defines the reliability and range of applicability of the final model. Any deficiencies that are present in the data will inevitably also lead to artifacts in models trained on it, a principle often colloquially stated as “garbage in, garbage out”216. As such, the reference data is one of the most important components of an ML-FF. The generation of datasets in computational chemistry and physics are challenges on their own. First of all, each reference point is a result of computationally expensive and often nontrivial calculations (see Section 2.1), which limits the amount of data that can be collected. Furthermore, the dimensionality of the configurational space of molecules, solids, or liquids is so vast that – except for trivial cases – it is not apparent how to identify the representative geometries in the ocean of possibilities. The optimal choice of reference data might even need to be adapted to the individual properties of the respective ML model that consumes it and/or its intended application. In the following, several strategies for sampling the PES and generating reference datasets are outlined (multiple of these approaches can be combined). Afterwards, problems that may occur due to insufficient sampling are highlighted and general remarks about the importance of a consistent reference dataset are given.

A good starting point to assemble the reference dataset is by sampling the PES using ab initio molecular dynamics (AIMD) simulations. Albeit expensive in terms of the amount of necessary reference calculations, this technique constitutes a straightforward way to explore configurational space. Here, the temperature of the simulation determines which regions of the PES and what energy ranges (according to the Boltzmann distribution) are explored (see Fig. 14). For example, if the aim is to construct an ML-FF for calculating the vibrational spectrum of ethanol at 300 K, generating the database at 500 K is a safe option since the subspace of configurations relevant at 300 K is contained in the resulting database (see Fig. 14A). Sampling at higher temperatures ensures that the model does not enter the extrapolation regime during production runs, which is practically guaranteed to happen when a lower temperature is used for sampling. In general, the resulting dataset will be biased towards lower energy regions of the PES, where the system spends most of the simulation time. For this reason, pure AIMD sampling is only advisable when the intended application of the final ML model involves MD simulations for equilibrium or close to equilibrium properties, where rare events do not play a major role. Examples of this are the study of vibrational spectra, minima population, or thermodynamic properties.

Constructing reliable reference datasets from AIMD simulations can be computationally expensive. While system size plays a major role, other phenomena, such as the presence of intramolecular interactions and fluxional groups, can also influence how quickly the PES is explored. Because of this, long simulation times may be required to visit all relevant regions. For example, generating $2\times 10^{5}$ conformations from AIMD using a relatively affordable level of theory (e.g. PBE+TS/DFT with a small basis set) can take between a few days to several weeks (depending on the size of the molecule). With higher levels of theory, the required computation time may increase to months, or, when highly accurate methods such as CCSD(T) are required, even become prohibitively long (several years).

To resolve this issue, a possible strategy is to sample the PES at a lower level of theory to generate a long trajectory that covers many regions on the PES. The collected dataset is then subsampled to generate a small, but representative set of geometries, which serve as input for performing single-point calculations at a higher level of theory (see Fig. 15). This strategy works best when the PES has a similar topology at both levels of theory, so it can be expected that configurations generated at the lower level are representative of configurations that would be visited in an AIMD simulation at the higher level (see the two-dimensional projections of the PES in Fig. 15). When the two PESs are topologically very different, e.g. when a semi-empirical method or even a conventional FF is used to generate the initial trajectory, it may happen that the relevant regions of the PES at the higher level of theory are not covered sufficiently. Then, when an ML-FF is trained on the collected dataset and used for running an MD simulation, the trajectory may enter the extrapolation regime and the model might give unphysical predictions. Thus, extra care should be taken when two very different levels of theory are used for sampling by proxy.

Another method to minimize the amount of expensive ab initio calculations is called adaptive sampling or on-the-fly ML218. Here, a preliminary ML-FF is trained on only a small initial set of reference data and then used to run an MD simulation. During the dynamics, additional conformations are collected whenever the model predictions become unreliable according to an uncertainty criterion. Then, new reference calculations are performed for the collected structures and the training of the ML model is continued or started from scratch on the augmented dataset. The process is repeated until no further unreliable regions can be discovered during MD simulations.

When following this strategy, the quality of the uncertainty estimate is crucial for an efficient sampling of the PES: If the estimate is overconfident, deviations from the reference PES might be missed. If the estimate is overly cautious, many redundant ab initio calculations have to be performed. There exist several ways to estimate the uncertainty of an ML-FF. For example, Bayesian methods learn a probability distribution over models, which enables straightforward uncertainty estimates (see the predictive variance of a Gaussian process, Eq. 9). For models where an explicit uncertainty estimate is not available, e.g. neural networks, a viable alternative is query-by-committee219, 204. Here, an ensemble of models is trained, for example on different subsets of the reference data and each starting from a different parameter initialization. Then, the discrepancy between their predictions can be used as uncertainty estimate. Query-by-committee has been successfully employed to sample PESs using neural networks for water dimers220, organic molecules71, 108 as well as across chemical compound space221. Other alternatives, for example using dropout222 as a Bayesian approximation223, could also be used.

Collecting data “on-the-fly” is even possible without uncertainty estimates. Instead, additional reference calculations are performed at fixed intervals during the MD simulation218, 224. This relies on the assumption that the probability of reaching the extrapolation regime of an ML model rises with increasing length of the MD trajectory. While performing ab initio calculations in regular intervals will discover all deviations of the model eventually, this variant of on-the-fly ML does not exploit any information about the already collected reference set and may thus lead to many redundant data points. More detailed reviews on uncertainty estimation and active sampling of PESs can be found in Refs. 225 and 226.

Similar to adaptive sampling, metadynamics sampling227, 228 uses a preliminary ML-FF to run MD simulations to find structures for which to run reference calculations. However, the dynamics are biased to increase the probability for visiting unexplored regions on the PES. This is achieved by placing “Gaussian bump functions” on the PES in regions that have already been visited, i.e. the potential energy of already known structures is artificially raised. It is possible to combine metadynamics with the uncertainty estimates used in adaptive sampling to only select the most relevant structures.

It is also possible to sample the PES without running any kind of MD simulation. In normal mode sampling193, the idea is to start from a minimum on the PES and generate distorted structures by randomly displacing atoms along the normal modes. They are the eigenvectors of the mass-weighted Hessian matrix obtained at the minimum position, i.e. a harmonic approximation of the molecular vibrations. From the associated force constants (related to the eigenvalues), the increase in potential energy for displacements along individual normal modes can be estimated. Since they are orthogonal to each other, it is straightforward to combine multiple random displacements along different normal modes such that the resulting structures are sampled from a Boltzmann distribution at a certain temperature. In other words, structures generated like this are drawn from the same distribution as if an “approximated PES” was sampled with a (sufficiently long) MD simulation. This approximated PES is equivalent to a Taylor expansion of the original PES around the minimum position, truncated after the quadratic term (the contribution of the linear term vanishes at extrema).

Structures generated from random normal mode sampling are not correlated, in contrast to those obtained from adjacent time steps in MD simulations, which makes this approach an efficient way to explore the PES. However, the disadvantage is that only regions close to minima can be sampled. Additionally, the harmonic approximation is only valid for small distortions, i.e. the larger the temperature, the more the sampled distribution diverges from the Boltzmann distribution on the true PES. Because of these limitations, it is best to combine normal mode sampling with other sampling methods, for example to generate an initial reference dataset, which is later expanded by adaptive sampling.

Because their extrapolation capabilities are limited, ML methods only give reliable predictions in regions where training data is present.229 When generating reference data, it is therefore important that all regions of the PES that may be relevant for a later study are sampled sufficiently. For example, when studying a reaction, the data should not only cover configurations corresponding to educt and product structures, but also the region around the transition state and along the transition pathway. When the reaction coordinate defining the transition process is already known, a straightforward way to generate the reference data would be to sample the transition path region. However, even when an ML model can reproduce the entire reference dataset with the required accuracy, it is still possible to run into issues when the model is used to study the reaction. If the rare transition process was not sampled sufficiently, it is not guaranteed that MD simulations with the ML-FF reproduce it correctly. The reference data may be restricted to a specific subset of molecular configurations along the transition pathway. Hence, the model can enter the extrapolation regime somewhere between the boundary states and the transition pathways generated by an MD simulation might be unreasonable. Another potential issue is that after passing the transition state region, typically, a large amount of potential energy is converted to internal motions such as bond vibrations. As a result, the effective temperature defined by the kinetic energy exceeds the ambient conditions by orders of magnitude. Even when using a thermostat in the simulation, thermal energy increases so rapidly that it may not be able to handle the increase in temperature immediately. As a consequence, the trajectory visits high-energy configurations, which may not be included in the reference data, and the model again has to extrapolate.

When ML-FFs enter the extrapolation regime, i.e. they are used to predict values outside the sampled regions of the PES, unphysical effects may be observed. Consider for example the dissociation of the O–H bond in the hydroxyl group of ethanol (Fig. 16). Here, different models were trained on data gathered from an MD simulation of ethanol at 500 K and used to predict how the energy changes when the O–H distance of the hydroxyl group is shortened or elongated to extreme values well outside the range sampled during the dynamics. In this example, while the sGDML model is able to accurately extrapolate to much shorter distances than are present in the training data, it still fails to predict the bond dissociation. The NNP models (PhysNet and SchNet) exhibit qualitatively wrong short-range behavior and spurious minima on the PES, which may trap trajectories during MD simulations. Because of these limited extrapolation capabilities, it is advisable to sample larger regions of the PES than are expected to be visited during MD simulations, so that there is a “buffer” and models never enter the unreliable extrapolation regime during production runs. For example, when an ML-FF is to be used for a study at a temperature of 300 K, the PES should be sampled around 500 K or higher.

Although it may appear trivial, it is crucial that all data used for training a model is internally consistent: A single level of theory (method and basis set) should be used to calculate the reference data. When multiple quantum chemical codes (or even different versions of the same code) are used for data generation, it should be checked that their output is numerically identical when given the same input geometry (if they do not then this will effectively manifest itself like noisy outputs, severely deteriorating the precision of the ML model). Further, many ab initio codes automatically re-orient the input geometry such that the principal moments of inertia are aligned with the $x$-, $y$- and $z$-axes, so extra care should be taken when forces or other orientation-dependent quantities (i.e. electric moments) are extracted to verify they are consistent with the input geometry. When some calculation settings need to be adapted for a subset of the data, e.g. for cases with difficult convergence, it is important to check that values computed with the modified settings are consistent with the rest of the data. Additionally, for training some ML models, it may be essential that atoms are ordered in a particular way throughout the data set. For example, the permutational symmetry of (s)GDML models is limited to the transformations recovered from the training set, whereas the NN models discussed in this review are fully agnostic with respect to atom indexing.

4 Data preparation

After the reference data is collected, it has to be prepared for the training procedure. This includes splitting the data into different subsets, which are reserved for separate purposes. Some models may also require that the data is preprocessed in some way before the training can start. In the following, important aspects of these preparation steps are highlighted.

Prior to training any ML model, it is necessary to split the reference data into disjoint subsets for training/validation and testing (see Section 2.2.3). While the training/validation set is used for fitting the model, the test set is only ever used after a model is trained to estimate its generalization error, i.e. to judge how well the model performs on unseen data.30, 158 It is very important to keep the two splits separate, as it is easily possible to achieve training errors that are several orders of magnitude lower than the true generalization error when the model is not properly regularized. Many models also feature hyperparameters, such as kernel widths, regularization terms or learning rates, that must be tuned by comparing several trained model variants on a third dataset used purely for validation (a subset of the training/validation set). Note that information from the validation set will still enter the model indirectly, i.e. it also participates in the training process. This is why a strict separation of the training/validation set from the test set is crucial. Undetected duplicates in the dataset can complicate splitting, as the contamination of the test set with training data (“data leakage”) might go unnoticed. In this case, the model is effectively trained on part of the test set and estimates of the generalization error might be too optimistic and unreliable. Such a scenario can occur even when no obvious mistakes were made, e.g. when the structures for a dataset are sampled by running a long MD simulation where snapshots are written very frequently. Structures collected from adjacent time steps may be highly correlated in this case and when splitting the data randomly into training and test sets, a large portion of both sets will be almost identical. In such a case, instead of using a random split, a better approach would be to use a time-split of the dataset,230 e.g. using the first 80% of the MD trajectory as the training/validation set and reserving the last 20% for testing.

Prior to training a model, the raw data is often processed in some way to improve the numerical stability of the ML algorithm. For example, a common practice is normalization, where inputs (or prediction targets) are scaled and shifted to lie in the range $-1\dots 1$ or to have a mean of zero and unit variance. The constants required for such transformations must never be extracted from the complete dataset. Instead, only the training set may be used to obtain this information.90, 230, 30 Otherwise, estimates of the generalization error on the test set may be overconfident (this is another form of data leakage). While normalization may be less common for the purpose of constructing ML-FFs, any “data-dependent“ transformation must be done carefully. For example, it may be desirable to subtract the mean energy of structures from the energy labels in order to obtain numbers with smaller absolute values (for numerical reasons). This mean energy should be calculated only from the structures in the training set.

If a model is trained using a hybrid loss that incorporates multiple interdependent properties, such as energy and forces, it is important to consider the effects of the normalization procedure on the functional relationship of those values. For example, multiplying the energy labels by a factor requires that the forces are treated in the same way, because the factor carries over to the derivative (scaling energies and forces by different factors would therefore introduce inconsistencies in the data). Also, while subtracting the mean value from energy labels is valid, it is not correct to add any constant to the force labels, because that would translate into a linear term in the energy domain (the energy is related to the forces through integration). Consequently, the consistency between both label types would be broken and an energy conserving model would be incapable of learning. Even when doing simple unit transformations, care should be taken not to introduce any inconsistencies. For example, when energy labels are converted from $E_{\rm h}$ to kcal mol-1 and atom coordinates from $a_{0}$ to Å, force labels have to also be converted to kcal mol-1 Å-1 so that all data is consistent. Depending on which code was used to obtain the reference data, it is even possible that units for some labels must be converted, because they may be given in different unit systems in the raw data (ab initio codes often report energy and forces in atomic units, whereas for coordinates, angstroms are popular).

5 Training the model

After the data has been collected and prepared, the next step is training the ML-FF. During the training process, the parameters of the model are tuned to minimize a loss function, which measures the discrepancy between the training data and the model predictions. In some cases, e.g. most kernel methods, the optimal solution can be found analytically. When this is not possible, e.g. when training neural networks, the parameters are typically optimized iteratively by gradient descent or a similar algorithm. Because standard gradient descent tends to converge very slowly, some authors have proposed to augment it with terms mimicking momentum231, 232 or adaptive step sizes233, 234. Not only training times, but also the achievable accuracy varies greatly with different optimization algorithms, so it is best to try different schemes (see Ref. 235 for an overview over different popular methods). A good default choice is the Adam optimizer236, which converges quickly and gives good results for many different NN architectures. The hyperparameters of a model (e.g. the number of layers or their width in the case of NNs) can also be selected in this step, albeit by checking the model performance on the validation set after training (instead of optimizing them directly). This section details the training process and highlights important points to consider, e.g. the choice of loss function or how to prevent overfitting of the model to the training data.

For regression tasks, a standard choice for the loss function is the mean squared error (MSE) given by $\mathcal{L}=\frac{1}{N}\sum_{i=1}^{N}(y_{i}-\hat{y}_{i})^{2}$, because it punishes outliers disproportionately. Here, the index $i$ runs over all $M$ samples of the training data, $y_{i}$ is the reference value for data point $i$ and $\hat{y}_{i}$ is the corresponding model prediction. When the MSE is used as loss function, it is implicitly assumed that any noise present in the reference data is distributed normally, which without additional information, is a sensible guess for most data. Further, the MSE loss allows finding the optimal parameters analytically (due to convexity) for linear ML algorithms, such as kernel ridge regression (see Eq. 10 in Section 2.2.1). However, the MSE is not necessarily the best choice for all cases. For example, to make the model less sensitive to outliers, a common alternative is to use a mean absolute error (MAE) loss given by $\mathcal{L}=\frac{1}{N}\sum_{i=1}^{N}\lvert y_{i}-\hat{y}_{i}\rvert$. Other functional forms, such as Huber loss237 or even an adaptive loss238, are also possible, provided they are a meaningful measure of model performance.

After deciding on the general form of the loss function, the question remains which labels $y$ to use as a reference. While the potential energy is an obvious choice, in classical MD, the PES is explored via integration of Newton’s second law of motion, which exclusively involves atomic forces. Since an important objective of ML-FFs is to reproduce the dynamical behavior of molecules in MD simulations as well as possible, it could even be argued that accurate force predictions should take priority over energy predictions in MD applications. However, since energy labels are usually available as a byproduct of force calculations, it seems reasonable to include both label types in the hope that this will help improve the overall prediction performance for both quantities. This gives rise to models based on hybrid loss functions that simultaneously penalize force $\mathbf{F}$ and energy $E$ training errors. Assuming an MSE loss, it generally takes the form

where the hyperparameter $\eta$ determines the relative weighting between both loss terms to account for differences in units, information content, and noise level of the label types. A bilateral reduction of both loss terms is only possible if the objectives are non-competing, i.e. when the optimal parameter set is equally effective across both tasks. For this to be true, it must hold that

at every training point $i$ (here, the relation $\mathbf{F}=-\nabla_{\mathbf{R}}E$ was substituted). Otherwise, the objectives $\mathcal{L}_{\mathbf{F}_{i}}$ and $\mathcal{L}_{E_{i}}$ are necessarily minimized by a different set of model parameters. Eq. 28 is only true in general when $\mathcal{L}_{E_{i}}=\mathcal{L}_{\mathbf{F}_{i}}=0$ for all $i$, which is not fulfilled in practice, because both labels may contain noise and they can usually not be fitted perfectly. A model trained using a hybrid loss (Eq. 27) will thus have to compromise between fulfilling both objectives on the training data, as opposed to joining energy and force labels for a performance gain on both. For this reason, the use of hybrid loss functions (or how to weight different contributions) warrants careful consideration depending on the intended application of the final model. Some models, e.g. (s)GDML (see 2.3.2), do not even include energy constraints in their loss function at all and are trained on forces only. The energy can still be recovered via integration, but it does not participate in the training procedure except for determining the integration constant. In the end, the ultimate measure of a model’s quality should not be how well it minimizes a particular loss function, but instead how well it is able to reproduce the experimental observables of interest. Also, it is important to keep in mind that the loss function measured on the training data is only a proxy for the true objective of any model, which is to generalize to unseen data. Compromising between the energy and force labels of the training data can even improve prediction accuracy for both label types on unseen data. For a more thorough discussion on the role of gradient reference data and how it can improve prediction performance, see Ref. 239, 240, 241.

Hyperparameters, such as kernel widths or the depth and width of a neural network, are typically optimized independently of the parameters that determine the model fit to the data: A hyperparameter configuration is chosen, the model is trained, and its performance is measured on the validation set. This process is repeated for as many trials as are affordable or until the desired accuracy is reached. Here it is crucial that no test data is used to measure model performance when tuning hyperparameters, so the ability to estimate the generalization error on the test set is not compromised. Choosing good values for the hyperparameter regimes requires some experience and intuition of the problem at hand. Fortunately, many models are quite robust and good default hyperparameters exist, which do not require any further tuning to arrive at good results. In other cases, hyperparameter tuning can be automated (for example via grid or random search149, 90, 242, 30, 164) and does not need to be performed manually. See also Section 2.2.3 for a more detailed discussion on tuning hyperparameters.

Because ML models contain many parameters (sometimes even more than the number of data points used for training), it is possible or even likely that they “overfit” to the training data. An overfitted model achieves low prediction errors on the training set, but performs significantly worse on unseen data (Fig. 17A). The aim of regularization methods is to prevent this unwanted effect by limiting or decreasing the complexity of a model.

When the loss function is minimized iteratively by gradient descent or similar algorithms, as is common practice for training NNs, one of the most simple methods to prevent overfitting is early stopping243): In the beginning of the training process, prediction errors typically decrease on both training and validation data. At some point however, because the validation set is not used to directly optimize parameters, the performance on the training data will continue to improve, whereas the loss measured on the validation set will stagnate at a constant value or even begin to increase again. This indicates that the model starts overfitting. Early stopping simply halts the training process as soon as the validation error converges (instead of waiting for convergence of the training error), see Fig. 17B. Early stopping also limits the size of the neural network weights and thus implicitly limits the complexity of the underlying function class. Similar to tuning hyperparameters, only the validation set, but never the test set, must be used for determining the stopping point.

Another method of regularization is the introduction of penalty terms to the loss function. Since overfitted models often are characterized by high variance in the prediction (see Fig. 17A), the idea is to penalize large model parameters. For example, $L^{2}$ regularization (adding the squared magnitude of parameters to the loss) shrinks the $L^{2}$-norm of the parameter vectors towards zero and prevents very large parameter values. On the other hand, $L^{1}$ regularization (adding the absolute values of parameters to the loss) shrinks their $L^{1}$-norm, i.e. it favors sparse parameter combinations. Typically, the regularization term is weighted by an additional hyperparameter $\lambda$ that determines its strength (like all hyperparameters, $\lambda$ has to be tuned on the validation set). Note that solving Eq. 10 to determine the parameters of a kernel method will result in an $L^{2}$-regularized model trained on the MSE loss function.

6 Using ML-FFs in production

The main motivation for training an ML-FF is to use it for some production task, such as running an MD simulation. Before doing so however, it is advisable to verify that it fulfills the accuracy requirements for its intended application. At this point in time, the test set becomes important: Since it was neither used directly nor indirectly during the training process, the data in the test set allows to estimate the performance of a model on truly unseen data, i.e. how well it generalizes. For this, it is common practice to compute summary errors on the test set, for example the mean absolute error (MAE) or root mean squared error (RMSE), as a measure of the overall accuracy of a model. In general, such a way of quantifying accuracy gives an overview of the ML model’s performance on the given dataset and provides a simple way to benchmark.

However, summary errors are biased towards the densely sampled regions of the PES, whereas much larger errors can be expected for less populated regions. Therefore, while summary errors measured on the test set are typically a good indicator for the quality of a model, they are not necessarily the best way to judge how well an ML-FF performs at its primal objective, namely capturing the relevant quantum interactions present in the original molecular system. In other words, performance measures evaluated on the test set should not be trusted blindly. They are only reliable when the test set is representative of the new data encountered during production tasks, i.e. when they are drawn from the same distribution. When a model has to extrapolate, it might give unreliable predictions, even when its performance on the test set is satisfactory. When in doubt, especially when an ML-FF is used for a different task than it was originally constructed for, it is better to collect a few new reference data points to verify that a model is still valid for its use case. Because of the generally limited extrapolation capabilities of ML models, results obtained from studies with ML-FFs should always be scrutinized more carefully than e.g. results obtained with conventional FFs. For example, it is advisable to randomly select a few trajectories and verify that the sampled structures look “physically sensible”, e.g. no extremely short or long bonds are present and atoms have no unusual valencies. Since the PES is a high-dimensional object, rare events, where a trajectory visits a part of configurational space that is not sampled in the reference data, are always possible, even when the PES was carefully sampled. If any questionable model predictions are found, it is advisable to double-check their accuracy with additional reference calculations.

7 Example code and software packages

While many modern ML-FFs are conceptually simple, their implementation is often not straightforward, involving many intricate details that can not be exhaustively covered in publications. Instead, those details are best conveyed by a reference implementation of the respective model. Publicly available well maintained codes allow to replicate numerical experiments and to build on top of existing models with minimal effort.

In this section, example code snippets for training and evaluating kernel- and NN-based ML-FFs with the sGDML164 (3.7.1) and SchNetPack244 (3.7.2) software packages are given. This is followed by a short description of other popular software packages for the construction of ML-FFs (3.7.3) as a first orientation for interested readers. Note that the list is not comprehensive and many other similar packages exist.

A reference implementation of the (s)GDML model is available as Python software package at http://www.sgdml.org 164. It includes a command-line interface that guides the user through the complete process of model creation and testing, in an effort to make this ML approach accessible to broad practitioners. Interfaces to the Atomic Simulation Environment (ASE)245 or i-PI246 make it straightforward to perform MD simulations, vibrational analyses, structure optimizations, nudged elastic band computations, and more.

To get started, only user-provided reference data is needed, specifically a set of Cartesian geometries with corresponding total energy and atomic-force labels. Force labels are necessary, because sGDML implements energy conservation as an explicit linear operator constraint by modeling the FF reconstruction as the transformation of an underlying energy model (see Section 2.3.2). The trained model will give predictions at the accuracy of the reference data and can be queried like any other FF.

The sGDML package uses a proprietary format for its datasets, but scripts to import and export to all file types supported by the ASE package245, which covers most popular standards, are included. To convert a , simply call

The most convenient way to reconstruct a FF is via the command line interface:

This command will automatically generate a fully trained and cross-validated model and save it to a file, i.e. model selection and hyperparameter tuning (see Section 2.2.3) are performed automatically. The parameters and specify the sample sizes for the training and validation subsets, respectively. All remaining points are reserved for testing. Each subset is sampled from the provided reference without overlap.

To use the trained model, the sGDML predictor is instantiated from the file generated above and energy and forces are queried for a given geometry (for example stored in an XYZ file ):

It is also possible to run MD simulations using ASE and the Calculator interface included with the sGDML package:

To run this script, a trained model () and an initial geometry () are needed. The resulting MD trajectory is stored in a file . For more details and applications examples, please visit the documentation at www.sgdml.org/doc/.

7.2 The SchNetPack package

SchNetPack244 is a toolbox for developing and applying deep neural networks to the atomistic modeling of molecules and materials available from https://schnetpack.readthedocs.io/. It offers access to models based on (weighted) atom-centered symmetry functions and the deep tensor neural network SchNet, which can be coupled to a wide range of output modules to predict potential energy surfaces and forces, as well as a growing number of other quantum-chemical properties. SchNetPack is designed to be readily extensible to other neural network potentials such as the DTNN 159 or PhysNet 108. It provides extensive functionality for training and deploying these models, including access to common benchmark datasets. It also provides an Atomic Simulation Environment (ASE)245 Calculator interface, which can be used for performing a wide variety of tasks implemented in ASE. Moreover, SchNetPack includes a fully functional MD suite, which can be used to perform efficient MD and PIMD simulations in different ensembles.

As it is based on the PyTorch deep learning framework247, SchNetPack models are highly efficient and can be applied to large datasets and across multiple GPUs. Combined with the modular design paradigm of the code package, these features also allow for a straightforward implementation and evaluation of new models. Similar to the sGDML package, the central commodity for training models in SchNetPack is a dataset containing the Cartesian geometries (including unit cells and periodic boundary conditions, if applicable) and atom types, as well as the target properties to be modeled (e.g. energies, forces, dipole moments, etc.). More information can be found in Ref. 244.

SchNetPack uses an adapted version of the ASE database format to handle reference data. The package provides several routines for preparing custom datasets, as well as a range of pre-constructed dataset classes for popular benchmarks (e.g. QM9248 and MD17105), which will automatically download and format the data. For example, molecular data from the MD17 dataset can be loaded via

where indicates the molecule for which data should be loaded (e.g. ethanol), while the second argument specifies where the data is stored locally.

SchNetPack also provides a utility script for converting data files in the extended XYZ format, which is able to handle a wide variety of properties, to the database format used internally. Conversion can be invoked with the command

where the arguments specify the file paths to the data file and database in SchNetPack format, respectively.

As for the sGDML package, training and evaluating ML models in SchNetPack can be performed via a command line interface. For example, a basic model can be trained with the script:

Here, [model_type] specifies which kind of NNP to use (wacsf for a descriptor-based NNP using wACSFs205, or schnet for the SchNet109 end-to-end NNP architecture) and [dataset_type] specifies either a preexisting dataset (e.g. qm9 or md17), or a custom dataset provided by the user. The next two arguments are the paths to the reference and the file the trained will be written to. The arguments and specify the sample sizes for the training and validation subsets, while the remaining points are reserved for testing. SchNetPack offers a wide range of additional settings to modify the training process (e.g. model composition, use of GPU, how different properties should be treated etc.), see https://schnetpack.readthedocs.io/.

Once a model has been trained, it can be evaluated in several different ways. The most basic method is to perform predictions via:

It is also possible to use the SchNetPack MD suite to perform various simulations with the trained model. Continuing the above example, a basic MD run can be carried out as:

Simulations can be further modified via hooks, which introduce temperature and pressure control, as well as various sampling schemes. Further documentation of the code package and usage tutorials can be found at https://schnetpack.readthedocs.io/.

7.3 Other software packages

AMP is a Python package designed to integrate closely with the Atomistic Simulation Environment245 (ASE) and aims to be as intuitive as possible. Its modular architecture allows many different combinations of structural descriptors and model types. The main idea of AMP is to construct ML-FFs on-demand, i.e. simulations are first started with an ab initio method and later switched to the ML-FF once the model is sufficiently accurate. The package is described in greater detail in Ref. 190 and on its official website https://amp.readthedocs.io/.

The Atomic Energy NETwork (ænet) package includes tools for constructing and applying neural network-based ML-FFs. It is written in Fortran 95/2003 and utilizes efficient BLAS (Basic Linear Algebra Subprograms) and LAPACK (Linear Algebra PACKage) routines for performing linear algebra. A Python interface is also included. More details can be found on https://github.com/atomisticnet/aenet/.

The DeePMD-kit is a package written in Python/C++ aiming to minimize the effort required to build deep NNPs with different structural descriptors. It is based on the TensorFlow deep learning framework249 and offers interfaces to the high-performance classical and path-integral MD packages LAMMPS250 and i-PI246. More details on the DeePMD-kit can be found in Ref. 251 or on https://github.com/deepmodeling/deepmd-kit/.

Dscribe is a Python package for transforming atomic structures into fixed-size numerical fingerprints.252 These descriptors can then be used as input for neural networks or kernel machines to construct ML-FFs. Supported representations include the standard Coulomb matrix28 and variants for the description of periodic systems119, ACSFs116, SOAP117, and MBTR169. More details can be found on the official website https://singroup.github.io/dscribe/ or in Ref. 252.

The neural network potential package (n2p2) allows to use existing parametrizations of Behler-Parinello NNPs to predict energies and forces (either with standalone tools or with the LAMMPS MD package250), but it also provides training tools for generating new potentials. It is mainly written in C++. For further details, refer to https://compphysvienna.github.io/n2p2/.

The PROPerty Prophet (or short: PROPhet) package uses neural networks to predict the relationship between chemical structure and material properties. As such, it can also be used to generate NN-based ML-FFs. It includes tools to automatically extract properties of interest from the output files of several ab initio codes and an interface to the LAMMPS MD package250. More details can be found on https://biklooost.github.io/PROPhet/.

QML is a toolkit for learning properties of molecules and solids written in Python.253 It supplies building blocks to construct efficient and accurate kernel-based ML models, such as different kernel functions and premade implementations of many different structural representations, e.g. Coulomb matrix28, SLATM130, and FCHL107. The package is primarily intended for the general prediction of chemical properties, but can also be used for the construction of ML-FFs. For further details, refer to the official website https://www.qmlcode.org or the github repository https://github.com/qmlcode/qml/.

The RKHS toolkit is mainly intended for constructing highly accurate and efficient PESs for studying scattering reactions of small molecules. As described in section 2.2.1, the evaluation of kernel-based methods scales linearly with the number of training points $M$ (see Eq. 2). By using special kernel functions and precomputed lookup tables, the RKHS toolkit allows to bring this cost down to $\mathcal{O}(\log N)$. However, it requires that the training data has grid structure, which limits its applicability to small systems, where it is meaningful to sample the PES by scanning a list of values for each internal coordinate. The implemented kernel functions also allow to encode physical knowledge about the long-range decay behavior of certain coordinates, which enables accurate extrapolation well beyond the range covered in the training data. A Fortran90 implementation of the toolkit can be downloaded from https://github.com/MMunibas/RKHS/ and the algorithmic details are described in Ref. 254.

The RuNNer Code was the first implementation of high-dimensional neural network potentials and the source code is freely available. Details on how to obtain access can be found on https://www.uni-goettingen.de/de/560580.html.

The TensorMol package allows to train NNPs that explictly account for electrostatic interactions. It is based on the TensorFlow deep learning framework249 and includes an interface to i-PI246 for performing path integral simulations. For further information, refer to Ref. 194 or https://github.com/jparkhill/TensorMol.

Physical and Chemical Insights from Machine Learned Force Fields

In nature, the atoms in chemical systems are in constant motion, giving rise to various configurations and reactive events. A large number of experimental observations are not based on a single molecule or atom, but instead on ensembles of various species subject to external conditions, such as temperature or pressure. Consequently, properties associated with individual structures are not sufficient to characterize macroscopic systems. One way to compute ensemble averages are molecular dynamics (MD) simulations, where the time evolution of a system is governed by the atomic forces derived from its associated potential energy surface (PES). From the ergodic hypothesis255 it is known that the expected value of an observable $A$ can also be obtained from the time average $\langle A\rangle_{t}=T^{-1}\sum_{t=1}^{T}{A_{t}}$, where $A_{t}$ is the value of $A$ corresponding to the structure at time step $t$ of the dynamics trajectory and $T$ is their total number. Of course, this relation is valid only when the dynamics is long enough to visit all configurations of the system accessible under the simulation conditions.

In order to obtain meaningful statistics with MD simulations, many thousands (or millions) of successive PES evaluations are necessary. Due to their high computational cost, accurate electronic structure PESs quickly become intractable for such simulations, which is why highly efficient classical force fields (FFs) are usually employed for running MD simulations. However, this efficiency comes at a cost: Conventional FFs completely neglect or misrepresent some potentially relevant contributions to the potential energy, such as polarization, charge transfer, or electronic effects, which limits their usefulness in modeling complex chemical phenomena. Machine-learned FFs (ML-FFs) offer a unique combination of computational efficiency and high accuracy, opening up tantalizing new possibilities in the simulation of the dynamics of molecules, surfaces, materials and condensed phases. They are able to model all chemical interactions – including those that are typically neglected by conventional FFs. The high accuracy of ML-FFs allows to obtain qualitatively different and novel insights, which would otherwise only be accessible from computationally infeasible ab initio MD (AIMD) simulations. In the following, some chemical insights made possible by ML-FFs, which could not have been obtained with conventional FFs, are highlighted in greater detail. A brief overview is given in Table 1. Note that the given examples represent only a tiny fraction of the published literature, an exhaustive list is beyond the scope of the current review. Interested readers can find further examples in other review articles, e.g. in Refs. 256, 257, 258.

A good example for the power of ML-FFs is a recent study of the dynamics of small molecules (malondialdehyde, ethanol, salicylic acid, paracetamol, aspirin) with atomic forces at CCSD(T) quality.70 AIMD simulations were run at 500 K at the PBE+TS/DFT level of theory166, 167 and the collected configurations randomly subsampled to calculate energies and forces at the CCSD(T) level of theory (reference data for aspirin was calculated at CCSD accuracy). For each molecule, an ML-FF was constructed from 1000 data points with the sGDML164 method (see Section 2.3.2) and used to run MD simulations at 300 K. Running simulations of this quality with ab initio methods is impossible, as they would require up to a billion times more computation time. Conventional FFs were shown to be no viable alternative to ML-FFs, as they do not adequately describe, or even completely neglect, effects which strongly influence the dynamics – and hence the properties – of the studied molecules (Fig. 18).

For example, in ethanol, the lone pairs of the oxygen atom interact with the partially positively charged hydrogen atoms of the methyl group. Due to this attraction, the configuration where both lone pairs are adjacent to a hydrogen atom is visited most frequently during a dynamics simulation. Any derived property, e.g. the Gibbs/Helmholtz free energy surface (FES) or the infrared spectrum, is only accurate when this effect is properly described. Conventional FFs do not account for lone pairs and are thus unable to predict the molecular properties correctly.

A similar effect can be observed in malondialdehyde. Here, the lone pairs of the two oxygen atoms strongly repel each other, which drives the dynamics away from configurations where they are close. While conventional FFs can crudely model electrostatic repulsion between the oxygen atoms with negative partial charges, the steric contributions from the overlap of the electron clouds is not described, causing a qualitatively different dynamics.

Paracetamol is another molecule where lone pairs influence the stability of specific configurations: The partially positively charged phenyl-hydrogen adjacent to the oxygen atom of the acetamide group interacts with its lone pairs and favors a specific dihedral angle. Additionally, the nitrogen atom of the acetamide group is $sp^{2}$ hybridized, which allows conjugation to the electrons in the phenyl system and leads to the planar geometry of paracetamol. When the nitrogen hybridization state is changed to $sp^{3}$, the energetically favorable interaction is broken and corresponding configurations are thus rarely visited during room temperature dynamics. However, at higher temperatures, the hybridization state may switch frequently – conventional FFs are unable to describe this.

Another important electronic effect can be observed in aspirin. Here, an occupied (lone pair) $n$ orbital of the carbonyl group overlaps with an unoccupied antibonding $\pi^{*}$ orbital in the ester group. This $n\rightarrow\pi^{*}$ interaction dictates the relative arrangement of these functional groups in the global minimum structure.213 The effect is even amplified during dynamics, since thermal fluctuations enhance the overlap.69

These and many other electronic effects, e.g. $n\rightarrow\sigma^{*}$ interactions, hyperconjugation, and Jahn-Teller distortions are captured automatically by ML-FFs. In contrast, including them in conventional FFs would require additional terms, whose functional form (and even which effects need to be modeled) are typically unknown a priori.

2 Thermodynamics

A typical application of classical FFs is the study of thermodynamic properties of bulk systems, such as enthalpies, entropies, and phase diagrams. However, their limited accuracy is a major obstacle for quantitative predictions, as small inaccuracies in the interaction of a few particles will inevitably lead to big discrepancies when studying many particles. A good example for this are van der Waals (vdW) interactions. They are weak contributions to the total potential energy for small molecules in gas phase, but they add up in large condensed systems and bulk materials and can strongly influence their properties and dynamics.280 While conventional FFs account for vdW interactions, they typically do so with a relatively crude model based on the Lennard-Jones potential14, which is insufficient for quantitative predictions in many cases. A prime example is water: It is the most studied liquid in literature and many different conventional FFs for water (some with additional special-purpose terms) have been proposed in the last decades, yet none of them is able to reproduce all experimentally measured properties of water in MD simulations281.

Here, ML-FFs offer a promising alternative. Morawietz et al. 259 trained a descriptor-based NNP on periodic configurations of liquid and crystalline water, for which reference data was calculated with different DFT functionals. MD simulations with the ML-FF revealed that the thermodynamic anomalies of water, such as its density maximum and negative volume of melting, are due to a delicate balance of weak vdW forces. The study was able to accurately predict experimentally measured radial distribution functions, as well as temperature dependent shear viscosities and diffusion coefficients. As ML-FFs are naturally able to describe bond breaking and formation, the study could even investigate proton transfer between different water molecules.

The ability to analyze thermodynamic properties of reactive events is a major advantage of ML-FFs over conventional methods. For example, a recent study investigated the Gibbs free energy of proton transfer in liquid water at a titanium oxide surface.260 A descriptor-based NNP was trained using reference data collected through an adaptive sampling approach and used to run MD simulations. The study revealed that a significant fraction of water molecules forms short-lived hydroxyl groups on the titanium oxide surface, which strongly influence its surface chemistry. Such insights are key to understanding phenomena such as surface functionalization and photocatalytic processes.

Another application where the flexibility of ML-FFs is a major advantage is the modeling of bulk materials. For example, Gaussian approximation potentials (GAPs, see Section 2.3.3) and NNPs have been constructed for elemental carbon261 and silicon262, 263, 264, 265. They allow to investigate a wide range of phenomena of liquid, crystalline, and amorphous solid phases, including defects and crack propagation. Modeling these effects accurately is only possible with ML-FFs or prohibitively expensive AIMD simulations. It is even possible to predict accurate phase diagrams of such systems with ML-FFs262, 263. Since this requires a correct model of bond formation and breaking, as well as changes of bonding patterns, such insights could not be obtained from conventional FFs.

3 Reactions

One of the most significant advantages of ML-FFs over conventional FFs is their natural ability to model chemical reactions. Even in cases where it is possible to construct special purpose classical FFs that are able to describe reactions, they are typically much less accurate than their ML-FF counterparts. For example, a recent study compared an ML-FF constructed with a message-passing NNP with two classical methods to obtain a reactive FF for the Cl–CH3–Br transformation.212 Here, the ML-FF achieved up to three orders of magnitude lower errors and yielded qualitatively and quantitatively different predictions for the Helmholtz free energy surface along the reaction path. It is therefore no surprise that one of the first fields where ML-FFs were employed with great success are reaction dynamics. Here, the chemical transformations associated with molecular collisions over short time and length scales are studied. These simulations offer detailed atomistic insights into the reaction mechanism, providing access to rate constants and scattering cross sections, as well as insights on how the molecular energy is distributed between different modes, all of which can be directly related to experiments. In order to yield quantitative predictions, sufficient statistics and highly accurate PESs are required, making them an excellent application for ML-FFs. Studies typically involve small molecular systems, which are treated at high levels of accuracy, such as the collision of N${}_{2}^{+}$ and Ar266, 267 or the Cl+CH3OH $\rightarrow$ HCl+CH3O/CH2OH reaction268. Typical conventional FFs require fixed bonding patterns and are thus intrinsically unsuited for studying chemical reactions. While there also exist reactive variants of classical FFs, they do not reach the accuracy of ML-FFs. For example, a recent study investigated the thermal activation of methane by MgO+ with a message-passing NNP (see Section 2.3.4) and a reactive classical FF.269 Here, the ML-FF achieved prediction errors up to two orders of magnitude lower than the classical variant compared to ab initio data. In addition, the disagreement between experimental rate constants and those predicted from MD simulations was lower by a factor of two with the ML-FF compared to the values obtained from the classical FF. The remaining discrepancy between prediction and experiment was further investigated and it was determined that the deviation was not due to inaccuracies of the ML-FF per se, but instead could be traced back to the multireference character of the transition state, i.e. problems with the ab initio reference data itself.

Even though it is possible to construct classical reactive FFs for specific reactions, there are cases where this is exceedingly difficult. A good example is a recent study where the photo-tautomerization reaction of acetaldehyde was investigated, which is speculated to be a major pathway for formic acid formation in the atmosphere270. After being photo-excited, acetaldehyde contains enough energy that it may not only tautomerize to ethenol, but also dissociate into carbon monoxide and methane, or into hydrogen and ethenone. An accurate description of all three possible reaction pathways with the same FF is extremely difficult to achieve with conventional methods. The NNP used for the study on the other hand was trained on MP2/aug-cc-pVTZ282, 283 reference data and allowed an unbiased description of all relevant processes at ab initio quality. Analyzing a total of 12000 individual trajectories, the study concluded that the formation of ethenol from photo-tautomerization of acetaldehyde is unlikely under atomospheric conditions. This insight could not have been obtained by running AIMD simulations in a reasonable time frame: The combined simulation time of 1 µs would amount to ten billion single point calculations (a time step of $\Delta t=0.1$ fs was used due to the large excitation energies). In contrast, less than 500k structures were used for training the ML-FF, i.e. the time spent for running ab initio calculations was reduced by more than five orders of magnitude by employing an ML-FF.

Due to the efficiency of ML-FFs, scattering simulations can now even be extended to involve larger organic molecules. For example, a study of the minimum dynamic path284 of Diels-Alder reactions of 1,3-dibromo-1,3-butadiene and maleic anhydride with an end-to-end NNP has revealed that molecular rotations are a major driving force for the formation of products271, an effect which had not been described previously in the literature for this type of reaction. ML-FFs can even be applied to reactions between molecules and surfaces. For example, a study by Liu et al. 272 investigated (reactive) HCl scattering on a gold surface using a descriptor-based NNP.

For a recent review on neural network-based PESs for small molecules and reactions, see Ref. 285.

4 Nuclear Quantum Effects

Predictive simulations of molecular systems and materials require not only highly accurate representations of the potential energy surface (PES), but also appropriate statistical sampling of the PES. While classical MD simulations are sufficient for this in some cases, the quantum nature of nuclei plays an important role in many systems. Nuclear quantum fluctuations are a fundamental phenomenon in nature resulting from Heisenberg’s uncertainty principle290, hence physical and chemical properties of molecular or biological systems, as well as nano- and bulk-materials, may be affected by them up to certain extent. In particular, light elements, such as protons and atoms in the first row of the periodic table, are prone to display nuclear quantum effects (NQEs) even at room temperature. Furthermore, materials or molecules formed by heavier atoms, but having strong bonds or being at low temperatures, exhibit significant NQEs291, 292, 293, 294, 295, 296, 297, 298.

Consequently, in order to generate predictive simulations of many physical properties, NQEs must be incorporated. A widely used methodology to perform quantum dynamics is path integral molecular dynamics (PIMD). This method is based on the isomorphism between a quantum particle and a classical harmonic ring polymer of $P$ beads (i.e. $P$ harmonically-coupled copies of the particle), where the equality holds for $P\to\infty$ (see Fig. 19A)82. In practice, convergence of thermodynamical properties can be achieved using only a small number of beads. For light atoms at room temperature for example, $P\sim$16–32 is often sufficient to converge mechanical properties164, 286, 288, 296. This number can be reduced even further by using more sophisticated thermodynamic estimators299.

Given that PIMD simulations require energies and forces for $P$ copies of the system of interest, it is infeasible to use ab initio methods to derive them in most cases. There are some exceptions: For example, PIMD simulations to study the IR spectrum of the porphycen molecule have been performed using DFT with the B3LYP functional and it was shown that the correct Helmholtz free energy and vibrational spectrum can only be recovered by considering NQEs (see Fig. 19B:3)288. Another example are PIMD simulations of the Zundel model at the CCSD level of theory to study the impact of NQEs on its structure and the 1H magnetic shielding tensor286. However, both of these studies required supercomputers to make the calculations possible in a reasonable time frame. On the other hand, ML-FFs can replicate the same results at a fraction of the computational cost, i.e. speed-ups by a factor of 105–107 (depending on the reference level of theory) can be achieved109, 69. This gain in computational efficiency makes it possible to run PIMD simulations for a wide range of systems and offers the chance to reveal new chemical and physical insights.

For example, Chmiela et al. 105 performed room temperature PIMD simulations of aspirin using a GDML model (see Section 3.7.1) trained on PBE+TS166, 167 reference data to investigate the paths followed between different minima on its PES. In a followup study, Chmiela et al. 69 compared free energies and vibrational density of states of a variety of medium-sized molecules obtained from PIMD simulations with a model trained on CCSD or CCSD(T) reference data to the same quantities obtained from a model trained on PBE+TS/DFT data. The authors found that even though the PESs at the two different levels of theory are very similar, tiny differences may still lead to largely different free energies. Additionally, it was shown that the experimentally determined populations for different conformations of ethanol can only be recovered from simulations when including NQEs (see Fig. 19C:2).

In another study, Schütt et al. 273 investigated the dynamics of C20 fullerene using a NNP trained on PBE+TS/DFT reference data. Here, including NQEs broadens the radial distribution function significantly, which also increases the molecular polarizability160. A change in the distribution of interatomic distances also influences electronic effects: A recent study of Ref. 70 (mentioned earlier in the paragraph on electronic effects) investigated NQEs in small organic molecules160. The study revealed that NQEs can dynamically strengthen molecular interactions by enhancing $n\rightarrow\pi^{*}$ donation through increasing orbital overlap, or by strengthening electrostatic interactions between neighboring charge densities (see Fig. 19C:3). Another interesting observed effect is a temporary change of bond orders, which can lead to emerging localized transient states of methyl rotors. The study also showed that vdW interactions are strengthened by NQEs: Since interatomic distances expand on average due to thermal and quantum dilations, the molecular polarizability is also increased (see Fig. 19C:3). Other observed implications of NQEs include “bonding” between hydroxyl groups and hindered rotor dynamics, which leads to molecular stiffening and smoother Helmholtz free energy surfaces.

ML-FFs also make it possible to go far beyond the system size accessible with standard electronic structure methods. In Ref. 274, a descriptor-based NNP was used to study the influence of NQEs on aqueous NaOH solutions of different concentrations ($\sim$1000 atoms). It could be shown that NQEs exert a subtle influence on the solvation structure in the Na+ environment and significantly increase the proton transfer rates and hence diffusion coefficients of the different species. The accuracy of the ML-FF also made it possible to identify error cancellation effects in the reference method, leading to artificially good agreement with experiment in the absence of NQEs.

5 Excited States

The Born–Oppenheimer approximation breaks down when modeling the dynamics of molecular excited states, which are essential for understanding photochemical processes. An extension to classical MD, which allows for the simulation of such phenomena, is quantum-classical surface hopping MD. In this approach, the excited state dynamics of a molecule are simulated by letting it evolve on a set of PESs associated with the different electronic states. To describe the distribution of the molecule between the different states, the effective PES governing the time evolution changes according to stochastic criteria, e.g. based on coupling terms between the relevant states. The correct quantum statistics are then recovered from multiple independent simulations. These simulations are computationally intensive, as they do not only require the computation of multiple PES, but also different coupling terms. This is further amplified by the need for a large number of trajectories in order to obtain reliable statistics. As such, quantum-classical surface hopping simulations can profit greatly from the efficiency and versatility of ML-FFs.

In Ref. 275, for example, the authors used descriptor-based NNPs to study the excited state dynamics of the methylene imine molecule, as well as regions close to the conical intersection between the singlet ground and excited states. It could be shown that the NNPs are able to recover the effective PES with high accuracy and allow for efficient simulations to estimate the state populations of the system. Here, the coupling between the different surfaces was computed based on the Zhu–Nakamura approximation300, which relies on the energy differences between states. More accurate quantum mechanical descriptions of the inter-state couplings rely on so-called non-adiabatic coupling vectors (NACs), which introduce several additional challenges from an ML perspective. First, NACs exhibit the same rotational equivariance as molecular forces. Second, they grow rapidly for states lying close in energy. And finally, as a quantity computed between different states, they are determined only up to an arbitrary phase. The latter property in particular complicates the construction of ML models, as the random nature of the phase factor needs to be compensated during training. Early works relied on a costly preprocessing of the reference data276. Ref. 277, however, demonstrated that the phase problem can be overcome by introducing phase-less loss functions during the training procedure. Using a modified end-to-end NNP to describe the excited state dynamics of the methylenimmonium cation, it could be shown that using such loss terms completely eliminates the need for a preprocessing step. In addition, the work modeled the NACs as derivatives of a proxy potential, thus accounting for their transformation under rotations of the molecule. The combination of these approaches not only made it possible to obtain accurate population statistics for the studied system, but could also greatly extend the time scales accessible by the simulation beyond the limits of conventional electronic structure approaches.

For a recent review on machine learning for electronically excited sates, see Ref. 301.

6 Spectroscopy

As stated at the beginning of this section, MD simulations are an excellent tool to model the temporal autocorrelation functions of various quantities, which can in turn be used to predict experimental observables, such as diffusion coefficients. These quantities need not be restricted to properties derived from the PES, but encompass other electronic properties such as dipole moments or polarizabilities. Access to the corresponding time autocorrelation functions enables the simulation of a wide range of molecular spectra, which can be directly related to experiment. The most prominent examples are infrared and Raman spectra derived from the autocorrelation functions of dipole moments and polarizabilities, respectively. Both types of vibrational spectra are of great practical interest, since they can be measured accurately via experiment and provide insights into the atomic structure of molecules and materials. However, these spectra can be subject to a series of complex quantum mechanical effects, such as vibrational anharmonicities. Hence, high level electronic structure treatments are required in order to obtain quantitatively accurate predictions of experimental results. Unfortunately, computing the required autocorrelation functions based purely on electronic structure calculations quickly grows prohibitively expensive, as simulations covering sufficient time scales are required in order to yield reliable spectra. In addition, if the influence of temperature or other phenomena should be studied in detail, a large number of such simulations is required. Recently, ML-FFs have emerged as invaluable tools for obtaining reliable molecular spectra. A growing number of ML-FFs now provide access to quantities beyond the PES, e.g. dipole moments or polarizabilities. As such, they offer the possibility to perform these simulations in only a fraction of the time required by an ab initio approach or even make them possible at all.

Ref. 71 demonstrates the potential inherent to ML-FFs based on the prediction of infrared spectra for organic molecules including the protonated alanine tripeptide. By combining a descriptor-based NNP model of the PES with a dipole moment model based on latent NN-predicted atomic charges, highly accurate infrared spectra could be obtained for all studied systems. The efficiency of such an approach was demonstrated based on an alkane containing more than 200 atoms, where it was possible to reduce a projected computation time of 9000 years with the original ab initio method to only a few days (including the reference calculations needed for training the ML models). Moreover, the high accuracy of the predictions made it possible to identify shortcomings of the original reference methods and study how they influence the infrared spectrum of the tripepdtide. A similar latent charge based approach was employed in Ref. 194 to model infrared spectra of various amino acids. This study could not only obtain accurate spectra, but also demonstrated that the latent charges predicted by the dipole model constitute a valid ML driven scheme for deriving atomic partial charges, which can be used to model long-range electrostatic interactions explicitly. This scheme has since been employed in many physically augmented models (e.g. TensorMol194 or PhysNet108).

In a similar manner, ML models capable of predicting polarizability tensors offer access to Raman spectra. Ref. 278 introduces a symmetry adapted approach for modeling polarizability tensors using Gaussian process regression (GPR) based on the SOAP117 kernel. The authors use this model to study the Raman spectra of paracetamol in gas phase and various molecular crystals and achieve excellent agreement with electronic structure methods in both cases. Not only is the proposed approach highly data efficient, requiring only a small number ($<$1000) of reference data, but it could also be shown that the resulting model is transferable between different polymorphic forms of the crystal. Ref. 279 models Raman spectra of liquid water using descriptor-based NNPs to predict molecular polarizabilities. The computational efficiency of the approach made it possible to obtain Raman spectra for a system containing 416 water molecules based on two nanosecond trajectories at DFT level accuracy, a feat which would be infeasible with the original reference method. As a consequence, the influence of temperature effects on the Raman spectra of water and heavy water could be studied in detail. The atomic resolution of the employed ML approach made it possible to decompose the simulated spectra into intramolecular and intermolecular contributions, offering insights into the mechanisms governing the temperature dependence of the different spectral features.

Challenges

Following the best practices outlined in the previous section, the current generation of ML-FFs is applicable to a wide range of problems in chemistry that involve small- to medium-sized systems. While this space of chemical compounds is already significant in size, the “dream scenario” of chemists and biologists referenced in the introduction can only be realized with access to larger system sizes. Not only does the number of stable structures increase exponentially with added atomistic degrees of freedom,302, 303 many interesting phenomena play out at nanoscale resolution, which is inaccessible to ML methods as of yet. This is because some steps involved in the construction of ML-FFs, like sampling the reference data, which are solvable at small scale, become seemingly insurmountable obstacles at larger scales due to unfeasible computing times. The complexity of interactions, e.g. the non-classical behavior of nuclei, as well as significant contributions from large fluctuations, increase the space of conformations that need to be learned. To further complicate things, the cost of accurate ab initio calculations increases steeply with expanding system size, limiting the amount of reference data that can be collected within a reasonable time frame. This also means that a growing number of atom correlations need to be represented by a model in order to capture the full scope of interactions present in the real system. Below, some considerations in reconciling the somewhat contradicting demands of scalability, transferability, data efficiency and accuracy in large-scale ML-FFs are outlined.

A fundamental challenge that must be faced by ab initio methods, conventional FFs, and ML models alike is the many-body problem. Most properties of a physical system are determined by the interaction of many particles, whether those are electrons or, on a higher abstraction level, atoms. In fact, the reason that ab initio calculations are expensive to obtain is due to the challenging computational scaling properties of high-dimensional many-body problems. As a result, the hierarchy of different levels of theory is directly defined by the level of correlation treatment in the respective wave function parametrization. Because the number of electronic degrees of freedom of a system is much higher than the number of atoms, the computational limitations of ab initio methods become evident very quickly, even for small systems. Atomistic approximations scale more favorably, because they need to correlate less particles, but they are subject to the same scaling laws. The only escape is to neglect some correlations in favor of a reduced problem size. Unfortunately, it is to date impossible to reliably determine which interactions can be removed with minimal impact, without compromising the full many-body solution. Thus, the ideal of a local model is in conflict with the very nature of many-body systems. Although it is possible to recover some effects such as non-local charge transfer by means of a charge equilibration scheme,304, 305 a general solution for this problem does not exist. While not fully justified from a physics perspective, assuming locality is still a useful inductive bias, which can help generalization and computational efficiency. It also helps when collecting reference data, as it implies that larger systems can be predicted using the information learned from smaller systems. Another assumption, which all ML-FFs discussed in this review make, is that the PES is smooth. This is a necessary requirement for most practical applications, since a non-smooth PES implies force discontinuities, which would lead to instabilities during MD simulations. Smoothness is also a requirement from the ML perspective, as only regular signals can be reconstructed from limited observations.

For most commonly used NNPs and many kernel-based ML-FFs, locality is built into the design explicitly through the introduction of a cutoff radius. The global interactions between atoms are modeled by accumulating individual local atomic contributions. In this “mean-field approximation”, the interaction of a particle with its surroundings is reduced to an effective one-body problem, i.e. an interaction of that particle with the average effect of its neighbors. As similar neighborhoods can be identified in different compounds across chemical space, these assumptions allow to build models from reference calculations of small molecules, which are transferable to much larger structures.306, 307 However, the lack of explicit higher-order terms comes at the cost of potentially loosing some important interaction effects, similar to the Hartree-Fock method and Kohn-Sham DFT in ab initio calculations.

On the other hand, some models (e.g. (s)GDML) capture global correlations in the sense that a single prediction is obtained for the whole structure. Of course, this relies on reference calculations that are accurate enough to contain the relevant information. Global interactions of large systems can not be accurately inferred from a training set of small molecules or molecular fragments, which is why reference calculations for the exact target structure are necessary. It can therefore become difficult to collect enough reference data for large structures. In addition, even models that use no atom-wise decomposition might still implicitly assume that interactions are local to some degree due to their chemical descriptor. For example, in (s)GDML models, systems are encoded as a vector of inverse pair-wise distances. Therefore, structural changes between distant atoms contribute less strongly to changes in the overall descriptor than proximal atoms.

While locality and smoothness are valid assumptions for the majority of chemical systems, there are pathological cases where they break down and ML models that rely on them perform poorly. As an example, consider cumulenes – hydrocarbons of the form C2+nH4 ($n\geq 0$) with $n+1$ cumulative double bonds. These molecules have a rigid linear geometry with the two terminal methylene groups forming an equilibrium dihedral angle of 0∘ (when $n$ is even) or 90∘ (when $n$ is odd). Rotating the dihedral angle out of its equilibrium position results in a sharp increase in potential energy even though the methylene groups may be separated by several angstroms when $n$ is large. This is due to the energetically favorable overlap of $\pi$-orbitals along the carbon chain (a highly non-local interaction), which is broken when the methylene groups are rotated against each other. Additionally, the potential energy exhibits a sharp “cusp” at the maximum energy (i.e. it is not smooth), because the ground state electronic configuration switches abruptly from one state to another (strictly speaking, multi-reference calculations would be necessary here). One-dimensional projections of the PESs predicted by ML-FFs along the rotation of the dihedral angle reveal several problems (Fig. 20). For example, all models predict smooth approximations by design, which is beneficial for running MD simulations, but results in large prediction errors around the cusp. Further, when the number of double bonds ($n+1$), i.e. the “non-locality” of relevant interactions, is increased, the quality of predictions decreases dramatically, until all models are unable to reproduce the energy profile.

Note that by design, NNPs relying on message-passing are unable to resolve information about the dihedral angle if information between hydrogen atoms on opposite ends of the molecule cannot be exchanged directly (i.e. $r_{\rm cut}$ is too small) and predict constant energies in this case. The same is true for descriptor-based NNPs, as fingerprints of chemical environments also only consider atoms up to a cutoff (see Eqs. 22 and 23). Any kernel method taking as input local structural descriptors relying on cutoff radii (e.g. SOAP117 or FCHL19107) will suffer from the same problems. Even when a “global” descriptor of chemical structure such as inverse pair-wise distances is chosen (e.g. Coulomb matrix28), changes in the dihedral angle between distant groups of atoms are not resolved sufficiently for accurate predictions (see sGDML model in Fig. 20). The only way to fix this problem in general is to drop the locality assumption completely, for example by including all $\binom{N}{4}$ possible dihedral angles in the structural descriptor (without introducing additional factors that decrease the weight of these features with increasing distance between atoms). However, due to the combinatorial explosion of the number of possible dihedral angles, this would lead to extremely large descriptors whenever the number of atoms $N$ is not very small. The resulting models would be slow to evaluate and require a lot of reference data to give robust predictions (to prevent them from entering the extrapolation regime). An expert choice, i.e. including only a single relevant dihedral angle in the descriptor, is a possible way around this issue, but requires prior knowledge of the problem at hand and goes somewhat against ML philosophy.

As a final remark, it should be mentioned that conventional FFs only include terms for dihedral angles between directly bonded atoms, so they are equally unable to predict the energy profiles of the larger cumulenes shown in Fig. 20. As such, relying on chemical locality is an assumption made by virtually all methods for approximating PESs and is not specific to just ML methods.

2 Transferability, scalability and long-range interactions

The concept of chemical locality discussed above also plays a central role in the transferabilty and scalability of ML models for atomistic systems. Transferability indicates how well models can adapt to compounds varying in their chemical composition, while scalability indicates how efficiently these models scale with respect to the size of systems modeled. Both concepts are closely related and inherently rooted in chemical locality. The assumption that interactions between atoms are local implies that similar structural motifs will give rise to comparable interactions and hence similar contributions to the properties of a molecule or material. In an ML context, chemical locality allows a model to reuse the information learned for different parts of a molecule for similar features in different systems. In this manner, a large atomistic system could in principle be assembled from smaller components like a jigsaw puzzle.306 The former aspect is crucial in order to make models transferable, while the latter allows for the development of architectures whose evaluation cost scales linearly with system size.

ML-FFs exploiting chemical locality offer several advantages compared to other models. If trained properly, they can be applied to systems of different size and composition. The training procedure benefits in a similar manner, as local models can be trained on structures containing different numbers of atoms. Moreover, it is also possible to use only fragments of the original system during construction of a model. This property is very attractive in situations where accurate reference computations for the whole system are infeasible due to system size and/or scaling of the computational method. Local chemical environments are also less diverse than complete structures, potentially reducing the need for extensive sampling and decreasing the chances that models enter the extrapolation regime in a production setting. In addition, local models scale linearly with system size, as interactions are limited to the cutoff radius and can be evaluated efficiently. In contrast, models without cutoffs are typically more limited in their practical applicability for extended systems. They always require reference computations to be performed for the whole system and, once trained, can only be reused for this particular molecule or material.

Despite these advantages, local ML models suffer from several inherent problems. In order to construct models which exploit locality, a chemical system needs to be partitioned in one way or another. This can for example be achieved by limiting interactions to terms involving only a certain number of atoms (similar to conventional FFs) or by restricting them to local atom-centered environments. These approximations place strong limitations on which kind of interactions can be described. As a result, local ML models have difficulty when dealing with the situations where non-local effects are important, such as strongly conjugated systems and excited states (see Section 5.1). For standard simulations, the presence of long-range interactions, such as electrostatic and dispersion effects, are much more common phenomena. These are particularly important for modeling extended systems, where ML models are typically believed to offer a significant advantage over more conventional FFs. Since the structure and dynamical behavior of such systems is influenced greatly by long-range interactions, ML models need to be able to account for them in a satisfying manner.

Recovering long-range effects necessitates a balancing act between physical accuracy and computational efficiency, as the scalability of local models hinges on there being a limited number of interactions which need to be evaluated. This feat is further complicated by the typical energy scales of these interactions, which are small compared to local contributions such as bond energies. For these reasons, it is not advisable to account for long-range interactions by simply increasing the size of local environments. While local models with sufficiently large cutoffs are able to learn the relevant effects in principle, it may require a disproportionately large amount of data to reach an acceptable level of accuracy for an interaction with a comparably simple functional form. The reason is that average gradients and curvature in different regions of the PES may differ by several orders of magnitude, which makes it difficult to achieve uniformly low prediction errors across all regions. Hence, an optimal description would require to employ different characteristic scales.

For illustration, consider the following toy examples: In the first variant, a Lennard-Jones (LJ) potential14 is separated into a region around its minimum, a repulsive short-range, and an attractive long-range part. The task is to learn each of the three regions with a separate model (see Fig. 21a). In the second variant, a single model is trained on all regions at once (see Fig. 21b). Here, all models are kernel-based and use a Gaussian kernel (Eq. 4). The kernel hyper-parameter $\gamma$ is optimized by a grid-search and cross-validation. Compared to the models trained on individual regions, the prediction errors of the model for all regions increase by around an order of magnitude. Further, it shows spurious oscillations between training points in the long-range region. When the optimal values of $\gamma$ for the different models are compared, the reason for failure when training on all regions at once becomes apparent: The optimal values of $\gamma$ are $198.88$, $75.47$, $0.08$ for the short-, middle-, and long-range models, respectively, which highlights the multi-scale nature of the PES. On the other hand, when training on all regions at once, the model necessarily has to compromise, which leads to an optimal value of $\gamma=22.12$. In this toy example, the multi-scale problem can be solved by switching from using $r$ as a structural descriptor to the more appropriate inverse distance $r^{-1}$ (Fig. 21c). Unfortunately, for realistic (high-dimensional) PESs with multiple minima, it can be difficult to find an appropriate descriptor to address the multi-scale nature of the PES, which leads to data-inefficient models. As a result, more training data is needed to reach an acceptable accuracy, which is problematic considering the computational cost of high-quality reference calculations.

One possibility of overcoming these limitations is by instead partitioning the energy into contributions modeled entirely via ML (short-range) and contributions described via explicit physical relations based on local quantities predicted via ML (long-range). A prime example for such an approach is the treatment of electrostatics, as was first introduced in Ref. 220. Here, an ML model is used to predict partial charges for each atom based on their local environment. These charges can then be used in standard Coulomb and Ewald summation to compute the long-range electrostatic energy of a system. While such schemes initially relied on point charge reference data obtained from (arbitrary) partitioning methods of the ab initio electron density (e.g. Hirshfeld charges309), they have since been extended to operate on charges derived from an ML model for dipole moments (a true quantum mechanical observable).71, 194, 108 Here, scalar partial charges $q_{i}$ are predicted for each atom $i$ and the molecular dipole moment is constructed as $\boldsymbol{\mu}=\sum_{i}q_{i}\mathbf{r}_{i}$, where $\mathbf{r}_{i}$ are the atomic positions (the predicted $q_{i}$ can be corrected to guarantee charge conservation108). The discrepancy between reference and predicted dipole moments is included in the loss function used for training the model (see Section 3.5) and the partial charges consequently derived in a purely data-driven manner.

Contrary to electrostatics, accounting for dispersion interactions is not as straightforward, because the exact physical form of dispersion interactions is still debated and a variety of approximate schemes have been proposed.280 In addition, dispersion corrections typically depend on coefficients computed from atomic polarizabilities as local properties. The corresponding quantum mechanical observable is the molecular polarizability tensor. In contrast to charges (scalars) derived from dipole moments (a vector quantity), predicting molecular polarizabilities requires rotationally equivariant ML models.310 Because of this, many ML approaches rely on the same empirical pair-wise dispersion potentials employed for correcting density functional theory computations.311, 312, 108

To summarize, local ML architectures are a promising approach towards transferable and scalable models, but they have a number of drawbacks which will still need to be addressed in the future. Promising alternative approaches to achieve transferability are ML models based directly on electronic structure methods, i.e. “semi-empirical ML”313, 314, 57 and models for electron density and Hamiltonians55. These approaches express fundamental quantum chemical quantities in a local representation, e.g. Hamiltonian matrix elements in an atomic orbital basis. Non-locality can then be introduced via the “correct” mathematical mechanism, e.g. matrix diagonalization in the case of Hamiltonians. This physically motivated structure allows such models to recover a wide range of interactions while still being transferable. They are also better suited to predict intensive properties of molecules (whose magnitude is independent of system size), where assuming additive atomic contributions is not valid. A downside of such models compared to conventional ML-FFs is the increased computational cost due to the additional matrix operations.

With respect to scalability, hybrid approaches similar to QM/MM25 might constitute valid alternatives to pure ML models. Although several orders of magnitude more efficient than electronic structure theory, even local ML models encounter problems when faced with systems containing tens of thousands of atoms. Compared to conventional FFs, the more complex functional form underlying ML-FFs leads to an increased computational cost. In such cases, partitioning the system into regions treated at different levels of approximation can lead to a significant speedup. ML models can for example be embedded into regions modeled by classical force fields, yielding ML/MM like simulation protocols. Restricting elaborate ML approaches to only a subset of a chemical systems would make it possible to employ more accurate approximations in a manner analogous to conventional QM/MM. For example, in Ref. 315, the authors study protein–ligand binding with a ML/MM approach: The ligand is described by an NN-based ML-FF and treated as if it was in gas phase. Coupling to the protein environment (described by a conventional FF) is achieved solely through non-bonded dispersion and electrostatic interactions. The disadvantage of such a simple embedding is that the “quantum region” cannot be polarized by the “classical region”. A more sophisticated embedding was recently proposed by Gastegger et al. 316. Here, the region described by the ML-FF is explicitly polarized by the electric field induced by surrounding point charges, i.e. the electric field is an additional input to the model. Alternative approaches, describe the effect of the classical environment by augmenting structural descriptors such as ACSFs by additional terms explicitly depending on the MM point charges.317 A similar approach is followed in Ref. 318, where the classical environment is described by auxiliary atom types.

Concluding remarks

The last decades have witnessed significant advances in statistical learning that allowed ML techniques to enter our daily lives, industrial practice and scientific research.

Classically, automation in industry and scientific fields relied on hand-crafted rules that represented human knowledge.319 Not only is the creation of rule-based systems laborious and may require to cover an excessive number of cases, it often leads to rigid structures that are unable to adapt well to new situations. Even worse, some concepts are difficult or impossible to formalize, such as human perception for image classification.

Modern statistical ML algorithms127, 98 such as deep learning94, 320, 321, 322 or kernel-based learning323, 324, 92, 90, 125 enable models that freely adapt to knowledge that is implicitly contained in datasets (in an abstract form) and thus offer a more robust way of solving problems than rule-based reasoning. For the field of molecular simulations, the potential of ML methods may help to bridge the accuracy-efficiency gap between first-principles electronic structure methods and conventional (rule-based) FFs. Bringing both fields together has raised many questions and still poses some fundamental challenges for new generations of ML-FFs. At this point in time, ML-FFs have already become a successful and practical tool in computational chemistry.

Starting from a broad perspective, this review has focused on the role of ML for constructing force fields and assessed what can be achieved with these new techniques at the current stage of development. This has been contrasted with problems that are (so far) beyond the reach of present methods. Illustrative examples of the relevant chemistry and ML concepts have been discussed to demonstrate the practical usefulness that modern ML techniques can bring to chemistry and physics. This includes an overview of the most important considerations behind the construction of modern ML-FFs, such as the incorporation of physical invariances, choice of ML algorithms, and loss functions. Special attention has been given to the topic of validating ML-FFs, which requires particular care in scientific applications.325 Furthermore, a comprehensive list of best practices, pitfalls, and challenges has been provided, which will serve as a useful guideline for practitioners standing on either side of this growing interdisciplinary field. These “tricks of the trade”156 are often assumed to be obvious and thus omitted from publications – here they have been deliberately spelled out to avoid unnecessary barriers to enter the field. Additionally, a small catalog of software tools that can enable and accelerate the implementation of ML-FFs has been provided as a pointer for readers wishing to adopt ML methods in their own research.

While routinely performing computational studies of condensed phase systems (e.g. proteins in solution) at the highest levels of theory is still beyond reach, ML methods have already made other “smaller dreams” a reality: Just a decade ago, it would have been unthinkable to study the dynamics of molecules like aspirin at coupled cluster accuracy. Today, a couple hundred ab initio reference calculations are enough to construct ML-FFs that reach this accuracy within a few tens of wavenumbers.326 In the past, even if suitable reference data was available, constructing accurate force fields was labor-intensive and required human effort and expertise. Nowadays, by virtue of automatic ML methods, the same task is as effortless as the push of a button. Thanks to the speed-ups offered by ML methods over conventional approaches, studies that previously required supercomputers to be feasible in a realistic time frame286, 288 can now be performed on a laptop computer109, 69.

In addition to enabling studies that were prohibitively expensive in the past, ML methods have also led to new chemical insights on systems that were thought to be already well understood. Even relatively small molecules were shown to display non-trivial electronic effects, influencing their dynamics and allowing a better understanding of experimental observations.160 Many other unknown chemical effects potentially wait to be discovered by studies now possible with ML-FFs. At the speed at which improvements to existing ML-FFs are published, it is not unreasonable to expect significant advances that will make similar studies possible for larger systems and help realize many more “dreams” in the near future.

Concluding, ML-FFs are a highly active line of research with many unexplored avenues and attractive applications in chemistry, with possibilities to contribute to a better understanding of fundamental quantum chemical properties and ample opportunity for novel theoretical, algorithmic and practical improvement. Given the success of this relatively young interdisciplinary field, it is to be expected that ML-FFs will become a fundamental part of modern computational chemistry.

OTU acknowledges funding from the Swiss National Science Foundation (Grant No. P2BSP2_188147). AT was supported by the European Research Council (ERC-CoG “BeStMo”). KRM was supported in part by the Institute of Information & Communications Technology Planning & Evaluation (IITP) grants funded by the Korea Government (No. 2017-0-00451, Development of BCI based Brain and Cognitive Computing Technology for Recognizing User’s Intentions using Deep Learning) and funded by the Korea Government (No. 2019-0-00079, Artificial Intelligence Graduate School Program, Korea University), and was partly supported by the German Ministry for Education and Research (BMBF) under Grants 01IS14013A-E, 01GQ1115, 01GQ0850, 01IS18025A, 031L0207D and 01IS18037A; the German Research Foundation (DFG) under Grant Math+, EXC 2046/1, Project ID 390685689. We would like to thank Stefan Ganscha for his valuable input to the manuscript. Correspondence to AT and KRM.

Biographies

Oliver T. Unke is an SNSF postdoctoral research fellow in the Machine Learning Group at Technische Universität Berlin. He received his Ph.D. in Chemistry from the University of Basel in 2019. His research has focused on developing methods for constructing accurate potential energy surfaces and their application in molecular dynamics simulations.

Stefan Chmiela is a senior researcher at the Berlin Institute for the Foundations of Learning and Data (BIFOLD). He received his Ph.D. from Technische Universität Berlin in 2019. His research interests include Hilbert space learning methods for applications in quantum chemistry, with particular focus on data efficiency and robustness.

Huziel E. Sauceda is a postdoctoral researcher in the Machine Learning Group at Technische Universität Berlin and part of the BASLEARN project at the same institution. He obtained his bachelor degree in Physics at the Universidad Autónoma de Sinaloa in Mexico, and his master and Ph.D. at the Institute of Physics of the Universidad Nacional Autónoma de México (UNAM) in Mexico City. Between 2016 and 2019, he was a postdoctoral researcher at the Fritz Haber Institute of the Max Planck Society in Berlin. His research interests include ab initio simulations, nuclear quantum effects, and thermodynamics of (nano)materials, as well as development and applications of machine learning methods to quantum chemistry and materials science.

Michael Gastegger is a postdoctoral researcher in the BASLEARN project of the Machine Learning Group at Technische Universität Berlin. He received his Ph.D. in Chemistry from the University of Vienna in Austria in 2017. His research interests include the development of machine learning methods for quantum chemistry and their application in simulations.

Igor Poltavsky is a senior researcher at the University of Luxembourg. He received his Ph.D. from B. Verkin Institute for Low Temperature Physics & Engineering in 2009. His research interests include statistical physics, imaginary-time path integral methods, nuclear quantum effects, ab initio simulations, and machine learning.

Kristof T. Schütt is a senior researcher at the Berlin Institute for the Foundations of Learning and Data (BIFOLD). He received his master’s degree in computer science in 2012 and his PhD at the machine learning group of Technische Universität Berlin in 2018. Until September 2020, he worked at the Audatic company developing neural networks for real-time speech enhancement. His research interests include interpretable neural networks, representation learning, generative models, and machine learning applications in quantum chemistry.

Alexandre Tkatchenko is a Professor of Theoretical Chemical Physics at the University of Luxembourg and Visiting Professor at Technische Universität Berlin. He obtained his bachelor degree in Computer Science and a Ph.D. in Physical Chemistry at the Universidad Autonoma Metropolitana in Mexico City. Between 2008 and 2010, he was an Alexander von Humboldt Fellow at the Fritz Haber Institute of the Max Planck Society in Berlin. Between 2011 and 2016, he led an independent research group at the same institute. Tkatchenko serves on editorial boards of two society journals: Physical Review Letters (APS) and Science Advances (AAAS). He received a number of awards, including elected Fellow of the American Physical Society, the 2020 Dirac Medal from WATOC, the Gerhard Ertl Young Investigator Award of the German Physical Society, and two flagship grants from the European Research Council: a Starting Grant in 2011 and a Consolidator Grant in 2017. His group pushes the boundaries of quantum mechanics, statistical mechanics, and machine learning to develop efficient methods to enable accurate modeling and obtain new insights into complex materials.

Klaus-Robert Müller has been a professor of computer science at Technische Universität Berlin since 2006; at the same time he is directing rsp. co-directing the Berlin Machine Learning Center and the Berlin Big Data Center. He studied physics in Karlsruhe from 1984 to 1989 and obtained his Ph.D. degree in computer science at Technische Universität Karlsruhe in 1992. After completing a postdoctoral position at GMD FIRST in Berlin, he was a research fellow at the University of Tokyo from 1994 to 1995. In 1995, he founded the Intelligent Data Analysis group at GMD-FIRST (later Fraunhofer FIRST) and directed it until 2008. From 1999 to 2006, he was a professor at the University of Potsdam. He was awarded the Olympus Prize for Pattern Recognition (1999), the SEL Alcatel Communication Award (2006), the Science Prize of Berlin by the Governing Mayor of Berlin (2014), the Vodafone Innovations Award (2017). In 2012, he was elected member of the German National Academy of Sciences-Leopoldina, in 2017 of the Berlin Brandenburg Academy of Sciences and also in 2017 external scientific member of the Max Planck Society. In 2019 and 2020 he became Highly Cited researcher in the cross-disciplinary area. His research interests are intelligent data analysis and Machine Learning in the sciences (Neuroscience (specifically Brain-Computer Interfaces), Physics, Chemistry) and in industry.

References