# NeuralNEB – Neural Networks can find Reaction Paths Fast

Mathias Schreiner1 Arghya Bhowmik2 Tejs Vegge2 Peter Bjørn Jørgensen2 Ole
Winther1,3,4
1 DTU Compute, Technical University of Denmark (DTU)
2 DTU Energy, Technical University of Denmark (DTU)
3 Department of Biology, University of Copenhagen (UCph)
4 Genomic Medicine, Copenhagen University Hospital, Rigshospitalet

###### Abstract

Quantum mechanical methods like Density Functional Theory (DFT) are used with
great success alongside efficient search algorithms for studying kinetics of
reactive systems. However, DFT is prohibitively expensive for large scale
exploration. Machine Learning (ML) models have turned out to be excellent
emulators of small molecule DFT calculations and could possibly replace DFT in
such tasks. For kinetics, success relies primarily on the models’ capability
to accurately predict the Potential Energy Surface (PES) around transition-
states and Minimal Energy Paths (MEPs). Previously this has not been possible
due to scarcity of relevant data in the literature. In this paper we train
state of the art equivariant Graph Neural Network (GNN)-based models on around
10.000 elementary reactions from the recently published Transition1x dataset.
We apply the models as potentials for the Nudged Elastic Band (NEB) algorithm
and achieve a Mean Average Error (MAE) of 0.13$\pm$0.03 eV on barrier energies
on unseen reactions. We compare the results against equivalent models trained
on QM9x and ANI1x. We also compare with and outperform Density Functional
based Tight Binding (DFTB) on both accuracy and required computational
resources. The implication is that ML models are now at a level where they,
given relevant data, aptly can be applied for studying reaction kinetics.

Figure 1:  Minimal Energy Paths (MEPs) found with Nudged Elastic Band (NEB)
applying the Graph Neural Network (GNN) architecture PaiNN trained on the
Transition1x dataset and Density Functional Theory (DFT) as potentials. The
MEPs are projected onto the planes in structural space, intersecting product,
reactant and transition-state of the converged MEPs. The PES has been
calculated on the planes in the vicinity of the MEPs with the respective
potential and is shown with colors and on the z-axis. The x and y-axes are
basis vectors describing the plane. The reaction involves a H-transfer coupled
with a C-C bond formation on C6H8. The reaction can be seen as a GIF by
following this link.

##  1 Introduction

ML models and especially GNNs 1, 2 have turned out to be potent emulators of
DFT potentials for small molecules 3, 4, 5, 6, 7, thanks to their remarkable
ability to find complex relations in high dimensional data. They have a
complexity-scaling orders of magnitudes lower than classic Quantum Mechanics
(QM) methods, but have in recent years achieved comparable accuracy 8, 9, 10,
11, 12. The capability of these models is manifested by their success in tasks
beyond simple prediction of molecular features such as structural optimization
or studying finite-temperature dynamical properties through molecular dynamics
13, 14. Despite their achievements, there has only been limited success in
applying ML-models as potentials for transition search algorithms. The
earliest work studied simple diatomic molecule dissociation and achieved
acceptable accuracy with tens of thousands of data points 15. Other works have
had success by limiting their scope to studying single or few reactions but
sacrificing the generality of the approach 16, 17, 18. Attempts to study
reactive systems with Gaussian Processs (GPs) 19 have been successful too, but
the GP is trained on the particular atomic system, sacrificing speed for
generality by requiring expensive DFT calculations at inference time.
Transition-states are notoriously hard to find as there is no well-defined
gradient on the PES to guide traditional optimization algorithms towards them.
A wealth of algorithms have been proposed to solve this problem – one is the
NEB 20 algorithm, which works by interpolating an initial path between
reactant and product and iteratively updating it to minimize energy by using
information about the PES. It shares a common bottleneck with other transition
search algorithms – the necessity to repeatedly evaluate energy and atomic
forces of molecular configurations, which is extremely costly, especially if
ab-initio or electron DFT calculations are used21.

Recent advances in ML have not alleviated the bottleneck as even modern Neural
Network (NN) architectures have not proved proficient potential approximators
for this type of application. The fault lies primarily with available data in
the literature rather than the models’ expressiveness 22. Most quantum
mechanical datasets are focused on molecular configurations in or near
equilibrium 23, 24, 25, 26. Without configurations on and around reaction
pathways in the training data, ML models cannot learn the interatomic
interactions that occur during chemical reactions and cannot reliably be
applied for transition-state search.

We compare ML models against DFTB, 27 a fast approximation to DFT that is
often used for fast screening of large quantities of configurations with an
acceptable trade-off between accuracy and speed, and our models outperform
DFTB with a factor five in accuracy and a factor 2.5 in CPU time.

In this work, we bridge generalization, speed, and accuracy for transition-
state search by applying Polarizable Atom interaction Neural Network (PaiNN)
models as surrogate potentials for DFT. We build on and showcase the utility
of our previous paper 28 where we released Transition1x, a dataset constituted
by DFT calculations for 10 million molecular configurations, all sampled
around reaction pathways from 10.000 elementary reactions. It is clear from
the results of this paper, that for precise modeling of transition-state
regions, and, consequently, transition states and barrier energies, hitherto
popular benchmark datasets have had insufficient relevant data. On the other
hand, training ML potentials on the Transition1x dataset allows for accurate
modeling of PESs in transition-state regions, underlining that relevant and
available data in the literature is as important as the efficiency of
available models.

