††thanks: equal contribution††thanks: equal contribution

# Learning Pair Potentials using Differentiable Simulations

Wujie Wang  Department of Materials Science and Engineering, Massachusetts
Institute of Technology. 77 Massachusetts Ave. Cambridge MA 02139 USA
Zhenghao Wu  Eduard-Zintl-Institut für Anorganische und Physikalische Chemie,
Technische Universität Darmstadt, Alarich-Weiss-Str. 8, 64287 Darmstadt,
Germany  Rafael Gómez-Bombarelli<EMAIL_ADDRESS>Department of Materials
Science and Engineering, Massachusetts Institute of Technology. 77
Massachusetts Ave. Cambridge MA 02139 USA

(August 28, 2024)

###### Abstract

Learning pair interactions from experimental or simulation data is of great
interest for molecular simulations. We propose a general stochastic method for
learning pair interactions from data using differentiable simulations
(DiffSim). DiffSim defines a loss function based on structural observables,
such as the radial distribution function, through molecular dynamics (MD)
simulations. The interaction potentials are then learned directly by
stochastic gradient descent, using backpropagation to calculate the gradient
of the structural loss metric with respect to the interaction potential
through the MD simulation. This gradient-based method is flexible and can be
configured to simulate and optimize multiple systems simultaneously. For
example, it is possible to simultaneously learn potentials for different
temperatures or for different compositions. We demonstrate the approach by
recovering simple pair potentials, such as Lennard-Jones systems, from radial
distribution functions. We find that DiffSim can be used to probe a wider
functional space of pair potentials compared to traditional methods like
Iterative Boltzmann Inversion. We show that our methods can be used to
simultaneously fit potentials for simulations at different compositions and
temperatures to improve the transferability of the learned potentials.

##  I Introduction

Molecular simulation has been an invaluable technique in various fields of
modern science, including chemistry, physics, and materials scienceVlachakis
_et al._ (2014); Blumberger (2015); Gartner and Jayaraman (2019); Fish,
Wagner, and Keten (2021). Ideally, an atomistic computer simulation would
combine the accurate description of the Born-Oppenheimer potential energy
surface (PES) with the treatment of nuclear motion, both at the quantum
mechanics level. However, this is computationally infeasible for most of the
relevant molecular sizes and time scales. It is typically necessary to
approximate the PES with simpler and faster surrogate functions, such as force
fields or interatomic potentials. Therefore, the development of accurate and
computationally efficient force fields has become one of the most important
topics in molecular simulationsMayo, Olafson, and Goddard (1990); Jorgensen,
Maxwell, and Tirado-Rives (1996); Pastor and MacKerell (2011).

Covalent bonds in organic molecules are typically described in force fields as
harmonic or polynomial expressions in short-range many-body terms (bonds,
angles, torsions), while numerous types of description have been proposed for
noncovalent interactionsRapaport (2004); Cisneros  _et al._ (2016). Of all
these variants, the pairwise potential is most commonly used in molecular
simulations, as it has been shown to be sufficient to capture significant
molecular mechanics, and is also the most computationally cost-effective way
to describe noncovalent interactions between particlesCisneros  _et al._
(2016). For example, rigid water models with pairwise additive PESs such as
SPC/EBerendsen, Grigera, and Straatsma (1987) and TIP4PJorgensen  _et al._
(1983) recover many thermodynamic properties, e.g., densities and heat
capacities, that are consistent with experimental measurements across a wide
range of temperatures. These models with pairwise additive PESs are efficient
in terms of computational cost and scalabilityGlaser  _et al._ (2015). It is
noted that the many-body expansion of PESs has been widely appreciated for its
role in accurately predicting behaviors or properties, such as phase
transitions of water, which the pairwise additive PES cannot predictCisneros
_et al._ (2016); Albaugh  _et al._ (2016). However, these potential
formulations that carry the many-body effect are usually not as efficient;
depending on the size of the system, they are usually 10 to 1000 times slower
to simulate than pairwise potentials.Plimpton and Thompson (2012)

In addition to all-atom potentials, the derivation of effective pair
potentials also plays an important role for the efficient simulation of
coarse-grained systems.Müller-Plathe (2002); Souza  _et al._ (2021); Dhamankar
and Webb (2021) Various variational frameworks have been developed to derive
coarse-grained potentials. Based on the idea of force matchingErcolessi and
Adams (1994), Izvekov and Voth developed a force-based coarse-graining method,
i.e., multiscale coarse-grained method (MS-CG)Izvekov and Voth (2005). The MS-
CG method aims to reproduce the many-body potential of the mean force (PMF) of
the fine-grained system (FG) in the mapped CG representation by matching the
force acting on the CG bead.Noid  _et al._ (2007) This framework can also be
extended to different architectures of force fields, e.g., machine learning
potentials, optimized with stochastic gradient descent. Zhang  _et al._
(2018); Wang  _et al._ (2019). However, it is a known issue that the
potentials learned from force matching do not necessarily recover target
structural distributions due to missing cross-correlations between degrees of
freedomRudzinski and Noid (2012) The use of machine learning potentials
naturally incorporates many-body correlation terms, but they often suffer from
instabilities in dynamics and often fail to generate robust simulation
trajectories. Wang  _et al._ (2019); Stocker  _et al._ (2022). Beyond force
matching, Shell et al. developed a series of coarse-graining approaches using
relative entropy to measure coarse-graining information lossChaimovich and
Shell (2011). The standard relative entropy framework uses the Newton-Raphson
algorithm as the solver, which, however, is computationally expensive for
optimizing models such as deep neural networks that have a large number of
parameters. Kocer, Ko, and Behler (2022) There is another class of structure-
based coarse-graining method which does not require atomistic position or
force data and directly aims at reproducing structural distributions, of which
iterative Boltzmann inversion (IBI) is the most notable.Reith, Pütz, and
Müller-Plathe (2003) IBI constructs potentials for a CG model by iteratively
correcting errors in simulated radial distribution function until
convergence.Reith, Pütz, and Müller-Plathe (2003) It has been widely employed
for coarse-graining various soft-matter systems such as organic liquidsQian
_et al._ (2008), ionic liquidsKarimi-Varzaneh  _et al._ (2010) and
polymersQian  _et al._ (2008). IBI does not necessarily require reference all-
atom simulations. However, its application is limited to tabulated pair
potentials and does not generalize to other forms of potential.

Most coarse-graining approaches only aim at reproducing statistics in a single
thermodynamic state, so the developed potentials are rarely transferable.
Several approaches have been proposed to improve the transferability of the
coarse-graining frameworks mentioned in previous sections. For example,
Mullinax and Noid proposed an extended ensemble method in the MS-CG
frameworkMullinax and Noid (2009), in which datasets from multiple equilibrium
ensembles are collected to develop CG potentials that reproduce the potential
of mean force of FG systems in these ensembles. This approach requires a
considerable amount of atomistic force data, which hinders its application.

CG models created by the IBI method can reproduce the structural correlations
of the FG model at a single thermodynamic state, but they do not produce
potentials that are transferable to a wide range of thermodynamic states with
different temperatures and compositions. Attempts have been made to derive
transferable CG potentials using multi-state data based on IBI, but these are
still challenging, for example, in finding appropriate weighting
functions.Moore, Iacovella, and McCabe (2014a)

In recent years, machine learning (ML) has emerged as an important tool for
molecular modeling. Much of the success has been based on the application of
differentiable programming, which enables the exact generation of gradient
programs. The process of generating the gradient program is termed Automatic
Differentiation (AD). Although emerging as a relatively new tool, AD has shown
great progress in many domains of computational science, such as quantum
computingLiao  _et al._ (2019) and fluid dynamicsSchenck and Fox (2018). For
the modeling of PES, AD has transformed the fitting of force fields from
first-principle calculations: once a forward program that transforms nuclear
coordinates to a scalar (energy) is constructed, the gradient (force) program
is automatically generated with AD. The generated force field is also
trainable with respect to the data, enabling parameterization of ML potentials
with near ab initio accuracy.Unke  _et al._ (2021)

