High Accuracy Uncertainty-Aware Interatomic Force Modeling with
  Equivariant Bayesian Neural Networks
Even though Bayesian neural networks offer a promising framework for modeling
uncertainty, active learning and incorporating prior physical knowledge, few
applications of them can be found in the context of interatomic force modeling.
One of the main challenges in their application to learning interatomic forces
is the lack of suitable Monte Carlo Markov chain sampling algorithms for the
posterior density, as the commonly used algorithms do not converge in a
practical amount of time for many of the state-of-the-art architectures. As a
response to this challenge, we introduce a new Monte Carlo Markov chain
sampling algorithm in this paper which can circumvent the problems of the
existing sampling methods. In addition, we introduce a new stochastic neural
network model based on the NequIP architecture and demonstrate that, when
combined with our novel sampling algorithm, we obtain predictions with
state-of-the-art accuracy as well as a good measure of uncertainty.
Tim Rensmeyer, Benjamin Craig, Denis Kramer, Oliver Niggemann
High Accuracy Uncertainty-Aware Interatomic Force Modeling with Equivariant Bayesian Neural Networks

###### Abstract

Even though Bayesian neural networks offer a promising framework for modeling uncertainty, active learning and incorporating prior physical knowledge, few applications of them can be found in the context of interatomic force modeling. One of the main challenges in their application to learning interatomic forces is the lack of suitable Monte Carlo Markov chain sampling algorithms for the posterior density, as the commonly used algorithms do not converge in a practical amount of time for many of the state-of-the-art architectures. As a response to this challenge, we introduce a new Monte Carlo Markov chain sampling algorithm in this paper which can circumvent the problems of the existing sampling methods. In addition, we introduce a new stochastic neural network model based on the NequIP architecture and demonstrate that, when combined with our novel sampling algorithm, we obtain predictions with state-of-the-art accuracy as well as a good measure of uncertainty.

## 1 Introduction

Despite the fact that the laws of quantum mechanics, which underly chemistry, were discovered almost a century ago, their application for the _ab initio_ prediction of many chemical phenomena remains a formidable task [1, 2]. This is in large parts a result of the sheer computational complexity involved in solving these equations numerically as well as the enormous phase space that needs to be considered[1, 2].

Molecular dynamics (MD), where the time evolution of molecular systems is investigated, is particularly affected by this challenge. While computational methods such as Density Functional Theory (DFT) have been developed that can calculate interatomic forces with very high accuracy, these methods are in general computationally expensive and scale badly with growing system sizes. Thus, it remains very difficult to model larger systems with high accuracy and a large enough time horizon.

To circumvent these difficulties, the learning of interatomic forces via machine learning models such as neural networks has become a highly active area of research [3]. However, the cost of generating training data from _ab initio_ simulations prevents the generation of large training sets commonly required to construct sufficiently predictive neural networks. Recent innovations in neural network designs have already made it possible to learn highly accurate interatomic force fields for a given compound with limited data by incorporating hard constraints in the form of symmetry properties and energy conservation into the model architecture [4, 5, 6, 7, 8, 9].

However, training models to have chemical accuracy for entire classes of compounds still requires large amounts of data [4, 5, 8, 9], owing in large part to the vastness of the space of possible atomic configurations. The creation of predictive and transferable models is, therefore, still a challenging task.

Additionally, the existing state-of-the-art models generally lack a measure of uncertainty in their prediction [3, 10], which makes it impossible to tell which of the predictions are reliable. This limits the applicability of active learning strategies, where uncertainty can be used to select useful training data from the large configuration space more efficiently [11, 12, 13, 14, 15] by only labeling outlier configurations where the model predicts a high likelihood of a large error.

Further, for a deployed model detected outliers can either be recomputed on-the-fly in DFT to ensure the accuracy of a neural network-based MD trajectory or used to retroactively discard trajectories that are likely to be physically inaccurate. All of these factors hinder their extension towards more general models, required for the prediction and optimization of complex (dynamic) properties of materials, where a vast amount of diverse atomic configurations need to be screened with an overall high accuracy and a good measure of uncertainty. The latter can be used to decide which of the samples that are predicted, should undergo additional investigation through _ab initio_ calculations or experiments and which are likely to be false positives. An uncertainty measure is also attractive to ensure high accuracy of modeling rare events that often govern longer time-scale dynamics, but are seriously under-represented in MD trajectories [14]. Configurations close to transition states are exponentially less likely to be sampled relative to near ground-state
configurations. But these states often have peculiar electronic configurations and bonding characteristics that do not feature in low energy configurations. Hence, an uncertainty measure would allow on-the-fly identification of configurations that require additional learning to bias active learning strategies towards relevant "new chemistry" such as transition states.

Furthermore, many existing state-of-the-art models are currently optimized in a somewhat intricate way on smaller datasets to improve performance, where artificial noise is introduced to the optimization process via a large learning rate and small batch size [4, 5, 6, 7, 9]. While this works well empirically, there is no established theoretical foundation for the utility of such optimization and consequently, there is no guide in choosing appropriate batch sizes and learning rates except trial and error. As a response to these challenges, we investigate here the viability of combining Bayesian Neural Networks (BNNs) with state-of-the-art neural networks for modeling interatomic forces, as BNNs have a robust measure of uncertainty and an inherent principled way of adding noise to the optimization process via Monte Carlo Markov Chain (MCMC) methods. Further, they offer a good framework for the incorporation and data-based updating of prior beliefs, such as approximate physics-guided interatomic force-fields, which may be utilized in the future to further improve data efficiency and generalization.