Reliable and fast analysis of reaction kinetics through ML will bring the
whole field of computational chemistry a considerable step closer to the
ultimate goal, a virtual laboratory, hyper-accelerating the discovery of
reaction mechanisms for synthesizing drugs and materials.

##  2 Methods

###  2.1 Nudged Elastic Band

NEB 20 is a method for finding MEP and transition-state given product and
reactant of a chemical reaction. It does so by iteratively nudging an
interpolated path between the reaction endpoints in the direction of the force
perpendicular to the path. Once the perpendicular force converges to zero, NEB
reports the maximal-energy configuration along the path as the transition-
state. The path is represented by an array of molecular configurations called
images, and there is no guarantee that, at convergence, the maximal energy
image corresponds to the maximal energy along the path. The maximum might lie
between two images. Climbing Image Nudged Elastic Band (CINEB) 29 addresses
this problem by letting the transition-state candidate (the maximal energy
image) further maximize its energy by following the gradient on the PES
parallel to the current path between iterations. If the current path has not
converged properly, the climbing image can pull the predicted MEP off the true
MEP and therefore, the path is first relaxed with regular NEB before turning
on CINEB. The MEP is considered converged once the maximal perpendicular force
on the path is below a threshold of 0.05 eVÅ-1. The spring constant between
images on the path is set to 0.1 eVÅ-2, and ten images are used to represent
the path.

###  2.2 Initial Path Generation

The endpoints of the reaction have to be minimized in their respective minima
before running NEB– otherwise the energetic difference between reactant and
transition-state cannot be evaluated properly. A configuration is considered
relaxed if the norm of the forces acting on it is below 0.01 eVÅ-1. Once the
endpoints have been minimized, the initial guess for the MEP is found by
running NEB with the Image Dependent Pair Potential (IDPP) 30 on a linearly
interpolated path between reactant and product. IDPP is an inexpensive
potential specifically designed to generate physically realistic MEP guesses
for NEB at an extremely low computational cost.

###  2.3 Optimizers

Reactants and products are relaxed using the BFGS 31 optimizer with
$\alpha=70$ and a maximal step size of 0.03 Å in configurational space. The
MEP is found with an optimizer 32 designed to reduce the computational cost of
transition-state search algorithms by applying an adaptive time step selection
algorithm with $\alpha=0.01$ and rtol $=0.1$, and a preconditioning scheme to
the PES given an estimate of its curvature.

##  3 Data

We train all models on ANI1x 24, QM9x 33, Transition1x 28. All datasets are
calculated with the 6-31G(d) 34 basis set and $\omega$B97x 35 functional which
has an accuracy comparable to the gold standard but expensive high-level
CCSD(T) 36 37 calculations. Given the compatibility of the datasets, it is
possible to train on either dataset alone or combinations of them to leverage
all of their strengths.

###  3.1 ANI1x

ANI1x 38 aims to provide varied data of off-equilibrium molecular
configurations by perturbing equilibrium configurations with pseudo molecular
dynamics. The data is collected through an active learning technique called
Query by Committee; an automated data diversification process that trains an
ensemble (committee) of models on a dataset and accepts or rejects new
proposed data based on the disagreement of models in the committee. The
assumption is that if the committee disagrees the data is sufficiently
different from what has already been learned, and the proposed data should be
included in the analysis. The procedure for proposing data and evaluating it
with the committee is cheap compared to the calculation of data using DFT. The
dataset is consecutively expanded by alternating between training committees
and adding new data points based on the committee uncertainty. In total, ANI1x
contains force and energy calculations for approximately 5 million
configurations.

###  3.2 Transition1x

We have recently published Transition1x 28, a dataset providing a collection
of molecular configurations on and along reaction paths for approximately
10.000 reactions. The reactions consist of up to 7 heavy atoms, including C,
N, and O. Transition events are rare, and it is not possible to collect
sufficient data in relevant regions by simple molecular dynamics if the
intention is to train NNs models to understand chemical reactions.
Transition1x addresses this problem by sampling molecular configurations
around reaction pathways proposed by NEB, using DFT as potential. The
procedure resulted in approximately 10 million DFT calculations that were
collected and saved during the process and constitute the dataset.
Transition1x is available through the repository
https://gitlab.com/matschreiner/Transition1x which includes data loaders and
scripts for downloading the dataset and generating ASE-database files.

###  3.3 QM9 and QM9x

QM9 33 is a dataset of 135k small organic molecules with various chemical
properties that has served as the benchmark for many existing ML methods for
quantum chemistry. All molecules in QM9 are in equilibrium. We have
recalculated QM9 with the 6-31G(d) basis set and $\omega$B97x functional to
make it compatible with Transition1x and ANI1x, and we refer to the
recalculated dataset as QM9x. Molecular configurations recalculated in the new
potential are not necessarily in equilibrium as the potential shifts when
changing functional and basis sets. QM9x is available through the repository
https://gitlab.com/matschreiner/QM9x which includes data loaders and scripts
for downloading the dataset and generating ASE-database files.

###  3.4 Models and Training