Differentiable programming is applicable to many controlled iterative
mathematical procedures in the physical sciences, ranging from learning
exchange correlation functionals for Density Functional Theory(DFT)
calculations from dataKasim and Vinko (2021); Kanungo, Zimmerman, and Gavini
(2019) or fine-tuning molecular basis setsTamayo-Mendoza  _et al._ (2018).
This requires differentiating through mathematical procedures such as fixed-
point iteration Blondel  _et al._ (2021) and eigendecompositionMagnus (1985).
Similarly, differentiation operations of molecular simulations, as solutions
to some ordinary differential equation (ODE), can also be defined.Chen  _et
al._ (2018) Recently, many applications have been proposed using
differentiable molecular dynamics (DiffSim). Wang, Axelrod, and Gómez-
Bombarelli (2020); Schoenholz and Cubuk (2020); Doerr  _et al._ (2021); Thaler
and Zavadlav (2021); Greener and Jones (2021) Specifically, recent work on
DiffSim suggests a differentiable top-down approach to learning interatomic
potentials directly based on macroscopic observations from experiments or
simulations. Compared to previously proposed machine learning potentials that
are trained from the bottom-up, DiffSim enables experimentally informed
parameterization of force fields. However, to improve fitting flexibility,
these methods utilize neural network (NN) potentials based on graph neural
networks, which have limited interpretability and are less scalable for large-
scale simulations. Another method of machine learning-assisted top-down
parameterization involves directly predicting force field parameters from
macroscopic structure correlation using a machine-learned function to bypass
the inverse problem.Moradzadeh and Aluru (2019) This approach requires
abundant labeled simulation data and is limited to a simple functional form of
Lennard-Jones type, which has limited capacity to encode complex structural
correlations.

In this work, we propose simulation-informed parameterization of pair
potentials using DiffSim. Compared to structure-based methods like IBI, our
approach can be used to optimize flexible functional forms and can be easily
configured to fit multiple systems simultaneously; compared to force-matching-
based methods, our method does not require atomistic position or force data
and aims at directly reproducing the target distributions. First, we introduce
our method and training protocols. For numerical examples, we demonstrate our
method on simple liquid systems where ground truth potentials are known. We
show that our method is good at obtaining a diverse set of pair potentials
that reproduce the target liquid structure. To demonstrate the flexibility of
the method, we show that the model can learn pair potentials for binary
mixture systems involving multiple interactions and compositions. Our results
indicate that fitting multiple states simultaneously helps to improve the
transferability of learned potentials. Furthermore, our method can be used to
learn coarse-grained potentials where there is no ground truth for pairwise
potentials. We show that the learning of interactions can be combined with
flexible functional forms that involve temperature dependence for improved
transferability.

##  II Method

###  II.1 Observable-based Differentiable Simulations

Figure 1: A schematic diagram for our learning framework

Differentiable programming is a programming paradigm in which the gradient of
a mathematical procedure can be automatically synthesized. In deep learning,
automatically computed gradients allow gradient-based optimization using
“backpropagation" to train deep neural networks. As universal approximators,
neural networks of various architectures have been used as a drop-in module
for optimization and learning. However, without the appropriate incorporation
of inductive biases, many NN models tend to overfit the data and often fail to
generalize.Battaglia  _et al._ (2018) With the development of differentiable
algorithms, many sampling and optimization procedures can now be directly
incorporated into the machine learning framework. These structures serve as
inductive biases to regularize learning, improve data efficiency, and improve
model interpretability. Baxter (2000) Machine learning of dynamical systems
and differentiable equations is a domain where differentiable sampling
operations are used as strong inductive biases. Han, Jentzen, and Weinan
(2018); Chen (2017); Brunton and Kutz (2022)

In molecular dynamics, the simulation of systems of many particles requires
integrating per-particle state variables in time following the governing
equations of motion. Although continuous, the computational solution of state
variable dynamics is solved in discrete time steps, involving the composition
of differentiable mathematical procedures. With AD, the computation of
observables from molecular simulations is, therefore, end-to-end
differentiable. Specifically, given a defined forward computational program,
the AD framework traces the computation backward and recursively applies chain
rules as vector Jacobian products. This requires the implementation of the
accompanying gradient function for every primitive operation used for the
forward computation. AD can be applied to arbitrarily complex numerical
programs, as long as the gradient of each elementary step can be computed.
Auto-differentiable code has been implemented in many modern numerical
computing packages. Abadi  _et al._ (2016); Frostig, Johnson, and Leary
(2018); Paszke  _et al._ (2019)

The evolution of atomic coordinates in MD is typically obtained by integrating
the equations of motion through a discretized update rule of the following
form:

$x_{t}=F_{\theta}(x_{t-1},v_{t-1},\eta_{t-1})dt+x_{t-1}$ (1)

where $x$ is the coordinate, $v$ is the velocity, $\eta$ is the bath variable,
$\theta$ represents the set of learning parameters, and $F$ is the discretized
update procedure. The simulation observables $O$ are usually extracted from
the trajectories. Examples of such observables include radial distribution
functions (RDF), velocity autocorrelation functions, etc. To learn or control
the simulated observables, a scalar loss function $L$ must be defined to
quantify the discrepancy between the simulated and target observables. The
loss function $L$ needs to be minimized using gradient descent with gradient
on $\theta$. We first consider the computation of the gradient on $x$, or the
adjoint state, at each time step:

$\displaystyle\frac{dL}{dx_{t}}$ $\displaystyle=\frac{\partial L}{\partial
x_{t}}+\frac{dL}{dx_{t+1}}\frac{\partial x_{t+1}}{\partial x_{t}}$ (2)
$\displaystyle=\frac{\partial L}{\partial
x_{t}}+\frac{dL}{dx_{t+1}}+\frac{dL}{dx_{t+1}}\frac{\partial
F_{\theta}}{\partial x_{t}}dt$

With $\frac{dL}{dxt_{t}}$ at each step $t$, the gradient of $\theta$ can be
accumulated using the chain rule:

$\frac{dL}{d\theta}=\sum_{t}^{T}\frac{dL}{dx_{t}}\frac{dx_{t}}{d\theta}=\sum_{t}^{T}\frac{dL}{dx_{t}}\frac{dF_{\theta}(x_{t-1},...)}{d\theta}dt$
(3)

Alternatively, one can treat eq. 3 as a dynamical system solved backward in
time.Chen  _et al._ (2018) With the gradient accumulated from MD simulations
for $T$ steps, gradient-based optimization can be applied to minimize the
observable optimization objective $L$ until convergence.

####  II.1.1 Differentiable Computation of Radial Distributions Functions

The goal of our method is to learn pair interactions from RDFs, which are
treated as our observable $O$. RDF is defined as:

$g(r)=\frac{V}{N^{2}4\pi r^{2}}p(r)$ (4)

where $V$ is the volume, $N$ is the number of particles and $r$ is the radial
distance from a target particle; $p(r)$ is the histogram of pair distances. To
propagate the gradient from $O$ to the model parameters, $p(r)$ must be
differentiable, while the operation to obtain a histogram is not
differentiable. To make the computation of $g(r)$ differentiable, we apply the
kernel density trick with Gaussian kernels to make the procedure for computing
histograms differentiable. Chen (2017) By definition,

$p(r)=\sum_{i\neq j}\frac{\delta(r_{ij}-r)}{\sum_{i\neq j}1}$ (5)

where $r_{ij}$ is the distance between the particle $i$ and $j$, and
$\delta(\cdot)$ is the Dirac Delta function. To estimate $\delta(\cdot)$ in a
differentiable way, we approximate the probability density of two particles
that are $\mu_{k}$ apart in distance as:

$p(\mu_{k})\approx\sum_{i\neq j}p_{k}(r_{ij})=\sum_{i\neq
j}\frac{e^{-(\mu_{k}-r_{ij})^{2}/d}}{\sum_{k^{\prime}=0}^{k=K-1}e^{-(\mu_{k^{\prime}}-r_{ij})^{2}/d}}$
(6)