However, we found that current off-the-shelf MCMC methods are unsuited for dealing with the vastly different gradient scales that different parameter groups exhibit in modern models. This results in an unpractically slow traversal through the parameter space to the point where no convergence is achievable in reasonable amounts of time, which explains their current lack of use. Hence, a new MCMC algorithm is needed for high-quality Bayesian uncertainty quantification with neural network-based interatomic force models.

The contributions of this paper are the following:

* We develop a new stochastic neural network model based on the NguuIP architecture [7] for uncertainty-aware interatomic force modeling with state-of-the-art accuracy.
* We provide a novel Bayesian MCMC algorithm that produces high-accuracy models with high-quality uncertainty estimates in a practical amount of time without relying on finetuned batch sizes and learning rates.
* We show that the current optimization methods based on high learning rates and small batch sizes already mimic Bayesian MCMC methods to a degree.
* We discuss shortcomings as well as possible further improvements and extensions of the proposed neural network model and sampling algorithm.

Notation commonly used in this paper is summarized in Table 1.

### Related Work

While BNNs to date have found very few applications in this context, some use cases exist in the literature [15, 16, 17, 18]. However, in [16, 17] Monte Carlo dropout was used for Bayesian inference. While this method can in some instances be interpreted as Bayesian variational inference [19], it yields a poor approximation of the true posterior distribution even for variational inference-based methods [19] and typically results in much poorer uncertainty quantification than MCMC methods [20]. While a much more accurate Bayesian algorithm was used in [15], the use case was very limited in scope to diatomic systems. Further, all of the above instances were limited to very small fully connected neural networks while we utilize a state-of-the-art graph neural network-based architecture. We only found one attempt to reconcile graph neural network-based architectures with sampling schemes based on simulating a stochastic process in the literature. This was done in a recently published preprint by Thaler et al. [18]. However, they failed to get competitive results while attempting to sample the Bayesian posterior and had to resort to sampling an artificially sharpened distribution. Even then they found that to reliably quantify uncertainty, they had to generate samples from multiple different stochastic processes while we generate samples from a single Markov chain.

## 2 The Base Model

### The Formal Setting

The aim of the neural network-based interatomic force model is to map a configuration of atoms \(\boldsymbol{x}=\{(\boldsymbol{r}_{1},z_{1}),...,(\boldsymbol{r}_{n},z_{n})\}\), where \(\boldsymbol{r}_{i}\) denotes the position of the nucleus of atom \(i\) and \(z_{i}\) refers to its nuclear charge, to the forces \(\{\boldsymbol{F}_{1},...,\boldsymbol{F}_{n}\}\) acting on the individual nuclei.

\begin{table}
\begin{tabular}{|c|c|} \hline
**Notation** & **Meaning** \\ \hline Bold case symbol & Vector-valued quantities \\ \(N(\boldsymbol{\mu},\Sigma)\) & Normal distribution \\ \(\boldsymbol{N}(\boldsymbol{\mu},\Sigma)\) & Random variable \(N(\boldsymbol{\mu},\Sigma)\) \\ \(\boldsymbol{F}\) & Potential energy \\ \(D\) & The dataset \\ \(l\) & The identity matrix \\ \(\boldsymbol{\theta}\) & Neural network parameters \\ \(\boldsymbol{x}\) & Independent variable \\ \(y\) & Dependent variable \\ \(\boldsymbol{y}|\boldsymbol{x}\) & \(\boldsymbol{y}\) conditioned on \(\boldsymbol{x}\) \\ \(r\) & Atomic coordinates \\ \(z\) & Nuclear charges \\ \(\boldsymbol{ab}\) & Element-wise product of \(\boldsymbol{a}\) \\  & and \(\boldsymbol{b}\) \\ \hline \end{tabular}
\end{table}
Table 1: Commonly used notation
These forces are then used to simulate the movement of the nuclei based on Newton's equations of motion via a solver for Ordinary Differential Equations (ODEs) (Figure 1).

One important aspect of this form of data is, that it is unstructured, meaning the order in which the atoms are enumerated is arbitrary. A suitable neural network model should therefore be equivariant under a reordering of the data.

Another important constraint is, that these forces are conservative and can therefore be derived as the negative gradients of a single potential energy surface

\[\boldsymbol{F}_{i}=-\nabla_{\boldsymbol{r}_{i}}E((\boldsymbol{r}_{1},z_{1}),...,(\boldsymbol{r}_{n},z_{n})).\]

The potential energy surface itself is invariant under any distance preserving transformation of the atomic coordinates. This was at first incorporated into neural network models by using only the interatomic distances \(r_{ij}\) and nuclear charges as input, which have these invariances themself. Even though it is in principle possible to extract directional information from the set of interatomic distances, in practice incorporating directional information explicitly can improve data efficiency and accuracy quite a lot [10] and has become a key feature of many state-of-the-art neural network models [4, 5, 6, 8, 7, 9].

### The NequIP Model

One of the most powerful neural network models for interatomic force modeling that currently exists is the NequIP model [7] (Figure 2). This model takes all the previous considerations into account and consists of a model base and projection layers. The model base maps an atomic configuration \(\{(\boldsymbol{r}_{1},z_{1}),...,(\boldsymbol{r}_{n},z_{n})\}\) to a set of high dimensional latent feature vectors \(\{\boldsymbol{v}_{1},...,\boldsymbol{v}_{n}\}\) that are invariant under distance-preserving transformations. The projection layers then map \(\{(\boldsymbol{v}_{1},z_{1}),...,(\boldsymbol{v}_{n},z_{n})\}\) to (virtual) atomic energies \(\{E_{1},...,E_{n}\}\).

The potential energy \(E\) is then calculated as the sum of the atomic energies. Finally, the forces acting on the nuclei are then calculated as the negative gradients of the potential energy with respect to the nuclear coordinates. The last step can be achieved straightforwardly using the automatic differentiation functionalities of modern deep-learning libraries.

