# Force Training Neural Network Potential Energy Surface Models

Christian Devereux  Yoona Yang  Carles Martí  Judit Zádor  Combustion Research
Facility, Sandia National Laboratories, Livermore, CA 94551, USA  Michael S.
Eldred  Sandia National Laboratories, Albuquerque, NM 87185, USA  Habib N.
Najm<EMAIL_ADDRESS>Sandia National Laboratories, Livermore, CA 94551, USA

###### Abstract

Machine learned chemical potentials have shown great promise as alternatives
to conventional computational chemistry methods to represent the potential
energy of a given atomic or molecular system as a function of its geometry.
However, such potentials are only as good as the data they are trained on, and
building a comprehensive training set can be a costly process. Therefore, it
is important to extract as much information from training data as possible
without further increasing the computational cost. One way to accomplish this
is by training on molecular forces in addition to energies. This allows for
three additional labels per atom within the molecule. Here we develop a neural
network potential energy surface for studying a hydrogen transfer reaction
between two conformers of ${\mathrm{C}}_{5}{\mathrm{H}}_{5}$. We show that,
for a much smaller training set, force training can greatly improve the
accuracy of the model compared to only training to energies. We also
demonstrate the importance of choosing the proper force to energy weight ratio
for the loss function to minimize the model test error.

##  1 Introduction

Neural networks 1 provide a functional representation with significant
expressive power 2, 3, 4, 5. Despite concerns about overfitting and lack of
reproducibility 6, the well-documented successes of neural networks (NNs) are
a significant spur for exploration of their utility in the natural sciences,
where unique challenges are present 7, 8, 9, 10, 11, 12, 13, 14, including in
particular utilization for potential energy surface (PES) representation in
chemistry 15, 16, 17, 18, 19, 20, 21, 22. There is, of course, a long
tradition of PESs using functional representations relying on specific
structure involving choices of basis functions 23, 24, 25, 26, 27, 28, 22.
Despite the many achievements of these methods, recent years have seen a
parade of successes of NNs and machine learning (ML) relative to this
established base, demonstrating superior performance and accuracy 29, 30, 31,
32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 22, 43. NNPESs offer multiple
advantages over classical force fields. The functions used to describe
chemical interactions within classical models can limit the types of molecules
that the model can describe well, and are rarely able to model chemical
reactions that involve breaking or formation of covalent bonds. NNPESs offer
much more flexible forms and are limited primarily by the amount and quality
of the data used to train them. Further, compared to quantum mechanical
methods, NNPESs can offer speed ups of several orders of magnitude.

A key ingredient of a successful NNPES representation is its satisfaction of
various necessary symmetries/invariances that are exhibited in molecular
structure. Moreover, a highly useful/desirable property is _extensibility_. In
other words, a PES representation built from some set of structures ought to
be useful for accurate description of the PES of other structures not in the
training set, but involving the same set of atoms, charge, and multiplicity.
One key element of a PES functional construction that governs its utility is
the set of variables used as a feature vector to represent molecular
structure. A ground state PES is in principle a function of internal
coordinates of all the atoms in the molecule and its charge, and to be sure,
many NN based studies rely directly on internal coordinates 44, 45, 46, 33,
47, 48. Other work also used internal coordinates with Gaussian processes
(GPs), rather than NNs 49. However, such representations lack extensibility,
as the PES representation for a different size molecule would have a different
dimension for its input space. Moreover, the dimensionality of this input
specification grows very fast with molecule size, and is not optimal, as
internal constraints on the coordinates imply the presence of lower
dimensional structure.

Lower dimensional representations of molecular structure are often referred to
as collective variables (CVs). The choice of symmetry-preserving, fixed-size,
extensible CVs, in which a low dimensional PES representation can be built in
a wide class of molecules, is challenging 50. There are available means to
discover such good coordinates for any given molecule using purely data-
centric methods 51, 52, 53. However, the enforcement of symmetries, fixed
feature size, and extensibility is yet lacking in these methods. There are
also recent constructions using local sub-networks, preserving symmetries 39,
54, however the extensibility of these methods beyond a limited database of
molecules has yet to be demonstrated.

One alternative to the exclusive use of data to construct CVs is to use
analytical CV formulations constructed based on chemical intuition and
experience, satisfying sought-after properties of fixed-size and symmetry
preservation, and to calibrate them based on data. There have been many
formulations along these lines, e.g. Coulomb matrix 55, 56, “bag-of-bonds” 57,
SMILES strings 58, 59, and others 60, 61, 62, 63, 16, 64, 65, 66, 67. In
particular, the symmetry functions (SFs) of Behler and Parrinello 68, have
found continued development 69, 70, 71, 72, and extensive successful
demonstrations in NN representations of both potential and free energy
surfaces 73, 74, 75, 76, 77, 78, 79, 31, 32, 80, 81. Notably, they have been
successfully demonstrated in gas phase systems, showing good accuracy in the
representation of quantum mechanical reaction probabilities for multiple
reactions 32. Behler & Parrinello SFs, with extensions including atomic-number
differentiated atomic environment vectors (AEVs), have been found to provide
extensibility in a NN setting for large databases of organic molecules 80. In
this work we employ these AEVs as atomic representations for our model.

For a potential energy surface to be useful for applications such as
optimization and dynamics it needs to be able to predict not only accurate
single point energies but also atomic forces. This can be achieved by training
solely on single point energies and using enough training data that the model
has good accuracy on force predictions. An alternate approach involves
training on atomic forces and integrating to get the energies 82. It has in
fact been demonstrated that, for an equal number of training points sampled
from near minima structures, the inclusion of atomic forces during the
training process can yield higher accuracy in both force and energy
predictions as compared to energy-only training 83. The inclusion of force
data improves model predictions away from the training data by ensuring good
representation of the energy gradient at the training data points. Christensen
and von Lilienfeld 83 showed that, when training to a variety of molecular
sizes and composition, the effects of force training on model performance for
potential energy prediction are diminished. However, these models still show a
significant improvement on the prediction of atomic forces when compared to
their non-force trained counterparts 83.

In the present study, we demonstrate the improvements brought about by force
training not only near minima structures but also along a reaction pathway
between two isomers. Accurate forces near transition states are of vital
importance for applications involving reactions as they are necessary to find
the correct transition state as well as the pathway connecting it to the
correct minima. Training on atomic forces has seen use in both materials 70,
84, 85, 86 and chemistry applications 87, 88, 89, 90, 83, 91, 92, 93, 94, 95,
96, 97, 98. Studies include training not just on forces and energies but also
on other properties such as dipoles 99. Generally, such multimodal training
can be advantageous when the various modalities/properties are intrinsically
linked (such as forces to energies through the first derivative), but can hurt
the accuracy of the model on specific quantities of interest (QoIs) if it is
forced to balance the importance of unrelated QoIs. For QoIs that depend on
the Hessian matrix, such as infrared spectra, force training can be useful as
QoIs are again linked through the derivative 95, 100, 41. In this work we
focus only on energies and forces of molecules.

Training on both energies and forces can result in greater memory use per data
point/structure, as well as increased training costs due the need for the
computation of the second derivative of the molecular energies during the
backpropagation process, but, depending on NN construction details, it need
not affect the predictive speed of the trained model. Further, the acquisition
of forces from many DFT calculations is not much more computationally costly
than the acquisition of energies. As might be expected, training to atomic
forces seems to be most impactful when a smaller number of training points are
available 90.

We explore the utility of joint energy+force training on the NNPES
representation of a portion of the PES of C5H5 spanning two wells connected by
a saddle point. The energy difference between the lowest of these wells and
the transition state ($\sim$61 kcal/mol above the lower more stable well) is
beyond the scope of most classical molecular dynamic simulations and cannot be
modeled properly by most classical force fields due to the breaking and
formation of a covalent bond, making it an excellent showcase of the utility
of NNPESs.

We present the relevant details of the construction, including the AEV
formulation, the NN construction, and the associated gradients/Jacobians
necessary for force-training. Using a computational database of C5H5
structures, which we generate using random mode sampling in the vicinity of
the three stationary points and along the intrinsic reaction coordinate
connecting them, we train the NNPES and examine its accuracy. In particular,
we examine the performance of energy+force versus energy-only training for
different training data sizes, and explore the effects of different loss-
function force-to-energy weight ratios during training to identify an optimal
balance between the two. This ratio was found after normalizing the
contribution of both the energy and forces in the loss function to prevent the
model from favoring one over the other due simply to the difference in
magnitude of the two modalities being learned.

We find that the inclusion of force training significantly reduces the test
error of the model on both energies and forces. For the largest data set, with
2048 training points, we find an order of magnitude difference between the L2
error of the energy only and energy+force trained models for both forces and
energies. Comparing the error of the energy+force trained and energy trained
models on molecules generated from normal mode sampling around the three
stationary points, we find that the force trained model shows lower error on
the majority of structures. We show how this increased accuracy on forces and
energy translates to increased accuracy in frequency and zero point energy
calculations, and that the force trained model predicts zero point energies
within chemical accuracy of the level of theory it was trained to.

In the following, we begin by outlining the problem formulation, including the
AEV design, energy and force prediction, as well as NN and loss function
construction. We then present highlights of the software implementation and
the training data generation. Finally, we proceed to presenting results and
discussing their implications.

##  2 Problem Formulation

We begin first with the definition of the feature vector, followed by the NN
construction.

###  2.1 Design of the Feature Vector