where $d=\mu_{k+1}-\mu_{k}$ is the size of the histogram bin and $K$ is the
total number of evenly spaced Gaussian kernels. The summation in the
denominator ensures that the histogram is normalized over
$[\mu_{0},\mu_{K-1}]$. This approximation provides a density estimate over
binned domains. Similar methods with triangular kernels have also been
proposed before.Ustinova and Lempitsky (2016) As $d\rightarrow 0$, the
Gaussian kernel function approximates $\delta(\cdot)$.

####  II.1.2 Design for Neural Pair Potentials

We propose a design for the neural network
$NN:\mathbb{R}\rightarrow\mathbb{R}$ with a parameter set $\theta$ as an
approximation function for pair potentials. We choose neural networks (NN) to
parameterize pair potentials, because NNs are universal approximatorsHornik,
Stinchcombe, and White (1989) and can be easily stacked in a model for end-to-
end training. The approach is also applicable to parameterizing fixed
functional forms such as Lennard-Jones or other choices. Given a distance
$r_{i}j$ within a cutoff distance $r_{\mathrm{cut}}$, we expand the distance
with a Gaussian basis to obtain a distance feature:

$e_{k}(r_{ij})=e^{-\frac{(\mu_{k}-r_{ij})^{2}}{d_{k}}}$ (7)

where the initial $d_{k}$ is learnable with the initial value set as
$r_{\mathrm{cut}}/K$ and $K$ is the number of Gaussian bases used to
parameterize $r_{ij}$. The $K$ dimensional distance feature is then fed into a
standard multilayer perceptron (MLP) with $n_{\mathrm{layer}}$ hidden layers
and $W$ hidden nodes per layer to parameterize a scalar energy output. We term
our neural pair potential $u_{\theta}$ with $\theta$ indicating the set of
learnable parameters that include $d_{k}$ and weights and biases in MLP. AD
can be applied to generate the force on an individual particle $i$ to be used
for simulation. The force program is also differentiable to receive gradient
signals to adjust the parameter set $\theta$. In addition to the MLP-based
pair potential, we also include an unlearnable prior pair potential
$U_{\mathrm{prior}}$ to guide the initial simulations to explore reasonable
regions of the configuration space. For $u_{\mathrm{prior}}$, we simply use a
repulsive potential of the form $\sum_{i,j}(\frac{\sigma}{r_{ij}})^{c}$ where
$c$ is the repulsive exponent, which is also treated as a hyperparameter. The
total pair potentials used for DiffSim take the form:

$U_{\theta}(x)=\sum_{m}\sum_{i,j\in\mathcal{P}_{m}}u^{m}_{\theta}(r_{ij})+u_{prior}(r_{ij})$
(8)

where $m$ indicates the type of pair interactions and $\mathcal{P}_{m}$
indicates the set of atom pairs that are governed by the pair potential of the
same type. In case faster convergence is needed, we also perform an extra pre-
training step to set the pair potentials with the direct Boltzmann
inversionTschöp  _et al._ (1998) so that
$u^{m}_{\theta}(r_{ij})+u_{prior}(r_{ij})\approx-kT\ln g^{m}(r)$ with
$g^{m}(r)$ being the target RDF.

####  II.1.3 Learning Protocols

In this section, we introduce our protocols for learning neural pair
potentials. The learning protocol follows the concept of learning-to-
simulate,Ruiz, Schulter, and Chandraker (2018) with the simulator performing
the coupled task of sampling data points and learning. We first initialize the
MD systems as crystal structures (FCCs or BCCs) at a target density. We then
apply DiffSim to evolve the system for $T$ steps and compute the observable of
interests. The examples used in this work are simulated in the canonical
ensemble (NVT) using the Nose-Hoover chain integrator.Martyna, Klein, and
Tuckerman (1992) To learn pair potentials, we are interested in learning a
$U_{\theta}$ that produces a target RDF $g_{\textnormal{ref.}}(r)$. We term
the RDF obtained from DiffSim as $g_{\textnormal{sim.}}(r)$. We then want to
minimize the L2 loss between targets and simulated RDFs:

$L=\sum_{m}\int_{r}(g^{m}_{\textnormal{sim.}}(r)-g^{m}_{\textnormal{ref.}}(r))^{2}dr$
(9)

where $m$ again indicates the type of pair potential. With DiffSim, we can
then obtain gradients on the parameters according to eq. 3. The simulation
explores the configuration space with an initial potential and accumulates
gradients based on the simulated trajectories.

Gradients can be fed into a gradient-based optimizer, such as AdamKingma and
Ba (2014), to update $U_{\theta}$. We simulate and update the pair potentials
for $N_{epoch}$ epochs until the learning converges. Because differentiable
simulations can be susceptible to exploding or vanishing gradients Ingraham
_et al._ (2018); Metz  _et al._ (2021), we tend not to evolve the system too
long before each gradient updates. In practice, we observe poor learning
outcomes with simulation steps that are more than 300 steps. We applied small
learning rates to perturb the simulation slightly between forward and backward
passes. The learning rate is also scheduled to decrease if it reaches a
plateau during learning. The decreasing learning rate facilitates fine-tuning
of the learned potentials when the sampled configurations produce an
observation close to the target. To report the equilibrium RDF, we simulate
our example systems sufficiently long to reflect the true RDF governed by the
learned pair interactions. In our experiments, we choose $T$ between 100 and
300 to ensure a stable gradient calculation and to provide sufficient time for
the system to evolve based on the updated force field.

###  II.2 Iterative Boltzmann Inversion

We compared our proposed learning protocol with iterative Boltzmann inversion
(IBI), a popular structure-based method to determine effective pair potentials
from RDF. Reith, Pütz, and Müller-Plathe (2003) IBI is the iterative version
of the Boltzmann inversion (BI).Tschöp  _et al._ (1998) BI is a simple method
based on the fact that the distribution of an independent degree of freedom
$x$ corresponds to a Boltzmann distribution in a canonical ensemble:

$p(x)\propto\exp(-\frac{u(x)}{k_{B}T})$ (10)

where $p(x)$ is the normalized distribution; $k_{B}$ and $T$ are the Boltzmann
constant and the temperature, respectively. The potential $u(x)$ can be
obtained by inverting the above equation:

$u(x)=-k_{B}T\log g_{\textnormal{ref.}}(x)$ (11)

The potential obtained is in the form of tabulated pair potentials. Such a
direct inversion ignores the indirect force contribution to the pair
statistics, and therefore introduces errors in reproducing the target RDF.
Tschöp  _et al._ (1998); Reith, Pütz, and Müller-Plathe (2003) As a natural
extension of BI, IBI iteratively corrects the potential using the difference
between the target distribution and the distribution obtained from the trial
potential energy function:

$u_{t+1}(x)=u_{t}(x)+\alpha\Delta u_{t}$ (12)

and

$\alpha\Delta
u_{t}=k_{B}T\ln\big{(}\frac{g_{t}(x)}{g_{\mathrm{\textnormal{ref.}}}(x)}\big{)}$
(13)

where $u_{t}$ is the effective potential at iteration $t$ and $\alpha$ is the
size of the update step. Here, $g_{t}(x)$ is the RDF simulated with $u_{t}$,
and $g_{\mathrm{\textnormal{ref.}}}(x)$ represents the RDF calculated from a
reference system. Since the IBI is only compatible with tabulated pair
potentials, it cannot be used to optimize a broader class of potentials such
as neural networks and pair potentials with explicit thermodynamic dependence,
which are used in this work. Practically, IBI also suffers from convergence
issues, as in some cases hundreds of iterations are required to converge.
Potestio (2013); Rosenberger, Hanke, and van der Vegt (2016) In our work, we
use IBI as a baseline comparison, and all CG potentials are obtained using the
Versatile Object-Oriented Toolkit for Coarse-Graining Applications (VOTCA).
Rühle  _et al._ (2009)

##  III Results

###  III.1 Learning Simple Pair Potentials

Figure 2: Transferability of DiffSim and IBI potentials fitted to reproduce
$g(r)$ of Lennard-Jones potentials at different densities. Each row contains
the results for the simulated g(r) at three densities (0.3, 0.5 and 0.7 from
left to right) with potential fitted at a particular density (light green
background, 0.3, 0.5 and 0.7 top to bottom) The last column reports the
learned potentials compared the ground truth. We include learned pair
potentials obtained from 10 independent runs.