## 3 The Bayesian Model

### Bayesian Neural Networks

The goal of classical neural network optimization is to find a single set of parameters, i.e. weights and biases, which models the underlying data distribution, from which the training data was sampled, well. In the Bayesian approach, in contrast, the parameters of the neural network are modeled probabilistically. To keep the notation simple, we use \(\boldsymbol{\theta}\) to denote a vector containing a complete set of neural network parameters. For BNNs it is assumed that some prior knowledge exists about what constitutes a good set of parameters, which is expressed in the form of a prior density \(p(\boldsymbol{\theta})\). This prior density then gets refined through the training data \(D=\{(\boldsymbol{x}_{1},\boldsymbol{y}_{1}),...,(\boldsymbol{x}_{l}, \boldsymbol{y}_{l})\}\) by invoking Bayes rule for calculating the posterior density:

\[p(\boldsymbol{\theta}|D)=\frac{p(D|\boldsymbol{\theta})p(\boldsymbol{\theta}) }{p(D)}=\frac{p(D|\boldsymbol{\theta})p(\boldsymbol{\theta})}{\int p(D| \boldsymbol{\theta})p(\boldsymbol{\theta})d\boldsymbol{\theta}}.\]

A prediction on a new data point \((\boldsymbol{x},\boldsymbol{y})\) can now be made via:

\[p(\boldsymbol{y}|\boldsymbol{x},D)=\int p(\boldsymbol{y}|\boldsymbol{x}, \boldsymbol{\theta})p(\boldsymbol{\theta}|D)d\theta\]

\[=\mathbb{E}_{p(\boldsymbol{\theta}|D)}\left[p(\boldsymbol{y}|\boldsymbol{x}, \boldsymbol{\theta})\right].\]

For BNNs this integral is almost always analytically intractable. However, by generating samples \(\boldsymbol{\theta}_{1},...,\boldsymbol{\theta}_{k}\) from \(p(\boldsymbol{\theta}|D)\) it can be estimated through the law of large numbers as:

\[p(\boldsymbol{y}|\boldsymbol{x},D)\approx\frac{1}{k}\sum_{i=1}^{k}p( \boldsymbol{y}|\boldsymbol{x},\boldsymbol{\theta}_{i}).\]

### The Stochastic Model

In order to achieve state-of-the-art accuracy combined with a good measure of predictive uncertainty in modeling interatomic forces, a new stochastic neural network model is required, which we will introduce in this section. To build a stochastic model of the data, we make the following assumptions on the conditional independence of the data:

\[p(D) =p(\boldsymbol{y}_{1},...,\boldsymbol{y}_{l}|\boldsymbol{x}_{1},...,\boldsymbol{x}_{l})p(\boldsymbol{x}_{1},...,\boldsymbol{x}_{l})\] \[=p(\boldsymbol{x}_{1},...,\boldsymbol{x}_{l})\Pi_{i=1}^{l}p( \boldsymbol{y}_{i}|\boldsymbol{x}_{i}).\]

Figure 1: An illustration of a molecular dynamics workflow with neural network predictions for the forces.
\(p(\mathbf{y}_{i}|\mathbf{x}_{i})\) will be inferred by a neural network as \(p(\mathbf{y}_{i}|\mathbf{x}_{i},\mathbf{\theta})\) while \(p(\mathbf{x}_{1},...,\mathbf{x}_{l})\) depends on the data generation process.

More specifically we model the conditional densities

\[p(\mathbf{y}|\mathbf{x},\mathbf{\theta})=p(E,\mathbf{F}_{1},...,\mathbf{F}_{n}|(\mathbf{r}_{1},z_{1}),...,(\mathbf{r}_{n},z_{n}),\mathbf{\theta})\]

as

\[\mathbf{y}|\mathbf{x},\mathbf{\theta}\sim N(\mu_{E }(\mathbf{\theta},\mathbf{x}),\sigma_{E}^{2}(\mathbf{\theta},\mathbf{x}))\]

\[\times\Pi_{i=1}^{n}N(\mathbf{\mu}_{Fi}(\mathbf{\theta}, \mathbf{x}),\sigma_{\mathbf{F}_{i}}^{2}(\mathbf{ \theta},\mathbf{x})I),\]

where \(\mathbf{y}|\mathbf{x},\mathbf{\theta}\) denotes \(y\) conditioned on \(x\) and \(\theta\), \(\mathbf{x}=\{(\mathbf{r}_{1},z_{1}),...,(\mathbf{r }_{n},z_{n})\}\), the means \(\mu_{E}(\mathbf{\theta},\mathbf{x})\) and \(\mathbf{\mu}_{Fi}(\mathbf{\theta},\mathbf{x})\) are the regular point predictions of the NequIp model, the standard deviations \(\sigma_{E}(\mathbf{\theta},\mathbf{x})\) and \(\sigma_{\mathbf{F}_{i}}(\mathbf{\theta},\mathbf{x})\) are predicted by two separate Multi-Layer Perceptrons (MLPs) from the outputs of the NequIp base \(\mathbf{v}_{1},...,\mathbf{v}_{n}\) and \(I\) denotes the identity matrix (Figure 3). Note that a single invariant standard deviation is calculated for all three force components to ensure equivariance of the predicted density under any distance-preserving transformation of the atomic coordinates.

We include the potential energy \(E\) here as an additional dependent variable because it can be useful in applications such as the determination of ground state configurations and is calculated as a byproduct of force calculations in DFT anyway. Lastly, we use the simple Gaussian mean field prior \(\mathbf{\theta}\sim N(\mathbf{0},I)\).

### Sampling from the Posterior Distribution

