<EMAIL_ADDRESS>
# Bayesian inference-driven model parameterization and model selection for
2CLJQ fluid models

Owen C. Madin  Department of Chemical & Biological Engineering, University of
Colorado Boulder, Boulder, CO 80309  Simon Boothroyd  Boothroyd Scientific
Consulting Ltd., 71-75 Shelton Street, London, Greater London, United Kingdom,
WC2H 9JQ  Richard A. Messerly  Theoretical Division, Los Alamos National
Laboratory, Los Alamos, NM 87545  John D. Chodera  Computational & Systems
Biology Program, Sloan Kettering Institute, Memorial Sloan Kettering Cancer
Center, New York, NY 10065  Josh Fass  Tri-Institutional PhD Program in
Computational Biology and Medicine, Weill Cornell Graduate School of Medical
Sciences, New York, NY 10065, USA. Current address: Relay Therapeutics,
Cambridge, MA 02139, USA  Michael R. Shirts  Department of Chemical &
Biological Engineering, University of Colorado Boulder, Boulder, CO 80309

###### Abstract

A high level of physical detail in a molecular model improves its ability to
perform high accuracy simulations, but can also significantly affect its
complexity and computational cost. In some situations, it is worthwhile to add
additional complexity to a model to capture properties of interest; in others,
additional complexity is unnecessary and can make simulations computationally
infeasible. In this work we demonstrate the use of Bayes factors for molecular
model selection, using Monte Carlo sampling techniques to evaluate the
evidence for different levels of complexity in the two-centered Lennard-Jones
+ quadrupole (2CLJQ) fluid model. Examining three levels of nested model
complexity, we demonstrate that the use of variable quadrupole and bond length
parameters in this model framework is justified only sometimes. We also
explore the effect of the Bayesian prior distribution on the Bayes factors, as
well as ways to propose meaningful prior distributions. This Bayesian Markov
Chain Monte Carlo (MCMC) process is enabled by the use of analytical surrogate
models that accurately approximate the physical properties of interest. This
work paves the way for further atomistic model selection work via Bayesian
inference and surrogate modeling.

##  1 Introduction

Parameterization of molecular force fields is a long-standing problem for the
classical molecular simulation community [1, 2, 3], as the choice of models
and their respective parameters are critical for achieving quantitative
accuracy in modeling of molecular interactions. The treatment of electrostatic
and dispersion-repulsion interactions, commonly referred to as _non-bonded
interactions_ , is a difficult parameterization task because these
interactions are described by relatively simple models that are crude
approximations of the underlying electronic behavior, requiring mapping
electronic probability distributions onto pairwise interactions between a
small number of points. Dispersion-repulsion interactions, commonly modeled
with a Lennard-Jones (LJ) potential are especially challenging to derive [4]
as values of the Lennard-Jones parameters are not directly extractable from
quantum mechanics (QM) calculations. Most force fields train these parameters
against macroscopic condensed phase physical properties, though there are some
attempts to obtain them with reference to QM calculations [5, 6].

Although most commonly employed classical force fields use electrostatic point
charges and pairwise LJ potentials [7, 8], other models for these interactions
exist. Atomic multipoles [9, 10], Drude oscillators [11], and fluctuating
charges [12] for electrostatic interactions and alternate functional forms [1,
13, 14, 15] for dispersion-repulsion interactions have been proposed over the
years. Even considering only LJ potentials, many different choices for atom
typing [16, 17] and combination rules [18, 19] exist. Choosing a potential
energy function with a more complex functional form can improve agreement with
experiment, however, increasing the complexity may also make the
parameterization more difficult due to an increased likelihood of overfitting
[20, 21].

This creates two distinct problems in designing force fields: (1) selecting
among discrete choices of models and (2) optimizing continuous choice of
parameters within models. It is common, although often expensive, to compare
the quantitative performance of different parameter sets within a single
model, but it is difficult to quantitatively compare fitness between models.

A common statistical formalism across for comparing the overall fitness of
discrete models is Bayesian inference [22], which allows users to evaluate
models in a way that automatically penalizes unnecessary complexity [23]. By
integrating over the marginal distribution of the parameters one can calculate
_Bayes factors_ which are interpreted as odds ratios between the separate
models [24]. The Bayesian framework combines specific information on a model’s
ability to reproduce target data, with more general prior knowledge based on
data collected with previous parameter sets, physical constraints or
invariances, or chemical intuition. In this way, the prior distribution
generalizes and influences the comparison between models. In any
parameterization process, the influence of the prior knowledge about the
system is always present; using a Bayesian approach makes the influence of any
prior information on the parameters explicit. Bayesian inference has
previously been used to incorporate uncertainty quantification [25, 26, 27,
28] and model selection [29, 27] into molecular models.

To compute Bayes factors, samples must be drawn from the model posterior
distributions which involves comparing model outputs to macroscopic
observables. Typically this is done through a MCMC sampling scheme in the
parameter space, since posteriors are generally non-analytical. Once
sufficient posterior samples are obtained, a range of statistical techniques
can be used to estimate the Bayes factors and make quantitative judgements on
the relative fitness of models. One established MCMC technique for cross-model
sampling is reversible jump Monte Carlo (RJMC) [30], which simultaneously
samples over several models and their respective parameter spaces. This
technique can be used to compare model posterior parameter distributions and
directly compare the support for each model. Bayes factors can also be
computed with _bridge sampling_ methods which compute the probability ratio
between a simple analytical reference distribution and the posterior
distribution in order to calculate model evidences. This method can be
enhanced by using intermediate probability distributions to bridge the gap
between the reference distribution and the posterior, increasing the overlap
between distributions being compared and thus improving the accuracy of the
calculations [31].

These statistical approaches for calculating Bayes factors have close
analogues in statistical physics techniques used to calculate free energy
differences in molecular simulations. Specifically, Bayes factors are
analogous to the ratio of partition coefficients of two models, but integrated
over the _parameter space_ of each model, rather than over atomic coordinates.
Bridge sampling techniques are analogous to free energy perturbation and
reweighting methods [32, 33], and RJMC is analogous to expanded ensemble
techniques or Hamilton replica exchange techniques [34], where the multiple
ensembles are different models rather than different potential functions or
temperatures.

Estimating physical observables with molecular simulations is often
computationally expensive, especially if some molecular degrees of freedom are
slow. To perform MCMC sampling in parameter space, observables must be
calculated many times. To avoid this cost and efficiently sample the posterior
distribution, surrogate models [35, 36] can be used to cheaply approximate the
response surface of the observables with respect to the force field
parameters.

To demonstrate the Bayesian model comparison strategy for atomistic models, we
have selected a problem for which a surrogate model is already well-defined.
The two-center Lennard-Jones + quadrupole (2CLJQ) model [37] is a simple
4-parameter fluid model that serves as a useful test bed for this problem.
Analytical surrogate models exist for predicting saturated densities
($\rho_{l}$), saturated vapor pressure ($P_{sat}$), and surface tension
($\gamma$) for a number of small nonpolar fluids. [38, 39]. A previous
parameterization study [40] of the model demonstrated the potential to
optimize the 2CLJQ model using surrogate models. In this study we use
quantitative evidence from Bayes factors to examine whether including an
adjustable bond length and/or quadrupole parameter are justified in this
model.

##  2 Methods

###  2.1 Molecular Model

We test our Bayesian model selection strategy on the 2CLJQ fluid model [37], a
simple 2-site fluid model that fits diatomic and other similar molecules [41,
42, 43, 44]. In particular we consider this model for a range of compounds, 4
diatomic ($\mathrm{O_{2},N_{2},Br_{2},F_{2}}$) and 4 “diatomic-like”
hydro/fluoro carbon compounds
($\mathrm{C_{2}H_{2}}$, $\mathrm{C_{2}H_{4}}$, $\mathrm{C_{2}H_{6}}$,
$\mathrm{C_{2}F_{4}}$), which are both well parameterized by the 2CLJQ model
[37]. Interactions within the 2CLJQ model are controlled by four parameters,
as illustrated in Figure 1: a Lennard-Jones $\sigma$ (nm) and $\epsilon$ (K),
a constant bond length $L$ (nm) between the two LJ sites, and a quadrupole
interaction strength parameter $Q(\mathrm{D}\cdot\mathrm{nm})$. The LJ
$\sigma$ represents the distance at which the potential energy between two
particles is equal to zero; the LJ $\epsilon$ represents the depth of the
potential well and the strength of the attraction between two particles. The
functional form of the molecular model is similar to the Lennard-Jones
potential, but adapted for a 2-center model and with a quadrupole interaction
at the model’s geometric center [37]:

$u_{2CLJQ}(\mathbf{r}_{ij},\mathbf{\omega}_{i},\mathbf{\omega}_{j},L,Q)=u_{2CLJ}(\mathbf{r}_{ij},\mathbf{\omega}_{i},\mathbf{\omega}_{j},L)+u_{Q}(\mathbf{r}_{ij},\mathbf{\omega}_{i},\mathbf{\omega}_{j},Q)$
(1)
$u_{2CLJ}(\mathbf{r}_{ij},\mathbf{\omega}_{i},\mathbf{\omega}_{j},L)=\sum_{a=1}^{2}\sum_{b=1}^{2}4\epsilon\left[\left(\frac{\sigma}{r_{ab}}\right)^{12}-\left(\frac{\sigma}{r_{ab}}\right)^{6}\right]$
(2)
$u_{Q}(\mathbf{r}_{ij},\mathbf{\omega}_{i},\mathbf{\omega}_{j},Q)=\frac{3}{4}\frac{Q^{2}}{|\mathbf{r}_{ij}|^{5}}f(\mathbf{\omega}_{i},\mathbf{\omega}_{j})$
(3)

Figure 1: The 2CLJQ model consists of two LJ sites parameterized by $\sigma$
and $\epsilon$, a bond length $L$, and a quadrupole interaction strength
parameter $Q$. In this example, an ethane ($\mathrm{C_{2}H_{6}}$) molecule is
parameterized as two identical sites, each consistent of a carbon atom and
three hydrogens.

####  2.1.1 Defining models of varying complexity

To investigate which parameters are justified to reproduce the physical
properties of interest, we split the model into three levels of complexity, as
shown in Table 1: united atom (UA), anisotropic united atom (AUA), and
anisotropic united atom + quadrupole (AUA+Q).

Model |  $\sigma$ (nm) |  $\epsilon$ (K) |  $L$ (nm) |  $Q$ ($\mathrm{D}\cdot\mathrm{nm}$)
---|---|---|---|---
UA | Variable | Variable | Fixed | Fixed (at 0)
AUA | Variable | Variable | Variable | Fixed (at 0)
AUA+Q | Variable | Variable | Variable | Variable
Table 1: Different 2CLJQ models considered in this study. While all models
share the same function form, some models have parameters/interactions fixed,
while others let them vary. UA: united atom; AUA: anisotropic united atom;
AUA+Q: anisotropic united atom + quadrupole.

The models here are nested—for example, the AUA model is the same as the AUA+Q
model, but with the quadrupole moment permanently set to zero. The UA model is
a subset of the AUA model, but with a fixed bond length chosen from literature
(typically taking an “experimental” value) [45]. This construction allows us
to test whether the variable quadrupole and bond length parameters are useful
in reproducing the chosen physical properties.

####  2.1.2 Surrogate Models for 2CLJQ output

Analytical surrogate models, which predict molecular properties as a function
of molecular parameters, were developed for saturated liquid density
($\rho_{l}$) and saturated vapor pressure ($P_{sat}$) by Stoll [38]; similar
surrogate models for surface tension ($\gamma$) were developed by Werth [39]
by fitting to 2CLJQ simulation results at a variety of parameter and
temperature conditions. Stöbener [40] previously optimized parameters for the
2CLJQ model by using surrogate models to drive fast optimization and a Pareto
optimization approach to choose parameter sets. Parameter set fitness in this
study was based on comparison to data from the Design Institute for Physical
Properties (DIPPR) correlations; here we train against NIST experimental data
for similar properties instead. In this work, we expand the idea of extensive
parameter sampling through surrogate models to include comparisons between
disparate models.

To apply this technique to an arbitrary force field, one will usually need to
construct such surrogate models for different properties. Common techniques to
build these surrogate models might include reweighting [46], Gaussian
processes [47], and machine learning methods [36]. Since the methods needed to
construct such a model depend substantially on the property and the parameters
of interest, that question is beyond the scope of this study; we focus only on
applying Bayesian inference given a surrogate model.

###  2.2 Bayesian inference

Bayesian inference characterizes the fitness of models and parameters by
combining specific information about a set of experiments with more general
prior information. The core of Bayesian inference is _Bayes’ rule_ :

$\underbrace{P(\theta|D,M)}_{\mathrm{Posterior}}\propto\underbrace{P(D|\theta,M)}_{\mathrm{Likelihood}}\underbrace{P(\theta|M)}_{\mathrm{Prior}}$
(4)

where $\theta$ represents a set of parameters belonging to a model $M$, and
$D$ represents a dataset (which can be compared to experimental data)
generated by the model $M$ and parameter set $\theta$. Bayes factors, computed
as the ratio of the model marginal likelihood ($P(D|M)$) between separate
models, facilitate comparisons between those models.

$B_{1/2}=\frac{P(D|M_{1})}{P(D|M_{2})}=\frac{\int_{\theta_{1}}P(D|\Theta_{1},M_{1})P(\theta_{1}|M_{1})\mathrm{d}\theta_{1}}{\int_{\theta_{2}}P(D|\theta_{2},M_{2})P(\theta_{2}|M_{2})\mathrm{d}\theta_{2}}=\frac{P(M_{1}|D)}{P(M_{2}|D)}\frac{P(M_{2})}{P(M_{1})}$
(5)

In order to compute the model evidence $P(D|M)$, one must integrate the
parameter posterior over all parameter space in the model (as shown in eq. 2).
For most molecular systems this Bayes factor integral cannot be calculated
analytically and must be estimated. Common interpretations for levels of
significance of Bayes Factor evidence were suggested by Kass and Raftery [24]
and are listed, with labels renamed for clarity, in Table 2.

Table 2:  The interpretation of Bayes factors for model 0 over model 1, $B_{01}$, due to Kass and Raftery [24]. Bayes factors provide a quantitative measure of the _model evidence_ for choosing one model over another.  $\ln{\left(B_{01}\right)}$ | $B_{01}$ | Evidence in favor of model 0
---|---|---
0 -1 | 1 - 3 | Inconclusive
1-3 | 3 - 20 | Significant
3-5 | 20 - 150 | Strong
$>$ 5 |  $>$ 150 | Very Strong

The advantage of using Bayesian inference to choose between models is that it
captures both model performance and model parsimony, striking a balance
between accuracy and generality. More precisely, more complex models are
penalized because their additional complexity allows them to make a wider
range of predictions. Unless this additional complexity makes the predictions
uniformly more accurate, many of these predictions will be extraneous to the
properties of interest. The lower the proportion of useful prediction to total
predictions, the more a model will be penalized.

###  2.3 Construction of Posterior

Since model evidence is based on the model posterior, one must consider the
choice of the prior distributions and likelihood function that form that
posterior carefully.

####  2.3.1 Priors enable model parsimony

The use of Bayesian priors allows us to consider model parsimony in the
evaluation of models. The prior distribution defines a multi-dimensional
probability landscape in parameter space containing all values of the prior to
be considered, weighted by our prior estimates of which values are most
likely. As the model’s complexity grows, so does the volume of this parameter
space. Since priors are normalized probability distributions and must
integrate to one, growing the parameter space volume _must_ reduce the
probability of any specific parameter set.

It is also important to emphasize that Bayesian model selection chooses the
best model _given the parameter space_. If parameter values that produce good
models are excluded from this space, or assigned very low probability, those
values will not factor into the calculated Bayes factors.

Figure 2:  Demonstration of how increased model complexity incurs a penalty.
In the left panel, with the model only having the set $\theta_{1}$ as variable
parameters, a higher proportion of parameter sets are high probability. When
the additional parameter $\theta_{2}$ is added, even though more high
probability parameter sets are available, they represent a lower proportion of
the total parameter space available to the model, which incurs a complexity
penalty in the Bayesian context.

In this way, demonstrated in figure 2 the prior distribution creates a
complexity penalty that the model must overcome with increased predictive
power in order to justify increased complexity. Setting a prior is crucial
because it encodes this penalty, which depends on the amount of initial
information about the parameter.

####  2.3.2 Choosing priors

In order to avoid excluding useful parameter space, or including too much
extraneous parameter space, we form the prior distribution using a training
sample method. In this method, data is split into a “training sample”, used to
inform the prior, and a “test sample” used to calculate the Bayes factor. This
method of including partial information in training samples is a relatively
simple method for calculating priors that are effective for each model and is
established in the statistical literature [48, 49].

In this case, we start from a wide, non-negative uniform distribution, and
simulate a posterior distribution (using a likelihood calculated as in section
2.3.4) with training samples of varying amounts of experimental data. We then
fit a simple analytical distribution (Gaussian for $\epsilon,\sigma,L$;
exponential or gamma for $Q$) to this _training sample posterior_. This
analytical distribution becomes the prior for calculation of Bayes factors
under a separate set of experimental data points. We set the priors using 3
levels of information in the training sample (low, medium, high; 3, 5, 8 data
points per property respectively). These data points are selected for each
property in the target criteria and are distributed across the available
temperature range, so they are not excluding information based on temperature
or property.

####  2.3.3 Evaluating the likelihood