Consider a configuration of $N$ atoms $\\{A_{1},\ldots,A_{N}\\}$. Let the
Cartesian coordinates of $A_{i}$ be $x_{i}\in\mathbb{R}^{3}$, and let
$\bm{x}=(x_{1},\ldots,x_{N})\in\mathbb{R}^{3N}$ be the vector of coordinates
of all $N$ atoms. Further, let the set of $n$ atom types in the system be
$T=\\{T_{1},\ldots,T_{n}\\}$, where the type of $A_{i}$ is defined as the
corresponding chemical element $\mathcal{T}_{i}:=\mathcal{T}(A_{i})$. For
convenience, we define the configuration as the $N$-tuple
$\mathcal{C}=(\mathcal{T}_{1},\ldots,\mathcal{T}_{N})$. Thus, _e.g._ , the
configuration $\mathcal{C}=({\mathrm{H}},{\mathrm{H}},{\mathrm{O}})$, involves
the set of atom types $T=\\{{\mathrm{O}},{\mathrm{H}}\\}$. Next, let the index
set of atoms of type $\tau$ in $\mathcal{C}$ be $S_{\tau}$, where

$S_{\tau}=\\{i\ |\ \mathcal{T}_{i}=\tau,\
i=1,\ldots,N\\},\quad\tau=T_{1},\ldots,T_{n}.$ (1)

Similarly, we define the index set of pair-wise atom indices, as

$\displaystyle S_{\tau,\kappa}$ $\displaystyle=$ $\displaystyle\\{(j,k)\ |\
\mathcal{T}_{j}=\tau,\ \mathcal{T}_{k}=\kappa,\ j=1,\ldots,N-1,\
k=j+1,\ldots,N\\},$ $\displaystyle\tau=T_{1},\ldots,T_{n};\
\kappa=\tau,\ldots,T_{n}.$

For example, for a system with
$T=\\{{\mathrm{C}},{\mathrm{H}},{\mathrm{O}}\\}$, and the two configurations
$({\mathrm{H}},{\mathrm{H}},{\mathrm{O}})$ and
$({\mathrm{C}},{\mathrm{O}},{\mathrm{O}})$, we have the index sets shown in
Table 1. Note that, with $n$ atom types, there are $n$ single-atom index sets
$S_{\tau}$, and $m=n(n+1)/2$ pair-wise index sets $S_{\tau,\kappa}$.

$\mathcal{C}$ | $S_{\mathrm{C}}$ | $S_{{\mathrm{H}}}$ | $S_{\mathrm{O}}$ | $S_{{\mathrm{C}},{\mathrm{C}}}$ | $S_{{\mathrm{C}},{\mathrm{H}}}$ | $S_{{\mathrm{C}},{\mathrm{O}}}$ | $S_{{\mathrm{H}},{\mathrm{H}}}$ | $S_{{\mathrm{H}},{\mathrm{O}}}$ | $S_{{\mathrm{O}},{\mathrm{O}}}$
---|---|---|---|---|---|---|---|---|---
(H,H,O) | {} | {1,2} | {3} | {} | {} | {} | {(1,2)} | {(1,3),(2,3)} | {}
(C,O,O) | {1} | {} | {2,3} | {} | {} | {(1,2),(1,3)} | {} | {} | {(2,3)}
Table 1: Sample configurations and associated index sets with
$T=\\{{\mathrm{C}},{\mathrm{H}},{\mathrm{O}}\\}$.

Following Smith _et al._ 80, we define the contribution to the potential
energy of the system due to $A_{i}$ as $E_{A_{i}}(\mathcal{C},\bm{x})$, such
that the total potential energy is

$E=\sum_{i=1}^{N}E_{A_{i}}(\mathcal{C},\bm{x}).$ (3)

We use the “atomic environment vector” (AEV) 80 as the feature vector that
summarizes the geometry of the system in the neighborhood of each atom. We
write the AEV of $A_{i}$ as $y_{i}(\mathcal{C},\bm{x})\in\mathbb{R}^{M}$, and
use one NN,
$\mathcal{N}_{\mathcal{\tau}}(y):\mathbb{R}^{M}\rightarrow\mathbb{R}$, for
each atom type $\tau$, such that

$E=\sum_{i=1}^{N}\mathcal{N}_{\mathcal{T}_{i}}(y_{i}(\mathcal{C},\bm{x})).$
(4)

The AEV involves both radial and angular geometry components following Behler
and Parinello 68. For $A_{i}$, we write the radial components of $y_{i}$ for
pairings with each available atom type separately, according to the $S_{\tau}$
index sets. Similarly, we write angular components of $y_{i}$, for groupings
with pairs of atom types separately according to the index sets
$S_{\tau,\kappa}$. This is done in each case by summing contributions of the
pairings/groupings within each set. The construction relies on summarizing
geometry information localized around each atom. A key element of the
localization involves the cutoff function

$f_{c}(R_{ij},R_{c})=\begin{cases}0.5\cos\Big{(}\frac{\pi
R_{ij}}{R_{c}}\Big{)}+0.5\quad&\mathrm{for}\ R_{ij}\leq R_{c}\\\
0&\mathrm{otherwise}\end{cases}$ (5)

where $R_{ij}=\|x_{i}-x_{j}\|$ is the Euclidean distance between
$(A_{i},A_{j})$, and $R_{c}$ is a cutoff radius. Then with
$\mu:=(\eta,\rho)\in\mathbb{R}^{2}$, we define the radial components of
$y_{i}$ for pairings with atom type $\tau$ as

$y^{\tau,\mu}_{i}=\sum_{j\in S_{\tau},j\neq
i}e^{-\eta(R_{ij}-\rho)^{2}}\,f_{c}(R_{ij},R_{c}^{r})$ (6)

with $\mu\in\\{\mu_{1},\ldots,\mu_{M_{r}}\\}$, and where $R_{c}^{r}$ is a
radial-SF cutoff radius. Further, with
$\nu:=(\xi,\gamma,\zeta,\alpha)\in\mathbb{R}^{4}$, we define the angular
components of $y_{i}$, for pairings with atom types $(\tau,\kappa)$, as

$y^{\tau,\kappa,\nu}_{i}=\sum_{(j,k)\in S_{\tau,\kappa},\,j\neq i,k\neq
i}(0.5+0.5\cos(\theta_{ijk}-\alpha))^{\zeta}\
e^{-\xi(0.5(R_{ij}+R_{ik})-\gamma)^{2}}\
f_{c}(R_{ij},R_{c}^{a})\,f_{c}(R_{ik},R_{c}^{a})$ (7)

where $\nu\in\\{\nu_{1},\ldots,\nu_{M_{a}}\\}$, $R_{c}^{a}$ is an angular-SF
cutoff radius, and $\theta_{ijk}\in[0,\pi]$ is the angle, centered on $A_{i}$,
between the two vectors $x_{ij}:=x_{j}-x_{i}$ and $x_{ik}:=x_{k}-x_{i}$, given
by

$\theta_{ijk}:=\arccos\frac{x_{ij}\cdot x_{ik}}{\|x_{ij}\|\|x_{ik}\|}.$ (8)

The full length of the AEV is $M=M_{r}\,n+M_{a}\,m$. With

$\displaystyle y_{i}^{\tau}$ $\displaystyle:=$
$\displaystyle(y_{i}^{\tau,\mu_{1}},\ldots,y_{i}^{\tau,\mu_{M_{r}}})\in\mathbb{R}^{M_{r}}$
(9) $\displaystyle y_{i}^{\tau,\kappa}$ $\displaystyle:=$
$\displaystyle(y_{i}^{\tau,\kappa,\nu_{1}},\ldots,y_{i}^{\tau,\kappa,\nu_{M_{a}}})\in\mathbb{R}^{M_{a}}$
(10)

we write the full AEV as

$y_{i}=(y^{\tau_{1}}_{i},\ldots,y^{\tau_{n}}_{i},y^{(\tau,\kappa)_{1}}_{i},\ldots,y^{(\tau,\kappa)_{m}}_{i})\in\mathbb{R}^{M}.$
(11)

Thus, _e.g._ , for $A_{i}$ in a system with atom types
$\\{{\mathrm{C}},{\mathrm{H}}\\}$,
$y_{i}=(y^{{\mathrm{C}}}_{i},y^{{\mathrm{H}}}_{i},y^{{\mathrm{C}},{\mathrm{C}}}_{i},y^{{\mathrm{C}},{\mathrm{H}}}_{i},y^{{\mathrm{H}},{\mathrm{H}}}_{i})\in\mathbb{R}^{M}$,
with $M=2M_{r}+3M_{a}$. Corresponding representative radial and angular
symmetry functions, whose superposition provides the above radial/angular AEV
components, are shown in Figs. 1, 2. Figure 1 illustrates how the radial
symmetry functions resolve the local environment around $A_{i}$ in a number of
radial shells, with decaying relevance at larger radii. Similarly Figure 2
shows the resolution of the angular environment around $A_{i}$ in terms of
angles formed at $A_{i}$ with a second and a third atom. Each frame
corresponds to a different location of the second atom, illustrating the
multiple radial/angular shells summarizing the localization of the third atom.

Figure 1: Radial symmetry functions around a reference atom $A_{i}$.

Figure 2: Angular symmetry functions around a reference atom $A_{i}$ at
$(0,0)$.

###  2.2 AEV Jacobian

The AEV provides the NN inputs for potential energy representation. This
construction is also usable for representation of forces, being the negative
spatial gradient of the energy. We present here the extension of the AEV
construction for force computations.