Generating high-quality sample weights from the posterior distribution is usually done via the simulation of a Markov chain which converges in distribution to the posterior. While several Markov chains methods have been constructed for neural network applications, we found them unsuitable for this use case. The main difficulty we encountered was, that the scale of the gradients with respect to the different sets of parameters corresponding to the different layers of the neural network vary by several orders of magnitude. This makes the use of a single step size for all parameters impossible, as it would cause the sets of parameters with smaller gradients to be frozen during the optimization. While modern neural network optimizers circumvent this problem through an adaptive step size for each parameter [21], this is only possible to a very limited degree for Markov chains [22] without changing the stationary distribution to which they converge. Furthermore, this typically requires the calculation of higher-order derivatives [22] which is computationally expensive. To deal with these challenges, we develop here a new approach to sampling the posterior distribution based on the Stochastic Gradient Hamiltonian Monte Carlo (SGHMC) algorithm introduced by Chen et al. [23]. In its basic form, the algorithm is given by the Markov chain

\[\Delta\mathbf{w}_{t}=-\left(\nabla_{\mathbf{ \theta}_{t}}u(\mathbf{\theta}_{t})+\mathbf{M}^{-1}\mathbf{C}\mathbf{w}_{t-1}\right)\Delta t+\sqrt{2\mathbf{C$ }\Delta t}\mbox{\boldmath$\mathcal{N}}_{t}(\mathbf{0},I),\] \[\Delta\mathbf{\theta}_{t}=\mathbf{M}^{-1}\mathbf{w}_{t}\Delta t.\]

where

\[u(\mathbf{\theta}):=-\ln p(\mathbf{\theta})-\sum_{i=1}^{l} \ln p(\mathbf{y}_{i}|\mathbf{x}_{i},\mbox{\boldmath$\theta$ }),\]

\(\mathbf{M}^{-1}\) and \(C\) are vectors containing strictly positive values, \(\Delta t<<1\) is the step size and \(\mathbf{\mathcal{N}}_{t}(\mathbf{0},I)\) denotes a random variable with a multivariate standard normal distribution. To keep the notation simple, we use the convention here, that all operations on vectors (multiplication, inversion, etc.) are to be taken element-wise. The auxiliary variable \(\mathbf{w}_{t}\) has the same dimension as \(\mathbf{\theta}_{t}\).

A physical intuition of SGHMC can be gained by noticing that these equations also describe the diffusion of particles with coordinates \(\theta\), potential energy \(u(\mathbf{\theta})\), momenta \(\mathbf{w}_{t}\) and masses \(M\) under friction. In this analogy, the vector \(C\) would represent the friction coefficients. From a machine learning perspective, \(\mathbf{w}_{t}\) represents a momentum term similar to the ones used in many modern neural network optimizers [21]. In fact, Chen et al. used some substitutions which lead to a Bayesian analog to stochastic gradient descent with momentum [23], but we will make somewhat different substitutions that lead to a natural way to include adaptive step sizes: \(\alpha=\Delta t\mathbf{M}^{-1}\mathbf{C}\), \(\gamma=(\Delta t)^{2}|D|\alpha^{-1}\) and \(\mathbf{v}_{t}=\gamma^{-1}\Delta t\mathbf{w}_{t}\) which yields:

\[\Delta\mathbf{v}_{t}=-\alpha\frac{1}{|D|}\nabla_{\mathbf{ \theta}_{t}}u(\mathbf{\theta}_{t})-\alpha\mathbf{v}_{t-1}+ \alpha\sqrt{\frac{2\mathbf{M}}{|D|\gamma}}\mathbf{\mathcal{N} }_{t}\left(0,I\right),\]

\[\Delta\mathbf{\theta}_{t}=\gamma\mathbf{M}^{-1}\mathbf{v}_{t}\]

where \(|D|\) is the size of the dataset \(D\).

This is up to the noise term \(\mathbf{\mathcal{N}}_{t}\left(0,I\right)\) equivalent to the updates of the Adam optimizer [24].

Figure 2: The computational graph of the NequIp model for predicting the potential energy \(E\) and atomic forces \(\mathbf{F}_{i}\) from the atomic numbers \(z_{i}\) and coordinates \(\mathbf{r}_{i}\).
Unfortunately, the mass term \(\mathbf{M}\) used in the standard Adam optimizer can vary a lot even during later stages of the training and this variation depends on the value of \(\mathbf{\theta}_{t}\) at the previous time steps. As a consequence, the resulting process \((\mathbf{\theta}_{t},\mathbf{v}_{t})\) would not even be a Markov chain, and the Bayesian posterior would most likely no longer be the distribution \(\mathbf{\theta}_{t}\) converges to [22]. For many neural network architectures, timely convergence to the posterior can be achieved by simply setting \(\mathbf{M}=I\). However, the vastly varying gradient scales of the different parameter groups in the model make this approach not feasible. As a solution, we introduce now a new adaptive step size method for the SGHMC algorithm which still converges to the posterior distribution without requiring the computation of higher-order derivatives. In order to achieve this, we set \(\mathbf{M}\) as the denominator of the AMSGrad algorithm [25] during the first phase of the optimization:

\[\mathbf{a}_{t}=(1-\beta)\frac{1}{|D|^{2}}\left(\nabla_{\mathbf{\theta}_{t}}u(\mathbf{ \theta}_{t})\right)\left(\nabla_{\mathbf{\theta}_{t}}u(\mathbf{\theta}_{t})\right)+ \beta\mathbf{a}_{t-1},\]

\[\mathbf{D}_{t}=\max(\mathbf{D}_{t-1},\mathbf{a}_{t}),\]

\[\mathbf{M}_{t}=\sqrt{\frac{\mathbf{D}_{t}}{1-\beta^{t}}}+\epsilon_{stability}.\]