As a first example, we apply our method to simple Lennard-Jones(LJ) systems at
different densities. We obtain target RDFs for LJ particles simulated at three
densities $\rho=(0.3,0.5,0.7)$ with temperature set at $1.0$, and we use the
obtained RDFs as the optimization target for our DiffSim protocol. We obtain
the target RDFs by simulating the systems for 500,000 steps after the initial
equilibration of 100,000 steps. We use a time step of 0.005. To test the
transferability of learned potentials at different system densities, we first
learn pair potentials at a single density and use the learned potential to
simulate the system with the other two densities that are not used for
learning. The learning results are presented in fig. 2, which also includes a
comparison with the potentials obtained using IBI. Our method shows an
accurate fitting of the RDFs for the dedicated densities, but the level of
transferability varies. The potential learned at the low density of $\rho=0.3$
demonstrates excellent transferability between different densities, and the
learned potential also agrees well with the ground-truth LJ potential. As the
density of the system used for fitting increases, the learned potentials show
high variability between replicate experiments, with overall poorer
transferability to systems that are not used for training. Many of the learned
potentials are dissimilar to the ground-truth potential, featuring concave-
down shapes. For a high-density system with $\rho=0.7$, neither IBI nor
DiffSim produces potentials that resemble the ground-truth LJ potential. The
learned potentials for $\rho=0.7$ also do not transfer well to $\rho=0.7,0.3$
with varied magnitudes of deviation.

To demonstrate that the observed trend is not specific to the LJ potential, we
tested our methods on the Modified Morse (MM) Potentials with different
softness at different densities Cheng and Yang (2007). The MM potential takes
the following form:

$MM(\rho,\psi)(r)=\frac{e^{2\rho(1-r^{\psi})}-2e^{\rho(1-r^{\psi}}-A)}{1+A}$
(14)

where $A=0$ for $\psi\leq 0$ and $A=e^{2\rho/\psi}-2e^{\rho/\psi}$ for
$\psi<0$. One can tune $(\rho,\psi)$ to generate a diverse set of pair
potentials. We tested our method for three different potentials
$MM(6.5,-0.45)$, $MM(5.5,0.44)$, and $MM(4.5,1.52)$ at four different
densities $\rho=(0.3,0.5,0.7,0.9)$ to infer pair potentials from target RDFs
in those systems. For quantitative evaluation, we calculate the average
deviation of learned potentials from the ground truth using the metric:

$\Delta=\int_{r_{0}}^{r_{c}}|U_{fit}(r)-U_{true}(r)|dr$ (15)

where we choose $r_{c}=2.5$ and $r_{0}=0.9$. fig. 4 shows that the learned
potentials deviate less from the ground truth as the density decreases, which
is consistent with the trend observed for LJ systems. This is also in line
with observations in the past literature.Moore, Iacovella, and McCabe (2014b)
For denser systems, we attribute the low sensitivities of the learned
potentials to dense local packing that downplays the effect of intermolecular
forces, largely due to the “entropic effect of packing". Chandler, Weeks, and
Andersen (1983); Wang, Stillinger, and Torquato (2020)

Figure 3: Comparison of fitted potentials with modified Morse potentials
MM(6.5, -0.45) (top), MM(5.5, 0.44) (middle), and MM(4.5, 1.52) (bottom). The
figure includes learned pair potentials obtained from 10 independent runs.
Figure 4: Left: Calculated mean deviation ($\Delta$) between the fitted
potentials and the ground-truth modified Morse potentials. For the three
potentials tested, the learned potentials manifest less deviation as density
decreases. Right: Modified Morse potentials used for learning.

###  III.2 Learning Pair Potentials for Binary Systems

Our DiffSim-based learning protocol is highly flexible. For example, it can be
modified to incorporate multiple learning targets to learn pair interactions
between different particle types. During training, this involves the
simultaneous optimization of different interaction functions between different
types of particles.

To demonstrate the efficacy of our method in this case, we construct a similar
example as presented in Rosenberger et al. Rosenberger and van der Vegt
(2019). The binary mixture system involves three types of LJ interactions with
the following parameters:
$\\{(\epsilon_{AA}:1.0,\sigma_{AA}:0.9),(\epsilon_{AB}:1.0,\sigma_{AB}:1.0),(\epsilon_{BB}:1.0,\sigma_{BB}:1.1)\\}$.
We simulate 500,000 steps at a temperature of $k_{B}T=1.0$ with a time step of
0.005 to generate the target RDFs for $x=(0.25,0.5,0.75)$, where $x$ is the
composition of particle A: $\frac{n_{A}}{n_{A}+n_{B}}$. We optimize based on
targets at $x=0.25$ and test the transferability of the learned potential at
$x=(0.5,0.75)$. During the simulation, we learn the three types of interaction
simultaneously by optimizing a joint loss of equally weighted loss terms.

$\displaystyle L$ $\displaystyle=L_{AA}+L_{AB}+L_{BB}$ (16)
$\displaystyle=\int_{r}(g_{\textnormal{sim.}}^{AA}(r)-g_{\textnormal{ref.}}^{AA}(r))^{2}dr+\int_{r}(g_{\textnormal{sim.}}^{AB}(r)-g_{\textnormal{ref.}}^{AB}(r))^{2}dr$
$\displaystyle+\int_{r}(g_{\textnormal{sim.}}^{BB}(r)-g_{\textnormal{ref.}}^{BB}(r))^{2}dr$

We optimize three potentials for 300 epochs using the Adam optimizer with a
learning rate of 0.001, with each epoch consisting of 250 MD steps. In fig. 5,
we show that DiffSim accurately reproduces the three target RDFs for $x=0.25$.
However, as the learned potential correctly recovers the repulsive part of the
target potential, the attractive part shapes are dissimilar to the ground
truth and feature a concave-down shape. As a transferability test, we apply
the learned potential at $x=0.25$ to systems with different compositions, that
is, $x=0.5$ and $x=0.75$, and observe poor transferability. Specifically, the
learned potentials overestimate the height of the first peak for the target
RDFs at $x=0.5$ and $x=0.75$. In this case, the learned potential overfits to
the composition state point on which it is trained.

For improved transferability, our model can also be configured to train
systems simultaneously in different states. This involves applying DiffSim to
simulate multiple systems and combining multiple gradient updates to the
potential arising from the different trajectories. In practice, we initialize
all three systems and simulate them in parallel with the same shared
potentials. Gradient updates are collected through eq. 2 and eq. 9 from each
trajectory, and the pair potentials of AA, AB, and BB are updated with the
combined three gradients from the three state points. For each training epoch,
three systems are simulated with gradients accumulated on the parameters of
the pair potential with the total loss as the sum of individual losses
evaluated for each system, i.e., $L_{tot}=L_{x=0.25}+L_{x=0.5}+L_{x=0.75}$.
Each loss takes the form described in eq. 16, giving us 9 optimization targets
in total, namely the RDF of AA, AB, and BB at the three concentrations. Our
training results in fig. 6 show that potentials can be learned to
simultaneously fit target RDFs at multiple systems and state points at the
same time. Compared to potentials solely fitted at $x=0.25$ (fig. 5), the
learned potentials in all three states show a better resemblance to ground-
truth LJ potentials for both repulsive and attractive portions.

Figure 5: We apply the pair potentials learned at $x=0.25$ (highlighted with
light green background) to systems with $x=(0.5,0.75)$. Unlike ground truth
potentials of the LJ type, the learned potentials feature concave-down shapes.
Figure 6: RDFs and pair potentials obtained from simultaneous learning of pair
potentials for systems at $x=0.25,x=0.5,x=0.75$. The learned potentials
resemble the ground truth potentials. The learned potentials accurately
recover the RDFs at all compositions.

###  III.3 Learning Transferable Coarse-Grained Potentials