Message Passing Neural Networks (MPNNs) 39 are a class of GNNs 1, 2 that build
their internal graph representation by running a series of message passing
steps. A single message passing step consists of two distinct operations: i)
Message Dispatching, each node computes a message given its state (and
possibly information about the edge connecting to – and the state of the
receiving node) and sends it to its neighbors. ii) State Update, incoming
messages are collected with an aggregation function, and are used to
simultaneously update the internal representation of all nodes. After the
message-passing phase, a readout function extracts the inner representation of
the nodes and computes a final feature vector of the graph for downstream
tasks. In the case of molecules, interesting properties are energy and forces
where conservative force fields can be computed via the back-propagation
algorithm as the negative gradient of the energy wrt. coordinates of the
atoms.

The PaiNN model 40 was used for all experiments – it is a GNN architecture
that implements rotationally equivariant representations for prediction of
tensorial properties of graph structures. We refer to the literature for
further details 40. A cut-off radius of 5 Å was used to generate the initial
molecular graph. All models have three message passing steps and 256 units in
each hidden layer, and are trained using the ADAM 41 optimizer with learning
rate $10^{-3}$ on training examples from QM9x, ANI1x, and Transition1x. A
batchsize of 75 was used for all datasets and a maximum of $10^{6}$ training
steps was allowed - however, models training on ANI1x and Transition1x reached
maximal scores on validation data after around $6\cdot 10^{5}$ steps.
Transition1x was stratified by reactions without attention to reaction-type.
500 reactions were set aside for testing such that no data from around any
test reaction was seen during training. 90% of the remaining data were used
for training and the last 10% were used for validation and early stopping.
ANI1x was stratified by chemical formula such that test, validation and
training sets consist of chemical formulas unique to that set. QM9x was split
randomly. In the case of QM9x and ANI1x, 80% of the data was used for
training, 10% for testing, and 10% was used for validation and early stopping.
In QM9x all configurations are unique as they are in distinct equilibria and
can therefore be split randomly. No attention was paid to the molecular
scaffold. For ANI1x, it is necessary to split on chemical formula to ensure
that configurations across splits are significantly different. Each chemical
formula contains similar configurations, since data is generated by randomly
perturbing identical initial configurations. In Transition1x the reactions
follow unique paths in chemical space, and thus configurations between
reactions are different. The Transition1x dataset forms a reaction network
where the product from one reaction is the reactant for the next. Therefore,
it is nontrivial to split it to get no shared reactants/products across train
validation and test sets. We believe that our splitting is still useful
because it achieves unique reaction paths and transition states between
datasets.

##  4 Results

| Barrier [ eV ] | NEB Convergence
---|---|---
| MAE | RMSE | Rate | Avg. CPU Time | Avg. Iterations
ANI1x | 0.43(1) | 0.56(3) | 71.9% | 36s | 130.91
T1x | 0.13(3) | 0.24(1) | 80.5% | 28s | 111.97
QM9x | 3.4(1) | 3.6(1) | 23.4% | 26s | 110.34
DFTB | 0.66 | 0.79 | 72.0% | 81s | 115.65
DFT | - | - | 84.1% | 12h14m43s | 100.74
Table 1:  Performance of various potentials used for Nudged Elastic Band (NEB) when compared to Density Functional Theory (DFT). ANI1x, Transition1x and QM9x indicate PaiNN models trained on the respective dataset. The Barrier column displays the Mean Average Error (MAE) and Root Mean Squared Error (RMSE) of barrier predictions, where the individual error is the difference between the barrier as predicted when using DFT as potential vs. using the surrogate potential. The convergence rate is the percentage of reactions that converged. Average CPU time is CPU time spent per reaction. Average iterations is the average number of Minimal Energy Path (MEP) updates before convergence.  Figure 2:  Histogram of barrier errors. The x-axis shows errors between reaction barriers calculated using Density Functional Theory (DFT) and surrogate potentials for Nudged Elastic Band (NEB). The x-axis has been truncated at +/- 1.5 eV error (see appendix for the full plot). The y-axis shows the frequency of each bin. Green, red and blue display results from PaiNN models trained on Transition1x, QM9x and ANI1x, respectively. Yellow displays results from Density Functional based Tight Binding (DFTB). The QM9x model has such a low convergence frequency, and general barrier error, that the model does not show in the plot.  Figure 3:  Comparison of reaction barriers found with Nudged Elastic Band (NEB) using Density Functional Theory (DFT) as potential on the x-axis vs. various surrogate potentials on the y-axis. Green, red and blue markers are PaiNN models trained on The Transition1x, QM9x, and ANI1x datasets respectively. Yellow is Density Functional based Tight Binding (DFTB). Points on the dashed line have been calculated perfectly. | MAE [ eV ] | RMSE [ eV ] | Systematic Error [ eV ]
---|---|---|---
ANI1x | 0.40 | 0.56 | 0.16
DFTB | 0.47 | 0.57 | -0.55
Transition1x | 0.17 | 0.30 | -0.07
Table 2:  MAE and RMSE of barrier errors found by PaiNN trained on
Transition1x and ANI1x and DFTB, after correcting for systematic error.

Table 1 shows the overall findings of the paper. Each row displays the
performance of a surrogate potential, where datasets in the leftmost column
indicate PaiNN models trained on the given dataset. The barrier error is the
difference in barrier heights found when applying DFT as potential for NEB
versus when applying the respective potential. As different initializations of
parameters in equivalent architectures result in variations in the trained
models, five models were trained on each dataset without changing training,
testing and validation splits. These models were used as potentials for NEB,
the table displays the mean and standard deviation of MAEs and RMSEs of each
model.