Here, \(\max(\cdot,\cdot)\) is the element-wise maximum, \(\epsilon_{stability}\) is a stability constant and \(\beta\) is a hyperparameter used in computing the running average \(\mathbf{a}_{t}\) and is typically set between 0.99 and 0.9999.

Because this mass term is still time- and path-dependent, \(\mathbf{\theta}_{t}\) can not, in general, be expected to converge exactly to the posterior distribution during this phase. However, because \(\mathbf{M}_{t}\) is not based on a running average of squared gradients, like most adaptive step size methods are, but instead on the maximum of such a running average, it typically changes very little during the later stages of the optimization. As a result, the process will already become fairly close in distribution to the Bayesian posterior during this stage of the optimization. Furthermore, because \(\mathbf{M}_{t}\) already remains almost constant after a while, we can keep it entirely constant after a certain amount of steps without causing instabilities, at which point the process \(\mathbf{\theta}_{t}\) becomes a regular SGHMC process which is known to converge to the Bayesian posterior for sufficiently small step sizes [23].

Because deep learning datasets are usually very large, evaluating \(\nabla_{\mathbf{\theta}}u(\mathbf{\theta})\) exactly is typically very time-consuming. Practical implementations instead estimate \(\nabla_{\mathbf{\theta}}u(\mathbf{\theta})\) on a smaller, randomly sampled subset \(\hat{D}\subset D\) via

\[\nabla_{\mathbf{\theta}}u(\mathbf{\theta})=\mathbf{\epsilon}(\mathbf{\theta})+\nabla_{\mathbf{ \theta}}\hat{u}(\mathbf{\theta})\]

\[:=\mathbf{\epsilon}(\mathbf{\theta})+\nabla_{\mathbf{\theta}}\left(-\ln p(\mathbf{\theta})- \frac{|D|}{|\hat{D}|}\sum_{(\mathbf{x},\mathbf{y})\in D}\ln p(\mathbf{y}|\mathbf{x},\mathbf{\theta })\right).\]

The resulting algorithm is summarized in Algorithm 1. Here \(\mathbf{\epsilon}(\mathbf{\theta})\) is the error in the estimation of \(\nabla_{\mathbf{\theta}}u(\mathbf{\theta})\) with \(\mathbb{E}\left[\mathbf{\epsilon}(\mathbf{\theta})\right]=\mathbf{0}\).

```
for\(t=1\)to\(T\)do  Sample minibatch \(\hat{D}\)  Evaluate \(\nabla_{\mathbf{\theta}_{t}}\hat{u}(\mathbf{\theta}_{t})\)  = \(\nabla_{\mathbf{\theta}_{t}}\left(-\ln p(\mathbf{\theta}_{t})-\frac{|D|}{|\hat{D}|} \sum_{(\mathbf{x},\mathbf{y})\in\hat{D}}\ln p(\mathbf{y}|\mathbf{x},\mathbf{\theta}_{t})\right)\) if\(t\leq t_{max}\)then \(\mathbf{a}_{t}=(1-\beta)\frac{1}{|D|^{2}}\left(\nabla_{\mathbf{\theta}_{t}}\hat{u}(\mathbf{ \theta}_{t})\right)\left(\nabla_{\mathbf{\theta}_{t}}\hat{u}(\mathbf{\theta}_{t}) \right)+\beta\mathbf{a}_{t-1}\) \(\mathbf{D}_{t}=\max(\mathbf{D}_{t-1},\mathbf{a}_{t})\) \(\mathbf{M}_{t}=\sqrt{\frac{\mathbf{D}_{t}}{1-\beta^{t}}}+\epsilon_{stability}\) else \(\mathbf{M}_{t}=\mathbf{M}_{t_{max}}\) endif \(\Delta\mathbf{v}_{t}=-\frac{\alpha}{|D|}\nabla_{\mathbf{\theta}_{t}}\hat{u}(\mathbf{\theta} _{t})-\alpha\mathbf{v}_{t-1}+\alpha\sqrt{\frac{2\mathbf{M}_{t}}{|D|\gamma}}\mathbf{N}_{t} \left(0,I\right)\) \(\Delta\mathbf{\theta}_{t}=\gamma\mathbf{M}_{t}^{-1}\mathbf{v}_{t}\) endfor
```

**Algorithm 1** The Stochastic Optimization Algorithm

The resulting algorithm is summarized in Algorithm 1. Here \(\mathbf{\epsilon}(\mathbf{\theta})\) is the error in the estimation of \(\nabla_{\mathbf{\theta}}u(\mathbf{\theta})\) with \(\mathbb{E}\left[\mathbf{\epsilon}(\mathbf{\theta})\right]=\mathbf{0}\).

As long as \(\frac{|\mathbf{\epsilon}(\mathbf{\theta})|}{|D|}<<\sqrt{\frac{2\mathbf{M}_{t}}{\gamma|D|}}\) it is clear that this additional noise does not have a large effect on the dynamic as the Gaussian noise term will dominate.