As the third example, we show that DiffSim can learn pair potentials for CG
simulations of water with a temperature-transferable functional form. For the
ground truth simulations, we simulate all-atom water molecules at $1g/cm^{3}$
at 288K, 338K and 388K with the SPC/E force field.Berendsen, Grigera, and
Straatsma (1987) We choose to map each water molecule onto the oxygen atom and
use the oxygen-oxygen RDF as our fitting target. We first run a single-state
fitting experiment at single temperatures and use the learned potentials to
simulate systems at other temperatures. For each system, we train the
potential using DiffSim for 500 epochs, with each epoch including 180
simulation steps with a time step of 2fs. For comparison, we also performed
the IBI at all three temperatures. fig. 7 shows the fitting and
transferability results. As DiffSim and IBI manage to accurately fit single
states, they manifest large RDF errors when using similar systems at other
temperatures. The learned potentials using IBI and DiffSim show similar
shapes, with an inner potential well at $2.7\textup{\AA}$, approximately where
the RDF of water has the first peak. Furthermore, the depth decreases with
temperature, showing a clear temperature dependence of the pair potential.
This observation agrees with the water potentials obtained in the past
literature Johnson, Head-Gordon, and Louis (2007); Pretti and Shell (2021),
suggesting that an explicit temperature dependence is required for the pair
potential to be transferable.

To incorporate the temperature-dependent nature of CG potentials, we made our
neural pair potentials temperature-dependent for improved transferability
across temperatures. With AD, our method enables more flexible forms of pair
potentials. Specifically, we optimize pair potentials with explicit
temperature dependence. Inspired by previous work studying the energy and
entropy decomposition of CG potentials,Kidder, Szukalo, and Noid (2021);
Pretti and Shell (2021) we propose the following form of neural pair
potential.

$u_{T}(r,T)=u_{1}(r)-k_{B}Tu_{2}(r)$ (17)

where $u_{1}$ and $u_{2}$ are two distinct neural pair potentials. $u_{1}(r)$
can be understood as the energy contribution to pair potentials, and
$u_{2}(r)$ is the entropy contribution. For training, we initialize three
coarse-grained water systems and simulate them at different temperatures. We
simulate and train the three systems simultaneously to optimize $u_{1}(r)$ and
$u_{2}(r)$ jointly. In fig. 7, we show the RDFs obtained and the learned pair
potentials at the three different temperatures used for learning. The learned
potentials at different temperatures show similar shapes. All three potentials
have a characteristic barrier around 3.2 $\AA$, and the change in temperature
affects their height. Compared to the potential obtained by fitting on single
states, the learned temperature-dependent pair potential shows excellent
transferability to all three temperatures used for optimization. To quantify
the deviation between the simulated RDF and the target RDF, we use the
following metric proposed by Pretti. et al.:

$D=4\pi\rho\int_{0}^{r_{c}}r^{2}(g_{(}r)-g_{\textnormal{ref.}}(r))^{2}dr$ (18)

In fig. 9, we compare the deviation of the simulated RDF obtained from
different training protocols. The result shows that fitting in a single state
can produce an accurate RDF in the state used for fitting, but it does not
transfer well to other states. The temperature-transferable potential produces
overall lower RDF deviations for all three states.

Figure 7: Learning coarse-grained water potential at single state points and
transfer to other state points. The results include learned pair potentials
and simulated RDFs from 20 independent runs with DiffSim. The learned
potentials using IBI is also included.  Figure 8: Simulated RDF and learned
transferable potential at 288K, 338K, 388K. The learned potential shows
temperature-dependent barrier height, with the height increasing with lowering
in temperature. Figure 9: Quantified RDF error ($D$) of learned potentials
simulated at different state points (288K, 338K, 388K).

##  IV Conclusion

In this work, we propose a flexible method based on DiffSim to directly learn
pair potentials. We demonstrate our methods in several computational examples.
First, we apply our method to learning simple pair potentials where the ground
truth is known. We show that our method recovers a diverse set of possible
pair potentials given a target RDF and, therefore, can be used to probe
sensitivities for fitting pair potentials. For a dense system, our results are
in line with previous reports, indicating that fitted potentials can be highly
insensitive. Wang, Stillinger, and Torquato (2020); Potestio (2013) In
contrast, the fitting for the dilute system shows a better sensitivity, with
the obtained potentials in reasonable agreement with the ground truth
potential. We further demonstrate the multistate fitting capabilities of the
method by applying it to binary mixtures that involve learning of multiple
pair interactions under multiple composition conditions. Learning involves the
incorporation of multiple loss functions in different RDF targets. Our
experiment shows that simultaneous learning on systems with different
compositions is possible and produces potentials that have better
transferability. Our method is also applied to learning the CG pair
potentials. We demonstrate the learning of water CG potentials at three
different state points using a single temperature-dependent CG potential that
is more transferable than a potential obtained from fitting a single state
point.

The proposed DiffSim protocol can be understood as a gradient-based extension
of the iterative Boltzmann inversion. Our method directly optimizes potentials
with the gradient signal obtained from simulations. It allows simultaneous
simulation and fitting for systems at different temperatures and compositions,
allowing direct parameterization of pair potentials to improve
transferability. The learning pipeline can be flexibly set up with the
incorporation of multiple optimization objectives and allows for the
optimization of a broader range of pair potentials, such as pair potentials
with explicit thermodynamic dependence. Our method opens up broader
possibilities for learning-based simulations: objectives can be extended to
include more observables such as stress tensors and velocity autocorrelation
functions to enable gradient-based top-down parameterization of force fields.

DiffSim learns potentials by directly optimizing unrolled simulation
operations that involve updates of positions, velocities, and forces. The
procedure directly learns to match an observation without atomistic data,
which requires extensive sampling that often requires active learning.
Compared to force-based optimization of machine learning potentials, our
method produces direct simulation feedback and guarantees the stable
production of desired observables. In comparison, relying on matching
atomistic forces, highly flexible machine learning force fields often suffer
from unstable dynamics, which is common even for models showing very low force
errors.Stocker  _et al._ (2022) In comparison, pair potentials constrain
interactions between pairs of particles and are therefore easier to interpret
and more efficient to simulate than machine learning potentials. Although
observable-informed and gradient-based, our learning protocol can also be
combined with force matching to incorporate top-down and bottom-up
parameterization into optimization.

In this work, we present a flexible machine learning method that uses
differentiable simulations to directly parameterize pair potentials from
radial distribution functions. We anticipate that the proposed framework
provides improved modeling flexibility for multiscale simulation of molecular
liquids.

## Acknowledgment

This work is supported by Toyota Research Institute.

## Author Declaration

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

Wujie Wang: Conceptualization (lead); Investigation (equal); Writing –
original draft (equal); Writing – review & editing (equal). Zhenghao Wu:
Conceptualization (supporting); Investigation (equal); Writing – original
draft (equal); Writing – review & editing (equal). Rafael Gómez-Bombarelli:
Conceptualization (supporting); Supervision; Funding acquisition; Writing –
review & editing (equal).

## Data Availability

The data and code that support the findings of this study are available in
https://github.com/torchmd/mdgrad.

##  V Appendix

###  V.1 Model and training hyper-parameters

Here, we detail the hyperparameters used in our experiments. For all training
runs, a learning rate scheduler is used to dynamically decrease the learning
when the pre-defined convergence criterion is met. The learning rate is
expected to decrease by half when the improvement in training loss is less
than $0.1\%$.

Lennard Jones and MM systems. The model is trained for 500 epochs, with each
epoch comprising 120 simulation steps. The learning rate is set at $0.002$.
Neural pair potentials are constructed with Gaussian-smeared distance inputs
of 100 bins ($K=100$) with $d=0.1$. The smeared inputs are fed into a
multilayer perceptron of 3 hidden layers with a hidden layer size of 128. The
activation function that we use is ELU.Clevert, Unterthiner, and Hochreiter
(2015) The prior potential $u_{prior}(r)$ that we use takes the form
$\epsilon(\frac{r}{\sigma})^{10}$ with $\sigma=0.9$ and $\epsilon=0.4$. During
simulations, the neighbor list is updated for every integration step.