The best models are trained on Transition1x, with the lowest MAE and RMSE and
the highest convergence ratio. The QM9x models have only seen data very close
to equilibrium and have not learned the structure of the PES between
equilibria which makes it unable to converge in most cases. In general DFT
performs the best in terms of convergence rate and average iterations run, but
it comes at a steep price, running almost a factor 1500 times slower than the
ML potentials. DFTB is the go-to fast potential, but the models trained on
Transition1x are roughly 2.5 times faster and five times as accurate.

Figure 1, on the frontpage, displays MEPs calculated with NEB using DFT and
PaiNN trained on Transition1x side by side. Each MEP is projected onto a plane
in configurational space spanned by the reaction’s transition-state, product,
and reactant. The x and y axes are basis vectors describing the plane in units
of Å, and the z-axis and color-coding show the atomization energy of
configurations in the plane in eV. Not only does PaiNN trained on the
Transition1x accurately calculate the barrier energy for the reaction, but it
also correctly identifies the plane spanned by the configurations defining the
reaction, and calculates an almost identical PES in the vicinity of the MEP.
Each MEP is projected from a high dimensional space onto the plane, and
therefore, only the atomization energy of equilibria and transition-states are
shown correctly in the plot. At these points, the MEP intersects with the
plane. The intermediate points have energies slightly shifted up the sides of
the energy valley. The MEP does not necessarily lie in the plane, and since
the MEP represents the energy valley, projecting it onto the plane, will
increase the energy. The $\times$ symbols on the surfaces are projections of
images predicted by NEB and the dashed lines connecting them are cubic spline
interpolations. The importance of accurate predictions in the vicinity of the
MEP is clear, as these calculations will guide the search for the transition-
state. The Transition1x model predicts smooth and well-behaved PESs resembling
DFT.

Figure 2 and 3 tell similar stories. Figure 2 is a histogram of barrier errors
where the error is the difference between activation energy found using the
surrogate potential and DFT. The Transition1x model is precise and accurate,
with a sharp peak around zero, whereas DFTB and ANI1x have wider spreads with
means below and above zero, respectively. The QM9x model is plotted on the
histogram, but due to high errors and low convergence, only a few calculated
barriers fall within an error of $\pm$1.5 eV, as shown in the figure. See
appendix for an equivalent figure without truncated x-axis.

Figure 3 compares activation energies found with DFT on the x-axis with those
found using various surrogate potentials on the y-axis. Each marker represents
a single reaction. Predictions from the model trained on Transition1x follow
the $x=y$ line with a MAE of only 0.13 eV. The QM9x model does not have a
proper representation of the transition-state regions as it has not seen that
type of data during training. Often, the QM9x model does not recognize nearby
initial equilibria as minima on the PES, and even before optimizing the MEP,
the reaction endpoints have dropped further on the PES to qualitatively
different endpoints which results in the model calculating the MEP for a
completely different reaction. The algorithm is not set up to detect this, and
as long as the reaction converges, it is included in the analysis. Even when
the QM9x model relaxes the endpoints of the reaction correctly, it either
finds low energy shortcuts in the faulty potential or does not converge, and
as a result the converged reactions are often only the energy difference
between reactant and product. The QM9x dataset was not designed with any type
of molecular dynamics or reaction kinetics in mind, and comparing it to ANI1x
and Transition1x for reaction path search is perhaps inappropriate. However,
given the ubiquity of QM9 in the literature, it is an important point to
convey, that new datasets are required for solving higher order problems in
computational chemistry. The Transition1x and ANI1x models drop in performance
above 5 eV. Data becomes scarcer at higher energies and consequently, models
are less accurate in high energy regions. DFTB and the ANI1x models have
systematic errors in their predictions. The ANI1x models are biased towards
high energies in the transition regions as they have not seen the low energy
valleys connecting equilibria. The DFTB potential systematically predicts
energies too low. In Table 2 the systematic errors are corrected based on the
training data. This leads to a lower test error for the ANI1x and DFTB, but a
higher test error for Transition1x underlining that Transition1x models are
already very accurate.

##  5 Discussion

To train models that can properly step in as surrogate potentials for DFT when
running NEB, it is necessary to have datasets with appropriate data in and
around transition-state regions. Finding reaction barriers with ML models and
NEB is a non-trivial test. ML models, and especially NNs, are known to perform
poorly for out of distribution tasks 42, 43. Table 3 illustrates this clearly
with results for training and testing ML models on various datasets.

Finding reaction barriers with NEB is a much more demanding test of the
models’ capabilities. When running NEB, the PES is swept by the path
connecting endpoints, and data encountered in the process can diverge wildly
from any data seen during testing and training. The model can get caught in
even a small region of high error, or it can be thrown off the correct MEP and
make it unable to converge altogether, so the model must be accurate across
the entire PES.

The reaction paths are represented by ten images in all reactions. A core
strength of NNs is their ability to utilize GPUs to evaluate multiple data
points at once, and in principle, NEB can be run with hundreds of images
instead of tens at little to no additional cost when using NNs as potentials.
We ran experiments with high density paths with the rest of the setup fixed
but saw no improvement in neither accuracy nor convergence speed. The
preconditioning scheme of the NEB optimizer relies on a sparsely populated
path. But this approach could possibly produce robust results by applying
other optimizers.