```
for\(t=1\)to\(T\)do  Sample minibatch \(\hat{D}\)  Evaluate \(\nabla_{\mathbf{\theta}_{t}}\hat{u}(\mathbf{\theta}_{t})\)  = \(\nabla_{\mathbf{\theta}_{t}}\left(-\ln p(\mathbf{\theta}_{t})-\frac{|D|}{|\hat{D}|} \sum_{(\mathbf{x},\mathbf{y})\in\hat{D}}\ln p(\mathbf{y}|\mathbf{x},\mathbf{\theta}_{t})\right)\) if\(t\leq t_{max}\)then \(\mathbf{a}_{t}=(1-\beta)\frac{1}{|D|^{2}}\left(\nabla_{\mathbf{\theta}_{t}}\hat{u}(\mathbf{\theta}_{t}) \right)\left(\nabla_{\mathbf{\theta}_{t}}\hat{u}(\mathbf{\theta}_{t})\right)+\beta\mathbf{a }_{t-1}\) \(\mathbf{D}_{t}=\max(\mathbf{D}_{t-1},\mathbf{a}_{t})\) \(\mathbf{M}_{t}=\sqrt{\frac{\mathbf{D}_{t}}{1-\beta^{t}}}+\epsilon_{stability}\) else \(\mathbf{M}_{t}=\mathbf{M}_{t_{max}}\) endif \(\Delta\mathbf{v}_{t}=-\frac{\alpha}{|D|}\nabla_{\mathbf{\theta}_{t}}\hat{u}(\mathbf{\theta}_ {t})-\alpha\mathbf{v}_{t-1}+\alpha\sqrt{\frac{2\mathbf{M}_{t}}{|D|\gamma}}\mathbf{N}_{t} \left(0,I\right)\) \(\Delta\mathbf{\theta}_{t}=\gamma\mathbf{M}_{t}^{-1}\mathbf{v}_{t}\) endfor
```

**Algorithm 2** The Stochastic Optimization Algorithm

The can always be achieved by choosing an appropriate
batch size and step size \(\gamma\).

However, the additional noise will become significant for large learning rates and small batch sizes.

In fact, even if we do not add the Gaussian noise to this process and subsequently just use the resulting AMS-Grad optimization, the injected noise at the commonly used large step sizes and small batch sizes will cause the optimization process to become a stochastic process which imitates the proposed algorithm to a degree, only with a noise term \(\frac{|\mathbf{w}(\mathbf{\theta})|}{|D|}\) that is not grounded on any theoretical foundation.

## 4 Empirical Evaluation

### The Benchmarks

To assess our model's accuracy and predicted uncertainty we utilize two different datasets. To analyze our model's performance with varying amounts of Monte Carlo samples and under domain shift, we used a dataset consisting of PEDOT polymers, which are conducting polymers, which due to their chemical stability and tuneable properties [26] have been used for a wide range of applications including sensors [27], supercapacitors [28], battery electrodes [29], bioelectronics, solar cells, electrochromic displays, electrochemical transistors [30], and spintronics [31]. The dataset consists of polymers of lengths 8, 12 and 16, where only the shorter chains were used for training (see appendix B for details). The total training set contained only 100 configurations, 50 of length 8 polymers and 50 of length 12 polymers. Equivalently a small validation set of size 30 was constructed.

The second dataset used is the MD17 dataset [32] consisting of long molecular dynamics trajectories of several small organic molecules calculated in DFT. Because this dataset had some numerical errors, a revised version RMD17 has recently been computed at higher accuracy [33] which we will utilize to evaluate our model's performance. Because dropout is the most common Bayesian method used for interatomic force modeling, we will compare our model's performance to a dropout-based version of the proposed stochastic model on this dataset (see appendix A.4 for details).

For each compound, we used 1000 randomly sampled configurations during training and the rest for testing. Of the 1000 configurations sampled, 30 were reserved as a small validation set, and the remaining ones comprised the actual training set. The same samples were used as training, validation and test sets for our proposed model and the dropout-based model.

### Evaluation Metrics

To measure prediction accuracy, we evaluate both the Mean Average Error in the Forces (MAE-F) and Energies (MAE-E) as well as the Root Mean Square Errors (RMSE-F, RMSE-E). In the evaluation, we use the expectation value under the estimated posterior density as a point prediction. I.e. \(\mathbf{F}_{pred}(\mathbf{x})=\frac{1}{k}\sum_{i=1}^{k}\mathbf{\mu}_{\mathbf{F}}(\mathbf{x},\mathbf{ \theta}_{i})\) and \(E_{pred}(\mathbf{x})=\frac{1}{k}\sum_{i=1}^{k}\mu_{E}(\mathbf{x},\mathbf{\theta}_{i})\). One metric used to evaluate and compare the predicted uncertainties is the Mean Log Likelihood (MLL). Since in the most common applications, the main task of the uncertainty measure is outlier detection, we further evaluate the ROC AUC scores for detecting force components with an error larger than 1 kcal/(mol\(\cdot\)A) on the basis of the variance of the predicted distribution of those force components. Because of the relatively small size of the PEDOT dataset, we only evaluate the outlier detection on the much larger RMD17 datasets.

Lastly, to evaluate if the proposed model is properly calibrated we utilize the Expected Calibration Er

Figure 4: Receiver Operating Characteristic curves for uncertainty-based detection of force components with a prediction error of at least 1 kcal/(mol\(\cdot\)Å) using k=8 Monte Carlo samples on the RMD17 datasets.
ror (ECE) [34] as well as a visual inspection of the predicted and observed error densities of the force components (Figures 4 and 5) (see appendix A.5 for details). To calculate an ECE, the predictions \(y_{i}\) are divided into \(m\) bins of equal width. The ECE is then calculated as

\[ECE=\sum_{i=1}^{m}f_{i}\left|f_{i}-e_{i}\right|\]

where \(f_{i}\) is the fraction of observed samples that fall into bin \(i\) and \(e_{i}\) is the predicted probability of a sample falling into bin \(i\). Because the ECE is highly dependent on the width of the bins, we divide the ECE by the bin width \(\delta\) to obtain a Normalized Expectation Calibration Error (NECE):

\[NECE=\sum_{i=1}^{m}\frac{f_{i}}{\delta}\left|f_{i}-e_{i}\right|\]

which is an estimate of

\[\mathbb{E}_{\rho(y)}\left[\left|\rho(y)-\rho_{predicted}(y)\right|\right].\]

### Results

**Results on the PEDOT Datasets**