We have so far expressed the 3D coordinates of $A_{j}$ as $x_{j}$. Here we
expressly write the formulation in terms of the 3 spatial coordinates, where
we define $x_{j}:=(x_{j1},x_{j2},x_{j3})$. The PES for a configuration with
$N$ atoms is the function $E(x_{1},\ldots,x_{N})$, being a function on a
$3N$-dimensional space, although rotational and translational potential energy
invariances imply that the true dependence is generally in a $3N-6$
dimensional space. DFT computations can typically provide both the energy and
its gradients with respect to components of $\bm{x}=(x_{1},\ldots,x_{N})$,
namely $\nabla_{x_{j}}E$, for $j=1,\ldots,N$. We write the gradient operator
as

$\nabla_{x_{j}}=\Big{(}\frac{\partial}{\partial
x_{j1}},\frac{\partial}{\partial x_{j2}},\frac{\partial}{\partial
x_{j3}}\Big{)},$ (12)

such that,

$\nabla_{x_{j}}E=\sum_{i=1}^{N}\nabla_{x_{j}}\mathcal{N}_{\tau_{i}}(y_{i}(\mathcal{C},\bm{x})),$
(13)

where, recalling that $y_{i}=(y_{i1},\ldots,y_{iM})$, we have

$\nabla_{x_{j}}\mathcal{N}_{\tau_{i}}(y_{i}(\mathcal{C},\bm{x}))=\sum_{r=1}^{M}\frac{\partial\mathcal{N}_{\tau_{i}}}{\partial
y_{ir}}\nabla_{x_{j}}y_{ir}.$ (14)

Both the NN and AEV gradients can be computed directly from their analytical
constructions. Notably, the AEV Jacobian

$J^{ij}_{rk}=\\{\partial y_{ir}/\partial x_{jk}\\}$ (15)

where $J^{ij}_{rk}$ is the derivative of the $A_{i}$ AEV component $r$ with
respect to the $A_{j}$ Cartesian coordinate $k$, is a key element of the
construction for force computations. The detailed formulation of the
analytical AEV Jacobian is omitted here, but is available elsewhere 101.

Recall that the AEV, given its particular atom-centered radial/angular
construction, is invariant to solid body translation or rotation of the
configuration. Further, the additive potential energy construction renders the
resultant energy prediction also invariant to permutations of same-type atoms.
The enforcement of these physical constraints on the energy representation by
construction is an important feature of the AEV-NN representation, ensuring
that NN training is constrained in a physical manifold. While the same
physical constraints hold for the force vectors in terms of their magnitudes
and relative orientations, this is not true for their Cartesian components,
specifically as concerns rotational and permutational invariance. This
dependence of force Cartesian components on rotation and permutation is indeed
captured in the AEV-NNPES representation. The NN-predicted force components
are computed as derivatives of the energy function as in Eq. 13, and thus of
the AEV as in Eq. 14, with respect to Cartesian coordinates. The NNPES
construction captures the rotational and permutational dependencies of the
force Cartesian components via the AEV Jacobian.

Figure 3: NNPES schematic for a molecule/configuration involving atom types
$({\mathrm{C}},{\mathrm{H}})$.

###  2.3 Neural Network Construction

Consider a fully connected feed forward NN with $L$ layers. Let the input be
the vector $\bm{u}_{0}\in\mathbb{R}^{n_{0}}$, and the output of layer
$\ell=1,\ldots,L$ be $\bm{u}_{\ell}\in\mathbb{R}^{n_{\ell}}$. With the
activation function defined as $f_{\ell}():\mathbb{R}\rightarrow\mathbb{R}$,
we write the output of each layer as
$\bm{u}_{\ell}=(u_{\ell,1},\ldots,u_{\ell,n_{\ell}})$, given by

$\bm{u}_{\ell}=\bm{f}_{\ell}(\bm{W}_{\ell}\bm{u}_{\ell-1}+\bm{b}_{\ell})$ (16)

where $\bm{f}_{\ell}(\bm{q}):=(f_{\ell}(q_{1}),\ldots,f_{\ell}(q_{n_{\ell}}))$
for $\bm{q}\in\mathbb{R}^{n_{\ell}}$,
$\bm{W}_{\ell}=\\{w_{\ell,ij}\\}\in\mathbb{R}^{n_{\ell}\times n_{\ell-1}}$ and
$\bm{b}_{\ell}\in\mathbb{R}^{n_{\ell}}$. Then, for any batch of $K$ input
vectors, we define the data input matrix
$\bm{U}_{0}=[\bm{u}_{0,1}\cdots\bm{u}_{0,K}]$ such that

$\bm{U}_{\ell}=\bm{F}_{\ell}(\bm{W}_{\ell}\bm{U}_{\ell-1}+\bm{B}_{\ell})$ (17)

so that
$\bm{B}_{\ell}=[\bm{b}_{\ell},\ldots,\bm{b}_{\ell}]\in\mathbb{R}^{n_{\ell}\times
K}$, and, for matrix $\bm{Q}\in\mathbb{R}^{n_{\ell}\times K}$, with column
vectors $\bm{q}_{1},\ldots,\bm{q}_{K}$, we have
$\bm{F}_{\ell}(\bm{Q}):=[\bm{f}_{\ell}(\bm{q}_{1})\cdots\bm{f}_{\ell}(\bm{q}_{K})].$
For our purposes, the output node (last layer $L$) is a scalar, thus
$n_{L}=1$, and we employ $f_{L}(u)\equiv u$.

As indicated above, we use a separate NN for each atom type, employing the
same above NN structure for each. With the AEV $y_{i}$, for $A_{i}$, provided
as input to the NN $\mathcal{N}_{\mathcal{T}_{i}}$ this construction provides
the corresponding NN output as $E_{A_{i}}$, to be used in Eq. 3. Considering
our present case with two atom types, Figure 3 illustrates the overall
computation of the potential energy for a molecule composed of
$\\{{\mathrm{C}},{\mathrm{H}}\\}$ atoms using the NNPES construction. Again,
the overall detailed formulation is presented elsewhere 101.

###  2.4 Loss Function

Given a dataset of $K$ configurations/structures where energies
$E^{D}=\\{E_{1},\ldots,E_{K}\\}$ and forces
$\bm{F}^{D}=\\{\bm{F}_{1},\ldots,\bm{F}_{K}\\}$ are known from quantum
chemistry computations, we quantify the accuracy of the NNPES using a loss
function that computes the norm of the error in NNPES predicted quantities
$(E,\bm{F})$ compared to their values in the database $(E^{D},\bm{F}^{D})$.
For configuration $k$ with $N_{k}$ atoms, energy $E_{k}$, and force vector
$\bm{F}_{k}:=-(\nabla_{x_{1}}E_{k},\ldots,\nabla_{x_{N_{k}}}E)=(F_{k,1},\ldots,F_{k,3N_{k}})\in\mathbb{R}^{3N_{k}}$,
we define the (mean for forces) square (L2) errors as

$\displaystyle\mathcal{E}^{E}_{k}$ $\displaystyle=$
$\displaystyle|E_{k}-E^{D}_{k}|^{2}$ (18) $\displaystyle\mathcal{E}^{F}_{k}$
$\displaystyle=$ $\displaystyle\|\bm{F}_{k}-\bm{F}^{D}_{k}\|^{2}/(3N_{k}).$
(19)

Further, implementing a weight $\lambda\in[0,1]$, and the normalizing scale
factors

$\displaystyle d_{E}$ $\displaystyle=$
$\displaystyle\max_{k}(E^{D}_{k})-\min_{k}(E^{D}_{k})$ (20) $\displaystyle
d_{F}$ $\displaystyle=$
$\displaystyle\max_{k,r}(F^{D}_{k,r})-\min_{k,r}(F^{D}_{k,r}),$ (21)

we define the aggregate mean-square loss as

$\mathcal{E}=\frac{(1-\lambda)}{Kd^{2}_{E}}\sum_{k=1}^{K}\mathcal{E}^{E}_{k}+\frac{\lambda}{Kd^{2}_{F}}\sum_{k=1}^{K}\mathcal{E}^{F}_{k}$
(22)

where division by $(d_{E},d_{F})$ provides normalization of the energy and
force magnitudes on $[0,1]$, and the weight $\lambda$ provides relative
weighting of energy and force errors. We define the ratio
$\lambda/(1-\lambda)$ as the force to energy ratio (F:E) and discuss its
effect on training and test errors in the results section.

##  3 Software Implementation

We used Pytorch 102 as a machine learning framework, in a Python context.
Given the costs associated with AEV and Jacobian computations, we implement
these (using the analytical Jacobian) in C++, relying on Pybind11 103 for
Python-C++ binding. We published this AEV code as an open source library
aevmod 104, 101. Although not directly relevant here, the use of the resulting
NNPES for geometry optimization requires the Hessian of the energy function,
and thus of the AEV. We implemented this Hessian similarly in C++ in aevmod,
relying on automatic differentiation using Sacado 105.

##  4 Data

For demonstration, we focus on DFT computations on the C5H5 PES using QChem
106. Energies and forces are evaluated at the $\omega$B97X-D/6-311++G(d,p)
level of theory for a portion of the 24-dimensional potential energy surface.
We consider a two-well system, where the two wells are connected via a saddle
point. The system is comprised of the lowest energy species on the C5H5 PES
(cyclopentadienyl), the saddle point corresponding to the reaction with the
lowest energy barrier cyclopentadienyl can undergo (1,2-H atom shift), and the
product resulting from it. The two wells and the associated reaction
(highlighted in red) are shown in Fig. 4.

Figure 4: PES of the reaction being studied (red) and other surrounding
reactions (black).

