[1]Rafael Gómez-Bombarelli

1]Department of Materials Science and Engineering, Massachusetts Institute of
Technology

2]Innovative Technology Laboratories, AGC Inc.

3]Computational and Systems Biology, Massachusetts Institute of Technology

# Single-model uncertainty quantification in neural network potentials does
not consistently outperform model ensembles

Aik Rui Tan  Shingo Urata  Samuel Goldman  Johannes C. B. Dietschreit
<EMAIL_ADDRESS>[  [  [

###### Abstract

Neural networks (NNs) often assign high confidence to their predictions, even
for points far out-of-distribution, making uncertainty quantification (UQ) a
challenge. When they are employed to model interatomic potentials in materials
systems, this problem leads to unphysical structures that disrupt simulations,
or to biased statistics and dynamics that do not reflect the true physics.
Differentiable UQ techniques can find new informative data and drive active
learning loops for robust potentials. However, a variety of UQ techniques,
including newly developed ones, exist for atomistic simulations and there are
no clear guidelines for which are most effective or suitable for a given case.
In this work, we examine multiple UQ schemes for improving the robustness of
NN interatomic potentials (NNIPs) through active learning. In particular, we
compare incumbent ensemble-based methods against strategies that use single,
deterministic NNs: mean-variance estimation, deep evidential regression, and
Gaussian mixture models. We explore three datasets ranging from in-domain
interpolative learning to more extrapolative out-of-domain generalization
challenges: rMD17, ammonia inversion, and bulk silica glass. Performance is
measured across multiple metrics relating model error to uncertainty. Our
experiments show that none of the methods consistently outperformed each other
across the various metrics. Ensembling remained better at generalization and
for NNIP robustness; MVE only proved effective for in-domain interpolation,
while GMM was better out-of-domain; and evidential regression, despite its
promise, was not the preferable alternative in any of the cases. More broadly,
cost-effective, single deterministic models cannot yet consistently match or
outperform ensembling for uncertainty quantification in NNIPs.

###### keywords:

Uncertainty quantification, neural network interatomic potentials, adversarial
sampling, single deterministic neural networks, neural network ensemble

##  1 Introduction

Over the last decade, neural networks (NN) have increasingly been deployed to
study complex materials systems. In addition to modeling Quantitative
Structure Property/Activity Relationship (QSPR/QSAR), NNs have been used
extensively to model interatomic potentials in atomistic simulations of
materials[1, 2, 3, 4, 5]. When employed to predict the potential energy
surfaces (PESes) of materials systems, NN interatomic potentials (NNIPs)
recover the accuracy of ab initio methods while being orders of magnitude
faster, enabling simulations of larger time- and length-scales[6, 7].
Furthermore, NNIPs have displayed sufficient flexibility to learn across
diverse chemical compositions and structures[8, 9, 10]. As a result, a wide
range of NNIP architectures has been developed[11, 12, 13, 14, 15, 16] and
used for applications in reactive processes[17, 18], protein design[19, 20],
solids[21, 22], solid-liquid interfaces[23], coarse-graining[24, 25], and
more. Nevertheless, NNIPs remain susceptible to making poor predictions in
extrapolative regimes. This often results in unphysical structures or
inaccurate thermodynamic predictions during atomistic simulations that can
compromise the scientific validity of the study[26, 27].

To maximize NNIP robustness and avoid distribution shift, the training data
should represent the same ensemble that the simulation will visit. However,
since high-quality ab initio calculations are computationally expensive,
acquiring new data points through exhaustive exploration of the chemical space
is intractable. Quantifying model uncertainty and active learning are thus key
to training robust NNIPs. A good uncertainty quantification (UQ) method should
estimate both the uncertainty arising from measurement noise (i.e., aleatoric
uncertainty), and the uncertainty in the prediction due to the model error
(i.e., epistemic uncertainty)[28]. Epistemic uncertainty originates from data
scarcity, model architecture limitations, and/or suboptimal model parameter
optimization. Since ab initio-calculated training data is free from aleatoric
error, NNIPs mostly suffer from epistemic uncertainties due to their strong
interpolative but weak extrapolative capabilities. Hence, it is important to
acquire new training data efficiently using a good UQ metric to facilitate
active learning and mitigate epistemic uncertainty[29, 30]. Out of the many UQ
schemes, Bayesian neural networks (BNNs), in which the weights and outputs are
treated as probability distributions, offer an inherent UQ scheme[31, 32].
However, BNNs suffer from scaling problems and are prohibitively expensive for
use as interatomic potentials for bulk materials systems. Another architecture
commonly used for UQ is the NN ensemble[33]. NN ensemble offers a
straightforward approach to estimating uncertainties since NN weight
initialization, training, and even hyperparameter choices are stochastic.
Although NN ensembles are straightforward and can yield good accuracy, the
computational cost associated with training and inference for an ensemble of
several independent NNs can be very high for large systems. Other uncertainty
estimation methods seek to address these bottlenecks by estimating uncertainty
with a single, deterministic NN. Such strategies are particularly interesting
due to the lower computational cost of training and inference for a single NN
model. Furthermore, the predictive uncertainty of single deterministic NNs can
be estimated from a single evaluation of the data, eliminating the need for
stochastic sampling to generate approximations of the underlying uncertainty
functions. UQ strategies using single deterministic NNs can be broadly grouped
into antecedent and succedent methods[34]. In the succedent scheme, the NN has
been trained and predictor uncertainty corresponding to the data set is
estimated from the feature space[35], latent space[36], or gradient matrices
of the trained NNs[37, 38]. One instance of a succedent method is introduced
by Janet et al.[36] for non-NNIP application and employs Euclidean norm
distance in latent features for uncertainty approximation. A more recent
development building upon this work involves fitting Gaussian mixture models
on the latent space of the NN[39]. In antecedent methods, on the other hand,
priors are placed on the input data[40, 41, 42] or conditions are placed on
the NN’s training objective[43]. The mean-variance estimation method, a
regression-prior scheme, is a commonly used antecedent scheme, where a
Gaussian prior distribution is placed on the input data and the NN is trained
to predict the mean and variance (also taken as the uncertainty)[40]. Another
antecedent method that has been introduced recently is the deep evidential
regression, where a high order prior distribution is placed on top of the
input data[41].

To the best of our knowledge, there exists no uncertainty quantification
method that significantly outperforms all others in the setting of NNIPs. In
fact, there seems to be a lack of comparison between UQ methods for this task.
Different from the prediction of physical properties (QSAR/QSPR), NNIPs both
necessitate and rely upon force-matching gradients in addition to fitting of
PESes, making UQ for NNIP applications more challenging and unique[44, 45]. In
this work, we study if UQ methods using single, deterministic NNs can
consistently outperform NN ensembles with respect to their capacity to produce
a strong ranking of uncertainties while having lower computational costs. In
particular, we use mean-variance estimation (MVE), deep evidential regression,
and Gaussian mixture models (GMM). We evaluate their performance in ranking
uncertainties using multiple metrics across three different data sets, and
analyze the improvement of the stability of molecular dynamics simulations via
active learning. In general, we find that ensemble-based methods consistently
perform well in terms of uncertainty rankings outside of the training data
domain and provide the most robust NNIPs. MVE has been shown to perform well
mostly in identifying data points within the training domain that corresponds
to high errors. Deep evidential regression, offers less accurate epistemic
uncertainty prediction, while GMM is more accurate and lightweight than deep
evidential regression and MVE.

##  2 Methods

Given a data set $\mathcal{D}$ containing a collection of tuples ${(x,y)}$, a
neural network takes the input, $x$ and learns a hypothesis function to
predict the corresponding output $\hat{y}$. In the case of an NNIP, $x$
contains the atomic numbers $\mathbf{Z}\in\mathbb{Z}^{n}_{+}$ and nuclear
coordinates $\mathbf{r}\in\mathbb{R}^{n\times 3}$ of a material system with n
atoms, while $\hat{y}$ denotes its potential energy value,
$\hat{E}\in\mathbb{R}$, and/or the energy-conserving forces,
$\hat{\mathbf{F}}_{i}$ acting on the atom i. The energy-conserving forces,
$\hat{\mathbf{F}}_{i}$ are calculated as the negative gradients of the
predicted potential energy with respect to the atomic coordinates, $r_{i}$

$\hat{\mathbf{F}}_{i}=-\frac{\partial\hat{E}}{\partial r_{i}}$ (1)

The Polarizable atom interaction Neural Network (PaiNN), an equivariant
message-passing NN was used to learn the interatomic potentials[12], given its
balance between speed and accuracy. Details of the architecture and training
are shown in Section I of the Supporting Information.

###  2.1 Uncertainty Quantification (UQ)

Four methods of uncertainty quantification (UQ) were compared: ensemble-based
uncertainty[33], mean-variance estimation (MVE)[40], deep evidential
regression[41], and Gaussian mixture models (GMM)[39]. Since forces are a
derivative of energy, they typically exhibit higher variability upon
overfitting or outside the training domain and thus describe epistemic
uncertainty better than the energy variance. Hence, the uncertainty estimates
were evaluated as $\sigma_{F}$ for all models[29].

Figure 1:  Illustration of the uncertainty quantification (UQ) methods used.
$x$ denotes the input to the neural networks (NNs), while $y$ is the predicted
property. In the case of NNIPs, $x$ generally represents the positions and
atomic numbers of the input structure, whereas $y$ is the energy and/or forces
of the corresponding input structure. In addition, red texts indicate the
variables used as the uncertainty estimates. (a) Multiple NNs trained in an
ensemble are made to predict the desired property of the same structure. The
mean and variance of the property can then be calculated, where higher
variation in the property implies higher uncertainty[33, 46]. (b) In the mean-
variance estimation (MVE) method, a Gaussian prior distribution is applied to
the input data and the NN is made to predict the mean and variance describing
the Gaussian distribution. A higher variance parameter indicates higher
uncertainty[40, 46]. (c) In the deep evidential regression method, an
evidential prior distribution is applied to the input data and the NN predicts
the desired property and the parameters to describe both the aleatoric and
epistemic uncertainties[41, 42]. (d) In the Gaussian mixture model (GMM)
method, the input data is assumed to be drawn from multiple Gaussian
distributions. The negative log-likelihood (NLL) function is calculated from a
fitted GMM on the learned feature vectors, $\xi_{x}$ of the structures. A
higher NLL value denotes higher uncertainty[39].

Ensemble-based models are among the most common, cost-effective and trusted
approaches to UQ in deep learning. Because there are various sources of
randomness in the training of NNs, from initializations to stochastic
optimization or hyperparameter choices, a distribution of predicted output
values can be obtained from the ensemble-based model that consists of $M$
distinct NNs[33, 46, 34, 29]. For an input system $x$, each neural network $m$
in the ensemble predicts the total energy, $\hat{E}_{m}(x)$ and atomic forces
$\hat{\mathbf{F}}_{m}(x)$ (calculated using Eq. 1). The predictions of the
ensemble are then taken to be the arithmetic mean of the predictions from each
individual NN

$\displaystyle\mu_{E}(x)$
$\displaystyle=\frac{1}{M}\sum^{M}_{m}\hat{E}_{m}(x),$ (2)

$\displaystyle\bm{\mu}_{F}(x)$
$\displaystyle=\frac{1}{M}\sum^{M}_{m}\hat{\mathbf{F}}_{m}(x)$ (3)

Lakshminarayanan et al. has proposed that variance of predictions from
individual NNs in an ensemble can be used as the uncertainty estimates since
the predicted values of NNs may become higher outside of the training data
domain[33]. In the application of NNIPs, the uncertainty estimates are
provided by variances of predicted energy and forces, which can be computed as

$\displaystyle\sigma_{E}^{2}(x)$
$\displaystyle=\frac{1}{M-1}\sum^{M}_{m}(\hat{E}_{m}(x)-\mu_{E}(x))^{2}$ (4)

$\displaystyle\sigma_{F}^{2}(x)$
$\displaystyle=\frac{1}{M-1}\sum^{M}_{m}\left[\frac{1}{3n}\sum_{i}^{n}(\hat{\mathbf{F}}_{i,m}(x)-\bm{\mu}_{F_{i}}(x))^{2}\right]$
(5)

Mean-variance estimation (MVE). As a function approximator, a standard NN does
not have the intrinsic capability to estimate the trustworthiness of its
prediction. By treating the training data set as Gaussian random variables
with variances, a NN can be constructed to predict the mean, variance pairs
that parametrize the Gaussian distributions $\mathcal{N}(\mu,\sigma^{2})$[40].
In other words, the NN predicts a distribution - itself parameterized by a
mean and a variance - as an output rather than a fixed point. This method is
called mean-variance estimation (MVE), in which the NN predicts the mean and
variance for each datum using maximum likelihood estimation.

The final layer of the NN is modified to predict both the mean, $\mu_{E_{i}}$
(gradients are subsequently computed to obtain $\mu_{F_{i}}$) and variance,
$\sigma_{F}^{2}$ for each datum, and a softplus activation is applied to
ensure that the variances are strictly positive. The predicted variances are
then taken as the uncertainty metric. Similar to the ensemble method, the
uncertainty estimates are described only by the variances of the atomic
forces.

Deep evidential regression. While MVE has been widely used for UQ, the
predicted variance used in lieu of a proper uncertainty fails to provide
estimates of epistemic uncertainties. This is due to the assumption underlying
MVE, where the predicted variance of each data point only describes the
Gaussian distribution over the data but not over the NN parameters. Thus, MVE
does not rigorously inform about distributions outside of the training data
regions[47]. To address this limitation, an emerging evidential regression
algorithm has been developed to directly handle representations of epistemic
uncertainty by parameterizing a probability distribution function over the
likelihood of the learned NN representations[41].

Instead of placing priors directly on the data set like MVE, evidential
regression places priors on the likelihood function of the data set. It is
assumed that the data set is drawn from an independent and identically
distributed (i.i.d.) Gaussian distribution ($\mathcal{N}(\cdot)$) with unknown
mean and variance ($\mu,\sigma^{2}$). Here we assume that the data is
represented by the atomic forces, $\mathbf{F}_{i}$. With this assumption, the
unknown mean can then be described as a Gaussian distribution while the
variance is drawn from an inverse-gamma distribution ($\Gamma(\cdot)$)[41].

$\begin{split}(\mathbf{F}_{1},\mathbf{F}_{2},...,\mathbf{F}_{N})&\sim\mathcal{N}(\mu_{F},\sigma^{2}_{F})\\\
\mu_{F}\sim\mathcal{N}(\gamma_{F},\sigma^{2}\nu^{-1}),\ \ &\
\sigma^{2}_{F}\sim\Gamma^{-1}(\alpha,\beta)\end{split}$ (6)

where $\gamma$, $\nu$, $\alpha$ and $\beta$ are predicted from the NN, and
$\gamma\in\mathbb{R}$, $\nu>0$, $\alpha>1$ and $\beta>0$. Consequently, this
higher-order evidential distribution is represented over the mean and variance
as

$p(\mu_{F},\sigma_{F}^{2}|\gamma_{F},\nu,\alpha,\beta)=\frac{\beta^{\alpha}\sqrt{\nu}}{\Gamma(\alpha)\sqrt{2\pi\sigma_{F}^{2}}}\left(\frac{1}{\sigma_{F}^{2}}\right)^{\alpha+1}\exp\left(-\frac{2\beta+\nu(\gamma_{F}-\mu_{F})^{2}}{2\sigma_{F}^{2}}\right)$
(7)

For detailed derivation, we refer to Ref. [47, 41]. In a similar fashion to
MVE, the final layer of the NN is modified to predict $\gamma_{E}$
($\gamma_{F}$ is then calculated from derivative shown in Eq. 1), $\nu$,
$\alpha$ and $\beta$. A softplus activation is applied to the 3 parameters
$\nu$, $\alpha$ and $\beta$ to ensure that they are strictly positive (with a
+1 for $\alpha$). Finally, the aleatoric, $\mathbb{E}[\sigma_{F}^{2}]$ and
epistemic uncertainty, $\text{Var}[\mu_{F}]$ can be calculated, respectively,
as

$\displaystyle\mathbb{E}[\sigma_{F}^{2}]$
$\displaystyle=\frac{\beta}{\alpha-1}$ (8) $\displaystyle\text{Var}[\mu_{F}]$
$\displaystyle=\frac{\beta}{\nu(\alpha-1)}.$ (9)

Gaussian mixture model (GMM). A GMM models the distribution of the data set
from a weighted mixture of $K$ multivariate Gaussians of dimension $D$[48].
Each $D$-variate Gaussian $k$ in the mixture has a mean vector,
$\mathbf{\mu}_{k}\in\mathbb{R}^{D\times 1}$ and a covariance matrix,
$\Sigma_{k}\in\mathbb{R}^{D\times D}$. The data distribution can then be
described as

$p(x|{\pi_{k},\mu_{k},\Sigma_{k}})=\sum_{k=1}^{K}\pi_{k}\mathcal{N}(\mu_{k},\Sigma_{k}),$
(10)

where $\pi_{k}$ are the mixture weights and satisfy the constraint
$\sum_{k=1}^{K}\pi_{k}=1$. To start, the latent features of each training set
datum $\xi_{\text{train}}\in\mathbb{R}^{D\times 1}$ are extracted from the
immediate layer of a trained NN before the per-atom energy prediction [39]. A
GMM is constructed on $\xi_{\text{train}}$ using the expectation-maximization
(EM) algorithm with full-rank covariance matrices and mean vectors
predetermined using k-means clustering. To determine the number of Gaussians
needed, we investigated the trade-off between the Bayesian Information
Criterion (BIC) and the silhouette score (see Figures S4 and S5). Using the
fitted GMM, we can obtain uncertainty of the test data by taking the latent
features, $\xi_{\text{test}}\in\mathbb{R}^{D\times 1}$ and evaluating the
negative log-likelihood function
$\text{NLL}(\xi_{\text{test}}|\xi_{\text{train}})$ for the test data. For all
our models, the dimension of latent features, $D$ were set to 128. A higher
$\text{NLL}(\xi_{\text{test}}|\xi_{\text{train}})$ indicates that the data is
“far” from the mean vectors of the Gaussians in the GMM and thus corresponds
to high uncertainty.

$\text{NLL}(\xi_{\text{test}}|\xi_{\text{train}})=-\log\left(\sum_{k=1}^{K}\pi_{k}\mathcal{N}(\xi_{\text{test}}|\mu_{k},\Sigma_{k})\right)$
(11)

###  2.2 Data Sets

Revised MD17 (rMD17). The rMD17 data set contains snapshots from long
molecular dynamics trajectories of ten small organic molecules[49, 50], in
which the energy and forces were calculated at the PBE/def2-SVP level of
electronic structure theory. For each molecule, there are 5 splits and each
split consists of 1000 training and 1000 testing data points. The data set was
obtained from https://doi.org/10.6084/m9.figshare.12672038.v3. Because the
test set is drawn randomly from the same trajectory as the training data,
rMD17 is more representative of an in-domain interpolative UQ challenge.

Ammonia. The ammonia training data set consists of 78 geometries, where the
energies and forces are calculated using the BP86-D3/def2-SVP level of theory
as implemented in ORCA. The geometries were generated using hessian-
displacement in the direction of normal mode vectors on initial molecular
conformers generated using RDKit with the MMFF94 force field[29]. The data set
was obtained from https://doi.org/10.24435/materialscloud:2w-6h. The test set
contains 200 geometries with energies ranging from 0 to 100 kcal/mol (2
geometries in each 1 kcal/mol bin), calculated using the same level theory.
The geometries were generated using adversarial sampling during trial runs of
this work, and are completely independent of the models and data set analyzed
in this work. As the test set comprises higher-energy geometries and the
training set encompasses low-energy structures near the ground state, the
ammonia dataset serves as a fundamental example of an out-of-domain,
extrapolative UQ challenge.

Silica glass. The silica glass (SiO2) training data set consists of 1590
structures, sampled from multiple long molecular dynamics trajectories
subjected to different conditions including temperature equilibration,
cooling, uniaxial tension, uniaxial compression, biaxial tension, biaxial
compression, and biaxial shear at different rates using the force-matching
potential described in Ref [51]. Each silica glass structure consists of 699
atoms (233 Si and 466 O atoms). In total, 1590 silica structures were sampled
from the molecular dynamic trajectories to form the train set and 402
structures were selected using adversarial sampling (obtained during
preliminary tests for this study, these samples were acquired separately from
the results presented here.[52]) to form the test set. DFT calculations were
performed on the structures using the Vienna Ab-initio Simulation Package
(VASP)[53, 54, 55, 56]. Details of the MD simulations[57, 58], DFT
calculations, and adversarial sampling of the silica structures are discussed
in Section II.A of the Supporting Information. Because of its high
configurational complexity - including high energy fracture geometries - and
low chemical complexity, our silica dataset represents a step-up in
generalization of UQ in out-of-domain extrapolative regimes.

###  2.3 Evaluation Metrics

Spearman’s rank correlation coefficient. Since UQ should measure how much the
model does not know about the data presented to it, we expect that a high true
error should correlate with a high predicted uncertainty. In other words, we
expect a monotonic relationship between the predicted uncertainty, $U$, and
the true error, $\epsilon$, for a good uncertainty estimator[46, 59]. Since we
are particularly interested in UQ in the context of active learning to improve
NNIPs where they are most erroneous, we use the Spearman’s rank correlation
coefficient to assess the reliability of the uncertainty estimator such that
the correlation between the rank variable of the predicted uncertainty,
$R_{U}$, and that of the true error, $R_{\epsilon}$, can be quantified. The
correlation coefficient can then be defined as

$\rho(U,\epsilon)=\frac{\text{cov}(R_{U},R_{\epsilon})}{\sigma_{R_{U}}\sigma_{R_{\epsilon}}}$
(12)

where $\text{cov}(R_{U},R_{\epsilon})$ is the covariance of the rank variables
while $\sigma_{R_{U}}$ and $\sigma_{R_{\epsilon}}$ are the standard deviations
of the rank variables. In the case of a perfect monotone function, a Spearman
correlation of +1 or -1 is expected, whereas a correlation of 0 indicates no
association between the ranks of the two variables. However, since this metric
only compares ranks of variables, it is still possible for an estimator
constantly predicting low uncertainties in the event of high errors to achieve
high Spearman coefficient. Hence, complementing this with another metric
allows us to understand the performance better.

Area under the receiver operating characteristic curve (ROC-AUC). Since
uncertainty estimates are expected to be high at high true error points, a
criterion to split uncertainties and errors into high and low values is
established[60, 39]. More precisely, an error threshold, $\epsilon_{c}$ is
specified such that data with true error higher than the error threshold,
$\epsilon>\epsilon_{c}$ are classified as “high error” points while data with
lower true error, $\epsilon\leq\epsilon_{c}$ are classified as “low error”
points. Similarly, a specified uncertainty threshold is applied to uncertainty
$U_{c}$. For a structure to be considered as a true positive (TP), both the
true error and estimated uncertainty have to be above their corresponding
thresholds ($(\epsilon>\epsilon_{c})\wedge(U>U_{c})$). A false positive (FP)
point is a structure with a true error below its threshold while the estimated
uncertainty is above the threshold
($(\epsilon\leq\epsilon_{c})\wedge(U>U_{c})$). Conversely, a true negative
(TN) point occurs when $(\epsilon\leq\epsilon_{c})\wedge(U\leq U_{c})$, and a
false negative (FN) when $(\epsilon>\epsilon_{c})\wedge(U\leq U_{c})$.

However, setting the threshold values requires manual optimization for each
uncertainty method, especially since the magnitude and distribution of
uncertainties for each method compared in this work vary extensively. Instead
of evaluating the true positive rate (TPR) and predictive positive value (PPV)
like in previous work[60, 39], we compute the area under curve (AUC) of the
receiver operating characteristic (ROC) curve such that the quality of the
uncertainty methods can be assessed independent of the discrimination
threshold[43]. Note that the threshold for “target” score, which in this case
is the threshold for true error, $\epsilon_{c}$, still has to be specified in
order for the ROC curve to be plotted. For all models, we set $\epsilon_{c}$
to be at the 20th percentile of the error distribution because we find that
varying the error percentile threshold does not affect the ROC-AUC score by a
significant amount (see Table S2). We also found that using system-dependent
absolute cutoff values give the same trend as using percentile cutoffs. The
more accurate the method is, the higher the ROC-AUC score. An ROC-AUC score of
0.5 indicates that the method has no discrimination power between high or low
uncertainty samples. A score lower than 0.5 implies that the method is better
at inverting the uncertainty prediction, predicting low uncertainty for high
error points.

Miscalibration area. Another way to evaluate the quality of the uncertainty
estimates is to quantify the calibration error of the models, which is a
numerical measure of the accuracy of calibration curve. A calibration curve
shows the relationship between the predicted and observed frequencies of data
points within different ranges. Tran et al. [61] proposed to compare how well
the expected fraction of errors for each data point falling within $z$
standard deviations of the mean follow the observed Gaussian random variables
constructed with the uncertainty estimates ($U(x)$) as the variances. In other
words, if the quality of uncertainty estimate is good, we would expect that
there will be around 68% of the errors falling within one standard deviation
of the mean of the Gaussian distributions constructed using the uncertainty
estimates[61, 46, 62, 59]. From this calibration curve, we can calculate the
area between this curve and the line of perfect calibration to provide a
single numerical summary of the calibration. This quantity is called the
miscalibration area. A perfect UQ method will show an exact agreement between
the observed and expected fractions and thus yield a miscalibration area of 0.

Calibrated negative log-likelihood function (cNLL). For observed (true)
errors, a negative log likelihood can be calculated as an evaluation metric
assuming that the errors are normally distributed with a mean of 0 and
variances given by the predicted uncertainties. However, given that some
methods do not directly produce uncertainty estimates equivalent to variances,
the uncertainty estimates can be calibrated such that they resemble variances
more closely[46, 36]. Specifically, the “variances” can be approximated as
$\hat{\sigma}^{2}(x)=aU(x)+b$, where the scalars can be $a$ and $b$ can be
estimated as values that minimize the negative log-likelihood function of
errors in the validation set, $\mathcal{D}_{val}$,

$a_{*},b_{*}=\operatorname*{argmin}_{a,b}\
\frac{1}{2}\sum_{x,\mathbf{F}\in\mathcal{D}_{val}}\ln{(2\pi)}+\ln{(aU(x)+b)}+\frac{(\hat{\mathbf{F}}-\mathbf{F})^{2}}{(aU(x)+b)}\
.$ (13)

With the calibrated scalars $a$ and $b$, the cNLL metric can then be computed
on the test set, $\mathcal{D}_{test}$, as

$\text{cNLL}=\frac{1}{2}\sum_{x,\mathbf{F}\in\mathcal{D}_{test}}\ln{(2\pi)}+\ln{(a_{*}U(x)+b_{*})}+\frac{(\hat{\mathbf{F}}-\mathbf{F})^{2}}{(a_{*}U(x)+b_{*})}\
.$ (14)

###  2.4 Adversarial Sampling

In addition to evaluating the reliability of an uncertainty estimator using
metrics like Spearman’s rank correlation coefficient and ROC-AUC, we would
like to assess whether new structures sampled using the estimated
uncertainties from each method help to increase robustness of the NNIPs in
different ways. To this end, we employed the adversarial sampling method
proposed by Schwalbe-Koda et al.[29] to sample new conformations that maximize
the estimated uncertainty such that these adversarial examples can be
incorporated into the retraining of the NNs in an active learning (AL) loop.
This would test reliability, smoothness, and utility of the UQ methods because
adversarial sampling requires the UQ method to systematically identify
underrepresented regions where the models are most uncertain about.
Effectively, improvement in robustness of NNIPs after active learning is a
better metric for uncertainty estimation methods since ideally a good
uncertainty estimator is able to acquire new training examples that improve
the model’s generalizability more than random sampling.

In the adversarial sampling strategy, a perturbation, $\delta$ is applied to
randomly chosen structures, $x$ from the initial training data. The
perturbation is then iteratively updated using gradient-ascent to maximize the
uncertainty of the perturbed structures, $U(x_{\delta})$, as defined below

$\displaystyle\max_{\delta}p(x_{\delta})U(x_{\delta}).$ (15)

More precisely, $U(x_{\delta})$ is taken to be $\sigma^{2}_{F}(x_{\delta})$
for the ensemble-based method (Eqs. 5) and MVE,
$\text{Var}[\mu_{F}(x_{\delta})]$ for evidential regression (Eq. 9), and
$\text{NLL}(x_{\delta}|x_{\text{train}})$ for GMM (Eq. 11). $p(x_{\delta})$ is
the thermodynamic probability of the perturbed structures based on the
distribution of data already available to the model

$\displaystyle
p(x_{\delta})=\frac{\exp\left(-\frac{\hat{E}(x_{\delta})}{k_{\mathrm{B}}T}\right)}{\sum_{x\in\text{train}}\exp\left(-\frac{E(x)}{k_{\mathrm{B}}T}\right)}$
(16)

where $k_{\mathrm{B}}$ is the Boltzmann constant, $T$ the absolute
temperature, and $\hat{E}(x_{\delta})$ in the numerator is the predicted
energy of the perturbed structure, while $E(x)$ in the denominator denotes
ground truth energy of structures in the initial training data set. Refer to
Section III in the Supporting Information for parameters used for adversarial
sampling.

###  2.5 Molecular Dynamics (MD) Simulations.

NN-based MD simulations were performed in the NVT ensemble with the Nosé-
Hoover thermostat. In the case of ammonia, 100 5 ps-long MD simulations were
run for four different temperatures (300 K, 500 K, 750 K, and 1000 K) and for
each UQ method at a timestep of 0.5 fs. The initial configurations for all
simulations were picked randomly from the initial training data set containing
78 geometries. The trajectories were considered unstable if the distance
between the atoms were closer than 0.75 Å or larger than 2.25 Å or if the
predicted energy was lower than the ground state energy, which is at a
reference point of 0 kcal/mol in this work. For simulations of silica, 10 1
ns-long trajectories at a temperature of 2500 K were performed for each UQ
methods at all generations with a time step of 0.25 fs. The initial structures
for all simulations were picked randomly from the structures sampled at 2500 K
equilibration simulation in the initial training data set. The trajectories
were considered unstable if the kinetic energy becomes 0, a null value or
greater than 10,000 kcal/mol (due to destabilization of the predicted
potential energy in a false minimum[27]), or if the distances between any
atoms are less than 1.0 Å. No upper bound was imposed on the atomic distances
because the simulations tend to fail before any stretching of bond occurs. All
simulations were performed using the Atomic Simulation Environment (ASE)
library in the Python language[63].

##  3 Results and Discussion

###  3.1 Revised MD17 Data Set (rMD17)

Figure 2:  (a) Hexbin plots showing (predicted) uncertainties versus squared
errors of atomic forces in all molecules of all 5-fold test sets from the
rMD17 data set for each considered UQ method. Note that since the
uncertainties (-NLL) for the GMM method contain negative values, all
uncertainties are scaled by the minimum value to yield positive quantities for
plotting of the log scale figure. (b) Total wall time used for training the
NNs and total GPU memory usage incurred by the NNs. (c) Bar chart showing the
prediction errors on the test set of the considered UQ models trained using
their corresponding loss functions. Each bar presents the overall average of
the metrics for all molecules in all 5-fold test sets of the rMD17 data set.
Prediction error is described by the mean absolute error (MAE) of the atomic
forces in the unit kcal/mol/Å. (d) Bar charts showing the evaluation metrics
for the considered UQ methods, including Spearman’s rank correlation
coefficient, ROC AUC score, miscalibration area and calibrated negative log-
likelihood (NLL) score. For (c) and (d), refer to Figure S2 for breakdown of
the metrics for each molecule.

As a starting point, we performed UQ on the rMD17 data set, which is readily
available online. We evaluated whether the predicted uncertainties are ranked
correctly with respect to the true error for all molecules in the data set
(Figure 2a). An immediate observation is that the uncertainty-error
distributions for the ensemble-based and mean-variance estimation (MVE)
methods follow a positively-correlated relationship desirable for good UQ
schemes. Similarly, the predicted uncertainty of Gaussian mixture model (GMM)
exhibits a relatively linear but slightly skewed distribution with respect to
the true errors. Evidential regression, on the other hand, shows a two-tiered
distribution that does not show good correlation with energy. In addition to
that, it is interesting to note that the predicted uncertainties by evidential
regression span eight orders of magnitudes, which is due to optimization of
the $\nu$, $\alpha$, and $\beta$ parameters during training such that they
approach 0. In terms of computational resources used for each method (Figure
2b), the training time and GPU memory usage for the ensemble-based method
containing 5 independent neural networks (NNs) are roughly 5 times greater
than the other single deterministic NN methods.

Since UQ is often utilized in the context of NNIPs to improve robustness for
atomistic simulations, it is critical to first assess the prediction accuracy
of each model. From Figure 2c, we can see that the ensemble-based method has a
lower overall error at inference than the other methods which all employ only
a single NN. Given that the choice of model architecture, input
representation, and training data for all UQ methods in the same splits are
the same, we can assume that all models possess the same model uncertainty
(also referred to as model bias), which is the component of epistemic
uncertainty associated with model form insufficiencies and lack of data
coverage. This suggests that the accuracies of models differ due to suboptimal
parameter optimization mainly arising from the difference in loss functions
used during training. It has been shown that ensemble learning removes
parametric uncertainty (also referred to as model variance) associated with
ambiguities in parameter optimization[34, 28, 64], and thus improves
generalization which reduces the overall error[29, 65, 66]. MVE, on the other
hand, shows the highest average test set error, possibly due to a harder-to-
optimize NLL loss function, which has been reported in Ref [67]. Since the
robustness of NNIPs largely depends on the accuracy of the models, the overall
error is an important criterion to consider while choosing UQ methods to be
incorporated into the active learning workflows. In terms of quality of
uncertainty estimates (Figure 2d), no single method consistently outperforms
the other. For the ensemble-based method, all the evaluated metrics are, by
comparison, in good ranges. Interestingly, MVE demonstrates good performance
in uncertainty ranking despite having low prediction accuracy. Evidential
regression, in contrast, outperforms the others in all evaluated metrics
except for the miscalibration area, possibly due to the bimodal distribution
in Figure 2a. Lastly, GMM shows the worst performance in all metrics although
within the error bar of the other approaches. We also found that the
performance of GMM is independent of the dimension of the latent feature
vector and the number of Gaussians (see Figure S4).

###  3.2 Ammonia

Figure 3:  (a) Plots showing predicted uncertainties versus squared errors of
atomic forces in the test data set for each considered UQ method trained only
on the initial training data set (subfigure (c) below). The box at the top
right corner of each subfigure shows the evaluation metrics of the methods,
namely Spearman’s rank correlation coefficient ($\rho$), ROC AUC score (AUC),
miscalibration area (Area), and calibrated negative log-likelihoods (cNLL).
(b) Energy barrier of nitrogen inversion calculated with NEB using DFT and the
considered UQ methods in generation 1. All NNs in the methods were trained on
the same initial training data. (c) Histogram showing distribution of energy
of geometries in the initial training data. (d) Fraction of stable MD
trajectories generated using the NNs of the UQ methods in generation 1 as
force field. (e) Energy barrier of nitrogen inversion calculated with NEB
using DFT and the UQ methods in generation 3. The NNs are trained on new
adversarial examples generated with their respective UQ method on top of the
initial training data. (f) Distribution of energy in training set after two
rounds of adversarial sampling is performed. (g) Fraction of stable MD
trajectories generated using the NNs of the UQ methods in generation 3 as
force field.

In this example, we tested the UQ methods by studying the process of nitrogen
inversion in ammonia. In the first generation, all NNs are trained on the same
initial training data that contains only 78 geometries with energies ranging
between 0 and 5 kcal/mol (Figure 3c). When tested on 200 conformations with
energies between 0 and 100 kcal/mol (2 conformations per 1 kcal/mol bin), it
can be clearly observed that the ensemble-based and GMM model produce good
monotonous relationship for the predicted uncertainty with respect to the true
error (Figure 3a). This can also be confirmed by examining the evaluation
metrics for ensemble and GMM in which all metrics outperform the other two
methods. Intriguingly, the performance of GMM on this data set is quite
different from the results for the rMD17 data set, where MVE and evidential
regression discernibly outperform GMM (Figure 2d). The most probable reason
for this is the difference in distribution of train and test sets. More
specifically, the energies of conformations in the train and test sets overlap
in the rMD17 data set (Figure S7), whereas for ammonia, the train set contains
only low energy conformations (0 - 5 kcal/mol) while the test set includes
conformations in both low and high energy ranges (0 - 100 kcal/mol). This
suggests that the ensemble-based and GMM methods are better at predicting
epistemic uncertainties outside of the training data domain. Evidential
regression is also good at ranking uncertainties with respect to the true
errors, but the distribution again shows a peculiar triangular trend. By
contrast, MVE underperforms compared to the other models in all evaluated
metrics and produces a two-tiered distributions. Comparing the performance of
MVE in this data set against the rMD17 data set, we can infer that MVE is
likely worse at predicting uncertainties for data outside of the training
domain.

When it comes to predicting the nitrogen inversion energy barrier, all the
methods seem to either under- or over-estimate the overall ground truth energy
curve (Figure 3b). This observation is reasonable since all models have only
seen low energy conformations in the initial (first generation) training data
(Figure 3c). When the different models are used as interatomic potentials for
MD simulations, NN-ensembles produced the highest fraction of stable MD
trajectories at 50%, in which there are no atomic dissociations or collisions
(Figure 3d). Evidential regression produces 42% of stable trajectories,
followed by GMM at 33% and lastly, MVE at 25%.

After performing adversarial sampling for two generations, 40 data points were
added to training of the NNIPs for each UQ method. It is important to note
that the new data points were generated independently by each UQ method and
differ from one method to the other. The energy distributions of the new
training data for each method, including the initial training data that
contains 78 geometries, are shown in Figure 3f. Looking closely at the energy
distributions, we can see that while the ensemble-based method are able to
drive some conformations to higher energies, the majority of the new data
points are still concentrated in the lower energy region around 0-30 kcal/mol.
This is also similar in the MVE method, but the number of new conformations
with energies in the 20-60 kcal/mol range is higher than in the ensemble-based
method. In contrast, evidential regression samples more new geometries with
energies in the range of 40-85 kcal/mol and below 20 kcal/mol. This may be due
to the diverging trend shown in Figure 3a, where the maximization of
uncertainty during the adversarial sampling pushes geometries to low and high
energy ranges in a tiered manner. Similar observations can be made for the
geometries sampled by the GMM method, but with the new geometries concentrated
more in the 70-90 kcal/mol. When attempting to predict the nitrogen inversion
barrier in ammonia, only the ensemble-based method improved after two rounds
of adversarial sampling and produced accurate predictions (Figure 3e).
However, after two rounds of adversarial sampling, the fraction of stable MD
trajectories produced by the ensemble-based model, MVE, and evidential
regression improved by roughly 30-40% points (Figure 3g). GMM, on the other
hand, achieved an improvement of 52% points, indicating that the quality of
uncertainty estimates improved robustness of the NNIP outside of the training
domain more than the other methods.

###  3.3 Silica Glass

Figure 4:  (a) Distributions of energy and forces in training and testing
sets. In the first generation (first row), all models are trained on the same
initial training set. In subsequent generations, NN models are trained on both
initial training data and new sampled data points using their corresponding UQ
methods. All models are finally evaluated on the test set, where the energy
and force magnitudes are considerably larger than the train sets. Note that
the x-axis scale of the test set energy distribution (last row first column)
is different from the other energy distributions. (b) Hexbin plots showing the
predicted uncertainties against squared errors of atomic forces from NNs
trained only on initial training data. (c) Average length of time in which no
atomic dissociation or collision occurs in the MD trajectories. Average is
taken from 10 trajectories. Item “ensemble (one network)” in legend describes
MD trajectories performed with predictions only from a single NN in the
ensemble. (d) Mean absolute error of predicted energy and forces by the NN
models on the test set.

Finally, we tested the UQ methods on silica glass, which has a long history
for interatomic potential fitting[68, 69, 70, 51, 71]. This is motivated by
the fact that silica glass is made up of large bulk structures and the
computational cost of performing quantum chemistry calculations on these bulk
systems is very high. NNs, which have increasingly been used for interatomic
potential fitting for vitreous silica[72, 73], are unfortunately still
unreliable outside of the training data domain. This calls for efficient
sampling methods to generate small sets of additional data to improve high-
error regions. These rely significantly on the quality of UQ methods, which
should additionally be lightweight and fast. We trained individual NNs and an
NN ensemble containing 4 networks on the same initial training data, and
evaluated the errors and uncertainties from the respective UQ method on the
test data set. Distributions of the initial training and testing data set are
shown in Figure 4a. The testing set contains structures with very high
energies and forces which have not been trained or validated on any of the
NNs, such that evaluation of the predicted uncertainties are guaranteed to be
epistemic uncertainties outside of training data domain. Evaluated true errors
and uncertainties, presented in Figure 4b, for the UQ methods show that all
methods are capable of ranking uncertainties relatively well with respect to
the true errors, even when the evaluated metrics deliver inconsistent results
as to which model is the best. Disputably, evidential regression seems to
perform the worse in uncertainty ranking for this data set considering that
the evaluated metrics, the order of magnitude of predicted uncertainties, and
the distribution of uncertainty-error are less favorable.

We then performed adversarial sampling using each of the UQ methods and
obtained new data points which are evaluated using DFT and incorporated into
re-training of the NNs in subsequent generations. Note that for NNs in each UQ
method, the training data thus becomes different starting in generation 2
since the new data points are sampled using different UQ methods (Figure 4a).
In generation 2, energies of the adversarial samples for the methods have
similar distributions with the exception of evidential regression, which
interestingly did not drive the adversarial samples to very high energies but
to energies lower than 0 kcal/mol (the reference energy). The reference energy
is chosen as the minimum energy of the initial training data, which consists
only of structures extracted from MD at 300 K or higher (refer to Figure S1),
and evidential regression instead sampled structures with energies lower than
those included in the training data. Ensemble-based method, on the other hand,
sampled structures that are within 0-5000 kcal/mol energy ranges, but contain
high atomic forces. While this is arguably the constraint constructed for
adversarial sampling (Equation 15), the effect of this constraint is most
apparent in the ensemble-based method, where these new data points
significantly increased the robustness of the NNIP. In the third generation,
GMM sampled new geometries with both energies and atomic forces much higher
than the other methods.

After training of NNs in each generation, 10 MD simulations are performed
using the NNs. The lengths of time where the simulations are stable are
recorded in Figure S15. Averages are plotted in Figure 4c. The figures show
that the simulations for generation 1 behaved unphysically early on (training
was conducted only on the initial data). This is potentially due to prediction
of false minima from the NN potentials that produce unphysical structures
during the simulations[27]. Additionally, the length of time of stable MD
trajectories differ despite the shared training data and model architecture.
This is consistent with the results on the ammonia data set, where the
robustness of NNIPs is directly influenced by the accuracy of the predictions,
represented as mean absolute error of the energy and atomic forces (Figure 4d
and Figure S14). In generation 2, the length of stable MD time for ensemble-
based method increases greatly alongside the accuracy of prediction,
potentially as a result of the inclusion of new adversarial data with very
high forces. A similar trend can be observed for GMM, but at a smaller
magnitude of improvement. However, in generation 3, there seems to be a
reversal of improvement for both ensemble-based and GMM methods even though
more data has been included and prediction errors on test set shown in Figure
4d have decreased. This suggests that adding more data does not always improve
robustness of NNs and that quality of new data is more crucial[74]. MVE shows
steady improvement but all MD simulations were unable to exceed 1 ps. In
contrast, evidential regression did not show improvement with inclusion of
more structures that were obtained through adversarial sampling (Figure S16).
Since the gradients of the different UQ measures could be different, we set
both the Boltzmann temperatures and learning rates in adversarial sampling to
be the same for all methods, but increased the number of iterations such that
the uncertainties could be maximized. For this data set, we can see that
ensemble-based method not only produces higher stability simulations, it also
samples higher quality uncertain data in a way that cannot be quantified using
the evaluation metrics.

###  3.4 General Remarks

In general, there seems to be a trade-off between the prediction accuracy,
confidence of UQ, and computational costs among the methods. For ensemble-
based methods, the prediction accuracy and UQ performance demonstrate that it
is more reliable than the other methods. However, the time and memory usage
incurred for training the model scales with the number of NNs used in the
ensemble and can be quite infeasible for extremely large data sets. MVE, on
the other hand, performs relatively well in ranking uncertainty with respect
to the true error and is computationally cheaper. The downside to this method
is the lower prediction accuracy achieved as compared to other methods even
with the same training parameters and conditions. Additionally, the predicted
uncertainty optimized using the loss function during training (Eq. S3)
intrinsically describes only the aleatoric uncertainty. For training data
generated via ab initio methods, the data set is relatively noise-less and
technically does not have meaningful aleatoric uncertainty, meaning MVE has
somehow co-opted epistemic uncertainty. In the case of evidential regression,
the method performs remarkably well in terms of uncertainty estimation but
achieves a slightly lower prediction accuracy than the ensemble-based method.
Nevertheless, it is important to note that training of the evidential
regression method often incurs numerical instability due to optimization of
the loss function (Eq. S6) that drives the parameters ($\nu$, $\alpha-1$, and
$\beta$) close to zero. As for GMM, the model is easy to implement,
computationally cheap, and can be used with any NN without compromising the
achieved prediction accuracy. The downsides to this method, however, are the
need for manual configuration of the number of Gaussians for each data set
(Figures S3 and S4) and the reliance on the accuracy of Gaussian mixture model
fitting.

##  4 Conclusions

The performance of the uncertainty quantification (UQ) methods vary broadly
between data sets based and evaluation metrics. The ensemble-based UQ method
performs consistently well across all data sets and metrics in this work. In
addition, due to the averaging effect, ensemble predictions typically show
lower prediction error albeit at a higher cost. However, other work has shown
that ensemble-based method can only be used to address epistemic uncertainty
stemming from parametric uncertainty (model variance), and is less reliable in
identifying data points outside of the training domain[28]. Contrary to these
findings, our study reveals that the ensemble-based method is capable of
identifying data points outside of the training domain when performing
adversarial sampling.

Mean-variance estimation (MVE), on the other hand, has shown inconsistent
performance across the data sets. We infer that MVE is less capable in
systematically predicting high uncertainty for underrepresented regions, but
is better at identifying regions within the training domain that correspond to
high errors. Deep evidential regression UQ has been shown to be decent but
does not outperform the ensemble-based method and tends to be numerically
unstable during parameter optimization. The Gaussian mixture model (GMM)
approach shows good performance in UQ, especially for out-of-domain frames,
and is more versatile since it is applicable, post-training, to any model
architecture, although not as well in the interpolative regime where MVE
seemed to excel.

Nevertheless, for applications in neural network interatomic potentials
(NNIPs), generalizability plays a very important role to ensure robustness of
the atomistic simulations. A model that seemingly estimates uncertainty well
and achieves good prediction accuracy but fails to achieve stable production
simulations from scratch or through active learning, is less helpful. The
unsolved challenge, thus is identifying a model that is able to rank
uncertainty well, while at the same time scales practically, has high
prediction accuracy, and produces stable atomistic simulations. UQ and
robustness in NNIPs may pose a unique challenge because forces are a
derivative function, which is easily accessible from the energy at training
time, but the relationship between force variance and energy variance is not
handled well by any of the methods explored.

Albeit cheaper than ensemble-based method, single deterministic NNs showed
lower generalizability, and thus failed to exceed the performance of the
ensemble-based method, especially in the active learning loops, since ensemble
learning lowers the parametric uncertainty. This is especially true for MVE
since it has been shown that joint optimization of the mean and variance
estimators lowers the stability of prediction[67]. When a single NN
architecture shows higher generalizability than that of a NN ensemble, then
usage of UQ methods such as deep evidential regression or GMM for the single
NNs would be justified.

More broadly, the lack of a one-size-fits-all solution, the high cost of the
more robust ensemble-based method, and the disappointing performance of
elsewhere-promising novel evidential approaches confirm that UQ in NNIPs is an
ongoing challenge for method development in AI for science.

##  5 Acknowledgements

This work was supported by funding from the AGC, Inc.

##  6 Declarations

  * •

Competing interest The authors declare no competing interests.

  * •

Data availability The molecular dynamic trajectories generated in the current
study have been deposited in the Materials Cloud Archive under accession code
https://doi.org/10.24435/materialscloud:55-sd[52].

  * •

Code Availability The code used to produce the results in this paper is
available at https://github.com/learningmatter-mit/UQ_singleNN.git under the
MIT license.

  * •

Author contributions R.G.-B and A.R.T. conceived the project. A.R.T. designed
and conducted the experiments, performed data analysis and wrote the
manuscript. S.U. performed molecular dynamics simulations and DFT calculations
on the initial silica training data set. S.G. and J.C.B.D. provided fruitful
advice and discussions. J.C.B.D. contributed to content organization and
manuscript review. R.G.-B. supervised the research, secured the funding and
contributed to manuscript writing.

## References

  *   * Butler et al. [2018] Butler, K.T., Davies, D.W., Cartwright, H., Isayev, O., Walsh, A.: Machine learning for molecular and materials science. Nature 559(7715), 547–555 (2018) https://doi.org/10.1038/s41586-018-0337-2
  * Schwalbe-Koda and Gómez-Bombarelli [2020] Schwalbe-Koda, D., Gómez-Bombarelli, R.: Generative Models for Automatic Chemical Design. In: Lecture Notes in Physics vol. 968, pp. 445–467. Springer, Cham (2020). https://doi.org/10.1007/978-3-030-40245-7_21 . http://link.springer.com/10.1007/978-3-030-40245-7_21
  * Kocer et al. [2022] Kocer, E., Ko, T.W., Behler, J.: Neural Network Potentials: A Concise Overview of Methods. Annual Review of Physical Chemistry 73(1), 163–186 (2022) https://doi.org/10.1146/annurev-physchem-082720-034254
  * Behler [2021] Behler, J.: Four Generations of High-Dimensional Neural Network Potentials. Chemical Reviews 121(16), 10037–10072 (2021) https://doi.org/10.1021/acs.chemrev.0c00868
  * Yang et al. [2019] Yang, K., Swanson, K., Jin, W., Coley, C.W., Eiden, P., Gao, H., Guzman-Perez, A., Hopper, T., Kelley, B., Mathea, M., Palmer, A., Settels, V., Jaakkola, T.S., Jensen, K.F., Barzilay, R.: Analyzing Learned Molecular Representations for Property Prediction. Journal of Chemical Information and Modeling 59(8), 3370–3388 (2019) https://doi.org/10.1021/acs.jcim.9b00237
  * Mueller et al. [2020] Mueller, T., Hernandez, A., Wang, C.: Machine learning for interatomic potential models. The Journal of Chemical Physics 152(5), 050902 (2020) https://doi.org/10.1063/1.5126336
  * Botu et al. [2017] Botu, V., Batra, R., Chapman, J., Ramprasad, R.: Machine learning force fields: Construction, validation, and outlook. Journal of Physical Chemistry C 121(1), 511–522 (2017) https://doi.org/10.1021/acs.jpcc.6b10908
  * Behler and Parrinello [2007] Behler, J., Parrinello, M.: Generalized Neural-Network Representation of High-Dimensional Potential-Energy Surfaces. Physical Review Letters 98(14), 146401 (2007) https://doi.org/10.1103/PhysRevLett.98.146401
  * Onat et al. [2018] Onat, B., Cubuk, E.D., Malone, B.D., Kaxiras, E.: Implanted neural network potentials: Application to Li-Si alloys. Physical Review B 97(9), 1–9 (2018) https://doi.org/10.1103/PhysRevB.97.094106
  * Gokcan and Isayev [2022] Gokcan, H., Isayev, O.: Learning molecular potentials with neural networks. WIREs Computational Molecular Science 12(2), 1–22 (2022) https://doi.org/10.1002/wcms.1564
  * Schütt et al. [2018] Schütt, K.T., Sauceda, H.E., Kindermans, P.-J.J., Tkatchenko, A., Müller, K.-R.R.: SchNet - A deep learning architecture for molecules and materials. Journal of Chemical Physics 148(24), 241722 (2018) https://doi.org/10.1063/1.5019779
  * Schütt et al. [2021] Schütt, K.T., Unke, O.T., Gastegger, M.: Equivariant message passing for the prediction of tensorial properties and molecular spectra. arXiv preprint arXiv:2102.03150 (2021)
  * Batzner et al. [2022] Batzner, S., Musaelian, A., Sun, L., Geiger, M., Mailoa, J.P., Kornbluth, M., Molinari, N., Smidt, T.E., Kozinsky, B.: E(3)-equivariant graph neural networks for data-efficient and accurate interatomic potentials. Nature Communications 13(1), 2453 (2022) https://doi.org/10.1038/s41467-022-29939-5
  * Zubatyuk et al. [2021] Zubatyuk, R., Smith, J.S., Nebgen, B.T., Tretiak, S., Isayev, O.: Teaching a neural network to attach and detach electrons from molecules. Nature Communications 12(1), 1–11 (2021) https://doi.org/10.1038/s41467-021-24904-0
  * Unke and Meuwly [2019] Unke, O.T., Meuwly, M.: PhysNet: A Neural Network for Predicting Energies, Forces, Dipole Moments, and Partial Charges. Journal of Chemical Theory and Computation 15(6), 3678–3693 (2019) https://doi.org/10.1021/acs.jctc.9b00181
  * Klicpera et al. [2020] Klicpera, J., Groß, J., Günnemann, S.: Directional message passing for molecular graphs. In: ICLR, pp. 1–13 (2020). http://arxiv.org/abs/2003.03123
  * Gastegger and Marquetand [2015] Gastegger, M., Marquetand, P.: High-Dimensional Neural Network Potentials for Organic Reactions and an Improved Training Algorithm. Journal of Chemical Theory and Computation 11(5), 2187–2198 (2015) https://doi.org/10.1021/acs.jctc.5b00211
  * Ang et al. [2021] Ang, S.J., Wang, W., Schwalbe-Koda, D., Axelrod, S., Gómez-Bombarelli, R.: Active learning accelerates ab initio molecular dynamics on reactive energy surfaces. Chem 7(3), 1–32 (2021) https://doi.org/10.1016/j.chempr.2020.12.009
  * Wang et al. [2020a] Wang, Z., Han, Y., Li, J., He, X.: Combining the Fragmentation Approach and Neural Network Potential Energy Surfaces of Fragments for Accurate Calculation of Protein Energy. J. Phys. Chem 2020, 3035 (2020) https://doi.org/10.1021/acs.jpcb.0c01370
  * Wang et al. [2020b] Wang, W., Axelrod, S., Gómez-Bombarelli, R.: Differentiable Molecular Simulations for Control and Learning. arXiv preprint arXiv: 2003.00868 (February) (2020)
  * Marchand et al. [2020] Marchand, D., Jain, A., Glensk, A., Curtin, W.A.: Machine learning for metallurgy I. A neural-network potential for Al-Cu. Physical Review Materials 4(10), 1–21 (2020) https://doi.org/10.1103/PhysRevMaterials.4.103601
  * Jakse et al. [2022] Jakse, N., Sandberg, J., Granz, L.F., Saliou, A., Jarry, P., Devijver, E., Voigtmann, T., Horbach, J., Meyer, A.: Machine learning interatomic potentials for aluminium: application to solidiﬁcation phenomena. Journal of Physics: Condensed Matter (2022) https://doi.org/10.1088/1361-648x/ac9d7d
  * Natarajan and Behler [2016] Natarajan, S.K., Behler, J.: Neural network molecular dynamics simulations of solid–liquid interfaces: water at low-index copper surfaces. Phys. Chem. Chem. Phys. 18(41), 28704–28725 (2016) https://doi.org/10.1039/C6CP05711J
  * Morawietz et al. [2016] Morawietz, T., Singraber, A., Dellago, C., Behler, J.: How van der waals interactions determine the unique properties of water. Proceedings of the National Academy of Sciences of the United States of America 113(30), 8368–8373 (2016) https://doi.org/10.1073/pnas.1602375113
  * Ruza et al. [2020] Ruza, J., Wang, W., Schwalbe-Koda, D., Axelrod, S., Harris, W.H., Gomez-Bombarelli, R., Gómez-Bombarelli, R.: Temperature-transferable coarse-graining of ionic liquids with dual graph convolutional neural networks. Journal of Chemical Physics 153(16), 164501 (2020) https://doi.org/10.1063/5.0022431
  * Fu et al. [2022] Fu, X., Wu, Z., Wang, W., Xie, T., Keten, S., Gomez-Bombarelli, R., Jaakkola, T.: Forces are not Enough: Benchmark and Critical Evaluation for Machine Learning Force Fields with Molecular Simulations (2022) https://doi.org/10.48550/arxiv.2210.07237
  * Morrow et al. [2022] Morrow, J.D., Gardner, J.L.A., Deringer, V.L.: How to validate machine-learned interatomic potentials. arXiv preprint arXiv: 2211.12484, 1–14 (2022)
  * Heid et al. [2023] Heid, E., Mcgill, C.J., Vermeire, F.H., Green, W.H.: Characterizing Uncertainty in Machine Learning for Chemistry. ChemRxiv Preprint. ChemrXiv:10.26434 (2023) https://doi.org/10.26434/chemrxiv-2023-00vcg
  * Schwalbe-Koda et al. [2021] Schwalbe-Koda, D., Tan, A.R., Gómez-Bombarelli, R.: Differentiable sampling of molecular geometries with uncertainty-based adversarial attacks. Nature Communications 12(1), 1–12 (2021) https://doi.org/10.1038/s41467-021-25342-8
  * Shuaibi et al. [2021] Shuaibi, M., Sivakumar, S., Chen, R.Q., Ulissi, Z.W.: Enabling robust offline active learning for machine learning potentials using simple physics-based priors. Machine Learning: Science and Technology 2(2), 025007 (2021) https://doi.org/10.1088/2632-2153/abcc44
  * Vandermause et al. [2020] Vandermause, J., Torrisi, S.B., Batzner, S., Xie, Y., Sun, L., Kolpak, A.M., Kozinsky, B.: On-the-fly active learning of interpretable Bayesian force fields for atomistic rare events. npj Computational Materials 6(1), 20 (2020) https://doi.org/10.1038/s41524-020-0283-z
  * Jinnouchi et al. [2019] Jinnouchi, R., Lahnsteiner, J., Karsai, F., Kresse, G., Bokdam, M.: Phase Transitions of Hybrid Perovskites Simulated by Machine-Learning Force Fields Trained on the Fly with Bayesian Inference. Physical Review Letters 122(22), 225701 (2019) https://doi.org/10.1103/PhysRevLett.122.225701
  * Lakshminarayanan et al. [2017] Lakshminarayanan, B., Pritzel, A., Blundell, C.: Simple and Scalable Predictive Uncertainty Estimation using Deep Ensembles. In: Guyon, I., Luxburg, U.V., Bengio, S., Wallach, H., Fergus, R., Vishwanathan, S., Garnett, R. (eds.) Advances in Neural Information Processing Systems, vol. 30, pp. 6402–6413. Curran Associates, Inc. (2017). https://proceedings.neurips.cc/paper/2017/file/9ef2ed4b7fd2c810847ffa5fa85bce38-Paper.pdf
  * Gawlikowski et al. [2021] Gawlikowski, J., Tassi, C.R.N., Ali, M., Lee, J., Humt, M., Feng, J., Kruspe, A., Triebel, R., Jung, P., Roscher, R., Shahzad, M., Yang, W., Bamler, R., Zhu, X.X.: A Survey of Uncertainty in Deep Neural Networks. arXiv preprint arXiv: 2107.03342, 1–41 (2021)
  * Nandy et al. [2018] Nandy, A., Duan, C., Janet, J.P., Gugler, S., Kulik, H.J.: Strategies and Software for Machine Learning Accelerated Discovery in Transition Metal Chemistry. Industrial & Engineering Chemistry Research 57(42), 13973–13986 (2018) https://doi.org/10.1021/acs.iecr.8b04015
  * Janet et al. [2019] Janet, J.P., Duan, C., Yang, T., Nandy, A., Kulik, H.J.: A quantitative uncertainty metric controls error in neural network-driven chemical discovery. Chemical Science 10(34), 7913–7922 (2019) https://doi.org/10.1039/c9sc02298h
  * Oberdiek et al. [2018] Oberdiek, P., Rottmann, M., Gottschalk, H.: Classification uncertainty of deep neural networks based on gradient information. Lecture Notes in Computer Science (including subseries Lecture Notes in Artificial Intelligence and Lecture Notes in Bioinformatics) 11081 LNAI, 113–125 (2018) https://doi.org/10.1007/978-3-319-99978-4_9
  * Lee and Alregib [2020] Lee, J., Alregib, G.: Gradients as a Measure of Uncertainty in Neural Networks. Proceedings - International Conference on Image Processing, ICIP 2020-Octob, 2416–2420 (2020) https://doi.org/10.1109/ICIP40778.2020.9190679
  * Zhu et al. [2022] Zhu, A., Batzner, S., Musaelian, A., Kozinsky, B.: Fast Uncertainty Estimates in Deep Learning Interatomic Potentials. arXiv preprint arXiv: 2211.09866, 1–38 (2022)
  * Nix and Weigend [1994] Nix, D.A., Weigend, A.S.: Estimating the mean and variance of the target probability distribution. IEEE International Conference on Neural Networks - Conference Proceedings 1, 55–60 (1994) https://doi.org/10.1109/icnn.1994.374138
  * Amini et al. [2020] Amini, A., Schwarting, W., Soleimany, A., Rus, D.: Deep evidential regression. Advances in Neural Information Processing Systems 2020-Decem(NeurIPS), 1–20 (2020)
  * Soleimany et al. [2021] Soleimany, A.P., Amini, A., Goldman, S., Rus, D., Bhatia, S.N., Coley, C.W.: Evidential Deep Learning for Guided Molecular Property Prediction and Discovery. ACS Central Science 7(8), 1356–1367 (2021) https://doi.org/10.1021/acscentsci.1c00546
  * van Amersfoort et al. [2020] Amersfoort, J., Smith, L., Teh, Y.W., Gal, Y.: Uncertainty Estimation Using a Single Deep Deterministic Neural Network. 37th International Conference on Machine Learning, ICML 2020 PartF16814, 9632–9642 (2020)
  * Behler [2016] Behler, J.: Perspective: Machine learning potentials for atomistic simulations. The Journal of Chemical Physics 145(17), 170901 (2016) https://doi.org/10.1063/1.4966192
  * Kulichenko et al. [2021] Kulichenko, M., Smith, J.S., Nebgen, B., Li, Y.W., Fedik, N., Boldyrev, A.I., Lubbers, N., Barros, K., Tretiak, S.: The Rise of Neural Networks for Materials and Chemical Dynamics. The Journal of Physical Chemistry Letters 12(26), 6227–6243 (2021) https://doi.org/10.1021/acs.jpclett.1c01357
  * Hirschfeld et al. [2020] Hirschfeld, L., Swanson, K., Yang, K., Barzilay, R., Coley, C.W.: Uncertainty Quantification Using Neural Networks for Molecular Property Prediction. arXiv preprint arXiv:2005.10036 60(8), 3770–3780 (2020) https://doi.org/10.1021/acs.jcim.0c00502
  * Gurevich and Stuke [2020] Gurevich, P., Stuke, H.: Gradient conjugate priors and multi-layer neural networks. Artificial Intelligence 278 (2020) https://doi.org/10.1016/j.artint.2019.103184
  * Reynolds [2009] Reynolds, D.: Gaussian Mixture Models BT - Encyclopedia of Biometrics, pp. 659–663. Springer, Boston, MA (2009). https://doi.org/10.1007/978-0-387-73003-5{_} . https://doi.org/10.1007/978-0-387-73003-5_196
  * Christensen and Anatole von Lilienfeld [2020] Christensen, A.S., Lilienfeld, O.: On the role of gradients for machine learning of molecular energies and forces. Machine Learning: Science and Technology 1(4) (2020) https://doi.org/10.1088/2632-2153/abba6f
  * Chmiela et al. [2017] Chmiela, S., Tkatchenko, A., Sauceda, H.E., Poltavsky, I., Schütt, K.T., Müller, K.-R.: Machine learning of accurate energy-conserving molecular force fields. Science Advances 3(5), 1603015 (2017) https://doi.org/10.1126/sciadv.1603015
  * Urata et al. [2021] Urata, S., Nakamura, N., Aiba, K., Tada, T., Hosono, H.: How fluorine minimizes density fluctuations of silica glass: Molecular dynamics study with machine-learning assisted force-matching potential. Materials & Design 197, 109210 (2021) https://doi.org/10.1016/j.matdes.2020.109210
  * Tan et al. [2023] Tan, A.R., Urata, S., Goldman, S., Gómez-Bombarelli, R.: Single-model uncertainty quantification in neural network potentials does not consistently outperform model ensembles. Materials Cloud Archive (2023) https://doi.org/10.24435/materialscloud:55-sd
  * Kresse and Furthmüller [1996a] Kresse, G., Furthmüller, J.: Efficiency of ab-initio total energy calculations for metals and semiconductors using a plane-wave basis set. Computational Materials Science 6(1), 15–50 (1996) https://doi.org/10.1016/0927-0256(96)00008-0
  * Kresse and Furthmüller [1996b] Kresse, G., Furthmüller, J.: Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set. Physical Review B 54(16), 11169–11186 (1996) https://doi.org/10.1103/PhysRevB.54.11169
  * Kresse and Joubert [1999] Kresse, G., Joubert, D.: From ultrasoft pseudopotentials to the projector augmented-wave method. Physical Review B 59(3), 1758–1775 (1999) https://doi.org/10.1103/PhysRevB.59.1758
  * Perdew et al. [1996] Perdew, J.P., Burke, K., Ernzerhof, M.: Generalized Gradient Approximation Made Simple. Physical Review Letters 77(18), 3865–3868 (1996) https://doi.org/10.1103/PhysRevLett.77.3865
  * Plimpton [1995] Plimpton, S.: Fast parallel algorithms for short-range molecular dynamics. Journal of computational physics 117(1), 1–19 (1995) https://doi.org/10.1006/JCPH.1995.1039
  * Evans [1983] Evans, D.J.: Computer “experiment” for nonlinear thermodynamics of Couette flow. The Journal of Chemical Physics 78(6), 3297–3302 (1983) https://doi.org/10.1063/1.445195
  * Varivoda et al. [2022] Varivoda, D., Dong, R., Omee, S.S., Hu, J.: Materials Property Prediction with Uncertainty Quantification: A Benchmark Study. arXiv preprint arXiv: 2211.02235 (Ml) (2022)
  * Kahle and Zipoli [2022] Kahle, L., Zipoli, F.: Quality of uncertainty estimates from neural network potential ensembles. Physical Review E 105(1), 1–11 (2022) https://doi.org/10.1103/PhysRevE.105.015311
  * Tran et al. [2020] Tran, K., Neiswanger, W., Yoon, J., Zhang, Q., Xing, E., Ulissi, Z.W.: Methods for comparing uncertainty quantifications for material property predictions. Machine Learning: Science and Technology 1(2), 025006 (2020) https://doi.org/10.1088/2632-2153/ab7e1a
  * Hu et al. [2022] Hu, Y., Musielewicz, J., Ulissi, Z.W., Medford, A.J.: Robust and scalable uncertainty estimation with conformal prediction for machine-learned interatomic potentials. Machine Learning: Science and Technology 3(4), 045028 (2022) https://doi.org/10.1088/2632-2153/aca7b1
  * Hjorth Larsen et al. [2017] Hjorth Larsen, A., Jørgen Mortensen, J., Blomqvist, J., Castelli, I.E., Christensen, R., Dułak, M., Friis, J., Groves, M.N., Hammer, B., Hargus, C., Hermes, E.D., Jennings, P.C., Bjerre Jensen, P., Kermode, J., Kitchin, J.R., Leonhard Kolsbjerg, E., Kubal, J., Kaasbjerg, K., Lysgaard, S., Bergmann Maronsson, J., Maxson, T., Olsen, T., Pastewka, L., Peterson, A., Rostgaard, C., Schiøtz, J., Schütt, O., Strange, M., Thygesen, K.S., Vegge, T., Vilhelmsen, L., Walter, M., Zeng, Z., Jacobsen, K.W.: The atomic simulation environment—a Python library for working with atoms. Journal of Physics: Condensed Matter 29(27), 273002 (2017) https://doi.org/10.1088/1361-648X/aa680e
  * Gabriel et al. [2021] Gabriel, J.J., Paulson, N.H., Duong, T.C., Tavazza, F., Becker, C.A., Chaudhuri, S., Stan, M.: Uncertainty Quantification in Atomistic Modeling of Metals and Its Effect on Mesoscale and Continuum Modeling: A Review. Jom 73(1), 149–163 (2021) https://doi.org/10.1007/s11837-020-04436-6
  * Tumer and Ghosh [1996] Tumer, K., Ghosh, J.: Error Correlation and Error Reduction in Ensemble Classifiers. Connection Science 8(3-4), 385–404 (1996) https://doi.org/10.1080/095400996116839
  * Mendes-Moreira et al. [2012] Mendes-Moreira, J., Soares, C., Jorge, A.M., De Sousa, J.F.: Ensemble approaches for regression: A survey. ACM Computing Surveys 45(1) (2012) https://doi.org/10.1145/2379776.2379786
  * Seitzer et al. [2022] Seitzer, M., Tavakoli, A., Antic, D., Martius, G.: On the Pitfalls of Heteroscedastic Uncertainty Estimation with Probabilistic Neural Networks, 1–24 (2022)
  * Beest et al. [1990] Beest, B.W.H.v., Kramer, G.J.J., Van Beest, B.W.H.H., Kramer, G.J.J., Van Santen, R.A.: Force fields for silicas and aluminophosphates based on ab initio calculations. Physical Review Letters 64(16), 1955–1958 (1990) https://doi.org/10.1103/PhysRevLett.64.1955
  * Tangney and Scandolo [2002] Tangney, P., Scandolo, S.: An ab initio parametrized interatomic force field for silica. Journal of Chemical Physics 117(19), 8898–8904 (2002) https://doi.org/10.1063/1.1513312
  * Tsuneyuki et al. [1988] Tsuneyuki, S., Tsukada, M., Aoki, H., Matsui, Y.: First-Principles Interatomic Potential of Silica Applied to Molecular Dynamics. Phys. Rev. Lett. 61(7), 869–872 (1988) https://doi.org/10.1103/PhysRevLett.61.869
  * Pedone et al. [2022] Pedone, A., Bertani, M., Brugnoli, L., Pallini, A.: Interatomic potentials for oxide glasses: Past , present , and future. Journal of Non-Crystalline Solids: X 15(August), 100115 (2022) https://doi.org/10.1016/j.nocx.2022.100115
  * Urata et al. [2021] Urata, S., Nakamura, N., Tada, T., Hosono, H.: Molecular dynamics study on the co‐doping effect of Al2O3 and fluorine to reduce Rayleigh scattering of silica glass. Journal of the American Ceramic Society 104(10), 5001–5015 (2021) https://doi.org/10.1111/jace.17774
  * Balyakin et al. [2020] Balyakin, I.A., Rempel, S.V., Ryltsev, R.E., Rempel, A.A.: Deep machine learning interatomic potential for liquid silica. Physical Review E 102(5), 52125 (2020) https://doi.org/10.1103/PhysRevE.102.052125
  * Koh and Liang [2017] Koh, P.W., Liang, P.: Understanding Black-box Predictions via Influence Functions. 34th International Conference on Machine Learning, ICML 2017 4, 2976–2987 (2017)

in 1,…,22 See pages \x of SupportingInformation.pdf