A clear application of this work is as a screening procedure for complex
reaction networks. Cheap methods, such as permuting bond order matrices, can
be used to automatically generate nodes for entire reaction networks. The
individual reactions can be screened fast using the method before
recalculating entire reaction networks with expensive methods. Usually this is
done with DFTB 27 but running NEB with NNs is faster and more accurate.

##  6 Conclusion

We have trained state-of-the-art GNN potentials on various datasets and used
them as surrogate potentials for DFT when running NEB for transition-state
search. A MAE of 0.13 eV and RMSE of 0.24 eV is achieved with the best model,
compared against running the same set up with DFT. The models converge 80.5%
of the time on unseen reactions. We show that expressive models alone are, by
no means, sufficient for solving complex tasks in quantum chemistry moving
forward, but just as much care has to be put into designing and generating
datasets. We tested 3 different datasets: ANI1x, QM9x and Transition1x and
only models trained on the latter could reliably solve the transition search
task.

Our results show that the future development of the field of ML for quantum
chemistry stands on two legs – the completeness of the available data, and the
expressiveness of the available models. Transition1x is a relatively simple
dataset that deals with four types of atoms. To apply the results of this
paper to general chemistry, larger datasets with more atom types should be
produced. Our results indicate that the machine learning approach scales: With
the right amount of the right data, accuracies at a sufficient level can be
achieved.

## Acknowledgements

The authors acknowledge support from the Novo Nordisk Foundation (SURE,
NNF19OC0057822) and the European Union’s Horizon 2020 research and innovation
program under Grant Agreement No. 957189 (BIG-MAP) and No. 957213
(BATTERY2030PLUS).

Ole Winther also receives support from Novo Nordisk Foundation through the
Center for Basic Machine Learning Research in Life Science (NNF20OC0062606)
and the Pioneer Centre for AI, DNRF grant number P1.

##  7 Code Availability

Code for training PaiNN models and running NEB is available through the
repository https://gitlab.com/matschreiner/neuralneb.