Likelihood evaluation requires a probabilistic model of observations in terms
of model parameters (referred to in statistical literature as a "forward
model"), as well as error model for comparing those observations with
experimental ones. In this case the forward model is composed of the surrogate
models for physical properties ($\rho_{l}$, $P_{sat}$, $\gamma$) as a function
of the model parameters ($\sigma,\epsilon,L,Q$) using the functional forms as
defined in Stoll et. al [38] (density, saturation pressure) and Werth et al.
[39] (surface tension). More details of the surrogate models are available in
Supporting Information, section 1.1.

Experimental measurements for the diatomic
($\mathrm{O_{2},N_{2},Br_{2},F_{2}}$) and hydro/fluorocarbon compounds
($\mathrm{C_{2}H_{2}}$, $\mathrm{C_{2}H_{4}}$, $\mathrm{C_{2}H_{6}}$,
$\mathrm{C_{2}F_{4}}$) for saturated liquid densities, saturation pressure,
and surface tension between 55% and 95% of critical temperature (some training
sets have slightly different temperature ranges due to data availability, full
information available in Supporting Information section 1.3) are taken from
the NIST ThermoData Engine database [50], with uncertainties assigned as a
typical class uncertainty for each compound and property. The use of class
uncertainties does not heavily impact the uncertainty as experimental
uncertainties are much smaller than correlation uncertainties. These surrogate
models are used to calculate the likelihood $P(D|\Theta)$ by estimating each
physical property at the temperatures corresponding to experimental physical
property data points and comparing those values to experimental with an error
model.

####  2.3.4 Error Model

We us a Gaussian error model to compare properties calculated from the
surrogate models to corresponding experimental properties. This error model is
chosen based on the assumption that the experimental value $x$ is the “true”
value and that the surrogate model produces a measurement of that value,
$\hat{x}$, with some error. We assume this error is Gaussian distributed, so
the probability of observing a measurement $\hat{x}$ is given by:

$P(\hat{x}|x)=\frac{1}{u_{tot}\sqrt{2\pi}}\exp\left(-\frac{1}{2}\left(\frac{x-\hat{x}}{u_{tot}}\right)^{2}\right)$
(6)

where the total uncertainty $u_{tot}$ is given by a sum of squares of the
average experimental uncertainty $u_{exp}$ and the surrogate model error
$u_{surr}$.

$u_{tot}^{2}=u_{exp}^{2}+u_{surr}^{2}$ (7)

This surrogate model is an analytical form fitted to several chemical species
at a range of thermodynamic conditions, so there is systematic uncertainty
when the model is applied to any specific compound. While the creators of the
surrogate models used in this study did not provide analytical correlations
for the uncertainty, they provided uncertainty estimates at difference
temperature and parameter values. We used these estimates to assign model
uncertainties $u_{surr}$ in our process. This model is piecewise, depending on
the temperature regime (defined as fraction of critical temperature) of the
system. Details of this error model are listed in the Supporting Information
section 1.2.

The full forms of the likelihood function (eq. 10) and prior distribution (eq.
11) are as follows:

$P(D|\Theta,M)=\left(\prod_{k=1}^{3}\left(\prod_{i=1}^{n}\frac{1}{u_{tot}\sqrt{2\pi}}\exp\left(-\frac{1}{2}\left(\frac{\vec{D_{i}}_{k}-\mathbf{f(\theta,T)}}{u_{tot}}\right)^{2}\right)\right)\right).$
(8)

The likelihood function evaluates agreement with experiment over a vector of
$i$ temperature data points from $k$ different properties.
$\mathbf{f(\theta,T)}$ is the output of the analytical surrogate model for the
$k^{th}$ physical property at the temperatures corresponding to the data
vector $\vec{D_{i}}_{k}$. It is important to note that for the saturation
pressure, the data vector $\vec{D_{i}}_{P_{sat}}$ is $\ln{P_{sat}}$ instead of
$P_{sat}$ due to the wide range of values at different temperatures. The prior
function is defined as:

$P(\Theta)=\left(\prod_{j=1}^{3}\frac{1}{\sigma_{j}\sqrt{2\pi}}\exp\left(-\frac{1}{2}\left(\frac{\theta_{j}-\mu_{j}}{\sigma_{j}}\right)^{2}\right)\right)\left(\frac{\beta^{\alpha}}{\Gamma(\alpha)}\theta_{4}^{\alpha-1}e^{-\beta\theta_{4}}\right).$
(9)

The first term of the total prior is the prior for
($\sigma=\theta_{1},\epsilon=\theta_{2},L=\theta_{3})$, and the second term is
the prior for $Q=\theta_{4}$, either an exponential ($\alpha=1,\beta$
determined by prior fitting) or gamma prior ($\alpha,\beta$ determined by
prior fitting). Distributions are fit to samples using the SciPy distributions
module and the fit functions. For compounds with quadrupole distributions
centered at values larger than zero, gamma priors are chosen; for compounds
with quadrupole support near zero, exponential distributions are chosen. For
the UA and AUA models, the prior terms for $L,Q$ (UA) or $Q$ (AUA) are set to
1, as the values are fixed. Combined, they form the posterior distribution as
described in eq. 1.

###  2.4 Sampling of Posteriors

For any method of calculating Bayes factors obtaining samples from the
parameter posterior distribution is essential. To draw samples from complex
distributions without simple closed-form distributions Markov Chain Monte
Carlo (MCMC) methods are used [51].

####  2.4.1 MCMC parameter proposals

Within a particular model, we propose moves using component-wise Metropolis
Hastings Monte Carlo [52], where at each step an individual variable parameter
is chosen with uniform probability and perturbed. This perturbation is done by
proposing a new value of the parameter based on a normal distribution centered
at the current value, with any proposed negative values rejected. The standard
deviation of this distribution is initially set to be 1/100 of the initial
value of the variable (or a minimum of 0.001 if the initial value is 0) and
then tuned to obtain a between model acceptance ratio between 20–50% during a
“burn-in” / tuning simulation ran for 1/5th the length of the production
simulation. After the burn-in period, the simulation proceeds with fixed
proposal distributions, as is required to satisfy detailed balance. Tuning
steps are discarded when computing model evidences.

####  2.4.2 Reversible jump Monte Carlo

We next turn to the more complicated question of choosing moves between
models. The reversible jump Monte Carlo (RJMC) technique [30] can be employed
to sample between multiple models with different numbers of parameters. RJMC
is an extension of traditional MCMC sampling from parameter space to combined
parameter _and_ model space. It allows for the simultaneous sampling of
parameter probability distributions from several models, as well as the
discrete model probability distribution of the collection of models. RJMC is a
powerful tool to simultaneously sample multiple models, but requires careful
construction in cases with dissimilar models. Good chain mixing requires (1)
the construction of high-quality mappings between the parameters in each model
and (2) efficient ways of proposing new values when the number of parameters
between models is uneven [53, 54, 55, 56].

For the implementation presented in this paper, proposals for variables not
present in current model, such as a quadrupole moment parameter when proposing
a move from a model without a quadrupole (UA, AUA) to the model with a
quadrupole (AUA+Q), are proposed independently from a Gaussian “proposal
distribution” that approximates the 1-D posterior distribution of that
parameter in the model they are moving to. These distributions are obtained by
fitting Gaussians to short simulations of the model posterior in order to
facilitate acceptances between inter-model moves.

####  2.4.3 Parameter mappings

To efficiently propose inter-model moves, one must define a _parameter mapping
function_ that transforms a point in one parameter space to another point in a
different parameter space (e.g. $\epsilon$ in the AUA model $\longrightarrow$
$\epsilon$ in the AUA+Q model). These mapping functions may range from a
simple direct mapping to a highly complex one, depending on the situation. The
form of the mapping function is often an important factor in model chain
mixing and convergence; as long as there is _some_ probability overlap
sampling between models is theoretically possible, but may be impossibly
inefficient. We explored several simple strategies for mapping common
parameters between models, described below and illustrated in figure 3;
although others exist, such as non-equilibrium candidate Monte Carlo [57].

  * •

Direct mapping: In the simplest cases, parameters can be taken directly from
one model and plugged into in the other model. This is only useful in the case
where the distribution of the parameter in the two models overlap
significantly.

  * •

Maximum a posteriori (MAP) mapping: Another possibility is to run a short MCMC
simulation of each model, and then map the highest probability (MAP) values of
each common parameter to each other, either by translation (additive mapping)
or rescaling (multiplicative mapping). This technique can improve overlap but
will ultimately fail when the shapes of the distributions differ
significantly, where only parameter sets near the MAP value will have good
overlap.

  * •

Gaussian affine mapping: One can also run short MCMC simulations of each
model, fit Gaussian distributions to each model posterior distribution, and
then map these distributions on top of each other using an affine map that
transforms the shape of one distribution to the shape of the other. This has
an advantage of capturing differently shaped posterior distributions, allowing
for more consistent quality mapping, especially away from the posterior
maximum, but can still fail if the distributions are multi-modal (or more
precisely, differently modal).

Figure 3: Mappings between distributions with poor overlap, visualized.  In
the first panel, there is enough overlap between the distributions to map
without any transformation. Second panel shows maximum a posteriori mapping,
improving overlap, but not matching shape. Third panel shows Gaussian affine
mapping, transforming location and matching shape.

####  2.4.4 Biasing Factors

Even with perfect mapping between the shape of the distributions between
models, inter-model moves will be rejected at a high rate if one model has
much stronger unnormalized evidence than the other. In this case a biasing
factor may be used to equalize the magnitude of the evidences. After achieving
improved sampling with the biased probabilities, one can easily unbias the
results if the bias is constant across the model. The biasing factors that
result in even sampling of models are the Bayes factors between models,
similar to how the weights leading to equal probability sampling in expanded
ensemble or simulated tempering simulations are the free energies of each
state. However, even approximate biasing factors can enable sampling such
model swaps are feasible and Bayes factors may be estimated.

###  2.5 Calculation of Bayes Factors

There are several approaches that can be used to calculate the model evidences
and Bayes factors including analysis of the reversible jump Monte Carlo model
posteriors [58] and bridge sampling [32] from the prior to the posterior.

####  2.5.1 Bridge sampling with intermediate states

One method for computing Bayes factors involves calculating the marginal
likelihoods (model normalizing constants) for each model $M$ individually, and
then calculating the Bayes factor as in eq. 5. The calculation of normalizing
constants using Monte Carlo methods is a well-known problem in statistical
inference and many methods are described in the literature [59, 60, 61, 33].

A technique common to statistical inference and statistical physics is bridge
sampling [32], which calculates the ratio of normalizing constants between two
distributions by utilizing information from Monte Carlo samples from both
distributions. The Bennett acceptance ratio (BAR) [62] is a specific case of
bridge sampling where $\alpha$ is chosen to be unbiased and minimize variance
of the estimator.

The BAR estimator is limited in cases when overlap between the two
distributions is poor. In this case, one can introduce intermediate
distributions designed to create a path between the target distributions with
sufficient overlap. In this case, which we call _bridge sampling with
intermediates_ , the ratio of normalizing constants can be estimated with the
unbiased, lowest-variance MBAR estimator [63], which extends BAR to data
collected over multiple states. The process for calculating log Bayes factors
from MBAR is analogous to the process of calculating free energies as
described by Shirts and Chodera [63], with the Boltzmann distributions
replaced with the prior $P(\theta|M)$, posterior $P(\theta|D,M)$, and
intermediate distributions. For a detailed discussion of how Bayes factors are
calculated with MBAR, see Supporting Information section 2.

Our strategy, described visually in figure 4, is to calculate the model
evidence (normalizing constant) of each model posterior $M$ by using MBAR to
calculate the ratio of normalizing constants between that model’s prior and
posterior. Since all priors used here are normalized analytical distributions,
the prior normalizing constant is known and we can calculate the posterior
normalizing constant in combination with the ration from MBAR. With this
information we can calculate the Bayes factors by comparing the posterior
normalizing constants between models.

Figure 4:  Bayes factor calculation scheme for 2CLJQ models. Bridge sampling
is used to calculate ratio of marginal likelihoods (normalizing constants)
between analytical reference distribution (with known normalizing constant)
and posterior distribution. Posterior distribution marginal likelihoods are
then calculated and used to compute Bayes factors.

Preliminary testing showed that the prior and posterior did not generally have
enough overlap to produce reliable estimates of normalizing constant ratios
with BAR. In order to improve overlap, we replace the prior in this sampling
process with a more informed auxiliary distribution. By running a short MCMC
simulation of the posterior distribution, we fit a simple multivariate normal
distribution to the posterior sample; with more complex posteriors, Gaussian
mixture models could be used, although we do not investigate them here. This
approach is not common in molecular simulation, with very complicated
distributions of samples as a functions of the coordinates, with the possible
exception of simple systems like the 2CLJQ model or crystals with near
harmonic behavior [64]. Testing also showed that intermediate distributions
were required to get sufficient overlap for bridge sampling. Distributions
$P_{int}$ are defined as:

$P_{int}=P_{ref}(\theta)^{f(\lambda)}\times P(\theta|D,M)^{f(1-\lambda)}$ (10)

where $P_{ref}$ is the analytical auxiliary distribution. Through testing we
determined that using one intermediate state between the reference
distribution and the posterior is enough to consistently achieve sufficient
overlap for Bayes factor calculations. Using intermediate distributions has an
analog in free energy simulations with $\lambda$-dynamics simulations [65],
and in staged methods of calculating host/guest binding affinities through
multiple windows.

####  2.5.2 Model distributions from RJMC

RJMC allows for the direct estimation of Bayes factors through two methods:
(1) estimating the ratio of normalizing constants by taking the sampling ratio
(SR) of the visits to each model and (2) estimating this ratio with what is
called in statistics the warped bridge sampling, using the BAR estimator. The
sampling ratio estimate is simply the ratio of the number of Monte Carlo
samples from each model; assuming each model’s parameter space is sampled
well, this will converge to the model Bayes factor. The warped bridge sampling
approach involves using bridge sampling on a set of samples from two models,
then using the proposal distributions from section 2.4.2 and parameter
mappings from section 2.4.3 to "warp" the distributions to improve overlap,
similar to mappings in configurational space. [66, 67].

####  2.5.3 Biased RJMC

As discussed in section 2.4.4, biasing factors can be used to equalize the
evidence between models and facilitate cross model jumps in cases where there
is a large difference in model evidences. In the case of our data, UA models
often have very low model evidence for this data set, compared to AUA or
AUA+Q. We can validate Bayes factor calculations as described in section 2.5.1
by performing RJMC with model evidences estimated from the importance sampling
method as biasing factors. If the biasing factors from the importance sampling
method in 2.5.1 are accurate, then biased RJMC should yield roughly equal
sampling between models. In practice, the sampling will not be perfectly equal
due to the stochastic nature of the simulation, but large deviations from
equal sampling would indicate that these biasing factors are not accurate.

###  2.6 Measuring Posterior Predictive Accuracy

To perform an benchmark of the models independent of the posterior
performance, we can calculate the _expected log pointwise predictive density_
for a new dataset (ELPPD), a metric introduced by Gelman _et al._ [68]. This
quantity is a sum of the _expected log predictive density_ (ELPD), which
evaluates the expectation value of a new data point’s ($\tilde{D}_{i}$) log
objective function over the previous posterior distribution $P(\theta|D)$. For
the objective function, we choose our standard likelihood function with the
new data point as the target, $\log P(\tilde{D}_{i}|\theta)$ .

$\mathrm{ELPPD(\tilde{D}|D)}=\sum_{i=1}^{n}\int_{\theta}\log
P(\tilde{D}_{i}|\theta)P(\theta|D)\mathrm{d}\theta$ (11)

In practice, we estimate this quantity by simulating the posterior of each
model, then drawing $k$ points at random from the posterior samples:

$\mathrm{ELPPD}_{k}(\tilde{D}|D)\approx\sum_{i=1}^{n}\sum_{j=1}^{k}\log
P(\tilde{D}_{i}|\theta_{j})\,;\,\,\,\theta_{i}\sim P(\theta|D)$ (12)

These quantities are useful in that they assess the model’s performance on
independent datasets and do not include the prior information, so they are a
more “traditional” metric to assess model performance.

###  2.7 Chosen methods for this investigation

Model priors are set using the training sample method described in section
2.3.3. To calculate Bayes factors, we used the approach of MBAR with three
intermediate states (section 2.5.1) is used, with three total $\lambda$
windows (analytical distribution, intermediate distribution, posterior
distribution) simulated for $4\times 10^{6}$ steps each. This method was
chosen over the RJMC method due to its more accurate Bayes factor
calculations. To validate these results, RJMC with biasing factors taken as
the model evidences from the MBAR calculation for $2\times 10^{6}$ steps. In
all cases, the biased RJMC calculations agree with the MBAR calculations.
ELPPD and average deviation calculations are taken from benchmarking posterior
simulations performed for $10^{6}$ steps for each model.

##  3 Results & Discussion

###  3.1 Effect of priors on Bayes factors

In order to compute meaningful Bayes factors between 2CLJQ models, we must
address the effect of priors on the model posteriors.

####  3.1.1 Illustration with uniform priors

One of the most important effects of the choice of prior distribution is that
the amount of _parameter uncertainty_ that the prior distribution encodes. In
the case of a Bayes factor between the model without a quadrupole (AUA) and
the model including a quadrupole (AUA+Q), different uniform priors over the
quadrupole have a strong effect on the Bayes Factor. This Bayes factor, shown
in figure 5 is calculated with a target of ethane (C2H6) $\rho_{l}$ and
$P_{sat}$ data.

Figure 5: Bayes factor comparison of AUA and AUA+Q models with varied
quadrupole priors. Left: Bayes factors and uncertainties between AUA and
AUA+Q, tested against C2H6 ($\rho_{L}$+$P_{sat}$) for an uniform prior
distribution on 2CLJQ quadrupole parameter $Q$ with several different
parameter values. Wider prior distributions over the quadrupole decreases
evidence in favor of AUA+Q model. Right: corresponding in color uniform prior
distributions, with sampled quadrupole posterior distribution superimposed
over priors.

As the uniform prior is widened, covering more parameter space, its
probability density decreases proportionally, which in turn, lowers the model
evidence for the AUA+Q model. This illustrates that a less certain prior over
a parameter (encoding more parameter uncertainty) will penalize a model
including that parameter. This simple example illustrates how wider uniform
priors penalize models with higher parameter uncertainty. This prior
uncertainty penalty is an intended part of Bayesian model selection. In Bayes
factors calculated in this work, we fit priors for each model in any given
comparison to the same set of training sample, so any differences in prior
uncertainty should be due to model differences rather than prior information.

####  3.1.2 Using more data in prior fitting increases prior certainty

To evaluate the effect of different amounts of prior information on the Bayes
factors between 2CLJQ models, we calculated Bayes factors with priors fit to
training samples with different amounts of experimental data. Our priors
become more certain as more training data is included. While there is no
objectively correct Bayes Factor due to its dependence on the prior, the prior
should effectively constrain the parameter space, so that extraneous parameter
values do not affect the model comparison. To ensure that the priors for the
2CLJQ model were effectively constraining the parameter space, we constructed
them with three levels of data in the training sample: n=3, 5, or 8 data
points per property respectively (referred to as low, medium, and high
information priors).

Figure 6:  Changes in prior distribution shape of AUA+Q model with different
amounts of data in training set. Lack of overlap between low data (n=3)
training sample and higher data training samples indicate that prior is not
well constrained with n=3 data points.

In this case, shown in figure 6, the overlap between the low information prior
(n=3 data points) and the medium and high information priors (n=5, n=8
respectively) is poor, indicating that the low information prior does not
constrain the prior in the same way as the other priors. In our Bayes factor
calculations, the high information prior is used, since it has excellent
overlap with the medium prior, and the increased amount of data in the
training sample increases our confidence that the prior represents parameters
that can model the whole range of experimental data points. We varied the
amount of information in this experiment to determine the effect of the amount
of data included on the prior and therefore the Bayes factor; in general, we
would use as much information as reasonably possible to inform the prior.

###  3.2 Results for $\rho_{l}$, $P_{sat}$ targets

For the $\rho_{l},P_{sat}$ target, the preferred model varies significantly
based on the molecule that it is applied to.

####  3.2.1 Diatomic molecules

For the set of diatomic molecules in figure 7, Bayes factor evidence favors
either the simplest UA model or the AUA model. For bromine, the AUA model
produces slightly different LJ parameters than the UA model, which yield
improved physical property accuracy for this target, but are incompatible with
the fixed bond length value of the UA model (sampling bond lengths of $\sim$
0.18 nm, rather than 0.23 nm, as in the UA model). While the value in the UA
model reflects a physical bond length, the anisotropic model’s assumption that
the interaction sites are not necessarily separated by the equilibrium bond
length produces better physical property estimates. This is similar to the
behavior of Oxygen, which also selects a different bond length then the
experimental value, but with weaker evidence for the AUA model.

Figure 7: Bayes factors between UA, AUA, AUA+Q models for diatomic molecules,
tested against $\rho_{l},P_{sat}$ data. Red bar indicates favored model, blue
bars indicate log Bayes factor (strength of evidence) against corresponding
models. Simpler models (UA/AUA) are favored in these cases.

####  3.2.2 Hydro/fluorocarbons

For the hydrocarbons/fluorocarbons in figure 10, the evidence generally points
towards the more complex models (AUA/AUA+Q), with the strength of the
experimental quadrupole moment generally dictating whether the $Q$ parameter
is required. C2H2 and C2F4 have the largest sampled quadrupole moments (Q
$\approx$ 0.5, 0.83), and overwhelming evidences in favor of the AUA+Q model.
C2H4 and C2H6 have lower sampled quadrupole moments (Q $\approx$ 0.41, <0.25),
and less evidence in favor of including this parameter; C2H4 and C2H6 have
significant evidence in favor of the AUA model.

Figure 8:  Top panel: Bayes factors between UA, AUA, AUA+Q models for
hydro/fluorocarbon molecules, tested against $\rho_{l},P_{sat}$ data. Red bar
indicates favored model, blue bars indicate log Bayes factor (strength of
evidence) between the favored model and the named model. Size of blue bar
indicates the strength of the evidence for the favored model (a larger bar
indicates stronger evidence for the favored model).

###  3.3 Results for $\rho_{l},P_{sat},\gamma$ targets

The addition of the surface tension target, tends to push the evidence in
favor of the more complex models, which could be expected with the need to
reproduce a very different type of molecular data.

####  3.3.1 Diatomics

Among diatomics in figure 9, both Br2 and N2 require more complex models, with
very strong evidence against the UA model, with weak evidence in favor of AUA
over AUA+Q for N2, and strong evidence in favor of AUA+Q in the case of Br2.
In the cases of F2 and O2, there is very little difference in model
performance (before accounting for parsimony) between UA and the more complex
models, leading Bayesian inference to favor UA due to parsimony. In the cases
of Br2 and N2, the more complex models (AUA and AUA+Q) correct a systematic
overprediction of both $P_{sat}$ and $\gamma$ compared to the UA model,
overcoming its parsimony advantage.

Figure 9:  Bayes factors between UA, AUA, AUA+Q models for diatomic molecules,
tested against $rho_{l},P_{sat},\gamma$ data. Red bar indicates favored model,
blue bars indicate log Bayes factor (strength of evidence) against
corresponding models. Addition of $\gamma$ data pushes model decision for Br2
and N2 towards more complex models (AUA/AUA+Q)

####  3.3.2 Hydro/fluorocarbons

For the hydro/fluorocarbons in figure 10 (C2F4 omitted due to lack of
experimental $\gamma$ data), C2H2 still has extremely strong evidence in favor
of AUA+Q, and the inclusion of $\gamma$ in the target yields weak evidence in
support of the AUA+Q model for C2H4. C2H6 has very weak evidence against the
inclusion of a quadrupole as the AUA and AUA+Q model sample very similar
values. All of these molecules have very strong evidence against the simplest
UA model.

Figure 10:  Bayes factors between UA, AUA, AUA+Q models for diatomic
molecules, tested against $rho_{l},P_{sat},\gamma$ data. Red bar indicates
favored model, blue bars indicate log Bayes factor (strength of evidence)
against corresponding models. C2F4 omitted due to lack of $\gamma$ data.

###  3.4 Parameter Correlations

An advantage of this Bayesian inference process is the information about
parameter probability distributions obtained from the MCMC process. This
information is valuable for examining parameter correlations and understanding
parameter sensitivity.

####  3.4.1 Correlation of Lennard-Jones parameters

One of notable parameter trends in this model, shown in figure 11 is the high
degree of correlation between the $\epsilon,\sigma,L$ that represent the
Lennard-Jones interactions. A degree of correlation between LJ $\epsilon$ and
$\sigma$ has been observed previously [46, 69], but the linear/planar nature
of the $(\epsilon,\sigma,L)$ probability surface suggests that the 2CLJQ model
has a number of degenerate/nearly degenerate parameter sets that do a
relatively good job of reproducing experimental data. Since the sampling is
based on the surrogate models, the correlation is probably overestimated
compared to the full model. The quadrupole parameter is also correlated with
the other parameters, but more weakly. Corner plots for all models and targets
are available in the Supporting Information, section 7.

Figure 11:  Parameter correlations shown in triangle plot of AUA+Q model’s 2-d
marginal parameter distributions for C2H4, $\rho_{l}+P_{sat}$ target. Plot
shows high degree of correlation between $\epsilon,\sigma,L$, and weaker
correlation with $Q$.

Notably, we do not find strong multi-modalities in any of these situations;
most have a single maximum of parameter probability. This is probably due to
the simplicity of the model and its surrogate models; we do not necessarily
expect this to be true for more complex models.

###  3.5 Benchmarking

Benchmarking based on ELPPD results was performed for all compounds with
enough experimental data available after prior fitting and Bayes factor
calculation. In general, measures of model performance based on Bayes factor
and likelihood-only ELPPD benchmarking agree, especially when the evidence in
favor of one model is strong. There are some situations where these methods
disagree; this is where the parsimonious nature of the Bayesian approach is
important.

Figure 12:  Average $\rho_{l}$ (left panel), $P_{sat}$ (right panel) %
deviation plots for C2H4. Parameter sets drawn from posterior probability
distribution, evaluated against separate benchmark data points (open points)
as well as points used in calculated Bayes factor (filled points). Data shows
AUA+Q model (yellow triangles) improves performance in reproducing $P_{sat}$
data at low temperatures, but Bayes factor evidence supports AUA model (blue
circles) due to parsimony.

This pattern is illustrated in figure 12, in the $\rho_{l},P_{sat}$ criteria
for C2H4. Although the performance of the AUA/AUA+Q models are very similar at
reproducing $rho_{l}$ (ELPPD average over test points is 1.78 for AUA vs 1.75
for AUA+Q), and $P_{sat}$ (0.81 for AUA vs 1.02 for AUA+Q), the Bayes factor
favors the AUA model due to the complexity penalty that the Bayesian method
assesses.

###  3.6 Challenges of implementation

A common challenge of using MCMC to sample points from a probability
distribution is ensuring good sampling and mixing within the MCMC chains. In
the posteriors sampled here, exhaustive sampling can be achieved with simple
MCMC techniques and long chains. In higher-dimension spaces with rougher
probability landscapes, directed sampling techniques like Hamiltonian Monte
Carlo [70] (HMC), Langevin Dynamics [71] or a No U-turn Sampler (NUTS) [72]
may be required to efficiently sample the space. Running multiple chains from
different starting points may also be important in situations with significant
multi-modality; advanced sampling techniques to “bridge the gap” between
disparate basins may also be helpful.

A more practical challenge to MCMC sampling is the implementation of surrogate
models used to sample parameter sets. For most parameterization problems,
implementing surrogate models will be much more difficult that for the 2CLJQ
model. For more complex models like biomolecular force fields, surrogate
models will not be general enough to capture model responses for large classes
of molecules or physical properties, so surrogate models will need to be
purpose-built for a process like this. This is a challenging task, but should
be possible by simulating properties of interest for multiple parameter sets,
and then using modeling techniques well suited to sparse data, such as
Gaussian processes (GP) [73, 74] regression. Befort et al. [47] recently had
success using GPs to build physical property surrogate models based on small
molecule force fields. These methods could be enhanced with reweighting
techniques like MBAR on simulated points to add gradient information to the
surrogate models. [46].

Another implementation challenge is selecting a technique to compute Bayes
factors from these MCMC chains. As discussed in section 2.4.2, RJMC, while
attractive due to its simplicity and simultaneous sampling, will be difficult
to implement in large model spaces, requiring much fine-tuning to facilitate
model jumps. The MBAR-based technique we use is more attractive in a large
model space, but requires some testing to ensure a good reference model, and
probability overlap between the reference model and the posterior.

##  4 Conclusions

In this work, we introduce Bayesian inference with surrogate modeling as a
paradigm for making model decisions for non-bonded molecular interactions. By
testing this strategy on the 2CLJQ model, we demonstrate that the technique
can choose between models in a way that balances model performance with
parsimony. We anticipate that this will be useful in force field
parameterization, where tradeoffs between model complexity and performance are
extremely common. We also assess the role the Bayesian prior has on the choice
of model, and offer insight into how these choices may be made in a force
field parameterization context.

The model decisions made in the Bayesian inference process vary widely due to
the compound and data targets considered, which highlights the variability of
the targets. Although this does not yield a blanket recommendation for a level
of model complexity, it identifies which compounds and properties require a
more complex model. This problem-specific parsimony is of particular interest
to the molecular modeling community, because the goal for molecular model
selection is choosing the simplest and computationally cheapest model that
describes the data of interest with sufficient experimental accuracy.

The most significant challenge of applying this technique to more complex
systems will be creating the surrogate models required for this process. For
biomolecular force fields, it is possible to build problem-specific Gaussian
process surrogate models built from simulation data and MBAR reweighting to
explore the parameter space of combinations of Lennard-Jones parameters,
though the sampling of parameter space required to construct them is a
significant challenge. This technique is also applicable to other common
dispersion-repulsion model selection problems such as LJ typing schemes, as
well as LJ combination rules. It is also potentially useful for choices
between fixed-charge electrostatic models, and could help identify where
virtual sites or polarizable sites could be most useful. Overall, this
technique appears to provide a systematic method of evaluating models based on
performance and parsimony applied to molecular simulations.

##  5 Data and Software Availability

Datasets and code used to generate the results shown in this study are
available from
https://github.com/SimonBoothroyd/bayesiantesting/tree/combined_run.

##  6 Author Contributions

Contributions based on CRediT taxonomy:

O.C.M.: Conceptualization, data curation, formal analysis, investigation,
methodology, software, visualization, original draft, review and editing
S.B.: Methodology, software, investigation, visualization, review and editing
R.A.M.: Conceptualization, data curation, resources, software, methodology,
review and editing
J.F.: Methodology, review and editing
J.D.C.: Resources, funding acquisition, review and editing
M.R.S.: Conceptualization, methodology, funding acquisition, resources,
project administration, supervision, review and editing.

##  7 Acknowledgements and Funding

We thank the Open Force Field Consortium for funding, including our industry
partners as listed at the Open Force Field website, and Molecular Sciences
Software Institute (MolSSI) for its support of the Open Force Field
Initiative. We gratefully acknowledge along with all current and former
members of the Open Force Field Initiative and the Open Force Field Scientific
Advisory Board. Research reported in this publication was in part supported by
National Institute of General Medical Sciences of the National Institutes of
Health under award number R01GM132386, specifically partial support of OCM,
MRS, and JDC. OCM, MRS, JDC, and JF acknowledge support from NSF CHE-1738975
for parts of the project. These findings are solely of the authors and do not
necessarily represent the offical views of the NIH or NSF.

##  8 Disclosures

The Chodera laboratory receives or has received funding from multiple sources,
including the National Institutes of Health, the National Science Foundation,
the Parker Institute for Cancer Immunotherapy, Relay Therapeutics, Entasis
Therapeutics, Silicon Therapeutics, EMD Serono (Merck KGaA), AstraZeneca, Vir
Biotechnology, XtalPi, Foresite Labs, the Molecular Sciences Software
Institute, the Starr Cancer Consortium, the Open Force Field Consortium, Cycle
for Survival, a Louis V. Gerstner Young Investigator Award, and the Sloan
Kettering Institute. A complete funding history for the Chodera lab can be
found at http://choderalab.org/funding. JDC is a current member of the
Scientific Advisory Boards of OpenEye Scientific Software, Interline, and
Redesign Science, and holds equity interests in Interline and Redesign
Science. MRS is an Open Science Fellow for Silicon Therapeutics. SB is a
director of Boothroyd Scientific Consulting Ltd.

## References

  * Dauber-Osguthorpe and Hagler [2019] Dauber-Osguthorpe P, Hagler AT.  Biomolecular Force Fields: Where Have We Been, Where Are We Now, Where Do We Need to Go and How Do We Get There?  J Comput Aided Mol Des. 2019 Feb; 33(2):133–203.  doi: 10.1007/s10822-018-0111-4.
  * Hagler [2019] Hagler AT.  Force Field Development Phase II: Relaxation of Physics-Based Criteria… or Inclusion of More Rigorous Physics into the Representation of Molecular Energetics.  J Comput Aided Mol Des. 2019 Feb; 33(2):205–264.  doi: 10.1007/s10822-018-0134-x.
  * Riniker [2018] Riniker S.  Fixed-Charge Atomistic Force Fields for Molecular Dynamics Simulations in the Condensed Phase: An Overview.  J Chem Inf Model. 2018 Mar; 58(3):565–578.  doi: 10.1021/acs.jcim.8b00042.
  * Monticelli and Salonen [2013] Monticelli L, Salonen E, editors.  Biomolecular Simulations: Methods and Protocols.  Methods in Molecular Biology, Humana Press; 2013.  https://www.springer.com/us/book/9781627030168.
  * Kantonen et al. [2020] Kantonen SM, Muddana HS, Schauperl M, Henriksen NM, Wang LP, Gilson MK.  Data-Driven Mapping of Gas-Phase Quantum Calculations to General Force Field Lennard-Jones Parameters.  J Chem Theory Comput. 2020 Feb; 16(2):1115–1127.  doi: 10.1021/acs.jctc.9b00713.
  * Tkatchenko and Scheffler [2009] Tkatchenko A, Scheffler M.  Accurate Molecular Van Der Waals Interactions from Ground-State Electron Density and Free-Atom Reference Data.  Phys Rev Lett. 2009 Feb; 102(7):073005.  doi: 10.1103/PhysRevLett.102.073005.
  * Wang et al. [2004] Wang J, Wolf RM, Caldwell JW, Kollman PA, Case DA.  Development and Testing of a General Amber Force Field.  J Comput Chem. 2004; 25(9):1157–1174.  doi: 10.1002/jcc.20035.
  * Banks et al. [2005] Banks JL, Beard HS, Cao Y, Cho AE, Damm W, Farid R, Felts AK, Halgren TA, Mainz DT, Maple JR, Murphy R, Philipp DM, Repasky MP, Zhang LY, Berne BJ, Friesner RA, Gallicchio E, Levy RM.  Integrated Modeling Program, Applied Chemical Theory (IMPACT).  J Comput Chem. 2005; 26(16):1752–1780.  doi: 10.1002/jcc.20292.
  * Ponder et al. [2010] Ponder JW, Wu C, Ren P, Pande VS, Chodera JD, Schnieders MJ, Haque I, Mobley DL, Lambrecht DS, DiStasio RA, Head-Gordon M, Clark GNI, Johnson ME, Head-Gordon T.  Current Status of the AMOEBA Polarizable Force Field.  J Phys Chem B. 2010 Mar; 114(8):2549–2564.  doi: 10.1021/jp910674d.
  * Jiao et al. [2008] Jiao D, Golubkov PA, Darden TA, Ren P.  Calculation of Protein–Ligand Binding Free Energy by Using a Polarizable Potential.  PNAS. 2008 Apr; 105(17):6290–6295.  doi: 10.1073/pnas.0711686105.
  * Lemkul and MacKerell [2017] Lemkul JA, MacKerell AD.  Polarizable Force Field for DNA Based on the Classical Drude Oscillator: I. Refinement Using Quantum Mechanical Base Stacking and Conformational Energetics.  J Chem Theory Comput. 2017 May; 13(5):2053–2071.  doi: 10.1021/acs.jctc.7b00067.
  * Patel and Brooks [2004] Patel S, Brooks CL.  CHARMM Fluctuating Charge Force Field for Proteins: I Parameterization and Application to Bulk Organic Liquid Simulations.  J Comput Chem. 2004; 25(1):1–16.  doi: 10.1002/jcc.10355.
  * Messerly et al. [2019] Messerly RA, Anderson MC, Razavi SM, Elliott JR.  Mie 16–6 Force Field Predicts Viscosity with Faster-than-Exponential Pressure Dependence for 2,2,4-Trimethylhexane.  Fluid Phase Equilibria. 2019 Sep; 495:76–85.  doi: 10.1016/j.fluid.2019.05.013.
  * Chiu et al. [2010] Chiu SW, Scott HL, Jakobsson E.  A Coarse-Grained Model Based on Morse Potential for Water and n-Alkanes.  J Chem Theory Comput. 2010 Mar; 6(3):851–863.  doi: 10.1021/ct900475p.
  * Halgren [1996] Halgren TA.  Merck Molecular Force Field. II. MMFF94 van Der Waals and Electrostatic Parameters for Intermolecular Interactions.  J Comput Chem. 1996; 17(5-6):520–552.  doi: 10.1002/(SICI)1096-987X(199604)17:5/6<520::AID-JCC2>3.0.CO;2-W.
  * Schauperl et al. [2020] Schauperl M, Kantonen SM, Wang LP, Gilson MK.  Data-Driven Analysis of the Number of Lennard–Jones Types Needed in a Force Field.  Commun Chem. 2020 Nov; 3(1):1–12.  doi: 10.1038/s42004-020-00395-w.
  * Boulanger et al. [2018] Boulanger E, Huang L, Rupakheti C, MacKerell AD, Roux B.  Optimized Lennard-Jones Parameters for Drug-Like Small Molecules.  J Chem Theory Comput. 2018 Jun; 14(6):3121–3131.  doi: 10.1021/acs.jctc.8b00172.
  * Halgren [1992] Halgren TA.  The Representation of van Der Waals (vdW) Interactions in Molecular Mechanics Force Fields: Potential Form, Combination Rules, and vdW Parameters.  J Am Chem Soc. 1992 Sep; 114(20):7827–7843.  doi: 10.1021/ja00046a032.
  * Waldman and Hagler [1993] Waldman M, Hagler AT.  New Combining Rules for Rare Gas van Der Waals Parameters.  J Comput Chem. 1993; 14(9):1077–1084.  doi: 10.1002/jcc.540140909.
  * Harrison et al. [2018] Harrison JA, Schall JD, Maskey S, Mikulski PT, Knippenberg MT, Morrow BH.  Review of Force Fields and Intermolecular Potentials Used in Atomistic Computational Materials Research.  Applied Physics Reviews. 2018 Aug; 5(3):031104.  doi: 10.1063/1.5020808.
  * Fröhlking et al. [2020] Fröhlking T, Bernetti M, Calonaci N, Bussi G.  Toward Empirical Force Fields That Match Experimental Observables.  J Chem Phys. 2020 Jun; 152(23):230902.  doi: 10.1063/5.0011346.
  * von Toussaint [2011] von Toussaint U.  Bayesian Inference in Physics.  Rev Mod Phys. 2011 Sep; 83(3):943–999.  doi: 10.1103/RevModPhys.83.943.
  * MacKay [????] MacKay DJC.  Information Theory, Inference, and Learning Algorithms;.  https://www.inference.org.uk/itprnn/book.pdf.
  * Kass and Raftery [1995] Kass RE, Raftery AE.  Bayes Factors.  Journal of the American Statistical Association. 1995 Jun; 90(430):773–795.  doi: 10.1080/01621459.1995.10476572.
  * Dutta et al. [2018] Dutta R, Brotzakis ZF, Mira A.  Bayesian Calibration of Force-Fields from Experimental Data: TIP4P Water.  J Chem Phys. 2018 Oct; 149(15):154110.  doi: 10.1063/1.5030950.
  * Angelikopoulos et al. [2012] Angelikopoulos P, Papadimitriou C, Koumoutsakos P.  Bayesian Uncertainty Quantification and Propagation in Molecular Dynamics Simulations: A High Performance Computing Framework.  J Chem Phys. 2012 Oct; 137(14):144103.  doi: 10.1063/1.4757266.
  * Farrell et al. [2015] Farrell K, Oden JT, Faghihi D.  A Bayesian Framework for Adaptive Selection, Calibration, and Validation of Coarse-Grained Models of Atomistic Systems.  Journal of Computational Physics. 2015 Aug; 295:189–208.  doi: 10.1016/j.jcp.2015.03.071.
  * Wu S. et al. [2016] Wu S , Angelikopoulos P , Papadimitriou C , Moser R , Koumoutsakos P .  A Hierarchical Bayesian Framework for Force Field Selection in Molecular Dynamics Simulations.  Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences. 2016 Feb; 374(2060):20150032.  doi: 10.1098/rsta.2015.0032.
  * Bacallado et al. [2009] Bacallado S, Chodera JD, Pande V.  Bayesian Comparison of Markov Models of Molecular Dynamics with Detailed Balance Constraint.  J Chem Phys. 2009 Jul; 131(4):045106.  doi: 10.1063/1.3192309.
  * Green [1995] Green PJ.  Reversible Jump Markov Chain Monte Carlo Computation and Bayesian Model Determination.  Biometrika. 1995; 82(4):711–32.
  * Neal [1998] Neal RM.  Annealed Importance Sampling.  arXiv:physics/9803008. 1998 Sep; http://arxiv.org/abs/physics/9803008.
  * Meng and Wong [1996] Meng XL, Wong WH.  Simulating Ratios of Normalizing Constants via a Simple Identity: A Theoretical Exploration.  Stat Sin. 1996; 6(4):831–860.  https://www.jstor.org/stable/24306045.
  * Gelman and Meng [1998] Gelman A, Meng XL.  Simulating Normalizing Constants: From Importance Sampling to Bridge Sampling to Path Sampling.  Stat Sci. 1998 May; 13(2):163–185.  doi: 10.1214/ss/1028905934.
  * Sugita and Okamoto [2000] Sugita Y, Okamoto Y.  Replica-Exchange Multicanonical Algorithm and Multicanonical Replica-Exchange Method for Simulating Systems with Rough Energy Landscape.  Chemical Physics Letters. 2000 Oct; 329(3):261–270.  doi: 10.1016/S0009-2614(00)00999-4.
  * Sidky et al. [2020] Sidky H, Chen W, L Ferguson A.  Molecular Latent Space Simulators.  Chem Sci. 2020; 11(35):9459–9467.  doi: 10.1039/D0SC03635H.
  * Kadupitiya et al. [2020] Kadupitiya JCS, Sun F, Fox G, Jadhao V.  Machine Learning Surrogates for Molecular Dynamics Simulations of Soft Materials.  Journal of Computational Science. 2020 Apr; 42:101107.  doi: 10.1016/j.jocs.2020.101107.
  * Vrabec et al. [2001] Vrabec J, Stoll J, Hasse H.  A Set of Molecular Models for Symmetric Quadrupolar Fluids.  J Phys Chem B. 2001 Dec; 105(48):12126–12133.  doi: 10.1021/jp012542o.
  * Stoll et al. [2009] Stoll J, Vrabec J, Hasse H.  Comprehensive Study of the Vapour-Liquid Equilibria of the Pure Two-Centre Lennard-Jones plus Pointquadrupole Fluid.  ArXiv09043679 Phys. 2009 Apr; http://arxiv.org/abs/0904.3679.
  * Werth et al. [2015] Werth S, Horsch M, Hasse H.  Surface Tension of the Two Center Lennard-Jones plus Quadrupole Model Fluid.  Fluid Phase Equilibria. 2015 Apr; 392:12–18.  doi: 10.1016/j.fluid.2015.02.003.
  * Stöbener et al. [2016] Stöbener K, Klein P, Horsch M, Küfer K, Hasse H.  Parametrization of Two-Center Lennard-Jones plus Point-Quadrupole Force Field Models by Multicriteria Optimization.  Fluid Phase Equilibria. 2016 Mar; 411:33–42.  doi: 10.1016/j.fluid.2015.11.028.
  * Stoll et al. [2003] Stoll J, Vrabec J, Hasse H.  A Set of Molecular Models for Carbon Monoxide and Halogenated Hydrocarbons.  J Chem Phys. 2003 Nov; 119(21):11396–11407.  doi: 10.1063/1.1623475.
  * Möller and Fischer [1994] Möller D, Fischer J.  Determination of an Effective Intermolecular Potential for Carbon Dioxide Using Vapour-Liquid Phase Equilibria from NpT + Test Particle Simulations.  Fluid Phase Equilibria. 1994 Sep; 100:35–61.  doi: 10.1016/0378-3812(94)80002-2.
  * Cheung and Powles [1976] Cheung PSY, Powles JG.  The Properties of Liquid Nitrogen.  Mol Phys. 1976 Nov; 32(5):1383–1405.  doi: 10.1080/00268977600102761.
  * Murthy et al. [1981] Murthy CS, Singer K, McDonald IR.  Interaction Site Models for Carbon Dioxide.  Mol Phys. 1981 Sep; 44(1):135–143.  doi: 10.1080/00268978100102331.
  * joh [2018] NIST Computational Chemistry Comparison and Benchmark Database; 2018.  https://cccbdb.nist.gov/.
  * Messerly et al. [2018] Messerly RA, Razavi SM, Shirts MR.  Configuration-Sampling-Based Surrogate Models for Rapid Parameterization of Non-Bonded Interactions.  J Chem Theory Comput. 2018 Jun; 14(6):3144–3162.  doi: 10.1021/acs.jctc.8b00223.
  * Befort et al. [2021] Befort BJ, DeFever RS, Tow GM, Dowling AW, Maginn EJ.  Machine Learning Directed Optimization of Classical Molecular Modeling Force Fields.  ArXiv210303208 Phys. 2021 Mar; http://arxiv.org/abs/2103.03208.
  * Berger and Pericchi [1996] Berger JO, Pericchi LR.  The Intrinsic Bayes Factor for Model Selection and Prediction.  J Am Stat Assoc. 1996; 91(433):109–122.  doi: 10.2307/2291387.
  * Berger [1990] Berger JO.  Robust Bayesian Analysis: Sensitivity to the Prior.  Journal of Statistical Planning and Inference. 1990 Jul; 25(3):303–328.  doi: 10.1016/0378-3758(90)90079-A.
  * Frenkel et al. [2005] Frenkel M, Chirico RD, Diky V, Yan X, Dong Q, Muzny C.  ThermoData Engine (TDE): Software Implementation of the Dynamic Data Evaluation Concept.  J Chem Inf Model. 2005 Jul; 45(4):816–838.  doi: 10.1021/ci050067b.
  * Chipman et al. [2001] Chipman H, George EI, McCulloch RE, Clyde M, Foster DP, Stine RA.  The Practical Implementation of Bayesian Model Selection.  Lect Notes-Monogr Ser. 2001; 38:65–134.  https://www.jstor.org/stable/4356164.
  * Johnson et al. [2013] Johnson AA, Jones GL, Neath RC.  Component-Wise Markov Chain Monte Carlo: Uniform and Geometric Ergodicity under Mixing and Composition.  Statist Sci. 2013 Aug; 28(3):360–375.  doi: 10.1214/13-STS423.
  * Karagiannis and Andrieu [2013] Karagiannis G, Andrieu C.  Annealed Importance Sampling Reversible Jump MCMC Algorithms.  J Comput Graph Stat. 2013 Jul; 22(3):623–648.  doi: 10.1080/10618600.2013.805651.
  * Brooks et al. [2003a] Brooks SP, Friel N, King R.  Classical Model Selection via Simulated Annealing.  J R Stat Soc Ser B Stat Methodol. 2003; 65(2):503–520.  doi: 10.1111/1467-9868.00399.
  * Brooks et al. [2003b] Brooks SP, Giudici P, Roberts GO.  Efficient Construction of Reversible Jump Markov Chain Monte Carlo Proposal Distributions.  J R Stat Soc Ser B Stat Methodol. 2003; 65(1):3–39.  doi: 10.1111/1467-9868.03711.
  * Hastie and Green [2012] Hastie DI, Green PJ.  Model Choice Using Reversible Jump Markov Chain Monte Carlo.  Stat Neerlandica. 2012; 66(3):309–338.  doi: 10.1111/j.1467-9574.2012.00516.x.
  * Nilmeier et al. [2011] Nilmeier JP, Crooks GE, Minh DDL, Chodera JD.  Nonequilibrium Candidate Monte Carlo Is an Efficient Tool for Equilibrium Simulation.  PNAS. 2011 Nov; 108(45):E1009–E1018.  doi: 10.1073/pnas.1106094108.
  * Bartolucci et al. [2006] Bartolucci F, Scaccia L, Mira A.  Efficient Bayes Factor Estimation from the Reversible Jump Output.  Biometrika. 2006 Mar; 93(1):41–52.  doi: 10.1093/biomet/93.1.41.
  * Geweke [1989] Geweke J.  Bayesian Inference in Econometric Models Using Monte Carlo Integration.  Econometrica. 1989; 57(6):1317–1339.  doi: 10.2307/1913710.
  * Neal [1993] Neal RM.  Probabilistic Inference Using Markov Chain Monte Carlo Methods; 1993\.
  * DiCiccio et al. [1997] DiCiccio TJ, Kass RE, Raftery A, Wasserman L.  Computing Bayes Factors by Combining Simulation and Asymptotic Approximations.  J Am Stat Assoc. 1997; 92(439):903–915.  doi: 10.2307/2965554.
  * Bennett [1976] Bennett CH.  Efficient Estimation of Free Energy Differences from Monte Carlo Data.  Journal of Computational Physics. 1976 Oct; 22(2):245–268.  doi: 10.1016/0021-9991(76)90078-4.
  * Shirts and Chodera [2008] Shirts MR, Chodera JD.  Statistically Optimal Analysis of Samples from Multiple Equilibrium States.  J Chem Phys. 2008 Sep; 129(12).  doi: 10.1063/1.2978177.
  * Abraham and Shirts [2018] Abraham NS, Shirts MR.  Thermal Gradient Approach for the Quasi-Harmonic Approximation and Its Application to Improved Treatment of Anisotropic Expansion.  J Chem Theory Comput. 2018 Nov; 14(11):5904–5919.  doi: 10.1021/acs.jctc.8b00460.
  * Knight and Brooks [2009] Knight JL, Brooks CL.  $\lambda$-Dynamics Free Energy Simulation Methods.  J Comput Chem. 2009; 30(11):1692–1700.  doi: 10.1002/jcc.21295.
  * Paliwal and Shirts [2013] Paliwal H, Shirts MR.  Multistate Reweighting and Configuration Mapping Together Accelerate the Efficiency of Thermodynamic Calculations as a Function of Molecular Geometry by Orders of Magnitude.  J Chem Phys. 2013 Apr; 138(15):154108.  doi: 10.1063/1.4801332.
  * Schieber and Shirts [2019] Schieber NP, Shirts MR.  Configurational Mapping Significantly Increases the Efficiency of Solid-Solid Phase Coexistence Calculations via Molecular Dynamics: Determining the FCC-HCP Coexistence Line of Lennard-Jones Particles.  J Chem Phys. 2019 Apr; 150(16):164112.  doi: 10.1063/1.5080431.
  * Gelman et al. [2014] Gelman A, Hwang J, Vehtari A.  Understanding Predictive Information Criteria for Bayesian Models.  Stat Comput. 2014 Nov; 24(6):997–1016.  doi: 10.1007/s11222-013-9416-2.
  * Shirts and Pande [2005] Shirts MR, Pande VS.  Solvation Free Energies of Amino Acid Side Chain Analogs for Common Molecular Mechanics Water Models.  J Chem Phys. 2005 Apr; 122(13):134508.  doi: 10.1063/1.1877132.
  * Neal [2012] Neal RM.  MCMC Using Hamiltonian Dynamics.  ArXiv12061901 Phys Stat. 2012 Jun; http://arxiv.org/abs/1206.1901.
  * Leimkuhler and Matthews [2013] Leimkuhler B, Matthews C.  Robust and Efficient Configurational Molecular Sampling via Langevin Dynamics.  J Chem Phys. 2013 May; 138(17):174102.  doi: 10.1063/1.4802990.
  * Hoffman and Gelman [2011] Hoffman MD, Gelman A.  The No-U-Turn Sampler: Adaptively Setting Path Lengths in Hamiltonian Monte Carlo.  ArXiv11114246 Cs Stat. 2011 Nov; http://arxiv.org/abs/1111.4246.
  * Liu et al. [2019] Liu H, Ong YS, Shen X, Cai J.  When Gaussian Process Meets Big Data: A Review of Scalable GPs.  ArXiv180701065 Cs Stat. 2019 Apr; http://arxiv.org/abs/1807.01065.
  * Oliver and Webster [1990] Oliver MA, Webster R.  Kriging: A Method of Interpolation for Geographical Information Systems.  Int J Geogr Inf Syst. 1990 Jul; 4(3):313–332.  doi: 10.1080/02693799008941549.