Binary mixtures. The model is trained for 500 epochs, with each epoch
comprising 120 simulation steps. The learning rate is set at $0.0003$. Three
neural pair potentials are initialized to represent the three types of pair
interaction involved in the system. For the computation of RDF loss, each pair
statistics is equally weighted. Each pair potential uses 100 bins ($K=100$) to
smear all distance inputs with $d=0.15$. The parameterization MLP consists of
3 hidden layers with each layer of 128 neurons. The activation function used
is SELU. Klambauer  _et al._ (2017) The prior potential $u_{prior}(r)$ we use
the form $\epsilon(\frac{r}{\sigma})^{6}$ with $\sigma=1.0$ and
$\epsilon=2.0$. Before training, a pre-training step is performed to
initialize neural network potentials with shapes similar to the Boltzmann
inverted potentials derived from the target RDF. This step is used to
facilitate faster convergence of training and sampling of relevant parts of
the configuration space.

Temperature-transferable water simulations.  The model is trained for 500
epochs, with each epoch comprising 190 simulation steps. The learning rate is
set at $0.00065$. We similarly apply 3-layer neural networks to represent
$u_{1}$ and $u_{2}$ with hidden layers of 115 neurons. The pair distances are
smeared with $K=120$ equally spaced Gaussian kernels with $d=0.15$. The
activation function used is ELU. Clevert, Unterthiner, and Hochreiter (2015).
$u_{prior}$ is used with $\sigma=40$ kcal/mol and $\epsilon=1.7\textup{\AA}$.
A pre-training procedure is also performed using the Boltzmann inverted RDF at
$298K$ as the target.

###  V.2 Supplementary results

We provide additional results in fig. 10 to show the RDFs generated by learned
pair potentials for different MM systems at different densities. fig. 10
provides supplmentary information for fig. 3. The simulations are carried out
for 100,000 steps with each step being 0.005.

Figure 10: Comparison between ground-truth RDFs and RDFs from simulations
using learned pair potentials for different MM potentials at different
densities. From top to bottom rows, the figures correspond to RDFs generated
by learned potentials for MM(6.5, -0.45) , MM(5.5, 0.44), and MM(4.5, 1.52)
systems respectively.

###  V.3 Molecular dynamics simulations of water

The FG model of water used in this work is an all-atom water model described
by SPC/E force fieldBerendsen, Grigera, and Straatsma (1987). The large-scale
atomic / molecular massively parallel simulator (LAMMPS)Plimpton (1995);
Thompson  _et al._ (2022) is used to perform MD simulations of systems with
$N=1500$ water molecules at constant density $\rho=1.0$ $\mathrm{g/cm^{3}}$
and at temperatures of $T=288K$, $T=338K$ and $T=388K$, with a time unit
$\delta t=1$ fs. We apply periodic boundary conditions in the x, y and z
directions, and the Nośe-Hoover thermostatNosé (1984); Hoover (1985) with a
damping time of $\tau=100$ fs to keep the temperature during the simulation.
The water systems are first equilibrated for $10$ ns. Subsequently, it was
subjected to a simulation run of $1$ ns to collect data every 0.1 ps to
calculate the reference radial distribution functions of oxygen atoms with bin
size 0.001 nm. In the CG representation, a water molecule is replaced by a
single bead placed on the oxygen atom. CG simulations are also performed in
LAMMPS. Plimpton (1995); Thompson  _et al._ (2022)

###  V.4 Protocol for Iterative Boltzmann Inversion

Lennard-Jones systems. The RDFs of the MD simulations with LJ potentials are
the target distribution fed to the IBI machinery. The cutoff radius used for
IBI optimization and simulations with IBI refined potential is $r_{cut}=2.0$
with space size $\delta r=0.01$ $\sigma$. For each state point, the MD
simulation performed under the canonical ensemble in each IBI iteration
consisted of $5\times 10^{5}$ equilibration steps and $5\times 10^{5}$
sampling steps for the calculation of the RDF. The entire IBI optimization
takes between 10 and 20 iterations to yield converged distributions that match
the target RDFs.

Water simulations. The RDFs from the SPC/E water simulation are obtained by
mapping the CG beads mapped to the oxygen atoms and are fed into the IBI
optimization protocol. The cutoff radii used for IBI optimization and CG
simulations are $r_{cut}=0.625$ nm with a space size $\delta r=0.0005$ nm. For
each state point, the CG-MD simulation performed in the canonical ensemble in
each IBI iteration consisted of $2\times 10^{6}$ equilibration steps and
$1\times 10^{6}$ sampling steps for the calculation of the RDF. The entire IBI
optimization takes between 10 and 20 iterations to yield converged
distributions that match the target RDFs.