## References

  * Bacciu et al. 2019 Davide Bacciu, Federico Errica, Alessio Micheli, and Marco Podda.  A gentle introduction to deep learning for graphs.  _Neural Networks_ , 129:203–221, 12 2019.  doi: 10.1016/j.neunet.2020.06.006.  URL http://arxiv.org/abs/1912.12693http://dx.doi.org/10.1016/j.neunet.2020.06.006.
  * Zhou et al. 2021 Jie Zhou, Ganqu Cui, Shengding Hu, Zhengyan Zhang, Cheng Yang, Zhiyuan Liu, Lifeng Wang, Changcheng Li, and Maosong Sun.  Graph neural networks: A review of methods and applications.  2021\.  doi: 10.1016/j.aiopen.2021.01.001.  URL https://doi.org/10.1016/j.aiopen.2021.01.001.
  * Faber et al. 2017a Felix A. Faber, Luke Hutchison, Bing Huang, Justin Gilmer, Samuel S. Schoenholz, George E. Dahl, Oriol Vinyals, Steven Kearnes, Patrick F. Riley, and O. Anatole Von Lilienfeld.  Prediction errors of molecular machine learning models lower than hybrid dft error.  _Journal of Chemical Theory and Computation_ , 13:5255–5264, 11 2017a.  ISSN 15499626.  doi: 10.1021/ACS.JCTC.7B00577/SUPPL_FILE/CT7B00577_SI_001.PDF.  URL https://pubs.acs.org/doi/abs/10.1021/acs.jctc.7b00577.
  * Westermayr et al. 2021 Julia Westermayr, Michael Gastegger, Kristof T. Schütt, and Reinhard J. Maurer.  Perspective on integrating machine learning into computational chemistry and materials science.  _The Journal of Chemical Physics_ , 154:230903, 6 2021.  ISSN 0021-9606.  doi: 10.1063/5.0047760.  URL https://aip.scitation.org/doi/abs/10.1063/5.0047760.
  * Campbell et al. 2020 Stuart I Campbell, Daniel B Allan, and Andi M Barbour.  Machine learning for the solution of the schrödinger equation.  _Machine Learning: Science and Technology_ , 1:013002, 4 2020.  ISSN 2632-2153.  doi: 10.1088/2632-2153/AB7D30.  URL https://iopscience.iop.org/article/10.1088/2632-2153/ab7d30https://iopscience.iop.org/article/10.1088/2632-2153/ab7d30/meta.
  * Behler and Parrinello 2007a Jörg Behler and Michele Parrinello.  Generalized neural-network representation of high-dimensional potential-energy surfaces.  _Physical review letters_ , 98(14):146401, 2007a.
  * Westermayr and Marquetand 2021 Julia Westermayr and Philipp Marquetand.  Machine learning for electronically excited states of molecules.  _Chemical Reviews_ , 121:9873–9926, 8 2021.  ISSN 15206890.  doi: 10.1021/ACS.CHEMREV.0C00749/ASSET/IMAGES/LARGE/CR0C00749_0011.JPEG.  URL https://pubs.acs.org/doi/full/10.1021/acs.chemrev.0c00749.
  * Reymond 2015 Jean Louis Reymond.  The chemical space project.  _Accounts of Chemical Research_ , 48:722–730, 3 2015.  ISSN 15204898.  doi: 10.1021/AR500432K/ASSET/IMAGES/LARGE/AR-2014-00432K_0003.JPEG.  URL https://pubs.acs.org/doi/full/10.1021/ar500432k.
  * Behler and Parrinello 2007b Jörg Behler and Michele Parrinello.  Generalized neural-network representation of high-dimensional potential-energy surfaces.  _Physical Review Letters_ , 98:146401, 4 2007b.  ISSN 00319007.  doi: 10.1103/PHYSREVLETT.98.146401/FIGURES/4/MEDIUM.  URL https://journals.aps.org/prl/abstract/10.1103/PhysRevLett.98.146401.
  * Behler 2015 Jörg Behler.  Constructing high-dimensional neural network potentials: A tutorial review.  _International Journal of Quantum Chemistry_ , 115:1032–1050, 8 2015.  ISSN 1097-461X.  doi: 10.1002/QUA.24890.  URL https://onlinelibrary.wiley.com/doi/full/10.1002/qua.24890https://onlinelibrary.wiley.com/doi/abs/10.1002/qua.24890https://onlinelibrary.wiley.com/doi/10.1002/qua.24890.
  * Faber et al. 2017b Felix A. Faber, Luke Hutchison, Bing Huang, Justin Gilmer, Samuel S. Schoenholz, George E. Dahl, Oriol Vinyals, Steven Kearnes, Patrick F. Riley, and O. Anatole Von Lilienfeld.  Prediction errors of molecular machine learning models lower than hybrid dft error.  _Journal of Chemical Theory and Computation_ , 13:5255–5264, 11 2017b.  ISSN 15499626.  doi: 10.1021/ACS.JCTC.7B00577/SUPPL_FILE/CT7B00577_SI_001.PDF.  URL https://pubs.acs.org/doi/abs/10.1021/acs.jctc.7b00577.
  * Gilmer et al. 2017a Justin Gilmer, Samuel S Schoenholz, Patrick F Riley, Oriol Vinyals, and George E Dahl.  Neural message passing for quantum chemistry.  2017a.
  * Kaappa et al. 2021 Sami Kaappa, Casper Larsen, and Karsten Wedel Jacobsen.  Atomic structure optimization with machine-learning enabled interpolation between chemical elements.  _Physical Review Letters_ , 127, 7 2021.  doi: 10.1103/PhysRevLett.127.166001.  URL http://arxiv.org/abs/2107.01055http://dx.doi.org/10.1103/PhysRevLett.127.166001.
  * Wang et al. 2020 Jiaqi Wang, Seungha Shin, and Sangkeun Lee.  Interatomic potential model development: Finite-temperature dynamics machine learning.  _Advanced Theory and Simulations_ , 3:1900210, 2 2020.  ISSN 2513-0390.  doi: 10.1002/ADTS.201900210.  URL https://onlinelibrary.wiley.com/doi/full/10.1002/adts.201900210https://onlinelibrary.wiley.com/doi/abs/10.1002/adts.201900210https://onlinelibrary.wiley.com/doi/10.1002/adts.201900210.
  * Malshe et al. 2007 M Malshe, LM Raff, MG Rockley, M Hagan, Paras M Agrawal, and R Komanduri.  Theoretical investigation of the dissociation dynamics of vibrationally excited vinyl bromide on an ab initio potential-energy surface obtained using modified novelty sampling and feedforward neural networks. ii. numerical application of the method.  _The Journal of chemical physics_ , 127(13):134105, 2007.
  * Lu et al. 2018 Xiaoxiao Lu, Qingyong Meng, Xingan Wang, Bina Fu, and Dong H Zhang.  Rate coefficients of the h+ h2o2→ h2+ ho2 reaction on an accurate fundamental invariant-neural network potential energy surface.  _The Journal of chemical physics_ , 149(17):174303, 2018.
  * Young et al. 2021 Tom A Young, Tristan Johnston-Wood, Volker L Deringer, and Fernanda Duarte.  A transferable active-learning strategy for reactive molecular force fields.  _Chemical science_ , 12(32):10944–10955, 2021\.
  * Manzhos and Carrington Jr 2020 Sergei Manzhos and Tucker Carrington Jr.  Neural network potential energy surfaces for small molecules and reactions.  _Chemical Reviews_ , 121(16):10187–10217, 2020\.
  * Koistinen et al. 2017 Olli Pekka Koistinen, Freyja B. Dagbjartsdóttir, Vilhjálmur Ásgeirsson, Aki Vehtari, and Hannes Jónsson.  Nudged elastic band calculations accelerated with gaussian process regression.  _Journal of Chemical Physics_ , 147, 2017.  ISSN 00219606.  doi: 10.1063/1.4986787.
  * Sheppard et al. 2008 Daniel Sheppard, Rye Terrell, and Graeme Henkelman.  Optimization methods for finding minimum energy paths.  _The Journal of Chemical Physics_ , 128:134106, 4 2008.  ISSN 0021-9606.  doi: 10.1063/1.2841941.  URL https://aip.scitation.org/doi/abs/10.1063/1.2841941.
  * Heinen et al. 2020 Stefan Heinen, Max Schwilk, Guido Falk von Rudorff, and O Anatole von Lilienfeld.  Machine learning the computational cost of quantum chemistry.  _Machine Learning: Science and Technology_ , 1(2):025002, 2020.
  * von Lilienfeld et al. 2020 O. Anatole von Lilienfeld, Klaus Robert Müller, and Alexandre Tkatchenko.  Exploring chemical compound space with quantum-based machine learning.  _Nature Reviews Chemistry 2020 4:7_ , 4:347–358, 6 2020\.  ISSN 2397-3358.  doi: 10.1038/s41570-020-0189-9.  URL https://www.nature.com/articles/s41570-020-0189-9.
  * Smith et al. 2017 Justin S. Smith, Olexandr Isayev, and Adrian E. Roitberg.  Ani-1, a data set of 20 million calculated off-equilibrium conformations for organic molecules.  _Scientific Data 2017 4:1_ , 4:1–8, 12 2017.  ISSN 2052-4463.  doi: 10.1038/sdata.2017.193.  URL https://www.nature.com/articles/sdata2017193.
  * Smith et al. 2020 Justin S. Smith, Roman Zubatyuk, Benjamin Nebgen, Nicholas Lubbers, Kipton Barros, Adrian E. Roitberg, Olexandr Isayev, and Sergei Tretiak.  The ani-1ccx and ani-1x data sets, coupled-cluster and density functional theory properties for molecules.  _Scientific Data 2020 7:1_ , 7:1–10, 5 2020.  ISSN 2052-4463.  doi: 10.1038/s41597-020-0473-z.  URL https://www.nature.com/articles/s41597-020-0473-z.
  * Fink and Raymond 2007 Tobias Fink and Jean Louis Raymond.  Virtual exploration of the chemical universe up to 11 atoms of c, n, o, f: Assembly of 26.4 million structures (110.9 million stereoisomers) and analysis for new ring systems, stereochemistry, physicochemical properties, compound classes, and drug discovery.  _Journal of Chemical Information and Modeling_ , 47:342–353, 2007.  ISSN 1549960X.  doi: 10.1021/CI600423U.
  * Fink et al. 2005 Tobias Fink, Heinz Bruggesser, and Jean Louis Reymond.  Virtual exploration of the small-molecule chemical universe below 160 daltons.  _Angewandte Chemie - International Edition_ , 44:1504–1508, 2 2005.  ISSN 14337851.  doi: 10.1002/ANIE.200462457.
  * Seifert and Joswig 2012 Gotthard Seifert and Jan Ole Joswig.  Density-functional tight binding—an approximate density-functional theory method.  _Wiley Interdisciplinary Reviews: Computational Molecular Science_ , 2:456–465, 5 2012.  ISSN 1759-0884.  doi: 10.1002/WCMS.1094.  URL https://onlinelibrary.wiley.com/doi/full/10.1002/wcms.1094https://onlinelibrary.wiley.com/doi/abs/10.1002/wcms.1094https://wires.onlinelibrary.wiley.com/doi/10.1002/wcms.1094.
  * Schreiner et al. 2022 Mathias Schreiner, Arghya Bhowmik, Tejs Vegge, Jonas Busk, and Ole Winther.  Transition1x.  6 2022.  doi: 10.6084/m9.figshare.19614657.v4.  URL https://figshare.com/articles/dataset/Transition1x/19614657.
  * Henkelman et al. 2000 Graeme Henkelman, Blas P. Uberuaga, and Hannes Jónsson.  A climbing image nudged elastic band method for finding saddle points and minimum energy paths.  _The Journal of Chemical Physics_ , 113:9901, 11 2000.  ISSN 0021-9606.  doi: 10.1063/1.1329672.  URL https://aip.scitation.org/doi/abs/10.1063/1.1329672.
  * Smidstrup et al. 2014 Søren Smidstrup, Andreas Pedersen, Kurt Stokbro, and Hannes Jónsson.  Improved initial guess for minimum energy path calculations.  _The Journal of Chemical Physics_ , 140:214106, 6 2014.  ISSN 0021-9606.  doi: 10.1063/1.4878664.  URL https://aip.scitation.org/doi/abs/10.1063/1.4878664.
  * Broyden 1970 C. G. Broyden.  The convergence of a class of double-rank minimization algorithms 1. general considerations.  _IMA Journal of Applied Mathematics_ , 6:76–90, 3 1970\.  ISSN 0272-4960.  doi: 10.1093/IMAMAT/6.1.76.  URL https://academic.oup.com/imamat/article/6/1/76/746016.
  * Makri et al. 2019 Stela Makri, Christoph Ortner, and James R. Kermode.  A preconditioning scheme for minimum energy path finding methods.  _The Journal of Chemical Physics_ , 150:094109, 3 2019.  ISSN 0021-9606.  doi: 10.1063/1.5064465.  URL https://aip.scitation.org/doi/abs/10.1063/1.5064465http://creativecommons.org/licenses/by/4.0/.
  * Ramakrishnan et al. 2014 Raghunathan Ramakrishnan, Pavlo O. Dral, Matthias Rupp, and O. Anatole Von Lilienfeld.  Quantum chemistry structures and properties of 134 kilo molecules.  _Scientific Data 2014 1:1_ , 1:1–7, 8 2014.  ISSN 2052-4463.  doi: 10.1038/sdata.2014.22.  URL https://www.nature.com/articles/sdata201422.
  * Ditchfield et al. 2003 R. Ditchfield, W. J. Hehre, and J. A. Pople.  Self-consistent molecular-orbital methods. ix. an extended gaussian-type basis for molecular-orbital studies of organic molecules.  _The Journal of Chemical Physics_ , 54:724, 9 2003.  ISSN 0021-9606.  doi: 10.1063/1.1674902.  URL https://aip.scitation.org/doi/abs/10.1063/1.1674902.
  * Chai and Head-Gordon 2008 Jeng Da Chai and Martin Head-Gordon.  Systematic optimization of long-range corrected hybrid density functionals.  _The Journal of Chemical Physics_ , 128:084106, 2 2008.  ISSN 0021-9606.  doi: 10.1063/1.2834918.  URL https://aip.scitation.org/doi/abs/10.1063/1.2834918.
  * Riley et al. 2010 Kevin E. Riley, Michal Pitončák, Petr Jurecčka, and Pavel Hobza.  Stabilization and structure calculations for noncovalent interactions in extended molecular systems based on wave function and density functional theories.  _Chemical Reviews_ , 110:5023–5063, 9 2010.  ISSN 00092665.  doi: 10.1021/CR1000173.
  * Thanthiriwatte et al. 2011 Kanchana S. Thanthiriwatte, Edward G. Hohenstein, Lori A. Burns, and C. David Sherrill.  Assessment of the performance of dft and dft-d methods for describing distance dependence of hydrogen-bonded interactions.  _Journal of Chemical Theory and Computation_ , 7:88–96, 1 2011.  ISSN 15499618.  doi: 10.1021/CT100469B.
  * Smith et al. 2018 Justin S. Smith, Ben Nebgen, Nicholas Lubbers, Olexandr Isayev, and Adrian E. Roitberg.  Less is more: Sampling chemical space with active learning.  _The Journal of Chemical Physics_ , 148:241733, 5 2018.  ISSN 0021-9606.  doi: 10.1063/1.5023802.  URL https://aip.scitation.org/doi/abs/10.1063/1.5023802.
  * Gilmer et al. 2017b Justin Gilmer, Samuel S Schoenholz, Patrick F Riley, Oriol Vinyals, and George E Dahl.  Neural message passing for quantum chemistry.  2017b.
  * Schütt et al. 2021 Kristof T Schütt, Sch¨ Schütt, Oliver T Unke, and Michael Gastegger.  Equivariant message passing for the prediction of tensorial properties and molecular spectra.  2021\.
  * 41 Diederik P Kingma and Jimmy Lei Ba.  Adam: A method for stochastic optimization.
  * Henriksson et al. 2021 Jens Henriksson, Christian Berger, Markus Borg, Lars Tornberg, Sankar Raman Sathyamoorthy, and Cristofer Englund.  Performance analysis of out-of-distribution detection on various trained neural networks.  2021\.  URL https://www.iso.org/deliverables-all.html.
  * Zhang et al. 2021 Lily H Zhang, Mark Goldstein, and Rajesh Ranganath.  Understanding failures in out-of-distribution detection with deep generative models.  2021\.