As can be seen in Table 2, the model achieves high accuracy in the force prediction on all three test sets with a decreasing accuracy for increasing chain lengths. A slight improvement in the accuracy is observed when going from one to eight Monte Carlo samples. The model has a much larger error in the auxiliary task of energy predictions when compared to the force predictions. This is not surprising though, as the PEDOT molecules are fairly large and have correspondingly large energy scales. Additionally, the PEDOT training set contains much fewer energy labels than force labels. Further, it can be seen from Table 3, that a single Monte Carlo sample does not yield good uncertainty estimates for the forces which vastly improves when using eight Monte Carlo samples. A similar trend is not evident for the uncertainty quantification of the energy predictions.

Even though no polymer chains of length 16 were included in the training set, we find that the model still achieves high accuracy and a good quantification of the uncertainty in the force predictions (Table 3) in the case of k=8 Monte Carlo samples and again much poorer uncertainty quantification for a single Monte Carlo sample. A complete histogram of predicted and actual force components can be found in appendix C. We observe some overconfidence even in the k=8 case (Figure 6).

**Influence of the Number of Monte Carlo Samples**

As can be seen in Table 3, both the calibration as

Figure 5: Observed and predicted error densities of the force components in kcal/(mol.Å) with k=8 Monte Carlo samples on the RMD17 datasets.

Figure 6: Observed and predicted error densities of the force components in kcal/(mol.Å) with k=8 Monte Carlo samples on the length 16 PEDOT dataset.
well as the MLL significantly improve with increasing sample size for the forces. Because these samples are generated from a single Markov chain, this demonstrates an efficient traversal of the parameter space by the proposed sampling algorithm. Further, it appears that the improvements in MLLs diminish exponentially with increasing sample size. No similar trend can be found for the auxiliary task of energy predictions. Overall we find \(k=8\) Monte Carlo samples to be a good tradeoff between quality of uncertainty quantification and computational complexity.

**Results on the RMD17 Datasets**

Using eight Monte Carlo samples for both models, the proposed model has comparable accuracy to the dropout model on many of the easier RMD17 datasets but significantly outperforms it on the aspirin and malonaldehyde datasets (Table 4). The achieved accuracies are very consistent with the original NequIP model on these datasets [7] (See appendix D for the results of a single Monte Carlo sample and the NequIP model on this dataset). Further, it consistently outperforms the dropout model in terms of mean log-likelihoods for both the forces and the energies. The only exception to this being the mean log-likelihoods of the uracil energies where both models achieve the same result.

A very large difference in the performance of the models is found in the outlier detection task where the proposed model consistently achieves ROC AUC scores much closer to the optimal score of 1 (Table 5). The benzene dataset was not included in this comparison because there were no instances of force prediction errors of the necessary scale for the proposed model. A complete plot of the receiver operating characteristic curves is given in Figure 4. As can be seen there, our proposed model is much more reliable at detecting outliers.

Lastly, as is evident from Table 6 and Figure 5 the resulting model is not accurately calibrated in the case of 8 Monte Carlo samples and has a tendency for overconfidenc

\begin{table}
\begin{tabular}{c c|c c c c|c c c} \hline \hline \multirow{2}{*}{Chain Length} & \multicolumn{5}{c|}{Proposed Model (k=1)} & \multicolumn{5}{c}{Proposed Model (k=8)} \\  & & MAE & RMSE & NICE & MIL & MAE & RMSE & NICE & Mul \\ \hline \multirow{4}{*}{16} & Energy & 15.86 & 32.22 & 0.009 & -4.57 & 14.31 & 22.87 & 0.006 & -4.33 \\  & Forces & 0.222 & 0.385 & 0.649 & -1.79 & 0.210 & 0.371 & 0.313 & -0.20 \\  & Energy & 5.40 & 6.79 & 0.023 & -3.70 & 5.11 & 6.24 & 0.022 & -3.68 \\  & Forces & 0.184 & 0.276 & 0.652 & -0.87 & 0.169 & 0.254 & 0.317 & 0.14 \\  & Energy & 4.15 & 5.27 & 0.044 & -3.59 & 5.37 & 6.54 & 0.045 & -3.65 \\
8 & forces & 0.151 & 0.221 & 0.619 & -0.284 & 0.139 & 0.203 & 0.215 & 0.44 \\ \hline \hline \end{tabular}
\end{table}
Table 2: Results on the PEDOT datasets. Force in kcal/(mol\(\cdot\)Å), energy in kcal/mol.

\begin{table}
\begin{tabular}{c c|c c c|c c} \hline \hline \multirow{2}{*}{
\begin{tabular}{c} Number of Samples k \\ \end{tabular} } & \multicolumn{5}{c}{Proposed Model} \\  & & MAE & RMSE & NICE & MILL \\ \hline \multirow{4}{*}{1} & Energy & 15.86 & 23.22 & 0.009 & -4.57 \\  & Forces & 0.222 & 0.385 & 0.649 & -1.79 \\  & Energy & 15.16 & 23.00 & 0.008 & -4.44 \\  & Forces & 0.219 & 0.381 & 0.473 & -0.87 \\  & Energy & 14.36 & 22.88 & 0.007 & -4.34 \\  & Forces & 0.216 & 0.378 & 0.376 & -0.44 \\  & Energy & 14.31 & 22.87 & 0.006 & -4.33 \\  & Forces & 0.210 & 0.371 & 0.313 & -0.20 \\
16 & Energy & 14.56 & 22.86 & 0.007 & -4.38 \\  & Forces & 0.201 & 0.356 & 0.254 & 0.00 \\ \hline \hline \end{tabular}
\end{table}
Table 3: Influence of the MC sample size on the length 16 PEDOT dataset. Force in kcal/(mol\(\cdot\)Å), energy in kcal/mol.