We train an ensemble of NNs using subsets of the DFT database. The training
and validation sets for each ensemble member were randomly sampled from the
same pool of structures containing a total of 10875 atomic configurations
(25$\times$435). The pool of structures used for training, validation, and
testing was generated from 25 anchor points, being the three stationary points
of the system and 22 points along the intrinsic reaction coordinate (IRC)
connecting them. From each of these 25 points, 434 additional structures were
generated using randomized normal mode sampling corresponding to $T=2000$ K.
This process involved stretching and compressing the structure at each anchor
point along its 24 vibrational modes. We modeled our normal mode sampling
based on the sampling method outlined in Smith _et al._ 80. For each mode $i$,
a uniformly distributed pseudo-random number ${\tilde{c}}_{i}$ is sampled on
the range $[0,1]$. The ${\tilde{c}}_{i}$ are then normalized by their sum and
scaled by another uniformly distributed random number $D$ on $[0,1]$, such
that

$c_{i}=\frac{D{\tilde{c}}_{i}}{\sum_{j=1}^{3N_{a}-6}{\tilde{c}}_{j}},\qquad
i=1,\ldots,3N_{a}-6.$ (23)

This scaling ensures that $\sum_{i=1}^{3N_{a}-6}c_{i}\equiv D\in[0,1]$, which
in turn ensures that structures will be randomly sampled at a range of
temperatures up to a given maximum corresponding to $D=1$. In this sense, $D$
serves as a measure of the “normalized temperature distance” from the anchor
point to the structure, where, for $D=0$ we have a structure at the anchor
point, and for $D=1$ we have a structure at the maximum temperature $T$ above
the anchor point.

However, we made an empirical correction to the normal mode sampling procedure
for low-frequency modes. The physical reason for this is that the small force
constants of these modes drive the sampling into regions well beyond the
validity of a harmonic approximation, yielding extremely high energy points or
leading to convergence failures in QChem due to clashing atoms. We designate a
frequency $f_{i}$ low if $f_{i}<150$ cm-1 and replace $\tilde{c}_{i}$ with
$\tilde{c}_{i,\mathrm{mod}}$ as

$\tilde{c}_{i,\mathrm{mod}}=\tilde{c}_{i}(f_{i}/150)^{4}.$ (24)

The displacement along mode $i$, $R_{i}$, is then given by

$R_{i}=\pm\sqrt{\frac{3c_{i}N_{a}k_{b}T}{K_{i}}},$ (25)

where $N_{a}$ is the number of atoms, $k_{b}$ is Boltzmann’s constant, and
$K_{i}$ is the force constant of the mode. The sign of $R_{i}$ is chosen
randomly from a Bernouli distribution with $p=0.5$ to ensure sampling of both
sides of the harmonic potential. The molecule is then perturbed along that
mode by its normalized normal mode coordinate scaled by $R_{i}$.

##  5 Results and discussion

For each case, we trained an ensemble of neural networks, each with different
random initial parameters and training/validation sets. All models were
trained with a 256-long AEV with the specific construction parameters outlined
in the supplementary information. For each atomic network, we used a fully
connected NN architecture of 256:128:64:64:1 with a Gaussian activation
function $f_{\ell}(q)=\exp(-q^{2})$ for each of the first three layers and a
liner function for the final layer. A test set of 59 structures was used to
assess the effects of force training and its accuracy compared to energy only
trained models. The test data is held apart, and thus not used in the NN
training/validation process. The same test data set is used in all cases. The
models trained only to energies are henceforth referred to as the energy
trained models and the models trained to both energies and forces as the force
trained models, for the sake of brevity.

###  5.1 Neural Network Training

For purposes of neural network training, we used a batch size of 8, and made
use of a randomized validation set of 64 structures for controlling the
learning rate and deciding when to stop training. Each network was trained
until its root mean square error (RMSE) on energy for the validation set
increased, checking every 1000 epochs. The ADAM algorithm 107 as implemented
in Pytorch 102, 108, was used as the optimizer. The initial learning rate was
set to $10^{-3}$. To determine when to update the learning rate, we compared
the energy RMSE for the validation and training sets averaged over 100 epochs
to those evaluated similarly 1000 epochs earlier. If the change in the average
validation RMSE increased while the average training RMSE decreased we
multiplied the current learning rate by a factor of 0.1 and continued. Once
the learning rate reached our chosen cutoff of $10^{-6}$, this process
continued until either the criterion to lower the learning rate was met again
or the model had been trained for 50000 epochs, at which point training was
stopped. Figures 5-6 show the evolution of the training/validation/test errors
of a single force-trained and energy-trained model with training epoch, each
trained to 512 structures. For the force trained model, the energy and forces
were given equal weights ($\lambda=0.5$ from equation 22) in the loss
function. For each model the energy error rapidly decreased in the early
epochs but fluctuated greatly before converging. A similar pattern can be seen
for the force error of the force trained model. The force error for the energy
trained model begins to decrease after the initial few hundred epochs then
fluctuates until converging at a similar level it started at. For both the
force and energy errors, the validation error for the model trained only to
energy reaches a minimum before increasing in subsequent epochs due to
overfitting, given the small size of the data. This does not happen for the
force trained models, consistent with the larger amount of information
available given the force data. For the subsequent analyses below, unless
otherwise stated, the models used to generate the results are those whose
parameters were saved at the point of lowest testing error during their
training.

Figure 5: Plot shows the evolution of the RMSE energy error (left panel) and
relative force error (right panel) in the training, testing, and validation
set for the energy trained model during training.

Figure 6: Plot shows the evolution of the RMSE energy error (left frame) and
relative force error (right frame) in the training, testing, and validation
set for the force trained model during training.

###  5.2 Weighting Energy and Force Errors in the Loss Function

Finding the optimal ratio of force to energy weight in the loss function is
important when training to forces, as the model needs to balance the
importance of the two components of the loss function. There are important
factors to consider in this regard, such as the difference in scale between
the two properties and the amount of data present for each. We discussed
earlier above the normalization of the loss for each task in the loss
function. This helps to balance the importance of the two quantities even if
the absolute scale of the values is different. To properly tune the weighting
between the energies and forces for our data set we tested a range of ratios
from a $10^{-2}$:1 F:E ratio to a $10^{5}$:1 F:E ratio using 512 structures in
the training set. These results are shown in Figure 7. We used the root mean
square (RMS) norm as a measure of the error of the energy predictions and the
relative RMS norm as a measure of the error in force prediction. We define the
relative RMS error in the forces as:

$\mathcal{E}_{\mathrm{F}}=\sqrt{\frac{\sum_{k}\|\bm{F}_{k}-\bm{F}^{D}_{k}\|^{2}}{\sum_{k}\|\bm{F}^{D}_{k}\|^{2}}}$
(26)

For each F:E ratio case we trained an ensemble of 3 NNs to provide a sense of
the scatter in the results. Both initial weights and training/validation data
were randomly sampled for each training ensemble member. We see that going to
low F:E ratios does not exhibit a reliable trend in the training error in
energy (TrE), with a possible increase in TrE as the F:E ratio is increased to
large values, if any trend is to be believed given the scatter. On the other
hand, the training error in force (TrF) exhibits a clear upward trend with
reduced F:E weight ratio, again without a clear trend in the high F:E limit.
Broadly speaking, this is expected, as in the low F:E limit, the model is
essentially trained only on energies, thus not providing much of a constraint
on the force error. For high F:E weight ratio, the model is trained largely on
the forces, with little constraint on the energies, so a possible rise in TrE
is conceivable. Clearly, the results indicate that the rise in TrF for low F:E
ratio is much more significant than any corresponding trend in TrE for high
F:E ratio. Looking at the test error is, we see that the test errors on both
energy and force increase significantly in the limit of very low F:E ratio.
This is a clear indication of poor generalizability in this limit,
highlighting the critical need for force training to ensure good test error,
i.e. good accuracy away from the training data. As for test error in the limit
of high F:E ratio, the trend, if any is to be inferred is slightly up in both
errors. Overall, we may say that a ratio of 100:1 F:E in the loss function
roughly yielded the lowest error on both the energies and forces on the test
set. For the training error no one ratio performed reliably best on both the
forces and energies for the range of ratios we used. Interestingly, even for
the ratios that most heavily favored forces the error in the energies was
always lower than the error in the forces.

Figure 7: Training and test errors in energy and forces for energy+force
trained models as a function of the F:E ratio. An ensemble of 3 trained models
is considered at each F:E ratio condition.

Generally, it is fair to say that the optimal weight ratio to use when
training for multiple tasks can vary greatly depending on the dataset used and
model architecture. A thorough search across many ratios can lead to lower
error and therefore a more useful model. Knowledge of the system being studied
and the properties of interest can help to narrow the range of values to
search but it is still important to perform a search as the factors
determining the optimal ratio vary from system to system. A model trained only
to near equilibrium points likely has a smaller range of forces and energies
in its training set than one trained to points along an IRC and this may
affect the ratio that should be used.

###  5.3 Dependence of Test Errors on Training Data Size

Figures (8,9) show the effect of force training on test error for a range of
training set sizes from 32 to 2048 in both energy and force predictions using
a F:E ratio of 100:1 in the loss function. We show the test error scatter for
an ensemble of 6 NNs, as well as the test error from the ensemble mean model.
In all cases every ensemble member of the force trained models outperformed
the energy only trained models. Each individual force trained network also
outperformed the ensemble average of the energy only trained models for the
same training set size for both forces and energies. The test error ratio
between the force and energy trained models increased with increasing training
set size. This is true for both the forces and energies. For the largest
training set, the ratio is about an order of magnitude. The difference in
accuracy of the force predictions is even more noticeable. The energy trained
models show little improvement in force prediction with increasing data set
size while the error in the force trained models decreases by almost an order
of magnitude. As the size of the dataset increases, the amount of data being
fitted in the force trained model increases much faster than that for the
energy trained model. In the case of these models, trained to C5H5, there are
31 times as many labels being used to train the force model compared to the
energy trained model (3 forces/atom $\times$ 10 atoms/structure + 1 energy,
per structure) so the difference between the amount of training data between
the models grows as $30N$ where $N$ is the number of structures in the
training set.