##  Appendix A PaiNN Performance On Test Data

Table displaying results of the models when training and testing on various
datasets. In all test set-ups the models that perform best, are models that
have been trained on training data from the corresponding dataset.

|  | Energy [eV] | Forces [eV/Å]
---|---|---|---
Trained on | Tested on | RMSE | MAE | RMSE | MAE
ANI1x |  | 0.04(1) | 0.02(0) | 0.04(0) | 0.01(0)
Transition1x | ANI1x | 0.35(2) | 0.22(1) | 0.34(2) | 0.08(0)
QM9x |  | 3.03(2) | 2.32(1) | 1.3(7) | 0.56(2)
ANI1x |  | 0.61(7) | 0.28(2) | 0.5(1) | 0.10(5)
Transition1x | Transition1x | 0.10(1) | 0.05(0) | 0.10(0) | 0.04(1)
QMx |  | 2.61(2) | 1.42(1) | 0.43(0) | 0.19(0)
ANI1x |  | 0.13(0) | 0.12(0) | 0.05(0) | 0.02(0)
Transition1x | QM9x | 0.12(0) | 0.07(1) | 0.07(0) | 0.04(0)
QM9x |  | 0.04(2) | 0.02(1) | 0.01(0) | 0.01(0)
Table 3:  Test results of PaiNN models trained on ANI1x, QM9x, Transition1x.
We report RMSE and MAE on energy and forces. Force error is the Euclidian
distance between the predicted and true force vector.

##  Appendix B Additional Figures

This section contains the unbounded version of Fig. 2 as well as additional
plots of MEPs and PESs comparing PaiNN trained on Transition1x with DFT.

Figure 4: Histogram of barrier errors. The x-axis shows errors between
reaction barriers calculated using Density Functional Theory (DFT) and
surrogate potentials for Nudged Elastic Band (NEB). The y-axis shows the
frequency of each bin. Green, red and blue display results from PaiNN models
trained on Transition1x, QM9x and ANI1x, respectively. Yellow displays results
from Density Functional based Tight Binding (DFTB). The QM9x model has such a
low convergence frequency, and general barrier error, that the model does not
show in the plot.

Figure 5: Reaction involving C5OH8. The reaction can be seen as a GIF by
following this link.

Figure 6: Reaction involving C3NCOH7. The reaction can be seen as a GIF by
following this link.

Figure 7: Reaction involving C3NCNH8. The reaction can be seen as a GIF by
following this link.

Figure 8: Reaction involving C3NC2OH9. The reaction can be seen as a GIF by
following this link.