\begin{table}
\begin{tabular}{c c|c c|c c c} \hline \hline \multirow{2}{*}{
\begin{tabular}{c} Molecule \\ \end{tabular} } & \multicolumn{5}{c|}{Dropout} & \multicolumn{2}{c}{Proposed Model} \\  & & MAE & RMSE & MILL & MAE & RMSE & NICE \\ \hline \multirow{4}{*}{Aspirin} & Energy & 0.072 & 0.096 & 0.88 & 0.045 & 0.067 & 1.41 \\  & Forces & 0.215 & 0.340 & -1.85 & 0.144 & 0.229 & 0.22 \\ \cline{2-6}  & Energy & 0.012 & 0.021 & 2.61 & 0.010 & 0.017 & 2.99 \\  & Forces & 0.078 & 0.147 & 0.74 & 0.064 & 0.115 & 1.03 \\  & Energy & 0.012 & 0.020 & 2.63 & 0.013 & 0.020 & 2.63 \\  & Forces & 0.082 & 0.141 & 0.71 & 0.085 & 0.138 & 0.83 \\  & Energy & 0.021 & 0.038 & 1.82 & 0.016 & 0.027 & 2.51 \\  & Forces & 0.145 & 0.259 & 0.18 & 0.099 & 0.170 & 0.49 \\ \hline \multirow{4}{*}{Salicylic acid} & Energy & 0.022 & 0.042 & 1.87 & 0.023 & 0.213 \\  & Forces & 0.105 & 0.199 & 0.14 & 0.109 & 0.132 & 0.63 \\ \cline{1-1}  & Energy & 0.009 & 0.018 & 2.55 & 0.009 & 0.012 & 2.74 \\ \cline{1-1}  & Forces & 0.044 & 0.060 & 1.319 & 0.044 & 0.060 & 1.63 \\ \cline{1-1}  & Energy & 0.011 & 0.015 & 2.63 & 0.010 & 0.014 & 2.79 \\ \cline{1-1}  & Forces & 0.029 & 0.084 & 1.11 & 0.051 & 0.082 & 1.43 \\ \cline{1-1}  & Benzene & Energy & 0.005 & 0.006 & 3.20 & 0.004 & 0.005 & 3.53 \\ \cline{1-1}  & Forces & 0.020 & 0.032 & -1.50 & 0.010 & 0.016 & 3.06 \\ \hline \hline \end{tabular}
\end{table}
Table 4: Mean results on the RMD17 dataset using \(k=8\) Monte Carlo samples. Force in kcal/(mol\(\cdot\)Å), energy in kcal/mol.

\begin{table}
\begin{tabular}{c|c|c} \hline \hline Molecule & Dropout & Proposed Model \\  & & Model \\ Aspirin & 0.800 & 0.952 \\ Ethanol & 0.897 & 0.993 \\ Uracil & 0.888 & 0.986 \\ Malonaldehyde & 0.824 & 0.983 \\ Salicylic acid & 0.809 & 0.977 \\ Naphthalene & 0.902 & 0.995 \\ Toluene & 0.894 & 0.995 \\ \hline \hline \end{tabular}
\end{table}
Table 5: ROC AUC scores for uncertainty-based detection of force components with a prediction error of at least 1 kcal/(mol\(\cdot\)Å) using k=8 Monte Carlo samples on the RMD17 datasets.
## 5 Conclusion and Outlook

As the results demonstrate, the proposed Bayesian neural network model achieves state-of-the-art accuracy while also achieving good uncertainty quantification. The proposed sampling algorithm appears to converge to the posterior distribution in a reasonable amount of time and maintains a fairly swift traversal through the parameter space afterward. Most importantly, the resulting model can reliably detect outliers while sampling parameters from the same Markov chain. This opens up new possibilities for Bayesian active learning procedures for learning interatomic forces. Further, it could enable the principled incorporation of existing physics-guided approximate interatomic force fields via the Bayesian prior distribution. This might potentially further improve both the data efficiency and generalizability of the model.

While the proposed model already achieves good results, some improvements could still be made.

While convergence to the posterior distribution becomes feasible with the proposed algorithm, in practice we find that it still takes about twice as many GPU hours to reach convergence when compared to the non-Bayesian optimization procedure described in [7]. However, while it was useful for demonstrative purposes of the convergence properties of the proposed sampler, faster convergence can most likely be achieved by simply reducing the injected gradient noise and the batch size during the initial phases of the optimization. Nevertheless, even without these measures, the generation of interatomic forces for the training set at _ab initio_ accuracy will likely be the computational bottleneck in most applications.

Further, the stochastic model could potentially still be improved by adequate incorporation of covariances between atoms and also between force components. However, both of these covariances are challenging to include. For covariance between force components, the main challenge is to maintain rotation equivariance of the predicted density. For the incorporation of interatomic covariances, a dense covariance matrix will very quickly become impractical for larger molecules, as the evaluation of the log-likelihoods becomes a computational bottleneck. While a sparse covariance matrix might be adequate due to the mostly local nature of interatomic forces, a way to parameterize sparse covariance matrices would be needed, which is not trivial and which we could not find in the existing literature.

Lastly, we found that the predicted uncertainties are not always properly calibrated and have a tendency for over-confidence. This might be a result of sampling from the same Markov chain, where the models can never be completely independently sampled from each other. This can lead to a reduced predictive variance between the sampled models and hence increased confidence. Here it might be beneficial to utilize two or more parallel chains to reduce the sample interdependency and explore the posterior distribution more efficiently or to recalibrate the uncertainties on a validation set.

## 6 Acknowledgement

This paper is funded by dtec.bw - Digitalization and Technology Research Center of the Bundeswehr which we gratefully acknowledge [project CoupleIT!]