Figure 8: Plot shows reduction in the RMS error in energy (kcal/mol) with
increasing training data size. It also highlights the two cases of energy-only
(ETrain) and energy+force (FTrain) training. Results illustrate the lower
energy test error in the FTrain case, by about an order of magnitude for 2048
training points.

Figure 9: Plot shows reduction in the RMS error in force (kcal/mol/A) with
increasing training data size. It also highlights the two cases of energy-only
(ETrain) and energy+force (FTrain) training. Results illustrate the lower
force test error in the FTrain case, by about an order of magnitude for 2048
training points.

Unfortunately, this extra data, as well as the need for the computation of the
second derivative of the energy during the backpropagation step, does lead to
increased computational training costs. However, when comparing the error of
the energy only trained models to the force trained models, it can be seen
that even when training to much larger datasets, the energy-only models
struggle to match the accuracy of the force trained models. For the energy
error the force trained models trained to only 64 structures show lower error
on average than the energy-only models trained to 2048 structures. For the
force error, even when training to the largest dataset tested the energy-only
model shows higher error than the force trained model trained to only the
smallest dataset tested. The fact that training to forces requires a much
smaller dataset to achieve the same accuracy as an energy-only trained model
can help offset the costs associated with training to forces. It is important,
of course, to still ensure that the training set used is diverse enough to
cover the region of chemical space being studied, but force training can help
the model learn the same region of chemical space with fewer data points.
Being able to achieve greater test accuracy with a smaller ensemble, using
force training rather than energy only can also help decrease training times
as fewer networks are required.

###  5.4 Error Dependence on Distance from Anchor Points

In addition to the 59 structures used in the previous test cases, we also
compared energy-only and energy+force trained models on structures generated
from normal mode sampling at 2000 K, using the same method described in the
data section above, centered around the three stationary points of the system
(well minima and the transition state between them). We used an ensemble of 6
neural nets for each model, each with different random initial parameters and
training/validation sets pulled from the same pool of molecules. For this
analysis, the models whose parameters were saved at the point of lowest
validation error lowest validation error were used. These models were trained
to 2048 structures and the force trained model used a 100:1 F:E ratio in its
loss function. For each stationary point, 20 structures were generated for
each vibrational mode for a total of 480 structures per stationary point. Some
of these structures failed to converge during the DFT calculations, such that
a total of 1423 points (out of 1440) were used for this test. We compared the
errors of the energy/force-trained models when predicting energies and forces
of these structures. For the energy we show the absolute difference in the
prediction of the models and DFT. For the forces we show the relative RMS
error as defined above in Eq. 26. The parity plots in Figures 10 and 11 show
the performance of each model on each point for energy and forces
respectively. They also show the normalized temperature distance ($D$ from Eq.
23) of each structure away from its anchor point. Points below the black
diagonal line are structures where the force trained model shows a lower error
than the energy trained model, while points above the diagonal exhibit the
reverse relationship.

Both models tend to show higher energy error on points further from the anchor
points. However, even for these far points, the force trained model performs
better on the majority of structures. The correlation between the model test
error and the distance from the anchor point is more obvious for the force
trained model than the energy trained model in Figure 11 but for both models
the highest error structures are structures that are far from the anchor
points. In Figure 10 there are 1188 structures below the black line and 235
above it. Additionally, the force trained model shows chemical accuracy on 82%
of the points, (1180 out of 1423). The energy trained model achieves this
accuracy on only 46% of the data (650 out of 1423) In Figure 11 there are 849
structures below the black line and 574 above it. The force trained model
shows a relative error as low as $10^{-2}$ on many structures, including the
anchor points, while the energy trained model does not achieve even an error
of less than of $10^{-1}$ on any of the tested structures. This means that
force trained model can be expected to show errors 1% of the magnitude of the
total forces on the anchor points while the energy trained model will be an
order of magnitude worse. Despite showing lower error on most of the points,
the highest errors in force and energy prediction for the force trained model
are worse than for the energy-trained model. However, these outliers represent
only a small percentage of the total number of points. For both the force
trained and energy trained models their worst points show places where the
model failed to learn and the results for each model are qualitatively
incorrect. The three anchor points (shown as stars) are also included in both
graphs as the darkest blue points with a distance of zero, although they were
also included in the training set for each model. The force trained model
outperforms the energy only trained model on each of these three points for
both energy and forces.

###  5.5 Frequency prediction

For a final test of the advantages of force training versus energy training
and to measure the usefulness of each model on a kinetically relevant
property, we computed the vibrational frequencies of each stationary point
using both models and compared the results to $\omega$B97X-D/6-311++G(d,p).
The same models used to generate figures 11 and 10 were used in this
comparison. For each model, we optimized the stationary points, starting from
the geometry predicted by $\omega$B97X-D/6-311++G(d,p), computed the
vibrational frequencies, and the average percent error of the frequencies
compared to $\omega$B97X-D/6-311++G(d,p). These results are shown in Figure
12. The force trained model shows excellent agreement with DFT across the
entire range of frequencies and has an average percent error of 4%. The energy
trained model underestimates most of the frequencies below 1000 cm-1 and
overestimates most of the ones above 1000 cm-1, showing an average percent
error of $30\%$. Both models correctly optimize the transition state to a
transition state, finding only one imaginary frequency, but the energy trained
model misses the value of this frequency by several hundred wavenumbers. The
zero point energies (ZPE) of the force trained model agree within 0.5 kcal/mol
with that of the DFT values, while the ZPE arising from the energy trained
model overestimates the ZPE by 6.9, 5.4 and 13.6 kcal/mol for the two wells
and the saddle point (respectively), causing a large error in the final,
chemically relevant reverse and forward barrier heights, despite the small
error in the electronic energy.

Figure 10: Energy absolute test error for mean 6-NN predictions at normal mode
samples. Plot shows mean prediction test errors resulting from energy+force
training versus those from energy training, with 2048 training points. Points
are colored according to their normalized distance from anchor points ($D$),
where normalization is done separately for each normal mode. The three anchor
points at distance 0 are shown as stars. The percentage of points with less
than 1.0 kcal/mol error from each model are also shown. 1423 points are
plotted in total.

Figure 11: Force relative RMS test error for mean 6-NN predictions at normal
mode samples. Plot shows mean prediction relative test errors resulting from
energy+force training versus those from energy training, with 2048 training
points. Points are colored according to their normalized distance from anchor
points, where normalization is done separately for each normal mode. The three
anchor points at distance 0 are shown as stars.

Figure 12: Frequencies predicted by the force trained(blue) and energy
trained(red) models versus frequencies predicted by
$\omega$B97X-D/6-311++G(d,p) for each stationary point.

##  6 Conclusions

Our results indicate that, with adequate training, a neural network potential
energy surface can provide a level of accuracy over reactive landscapes that
is comparable to QM methods, but at a much lower predictive computational
cost. Of course, the generation of a training set that adequately describes
the system of interest can still be a difficult and costly process. Training
to forces in addition to energy can greatly decrease the amount of training
data needed without sacrificing model accuracy. Training to both energies and
forces increased the accuracy of the model for both energy and force
prediction. A force trained model trained on a small set of structures shows
lower error than a model trained to an order of magnitude more structures
without force data. However, the trained model accuracy does depend on the
relative weights given to the forces versus the energies in the loss function.
The models trained only to energy show slower improvement in force prediction
with increasing training set size as compared to force trained models. Force
training also results in greater accuracy when predicting kinetically relevant
properties than training to energy alone. The force trained model shows
excellent correlation with DFT when predicting vibrational frequencies, and
its zero point energy predictions are more than an order of magnitude lower in
error than the energy trained model.

We believe force training can allow for larger and more complex systems to be
studied using NNPESs, systems which would otherwise be too costly to generate
enough training data for to achieve requisite accuracy from energy-only
training. Further, the cost of running DFT calculations grows as
$\mathcal{O}(N^{3})$, making the creation of large training sets for large
molecules a costly process. Force training allows the same level of accuracy
to be achieved with a much smaller training set so fewer of these costly
calculations are required.

##  7 Acknowledgement

This work was supported by the US Department of Energy (DOE), Office of Basic
Energy Sciences (BES) Division of Chemical Sciences, Geosciences, and
Biosciences, under contract number 022187. This research used resources of the
National Energy Research Scientific Computing Center (NERSC), a DOE Office of
Science User Facility supported by the Office of Science of the U.S. DOE under
Contract No. DE-AC02-05CH11231 using NERSC Award No. BES-ERCAP0021664. This
article has been authored by employees of National Technology & Engineering
Solutions of Sandia, LLC under contract no. DE-NA0003525 with the U.S.
Department of Energy (DOE). The employees co-own right, title, and interest in
and to the article and are responsible for its contents. The United States
Government retains and the publisher, by accepting the article for
publication, acknowledges that the United States Government retains a
nonexclusive, paid-up, irrevocable, worldwide license to publish or reproduce
the published form of this article or allow others to do so, for United States
Government purposes. The DOE will provide public access to these results of
federally sponsored research in accordance with the DOE Public Access Plan
https://www.energy.gov/downloads/doe-public-access-plan.

We specify the various parameters of the {C,H} AEV construction as follows. We
use

  * •

Cutoff radius for radial SF: $R_{c}^{r}=4.6$Å

  * •

Cutoff radius for angular SF: $R_{c}^{a}=3.1$Å

  * •