## References

  * Vlachakis  _et al._ (2014) D. Vlachakis, E. Bencurova, N. Papangelopoulos, and S. Kossida, “Current state-of-the-art molecular dynamics methods and applications,” in  _Adv. Protein Chem. Struct. Biol._, Vol. 94 (Elsevier, 2014) pp. 269–313.
  * Blumberger (2015) J. Blumberger, “Recent advances in the theory and molecular simulation of biological electron transfer reactions,” Chem. Rev. 115, 11191–11238 (2015).
  * Gartner and Jayaraman (2019) T. E. Gartner and A. Jayaraman, “Modeling and Simulations of Polymers: A Roadmap,” Macromolecules 52, 755–786 (2019).
  * Fish, Wagner, and Keten (2021) J. Fish, G. J. Wagner, and S. Keten, “Mesoscopic and multiscale modelling in materials,” Nat. Mater. 20, 774–786 (2021).
  * Mayo, Olafson, and Goddard (1990) S. L. Mayo, B. D. Olafson, and W. A. Goddard, “DREIDING: a generic force field for molecular simulations,” J. Phys. Chem. 94, 8897–8909 (1990).
  * Jorgensen, Maxwell, and Tirado-Rives (1996) W. L. Jorgensen, D. S. Maxwell, and J. Tirado-Rives, “Development and Testing of the OPLS All-Atom Force Field on Conformational Energetics and Properties of Organic Liquids,” J. Am. Chem. Soc. 118, 11225–11236 (1996).
  * Pastor and MacKerell (2011) R. W. Pastor and A. D. MacKerell, “Development of the CHARMM Force Field for Lipids,” J. Phys. Chem. Lett. 2, 1526–1532 (2011).
  * Rapaport (2004) D. C. Rapaport,  _The art of molecular dynamics simulation_ (Cambridge University Press, Cambridge, UK; New York, NY, 2004) oCLC: 928698870.
  * Cisneros  _et al._ (2016) G. A. Cisneros, K. T. Wikfeldt, L. Ojamäe, J. Lu, Y. Xu, H. Torabifard, A. P. Bartók, G. Csányi, V. Molinero, and F. Paesani, “Modeling Molecular Interactions in Water: From Pairwise to Many-Body Potential Energy Functions,” Chem. Rev. 116, 7501–7528 (2016).
  * Berendsen, Grigera, and Straatsma (1987) H. J. C. Berendsen, J. R. Grigera, and T. P. Straatsma, “The missing term in effective pair potentials,” J. Phys. Chem. 91, 6269–6271 (1987).
  * Jorgensen  _et al._ (1983) W. L. Jorgensen, J. Chandrasekhar, J. D. Madura, R. W. Impey, and M. L. Klein, “Comparison of simple potential functions for simulating liquid water,” J. Chem. Phys. 79, 926–935 (1983).
  * Glaser  _et al._ (2015) J. Glaser, T. D. Nguyen, J. A. Anderson, P. Lui, F. Spiga, J. A. Millan, D. C. Morse, and S. C. Glotzer, “Strong scaling of general-purpose molecular dynamics simulations on GPUs,” Comput. Phys. Commun. 192, 97–107 (2015).
  * Albaugh  _et al._ (2016) A. Albaugh, H. A. Boateng, R. T. Bradshaw, O. N. Demerdash, J. Dziedzic, Y. Mao, D. T. Margul, J. Swails, Q. Zeng, D. A. Case, P. Eastman, L.-P. Wang, J. W. Essex, M. Head-Gordon, V. S. Pande, J. W. Ponder, Y. Shao, C.-K. Skylaris, I. T. Todorov, M. E. Tuckerman, and T. Head-Gordon, “Advanced Potential Energy Surfaces for Molecular Simulation,” J. Phys. Chem. B 120, 9811–9832 (2016).
  * Plimpton and Thompson (2012) S. J. Plimpton and A. P. Thompson, “Computational aspects of many-body potentials,” MRS Bull. 37, 513–521 (2012).
  * Müller-Plathe (2002) F. Müller-Plathe, “Coarse-Graining in Polymer Simulation: From the Atomistic to the Mesoscopic Scale and Back,” ChemPhysChem 3, 754–769 (2002).
  * Souza  _et al._ (2021) P. C. T. Souza, R. Alessandri, J. Barnoud, S. Thallmair, I. Faustino, F. Grünewald, I. Patmanidis, H. Abdizadeh, B. M. H. Bruininks, T. A. Wassenaar, P. C. Kroon, J. Melcr, V. Nieto, V. Corradi, H. M. Khan, J. Domański, M. Javanainen, H. Martinez-Seara, N. Reuter, R. B. Best, I. Vattulainen, L. Monticelli, X. Periole, D. P. Tieleman, A. H. de Vries, and S. J. Marrink, “Martini 3: a general purpose force field for coarse-grained molecular dynamics,” Nat. Methods 18, 382–388 (2021).
  * Dhamankar and Webb (2021) S. Dhamankar and M. A. Webb, “Chemically specific coarse-graining of polymers: Methods and prospects,” J. Polym. Sci. 59, 2613–2643 (2021).
  * Ercolessi and Adams (1994) F. Ercolessi and J. B. Adams, “Interatomic Potentials from First-Principles Calculations: The Force-Matching Method,” Europhys. Lett. 26, 583–588 (1994).
  * Izvekov and Voth (2005) S. Izvekov and G. A. Voth, “A Multiscale Coarse-Graining Method for Biomolecular Systems,” J. Phys. Chem. B 109, 2469–2473 (2005).
  * Noid  _et al._ (2007) W. G. Noid, J.-W. Chu, G. S. Ayton, and G. A. Voth, “Multiscale Coarse-Graining and Structural Correlations: Connections to Liquid-State Theory,” J. Phys. Chem. B 111, 4116–4127 (2007).
  * Zhang  _et al._ (2018) L. Zhang, J. Han, H. Wang, R. Car, and W. E, “Deepcg: Constructing coarse-grained models via deep neural networks,” The Journal of chemical physics 149, 034101 (2018).
  * Wang  _et al._ (2019) J. Wang, S. Olsson, C. Wehmeyer, A. Pérez, N. E. Charron, G. de Fabritiis, F. Noé, and C. Clementi, “Machine Learning of Coarse-Grained Molecular Dynamics Force Fields,” ACS Cent. Sci. , acscentsci.8b00913 (2019).
  * Rudzinski and Noid (2012) J. F. Rudzinski and W. G. Noid, “The role of many-body correlations in determining potentials for coarse-grained models of equilibrium structure,” The Journal of Physical Chemistry B 116, 8621–8635 (2012).
  * Stocker  _et al._ (2022) S. Stocker, J. Gasteiger, F. Becker, S. Günnemann, and J. Margraf, “How robust are modern graph neural network potentials in long and hot molecular dynamics simulations?”  (2022), 10.26434/chemrxiv-2022-mc4gb.
  * Chaimovich and Shell (2011) A. Chaimovich and M. S. Shell, “Coarse-graining errors and numerical optimization using a relative entropy framework,” J. Chem. Phys. 134, 094112 (2011).
  * Kocer, Ko, and Behler (2022) E. Kocer, T. W. Ko, and J. Behler, “Neural Network Potentials: A Concise Overview of Methods,” Annu. Rev. Phys. Chem. 73, annurev–physchem–082720–034254 (2022).
  * Reith, Pütz, and Müller-Plathe (2003) D. Reith, M. Pütz, and F. Müller-Plathe, “Deriving effective mesoscale potentials from atomistic simulations: Mesoscale Potentials from Atomistic Simulations,” J. Comput. Chem. 24, 1624–1636 (2003).
  * Qian  _et al._ (2008) H.-J. Qian, P. Carbone, X. Chen, H. A. Karimi-Varzaneh, C. C. Liew, and F. Müller-Plathe, “Temperature-Transferable Coarse-Grained Potentials for Ethylbenzene, Polystyrene, and Their Mixtures,” Macromolecules 41, 9919–9929 (2008).
  * Karimi-Varzaneh  _et al._ (2010) H. A. Karimi-Varzaneh, F. Müller-Plathe, S. Balasubramanian, and P. Carbone, “Studying long-time dynamics of imidazolium-based ionic liquids with a systematically coarse-grained model,” Phys. Chem. Chem. Phys. 12, 4714 (2010).
  * Mullinax and Noid (2009) J. W. Mullinax and W. G. Noid, “Extended ensemble approach for deriving transferable coarse-grained potentials,” J. Chem. Phys. 131, 104110 (2009).
  * Moore, Iacovella, and McCabe (2014a) T. C. Moore, C. R. Iacovella, and C. McCabe, “Derivation of coarse-grained potentials via multistate iterative Boltzmann inversion,” J. Chem. Phys. 140, 224104 (2014a).
  * Liao  _et al._ (2019) H.-J. Liao, J.-G. Liu, L. Wang, and T. Xiang, “Differentiable programming tensor networks,” Phys. Rev. X. 9, 031041 (2019).
  * Schenck and Fox (2018) C. Schenck and D. Fox, “Spnets: Differentiable fluid dynamics for deep neural networks,” in  _Conference on Robot Learning_ (PMLR, 2018) pp. 317–335.
  * Unke  _et al._ (2021) O. T. Unke, S. Chmiela, H. E. Sauceda, M. Gastegger, I. Poltavsky, K. T. Schutt, A. Tkatchenko, and K.-R. Muller, “Machine learning force fields,” Chem. Rev. 121, 10142–10186 (2021).
  * Kasim and Vinko (2021) M. F. Kasim and S. M. Vinko, “Learning the exchange-correlation functional from nature with fully differentiable density functional theory,” Phys. Rev. Lett. 127, 126403 (2021).
  * Kanungo, Zimmerman, and Gavini (2019) B. Kanungo, P. M. Zimmerman, and V. Gavini, “Exact exchange-correlation potentials from ground-state electron densities,” Nat. Commun. 10, 1–9 (2019).
  * Tamayo-Mendoza  _et al._ (2018) T. Tamayo-Mendoza, C. Kreisbeck, R. Lindh, and A. Aspuru-Guzik, “Automatic differentiation in quantum chemistry with applications to fully variational hartree–fock,” ACS Cent. Sci. 4, 559–566 (2018).
  * Blondel  _et al._ (2021) M. Blondel, Q. Berthet, M. Cuturi, R. Frostig, S. Hoyer, F. Llinares-López, F. Pedregosa, and J.-P. Vert, “Efficient and modular implicit differentiation,” arXiv preprint arXiv:2105.15183 (2021).
  * Magnus (1985) J. R. Magnus, “On differentiating eigenvalues and eigenvectors,” Econometric theory 1, 179–191 (1985).
  * Chen  _et al._ (2018) R. T. Chen, Y. Rubanova, J. Bettencourt, and D. K. Duvenaud, “Neural ordinary differential equations,” Adv. Neural Inf. Process. Syst. 31 (2018).
  * Wang, Axelrod, and Gómez-Bombarelli (2020) W. Wang, S. Axelrod, and R. Gómez-Bombarelli, “Differentiable molecular simulations for control and learning,” arXiv preprint arXiv:2003.00868 (2020).
  * Schoenholz and Cubuk (2020) S. Schoenholz and E. D. Cubuk, “Jax md: a framework for differentiable physics,” Adv. Neural Inf. Process. Syst. 33, 11428–11441 (2020).
  * Doerr  _et al._ (2021) S. Doerr, M. Majewski, A. Pérez, A. Kramer, C. Clementi, F. Noe, T. Giorgino, and G. De Fabritiis, “Torchmd: A deep learning framework for molecular simulations,” J. Chem. Theory Comput. 17, 2355–2363 (2021).
  * Thaler and Zavadlav (2021) S. Thaler and J. Zavadlav, “Learning neural network potentials from experimental data via differentiable trajectory reweighting,” Nat. Commun. 12, 1–10 (2021).
  * Greener and Jones (2021) J. G. Greener and D. T. Jones, “Differentiable molecular simulation can learn all the parameters in a coarse-grained force field for proteins,” PloS one 16, e0256990 (2021).
  * Moradzadeh and Aluru (2019) A. Moradzadeh and N. R. Aluru, “Transfer-learning-based coarse-graining method for simple fluids: Toward deep inverse liquid-state theory,” The journal of physical chemistry letters 10, 1242–1250 (2019).
  * Battaglia  _et al._ (2018) P. W. Battaglia, J. B. Hamrick, V. Bapst, A. Sanchez-Gonzalez, V. Zambaldi, M. Malinowski, A. Tacchetti, D. Raposo, A. Santoro, R. Faulkner, _et al._ , “Relational inductive biases, deep learning, and graph networks,” arXiv preprint arXiv:1806.01261 (2018).
  * Baxter (2000) J. Baxter, “A model of inductive bias learning,” J. Artif. Intell. Res. 12, 149–198 (2000).
  * Han, Jentzen, and Weinan (2018) J. Han, A. Jentzen, and E. Weinan, “Solving high-dimensional partial differential equations using deep learning,” Proc. Natl. Acad. Sci. 115, 8505–8510 (2018).
  * Chen (2017) Y.-C. Chen, “A tutorial on kernel density estimation and recent advances,” Biostatistics & Epidemiology 1, 161–187 (2017).
  * Brunton and Kutz (2022) S. L. Brunton and J. N. Kutz,  _Data-driven science and engineering: Machine learning, dynamical systems, and control_ (Cambridge University Press, 2022).
  * Abadi  _et al._ (2016) M. Abadi, P. Barham, J. Chen, Z. Chen, A. Davis, J. Dean, M. Devin, S. Ghemawat, G. Irving, M. Isard, _et al._ , “$\\{$TensorFlow$\\}$: A system for $\\{$Large-Scale$\\}$ machine learning,” in  _12th USENIX symposium on operating systems design and implementation (OSDI 16)_ (2016) pp. 265–283.
  * Frostig, Johnson, and Leary (2018) R. Frostig, M. J. Johnson, and C. Leary, “Compiling machine learning programs via high-level tracing,” Systems for Machine Learning , 23–24 (2018).
  * Paszke  _et al._ (2019) A. Paszke, S. Gross, F. Massa, A. Lerer, J. Bradbury, G. Chanan, T. Killeen, Z. Lin, N. Gimelshein, L. Antiga, _et al._ , “Pytorch: An imperative style, high-performance deep learning library,” Adv. Neural Inf. Process. Syst. 32 (2019).
  * Ustinova and Lempitsky (2016) E. Ustinova and V. Lempitsky, “Learning deep embeddings with histogram loss,” Adv. Neural Inf. Process. Syst. 29 (2016).
  * Hornik, Stinchcombe, and White (1989) K. Hornik, M. Stinchcombe, and H. White, “Multilayer feedforward networks are universal approximators,” Neural networks 2, 359–366 (1989).
  * Tschöp  _et al._ (1998) W. Tschöp, K. Kremer, J. Batoulis, T. Bürger, and O. Hahn, “Simulation of polymer melts. i. coarse-graining procedure for polycarbonates,” Acta Polymerica 49, 61–74 (1998).
  * Ruiz, Schulter, and Chandraker (2018) N. Ruiz, S. Schulter, and M. Chandraker, “Learning to simulate,” arXiv preprint arXiv:1810.02513 (2018).
  * Martyna, Klein, and Tuckerman (1992) G. J. Martyna, M. L. Klein, and M. Tuckerman, “Nosé–hoover chains: The canonical ensemble via continuous dynamics,” J. Chem. Phys. 97, 2635–2643 (1992).
  * Kingma and Ba (2014) D. P. Kingma and J. Ba, “Adam: A method for stochastic optimization,” arXiv preprint arXiv:1412.6980 (2014).
  * Ingraham  _et al._ (2018) J. Ingraham, A. Riesselman, C. Sander, and D. Marks, “Learning protein structure with a differentiable simulator,” in  _International Conference on Learning Representations_ (2018).
  * Metz  _et al._ (2021) L. Metz, C. D. Freeman, S. S. Schoenholz, and T. Kachman, “Gradients are not all you need,” arXiv preprint arXiv:2111.05803 (2021).
  * Potestio (2013) R. Potestio, “Is henderson’s theorem practically useful?” JUnQ 3, 13–15 (2013).
  * Rosenberger, Hanke, and van der Vegt (2016) D. Rosenberger, M. Hanke, and N. F. van der Vegt, “Comparison of iterative inverse coarse-graining methods,” Eur. Phys. J. Spec. Top. 225, 1323–1345 (2016).
  * Rühle  _et al._ (2009) V. Rühle, C. Junghans, A. Lukyanov, K. Kremer, and D. Andrienko, “Versatile Object-Oriented Toolkit for Coarse-Graining Applications,” J. Chem. Theory Comput. 5, 3211–3223 (2009).
  * Cheng and Yang (2007) L. Cheng and J. Yang, “Modified morse potential for unification of the pair interactions,” J. Chem. Phys. 127, 124104 (2007).
  * Moore, Iacovella, and McCabe (2014b) T. C. Moore, C. R. Iacovella, and C. McCabe, “Derivation of coarse-grained potentials via multistate iterative boltzmann inversion,” J. Chem. Phys. 140, 06B606_1 (2014b).
  * Chandler, Weeks, and Andersen (1983) D. Chandler, J. D. Weeks, and H. C. Andersen, “Van der waals picture of liquids, solids, and phase transformations,” Science 220, 787–794 (1983).
  * Wang, Stillinger, and Torquato (2020) H. Wang, F. H. Stillinger, and S. Torquato, “Sensitivity of pair statistics on pair potentials in many-body systems,” J. Chem. Phys. 153, 124106 (2020).
  * Rosenberger and van der Vegt (2019) D. Rosenberger and N. F. van der Vegt, “Relative entropy indicates an ideal concentration for structure-based coarse graining of binary mixtures,” Phys. Rev. E 99, 053308 (2019).
  * Johnson, Head-Gordon, and Louis (2007) M. E. Johnson, T. Head-Gordon, and A. A. Louis, “Representability problems for coarse-grained water potentials,” J. Chem. Phys. 126, 144509 (2007).
  * Pretti and Shell (2021) E. Pretti and M. S. Shell, “A microcanonical approach to temperature-transferable coarse-grained models using the relative entropy,” J. Chem. Phys. 155, 094102 (2021).
  * Kidder, Szukalo, and Noid (2021) K. M. Kidder, R. J. Szukalo, and W. Noid, “Energetic and entropic considerations for coarse-graining,” Eur. Phys. J. B 94, 1–29 (2021).
  * Clevert, Unterthiner, and Hochreiter (2015) D.-A. Clevert, T. Unterthiner, and S. Hochreiter, “Fast and accurate deep network learning by exponential linear units (elus),” arXiv preprint arXiv:1511.07289 (2015).
  * Klambauer  _et al._ (2017) G. Klambauer, T. Unterthiner, A. Mayr, and S. Hochreiter, “Self-normalizing neural networks,” Advances in neural information processing systems 30 (2017).
  * Plimpton (1995) S. Plimpton, “Fast Parallel Algorithms for Short-Range Molecular Dynamics,” J. Comput. Phys. 117, 1–19 (1995).
  * Thompson  _et al._ (2022) A. P. Thompson, H. M. Aktulga, R. Berger, D. S. Bolintineanu, W. M. Brown, P. S. Crozier, P. J. in’t Veld, A. Kohlmeyer, S. G. Moore, T. D. Nguyen, _et al._ , “Lammps-a flexible simulation tool for particle-based materials modeling at the atomic, meso, and continuum scales,” Computer Physics Communications 271, 108171 (2022).
  * Nosé (1984) S. Nosé, “A unified formulation of the constant temperature molecular dynamics methods,” J. Chem. Phys. 81, 511–519 (1984).
  * Hoover (1985) W. G. Hoover, “Canonical dynamics: Equilibrium phase-space distributions,” Phys. Rev. A 31, 1695–1697 (1985).