Discretization of $\eta$ uses $n_{\eta}=1$, $\eta=1.0/\delta_{r}^{2}$, with
$\delta_{r}=2h_{r}/3$, $h_{r}=R^{r}_{c}/n_{\rho}$

  * •

Discretization of $\rho$ uses $n_{\rho}=32$, $\rho_{i}=(i-\frac{1}{2})h_{r}$,
$i=1,\ldots,n_{\rho}$

  * •

We have $M_{r}=n_{\eta}n_{\rho}=32$, where
$(\eta,\rho)\in\mathbb{R}^{n_{\eta}}\times\mathbb{R}^{n_{\rho}}$

  * •

Discretization of $\xi$ uses $n_{\xi}=1$, $\xi=1.0/\delta_{a}^{2}$,
$\delta_{a}=2h_{a}/3$, $h_{a}=R^{a}_{c}/n_{\gamma}$

  * •

Discretization of $\gamma$ uses $n_{\gamma}=8$,
$\gamma_{i}=(i-\frac{1}{2})h_{a}$, $i=1,\ldots,n_{\gamma}$

  * •

Discretization of $\zeta$ uses $n_{\zeta}=1$, $\zeta=8$

  * •

Discretization of $\alpha$ uses $n_{\alpha}=8$, $h=\pi/(n_{\alpha}-1)$,
$\alpha_{i}=(i-1)h$, $i=1,\ldots,n_{\alpha}$

  * •

We have $M_{a}=n_{\xi}n_{\gamma}n_{\zeta}n_{\alpha}=64$, where
$(\xi,\gamma,\zeta,\alpha)\in\mathbb{R}^{n_{\xi}}\times\mathbb{R}^{n_{\gamma}}\times\mathbb{R}^{n_{\zeta}}\times\mathbb{R}^{n_{\alpha}}$

  * •

$M=2M_{r}+3M_{a}=256$

## References

  * Haykin 1998 Haykin, S. _Neural Networks: A Comprehensive Foundation_ , 2nd ed.; Prentice Hall PTR: Upper Saddle River, NJ, USA, 1998
  * Chen and Lin 2014 Chen, X.-W.; Lin, X. Big Data Deep Learning: Challenges and Perspectives. _IEEE Access_ 2014, _2_ , 514–525
  * Schmidhuber 2015 Schmidhuber, J. Deep learning in neural networks: An overview. _Neural Networks_ 2015, _61_ , 85–117
  * Jordan and Mitchell 2015 Jordan, M. I.; Mitchell, T. M. Machine learning: trends, perspectives, and prospects. _Science_ 2015, _349_ , 255–260
  * Sze et al. 2017 Sze, V.; Chen, Y.-H.; Yang, T.-J.; Emer, J. S. Efficient Processing of Deep Neural Networks: A Tutorial and Survey. _Proceedings of the IEEE_ 2017, _105_ , 2295–2329
  * Hutson 2018 Hutson, M. Artificial intelligence faces reproducibility crisis. _Science_ 2018, _359_ , 725–726
  * Hartana and Richards 1993 Hartana, R. K.; Richards, G. G. Constrained neural network-based identification of harmonic sources. _IEEE Transactions on Industry Applications_ 1993, _29_ , 202–208
  * Stewart and Ermon 2017 Stewart, R.; Ermon, S. Label-Free Supervision of Neural Networks with Physics and Domain Knowledge. AAAI. 2017; pp 2576–2582
  * Bouzerdoum and Pattison 1993 Bouzerdoum, A.; Pattison, T. R. Neural network for quadratic optimization with bound constraints. _IEEE transactions on neural networks_ 1993, _4_ , 293–304
  * Xia et al. 2005 Xia, Y. S.; Feng, G.; Wang, J. A primal-dual neural network for online resolving constrained kinematic redundancy in robot motion control. _IEEE Transactions on Systems, Man, and Cybernetics, Part B (Cybernetics)_ 2005, _35_ , 54–64
  * Rudd and Ferrari 2015 Rudd, K.; Ferrari, S. A constrained integration (CINT) approach to solving partial differential equations using artificial neural networks. _Neurocomputing_ 2015, _155_ , 277 – 285
  * Ling et al. 2016 Ling, J.; Jones, R.; Templeton, J. Machine learning strategies for systems with invariance properties. _Journal of Computational Physics_ 2016, _318_ , 22–35
  * Hennigh 2017 Hennigh, O. Lat-Net: Compressing Lattice Boltzmann Flow Simulations using Deep Neural Networks. _ArXiv e-prints_ 2017,
  * Baker et al. 2019 Baker, N.; Alexander, F.; Bremer, T.; Hagberg, A.; Kevrekidis, Y.; Najm, H.; Parashar, M.; Patra, A.; Sethian, J.; Wild, S.; Willcox, K. Workshop Report on Basic Research Needs for Scientific Machine Learning: Core Technologies for Artificial Intelligence. 2019,  https://www.osti.gov/servlets/purl/1478744
  * Ramakrishnan et al. 2015 Ramakrishnan, R.; Dral, P. O.; Rupp, M.; von Lilienfeld, O. A. Big Data Meets Quantum Chemistry Approximations: The $\Delta$-Machine Learning Approach. _Journal of Chemical Theory and Computation_ 2015, _11_ , 2087–2096, PMID: 26574412
  * Wei et al. 2016 Wei, J. N.; Duvenaud, D.; Aspuru-Guzik, A. Neural Networks for the Prediction of Organic Chemistry Reactions. _ACS Central Science_ 2016, _2_ , 725–732
  * Khorshidi and Peterson 2016 Khorshidi, A.; Peterson, A. A. Amp: A modular approach to machine learning in atomistic simulations. _Computer Physics Communications_ 2016, _207_ , 310–324
  * Liu et al. 2017 Liu, B.; Ramsundar, B.; Kawthekar, P.; Shi, J.; Gomes, J.; Luu Nguyen, Q.; Ho, S.; Sloane, J.; Wender, P.; Pande, V. Retrosynthetic Reaction Prediction Using Neural Sequence-to-Sequence Models. _ACS Central Science_ 2017, _3_ , 1103–1113
  * Liu et al. 2018 Liu, Q.; Zhou, X.; Zhou, L.; Zhang, Y.; Luo, X.; Guo, H.; Jiang, B. Constructing High-Dimensional Neural Network Potential Energy Surfaces for Gas–Surface Scattering and Reactions. _The Journal of Physical Chemistry C_ 2018, _122_ , 1761–1769
  * Ferguson 2018 Ferguson, A. L. ACS Central Science Virtual Issue on Machine Learning. _ACS Central Science_ 2018, _4_ , 938–941
  * Amabilino et al. 2019 Amabilino, S.; Bratholm, L. A.; Bennie, S. J.; Vaucher, A. C.; Reiher, M.; Glowacki, D. R. Training Neural Nets to Learn Reactive Potential Energy Surfaces using Interactive Quantum Chemistry in Virtual Reality. _The Journal of Physical Chemistry A_ 2019, _0_ , null, PMID: 30892040
  * Unke et al. 2021 Unke, O. T.; Chmiela, S.; Sauceda, H. E.; Gastegger, M.; Poltavsky, I.; Schütt, K. T.; Tkatchenko, A.; Müller, K.-R. Machine Learning Force Fields. _Chemical Reviews_ 2021, _121_ , 10142–10186, PMID: 33705118
  * Sun 1998 Sun, H. . _J. Phys. Chem. B_ 1998, _102_ , 7338
  * Kirschner et al. 2008 Kirschner, K.; Yongye, A.; Tschampel, S.; no, J. G.-O.; Daniels, C.; Foley, B.; Woods, R. . _J. Comput. Chem._ 2008, _29_ , 622
  * Huang and MacKerell Jr 2013 Huang, J.; MacKerell Jr, A. D. CHARMM36 all-atom additive protein force field: Validation based on comparison to NMR data. _Journal of Computational Chemistry_ 2013, _34_ , 2135–2145
  * Maier et al. 2015 Maier, J.; Martinez, C.; Kasavajhala, K.; Wickstrom, L.; Hauser, K.; Simmerling, C. . _J. Chem. Theory Comput._ 2015, _11_ , 3696
  * Chen et al. 2017 Chen, C.; Deng, Z.; Tran, R.; Tang, H.; Chu, I.-H.; Ong, S. Accurate Force Field for Molybdenum by Machine Learning Large Materials Data. _Phys. Rev. Materials_ 2017, _1_ , 043603
  * Wood and Thompson 2017 Wood, M. A.; Thompson, A. P. https://arxiv.org/abs/1702.07042
  * Behler et al. 2008 Behler, J.; Martoňák, R.; Donadio, D.; Parrinello, M. Pressure-induced phase transitions in silicon studied by neural network-based metadynamics simulations. _physica status solidi (b)_ 2008, _245_ , 2618–2629
  * Handley and Popelier 2010 Handley, C. M.; Popelier, P. L. A. Potential Energy Surfaces Fitted by Artificial Neural Networks. _The Journal of Physical Chemistry A_ 2010, _114_ , 3371–3383
  * Kondati Natarajan et al. 2015 Kondati Natarajan, S.; Morawietz, T.; Behler, J. Representing the potential-energy surface of protonated water clusters by high-dimensional neural network potentials. _Phys. Chem. Chem. Phys._ 2015, _17_ , 8356–8371
  * Kolb et al. 2016 Kolb, B.; Zhao, B.; Li, J.; Jiang, B.; Guo, H. Permutation invariant potential energy surfaces for polyatomic reactions using atomistic neural networks. _The Journal of Chemical Physics_ 2016, _144_ , 224103
  * Ho et al. 2016 Ho, T. H.; Pham-Tran, N.-N.; Kawazoe, Y.; Le, H. M. Ab Initio Investigation of O–H Dissociation from the Al–OH2 Complex Using Molecular Dynamics and Neural Network Fitting. _The Journal of Physical Chemistry A_ 2016, _120_ , 346–355, PMID: 26741404
  * Hellström and Behler 2017 Hellström, M.; Behler, J. Structure of aqueous NaOH solutions: insights from neural-network-based molecular dynamics simulations. _Phys. Chem. Chem. Phys._ 2017, _19_ , 82–96
  * Behler 2017 Behler, J. First Principles Neural Network Potentials for Reactive Simulations of Large Molecular and Condensed Systems. _Angewandte Chemie International Edition_ 2017, _56_ , 12828–12840
  * Botu et al. 2017 Botu, V.; Batra, R.; Chapman, J.; Ramprasad, R. Machine Learning Force Fields: Construction, Validation, and Outlook. _The Journal of Physical Chemistry C_ 2017, _121_ , 511–522
  * Yao et al. 2017 Yao, K.; Herr, J. E.; Brown, S. N.; Parkhill, J. Intrinsic Bond Energies from a Bonds-in-Molecules Neural Network. _The Journal of Physical Chemistry Letters_ 2017, _8_ , 2689–2694, PMID: 28573865
  * Pietrucci 2017 Pietrucci, F. Strategies for the exploration of free energy landscapes: Unity in diversity and challenges ahead. _Reviews in Physics_ 2017, _2_ , 32 – 45
  * Schütt et al. 2017 Schütt, K. T.; Arbabzadah, F.; Chmiela, S.; Müller, K. R.; Tkatchenko, A. Quantum-chemical insights from deep tensor neural networks. _Nature Communications_ 2017, _8_ , 13890
  * Schütt et al. 2017 Schütt, K. T.; Kindermans, P.-J.; Sauceda, H. E.; Chmiela, S.; Tkatchenko, A.; Müller, K.-R. SchNet: A Continuous-filter Convolutional Neural Network for Modeling Quantum Interactions. Proceedings of the 31st International Conference on Neural Information Processing Systems. USA, 2017; pp 992–1002
  * Yao et al. 2018 Yao, K.; Herr, J. E.; Toth, D. W.; Mckintyre, R.; Parkhill, J. The TensorMol-0.1 model chemistry: a neural network augmented with long-range physics. _Chem. Sci._ 2018, _9_ , 2261–2269
  * Lubbers et al. 2018 Lubbers, N.; Smith, J. S.; Barros, K. Hierarchical modeling of molecular energies using a deep neural network. _The Journal of Chemical Physics_ 2018, _148_ , 241715
  * Jiang et al. 2020 Jiang, B.; Li, J.; Guo, H. High-Fidelity Potential Energy Surfaces for Gas-Phase and Gas–Surface Scattering Processes from Machine Learning. _The Journal of Physical Chemistry Letters_ 2020, _11_ , 5120–5131
  * Blank et al. 1995 Blank, T.; Brown, S.; Calhoun, A.; Doren, D. Neural network models of potential energy surfaces. _J. Chem. Phys._ 1995, _103_ , 4129
  * Gassner et al. 1998 Gassner, H.; Probst, M.; Lauenstein, A.; Hermansson, K. Representation of Intermolecular Potential Functions by Neural Networks. _The Journal of Physical Chemistry A_ 1998, _102_ , 4596–4605
  * Lorenz et al. 2004 Lorenz, S.; Groß, A.; Scheffler, M. Representing high-dimensional potential-energy surfaces for reactions at surfaces by neural networks. _Chemical Physics Letters_ 2004, _395_ , 210 – 215
  * Rodriguez et al. 2018 Rodriguez, A.; d’Errico, M.; Facco, E.; Laio, A. Computing the Free Energy without Collective Variables. _Journal of Chemical Theory and Computation_ 2018, _14_ , 1206–1215
  * Schneider et al. 2017 Schneider, E.; Dai, L.; Topper, R. Q.; Drechsel-Grau, C.; Tuckerman, M. E. Stochastic Neural Network Approach for Learning High-Dimensional Free Energy Surfaces. _Phys. Rev. Lett._ 2017, _119_ , 150601
  * Stecher et al. 2014 Stecher, T.; Bernstein, N.; Csányi, G. Free Energy Surface Reconstruction from Umbrella Samples Using Gaussian Process Regression. _Journal of Chemical Theory and Computation_ 2014, _10_ , 4079–4097
  * Chen et al. 2015 Chen, M.; Yu, T.-Q.; Tuckerman, M. Locating landmarks on high-dimensional free energy surfaces. _PNAS_ 2015, _112_ , 3235–3240
  * Zheng et al. 2013 Zheng, W.; Rohrdanz, M. A.; Clementi, C. Rapid Exploration of Configuration Space with Diffusion-Map-Directed Molecular Dynamics. _The Journal of Physical Chemistry B_ 2013, _117_ , 12769–12776
  * Chiavazzo et al. 2017 Chiavazzo, E.; Covino, R.; Coifman, R.; Gear, C.; Georgiou, A.; Hummer, G.; Kevrekidis, I. Intrinsic map dynamics exploration for uncharted effective free-energy landscapes. _PNAS_ 2017, _114_ , E5494–E5503
  * Georgiou et al. 2017 Georgiou, A. S.; Bello-Rivas, J. M.; Gear, C. W.; Wu, H.-T.; Chiavazzo, E.; Kevrekidis, I. G. An Exploration Algorithm for Stochastic Simulators Driven by Energy Gradients. _Entropy_ 2017, _19_ , 294
  * Han et al. 2018 Han, J.; Zhang, L.; Car, R.; E, W. Deep Potential: A General Representation of a Many-Body Potential Energy Surface. _Commun. Comput. Phys._ 2018, _23_ , 629–639
  * Rupp et al. 2012 Rupp, M.; Tkatchenko, A.; Müller, K.-R.; von Lilienfeld, O. A. Fast and Accurate Modeling of Molecular Atomization Energies with Machine Learning. _Phys. Rev. Lett._ 2012, _108_ , 058301
  * Hansen et al. 2013 Hansen, K.; Montavon, G.; Biegler, F.; Fazli, S.; Rupp, M.; Scheffler, M.; von Lilienfeld, O. A.; Tkatchenko, A.; Müller, K.-R. Assessment and Validation of Machine Learning Methods for Predicting Molecular Atomization Energies. _Journal of Chemical Theory and Computation_ 2013, _9_ , 3404–3419
  * Hansen et al. 2015 Hansen, K.; Biegler, F.; Ramakrishnan, R.; Pronobis, W.; von Lilienfeld, O. A.; Müller, K.-R.; Tkatchenko, A. Machine Learning Predictions of Molecular Properties: Accurate Many-Body Potentials and Nonlocality in Chemical Space. _The Journal of Physical Chemistry Letters_ 2015, _6_ , 2326–2331, PMID: 26113956
  * Weininger et al. 1989 Weininger, D.; Weininger, A.; Weininger, J. L. SMILES. 2. Algorithm for generation of unique SMILES notation. _Journal of Chemical Information and Computer Sciences_ 1989, _29_ , 97–101
  * Sanchez-Lengeling and Aspuru-Guzik 2018 Sanchez-Lengeling, B.; Aspuru-Guzik, A. Inverse molecular design using machine learning: Generative models for matter engineering. _Science_ 2018, _361_ , 360–365
  * Bartók et al. 2010 Bartók, A. P.; Payne, M. C.; Kondor, R.; Csányi, G. Gaussian Approximation Potentials: The Accuracy of Quantum Mechanics, without the Electrons. _Phys. Rev. Lett._ 2010, _104_ , 136403
  * Rogers and Hahn 2010 Rogers, D.; Hahn, M. Extended-Connectivity Fingerprints. _Journal of Chemical Information and Modeling_ 2010, _50_ , 742–754, PMID: 20426451
  * Bartók et al. 2013 Bartók, A. P.; Kondor, R.; Csányi, G. On representing chemical environments. _Phys. Rev. B_ 2013, _87_ , 184115
  * Duvenaud et al. 2015 Duvenaud, D. K.; Maclaurin, D.; Iparraguirre, J.; Bombarell, R.; Hirzel, T.; Aspuru-Guzik, A.; Adams, R. P. In _Advances in Neural Information Processing Systems 28_ ; Cortes, C., Lawrence, N. D., Lee, D. D., Sugiyama, M., Garnett, R., Eds.; Curran Associates, Inc., 2015; pp 2224–2232
  * Pham et al. 2018 Pham, T.-L.; Nguyen, N.-D.; Nguyen, V.-D.; Kino, H.; Miyake, T.; Dam, H.-C. Learning structure-property relationship in crystalline materials: A study of lanthanide–transition metal alloys. _The Journal of Chemical Physics_ 2018, _148_ , 204106
  * Hirn et al. 2017 Hirn, M.; Mallat, S.; Poilvert, N. Wavelet Scattering Regression of Quantum Chemical Energies. _Multiscale Model. Simul._ 2017, _15_ , 827–863
  * Sidky and Whitmer 2018 Sidky, H.; Whitmer, J. K. Learning free energy landscapes using artificial neural networks. _The Journal of Chemical Physics_ 2018, _148_ , 104111
  * Tang et al. 2018 Tang, Y.-H.; Zhang, D.; Karniadakis, G. E. An atomistic fingerprint algorithm for learning ab initio molecular force fields. _The Journal of Chemical Physics_ 2018, _148_ , 034101
  * Behler and Parrinello 2007 Behler, J.; Parrinello, M. Generalized Neural-Network Representation of High-Dimensional Potential-Energy Surfaces. _Physical Review Letters_ 2007, _98_ , 146401
  * Behler 2011 Behler, J. Atom-centered symmetry functions for constructing high-dimensional neural network potentials. _J. Chem. Phys._ 2011, _134_ , 074106
  * Behler 2015 Behler, J. Constructing High-Dimensional Neural Network Potentials: A Tutorial Review. _Int. J. Quant. Chem._ 2015, _115_ , 1032–1050
  * Behler 2017 Behler, J. Neural network potential-energy surfaces for atomistic simulations. _Chem. Modell._ 2017, _7_ , 1–41
  * Jose et al. 2012 Jose, K. V. J.; Artrith, N.; Behler, J. Construction of high-dimensional neural network potentials using environment-dependent atom pairs. _The Journal of Chemical Physics_ 2012, _136_ , 194111
  * Cubuk et al. 2017 Cubuk, E.; Malone, B.; Onat, B.; Waterland, A.; Kaxiras, E. Representations in neural network based empirical potentials. _J. Chem. Phys._ 2017, _147_ , 024104
  * Khaliullin et al. 2011 Khaliullin, R. Z.; Eshet, H.; Kühne, T. D.; Behler, J.; Parrinello, M. Nucleation mechanism for the direct graphite-to-diamond phase transition. _Nature Materials_ 2011, _10_ , 693
  * Behler et al. 2008 Behler, J.; Martoňák, R.; Donadio, D.; Parrinello, M. Metadynamics Simulations of the High-Pressure Phases of Silicon Employing a High-Dimensional Neural Network Potential. _Phys. Rev. Lett._ 2008, _100_ , 185501
  * Artrith and Behler 2012 Artrith, N.; Behler, J. High-dimensional neural network potentials for metal surfaces: A prototype study for copper. _Phys. Rev. B_ 2012, _85_ , 045439
  * Artrith and Kolpak 2014 Artrith, N.; Kolpak, A. M. Understanding the Composition and Activity of Electrocatalytic Nanoalloys in Aqueous Solvents: A Combination of DFT and Accurate Neural Network Potentials. _Nano Letters_ 2014, _14_ , 2670–2676, PMID: 24742028
  * Galvelis and Sugita 2017 Galvelis, R.; Sugita, Y. Neural Network and Nearest Neighbor Algorithms for Enhancing Sampling of Molecular Dynamics. _Journal of Chemical Theory and Computation_ 2017, _13_ , 2489–2500
  * Onat et al. 2018 Onat, B.; Cubuk, E. D.; Malone, B. D.; Kaxiras, E. Implanted neural network potentials: Application to Li-Si alloys. _Phys. Rev. B_ 2018, _97_ , 094106
  * Smith et al. 2017 Smith, J. S.; Isayev, O.; Roitberg, A. E. ANI-1: an extensible neural network potential with DFT accuracy at force field computational cost. _Chem. Sci._ 2017, _8_ , 3192–3203
  * Smith et al. 2018 Smith, J. S.; Nebgen, B.; Lubbers, N.; Isayev, O.; Roitberg, A. E. Less is more: Sampling chemical space with active learning. _The Journal of Chemical Physics_ 2018, _148_ , 241733
  * Chmiela et al. 2017 Chmiela, S.; Tkatchenko, A.; Sauceda, H. E.; Poltavsky, I.; Schütt, K. T.; Müller, K.-r.; Poltavsky, I.; Sch, K. T. Machine learning of accurate energy-conserving molecular force fields. _Science Advances_ 2017, _3_ , e1603015
  * Christensen and von Lilienfeld 2020 Christensen, A. S.; von Lilienfeld, O. A. On the role of gradients for machine learning of molecular energies and forces. _Machine Learning: Science and Technology_ 2020, _1_ , 045018
  * Bartók and Csányi 2015 Bartók, A. P.; Csányi, G. Gaussian approximation potentials: A brief tutorial introduction. _International Journal of Quantum Chemistry_ 2015, _115_ , 1051–1057
  * Huang et al. 2017 Huang, S.-D.; Shang, C.; Zhang, X.-J.; Liu, Z.-P. Material discovery by combining stochastic surface walking global optimization with a neural network. _Chem. Sci._ 2017, _8_ , 6327–6337
  * Cooper et al. 2020 Cooper, A. M.; Kästner, J.; Urban, A.; Artrith, N. Efficient training of ANN potentials by including atomic forces via Taylor expansion and application to water and a transition-metal oxide. _npj Computational Materials_ 2020, _6_ , 54
  * Pukrittayakamee et al. 2009 Pukrittayakamee, A.; Malshe, M.; Hagan, M.; Raff, L. M.; Narulkar, R.; Bukkapatnum, S.; Komanduri, R. Simultaneous fitting of a potential-energy surface and its corresponding force fields using feedforward neural networks. _The Journal of Chemical Physics_ 2009, _130_ , 134101
  * Nguyen and Le 2012 Nguyen, H. T. T.; Le, H. M. Modified Feed-Forward Neural Network Structures and Combined-Function-Derivative Approximations Incorporating Exchange Symmetry for Potential Energy Surface Fitting. _The Journal of Physical Chemistry A_ 2012, _116_ , 4629–4638, PMID: 22548349
  * Nguyen-Truong and Le 2015 Nguyen-Truong, H. T.; Le, H. M. An implementation of the Levenberg–Marquardt algorithm for simultaneous-energy-gradient fitting using two-layer feed-forward neural networks. _Chemical Physics Letters_ 2015, _629_ , 40–45
  * Nandi et al. 2019 Nandi, A.; Qu, C.; Bowman, J. M. Using Gradients in Permutationally Invariant Polynomial Potential Fitting: A Demonstration for CH4 Using as Few as 100 Configurations. _Journal of Chemical Theory and Computation_ 2019, _15_ , 2826–2835, PMID: 30896950
  * Conte et al. 2020 Conte, R.; Qu, C.; Houston, P. L.; Bowman, J. M. Efficient Generation of Permutationally Invariant Potential Energy Surfaces for Large Molecules. _Journal of Chemical Theory and Computation_ 2020, _16_ , 3264–3272, PMID: 32212729
  * Bowman et al. 2022 Bowman, J. M.; Qu, C.; Conte, R.; Nandi, A.; Houston, P. L.; Yu, Q. The MD17 datasets from the perspective of datasets for gas-phase “small” molecule potentials. _J. Chem. Phys_ 2022, _156_ , 240901
  * Zaverkin et al. 2023 Zaverkin, V.; Holzmüller, D.; Bonfirraro, L.; Kästner, J. Transfer learning for chemically accurate interatomic neural network potentials. _Phys. Chem. Chem. Phys._ 2023, _25_ , 5383–5396
  * Singraber et al. 2019 Singraber, A.; Morawietz, T.; Behler, J.; Dellago, C. Parallel Multistream Training of High-Dimensional Neural Network Potentials. _Journal of Chemical Theory and Computation_ 2019, _15_ , 3075–3092
  * Gastegger et al. 2017 Gastegger, M.; Behler, J.; Marquetand, P. Machine learning molecular dynamics for the simulation of infrared spectra. _Chem. Sci._ 2017, _8_ , 6924–6935
  * Nandi et al. 2019 Nandi, A.; Qu, C.; Bowman, J. M. Full and fragmented permutationally invariant polynomial potential energy surfaces for trans and cis N-methyl acetamide and isomerization saddle points. _J. Chem. Phys_ 2019, _151_ , 084306
  * Houston et al. 2020 Houston, P.; Conte, R.; Qu, C.; Bowman, J. M. Permutationally invariant polynomial potential energy surfaces for tropolone and H and D atom tunneling dynamics. _J. Chem. Phys_ 2020, _153_ , 024107
  * Chen and Goldsmith 2020 Chen, X.; Goldsmith, C. F. Accelerating Variational Transition State Theory via Artificial Neural Networks. _The Journal of Physical Chemistry A_ 2020, _124_ , 1038–1046
  * Unke et al. 2021 Unke, O. T.; Chmiela, S.; Gastegger, M.; Schütt, K. T.; Sauceda, H. E.; Müller, K.-R. SpookyNet: Learning force fields with electronic degrees of freedom and nonlocal effects. _Nature Communications_ 2021, _12_ , 7273
  * Herr et al. 2018 Herr, J. E.; Yao, K.; McIntyre, R.; Toth, D. W.; Parkhill, J. Metadynamics for training neural network model chemistries: A competitive assessment. _The Journal of Chemical Physics_ 2018, _148_ , 241710
  * Najm and Yang 2021 Najm, H. N.; Yang, Y. _AEVmod – Atomic Environment Vector Module Documentation_ ; 2021; Technical Report SAND2021–9473 https://doi.org/10.2172/1817835
  * Paszke et al. 2019 Paszke, A. et al. _Advances in Neural Information Processing Systems 32_ ; Curran Associates, Inc., 2019; pp 8024–8035
  * Jakob et al. 2017 Jakob, W.; Rhinelander, J.; Moldovan, D. pybind11 – Seamless operability between C++11 and Python. 2017; https://github.com/pybind/pybind11
  * Najm and Yang 2021 Najm, H. N.; Yang, Y. aevmod. 2021; https://github.com/sandialabs/aevmod
  * Phipps 2022 Phipps, E. Sacado. 2022; https://trilinos.github.io/sacado
  * 106 https://www.q-chem.com
  * Kingma and Ba 2017 Kingma, D. P.; Ba, J. Adam: A Method for Stochastic Optimization. 2017; https://doi.org/10.48550/arXiv.1412.6980
  * 108 https://pytorch.org/docs/stable/generated/torch.optim.Adam.html