Home | About | Journals | Submit | Contact Us | Français |

**|**NIST Author Manuscripts**|**PMC5087108

Formats

Article sections

- Abstract
- I. INTRODUCTION
- II. CALIBRATION WORKFLOW
- III. DEPENDENCE ON UNDERLYING ASSUMPTIONS AND ACCELERATED COARSE-GRAINING
- IV. DISCUSSION
- References

Authors

Related links

J Chem Phys. Author manuscript; available in PMC 2017 April 21.

Published in final edited form as:

PMCID: PMC5087108

NIHMSID: NIHMS819324

National Institute of Standards and Technology

The publisher's final edited version of this article is available at J Chem Phys

Generating and calibrating forces that are transferable across a range of state-points remains a challenging task in coarse-grained (CG) molecular dynamics (MD). In this work, we present a coarse-graining *workflow*, inspired by ideas from uncertainty quantification and numerical analysis, to address this problem. The key idea behind our approach is to introduce a Bayesian correction algorithm that uses functional derivatives of CG simulations to rapidly and inexpensively recalibrate initial estimates *f*_{0} of forces anchored by standard methods such as Force-Matching (FM). Taking density-temperature relationships as a running example, we demonstrate that this algorithm, in concert with various interpolation schemes, can be used to efficiently compute physically reasonable force curves on a fine grid of state-points. Importantly, we show that our work flow is robust to several choices available to the modeler, including the interpolation schemes and tools used to construct *f*_{0}. In a related vein, we also demonstrate that our approach can speed up coarse-graining by reducing the number of atomistic simulations needed as inputs to standard methods for generating CG forces.

In recent years, the high-throughput capabilities of molecular dynamics (MD) simulations have attracted the interest of companies involved in materials research and development.^{1–3} In such settings, these tools are becoming an integral component of the decision-making process, given their ability to rapidly search design spaces for novel systems that warrant further experimental study. While it is generally agreed that this procedure can significantly reduce material development costs, a key theme has repeatedly emerged from practical implementations thereof: in order to be informative for commercial applications, computational models require techniques that can rapidly calibrate inputs and assess uncertainties in simulated predictions.^{2–4}

In the context of coarse-grained (CG) MD, this theme has proven to be especially germane. In particular, high-throughput applications generally require the ability to rapidly customize the level of coarse-graining for arbitrary systems across a range of state points (e.g. temperatures); *inexpensive and robust* calibration techniques are critical for enabling this degree of exibility. Often, the key input to be determined is the CG interparticle forcefield *f*(*r, q*) as a function of separation *r* and state-point *q*, and many techniques have been developed to compute these functions from thermodynamic considerations of atomistic systems.^{5–13} However, such methods typically draw information from a limited number of state-points, and the resulting forces can yield incorrect predictions for even small changes away from the calibration conditions (see Fig. 1). This has led several authors to consider the following questions:^{14–25}

Comparison of coarse-grained density-temperature relation (solid line) and atomistic results (
) for polystyrene. To generate atomistic data, we simulated a roughly 8000 atom system for 10 ns at anchor temperatures *T* = 100 K, 200 K,..., 800 K and averaged **...**

- Is there a general and inexpensive method for calibrating f(r, q) so that a CG simulation yields the same predictions as a corresponding atomistic simulation on a large interval of state-points?
- Can uncertainty quantification (UQ) be used to speed up the CG process?

In this paper, we propose a workflow that addresses these questions through a combination of techniques inspired by Bayesian probability, uncertainty quantification, and numerical analysis. This main steps in this work flow amount to:

- building an analytical, target model (
*q*) to describe the predictions that the CG simulations should reproduce; - postulating an initial or guess set of forces
*f*_{0}(*r, q*); - iteratively correcting
*f*_{0}via a pointwise (in*q*) but inexpensive Bayesian algorithm that- estimates functional derivatives $\frac{\delta \mathbb{S}}{\delta f}[{f}_{0},q]$ of the CG prediction [
*f*_{0}*, q*], and - uses these derivatives to correct
*f*_{0}by comparing CG simulations to the target.

Importantly, this work flow overcomes the difficulties of generating transferable forces via the low cost of Step III, which can quickly calibrate *f*_{0} constructed by non-robust but cheap interpolations in state-space. Moreover, this approach is meant to complement existing methods for coarse-graining insofar as it can use their outputs to inform these initial guesses. In contrast to previous works, our procedure also yields a probability *P*[*f; q*] expressing our belief that a given *f* is a physically plausible and well-calibrated input to the CG simulations.^{26} By means of this *P*, we can compute an average (or most-plausible) estimate (*r, q*) for which the difference [ (*r, q*)*, q*]*–* (*q*) can be rendered statistically insignificant by iteratively applying Step III.

A central idea that facilitates our analysis is a spectral decomposition of the forces via

$$f(r,q)=\sum _{m=1}^{M}{w}_{m}(q){\varphi}_{m}(r)=\mathit{w}\xb7\mathit{\varphi},$$

(1)

where *ϕ _{m}*(

$$\frac{\delta \mathbb{S}}{\delta f}\approx \frac{\mathbb{S}[\mathit{w}+\epsilon \mathrm{\Delta}\mathit{w}]-\mathbb{S}[\mathit{w}]}{\epsilon},\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\epsilon \ll 1$$

(2)

(for a given direction Δ** w**) to alter the weights in systematic and meaningful ways.

As we have suggested, though, a key difficulty repeatedly appears: it is generally too expensive to realize any simulated quantity at a dense set of state-points. The problem is often most pronounced when generating individual *f*_{0} by CG methods such as Force Matching (FM), which can take thousands of CPU hours; similar problems arise when estimating functional derivatives for large *M*. As a result, a frequent and unavoidable task in our analysis is to interpolate weights and derivatives between so-called anchor points in the variable *q. Importantly, we desire that these interpolations yield an estimate of their uncertainties in order to avoid overstating our confidence in *(*r, q*). In what follows, we make extensive use of Gaussian-process regression (GPR). We refer the reader to Refs. ^{29}^{–}^{31} for excellent reviews of this subject.

The problem of *non-transferable* force fields has been considered by several authors; see, for example, Refs. ^{7}^{, }^{9}^{, }^{12}^{, }^{14}^{–}^{18}^{, }^{21}^{–}^{25}^{, }^{32}^{–}^{36} and Sec. IV A. In many instances, these works have been successful in computing interparticle potentials that correctly predict material properties across a range of state points. To the best of our knowledge, however, these methods are limited in at least one of the following ways: they (i) cannot be calibrated to arbitrary reference data; (ii) are tailored for specific systems and/or material properties; (iii) lack a systematic correction algorithm to update force curves that yield incorrect predictions. Here our goal is to propose a generic, inexpensive, and *robust* calibration method whose usefulness does not strongly depend on either the system at hand or the details of coarse-graining. By making uncertainty quantification an integral part of the workflow, we also emphasize the *subjective* nature of coarse-graining, thereby uncovering decision-making elements of the problem. Within this framework, the CG methods advocated by the aforementioned references can be incorporated to inform such decisions in a systematic way.

We mention a few limitations of our analysis. For one, estimates of (*r, q*) only yield CG predictions that are as good as our knowledge of the target properties in the first place. Given that the latter are never known with infinite precision, there is inherent uncertainty in our assessment of what constitutes a calibrated force curve. Moreover, the Bayesian algorithm may not necessarily calibrate *f* against all properties of interest, even though the procedure is meant to be general in nature. This question is difficult to address insofar as its answer depends on how smoothly the CG predictions change relative to *w _{m}*. We expand on these and further limitations in the following sections.

The rest of the paper is organized as follows. In Sec. II, we present the main ideas and key equations behind our calibration workflow. We begin by assuming an initial set of force curves and show how corrections are computed. In Sec. III, we consider stability of the correction algorithm to changes in the underlying inputs, e.g. the interpolation methods used in Sec. II. Section IV frames our analysis in the context of past works and more fully considers extensions and limitations of our approach. Throughout, we illustrate aspects of the analysis via a running example of calibrating density-temperature relations for CG polystyrene. Appendices A–D provide relevant information on simulation techniques, FM and the coarse-grained system, building target models, and interpolation procedures.

We review notation and other conventions used throughout this work.

- The capital letters
*N, M, Q, S*denote positive integers. - The letter
*q*denotes a state-point such as temperature or pressure. In our running examples throughout, we will take*q*to be the temperature*T*. - Non-bold letters are scalars
*if no argument or a non-bold (scalar) argument is indicated.*Examples include*q*,*w*, and*w*(*q*). - Bold letters
**q**= (*q*_{1},…,*q*) and_{Q}= (*w**w*_{1}*,*…*,w*) = (_{M}*w*_{1}(*q*) ,…,*w*(_{M}*q*)) =(*w**q*) denote vectors. In particular,**q**is a vector of state-points, while(*w**q*) = (*w*_{1}(*q*) ,…,*w*(_{M}*q*)) denotes the weights in Eq. (1) evaluated at a fixed*q*. - In contrast to normal conventions, we allow a scalar function of a scalar argument to also take a vector argument. In such cases, this denotes the vector whose components are the scalar function evaluated at each element of the vector. That is,
*w*(**q**) = (w(*q*_{1}) ,…,*w*(*q*)) if_{Q}**q**= (*q*_{1},…,*q*)._{Q} - We use
*k*(*q, q*′) to denote a scalar function of two variables. If we replace the arguments with*Q*-dimensional vectors**q**and**q**′, the object*k*(**q***,***q**′) is a*Q × Q*matrix$$k(\mathbf{q},{\mathbf{q}}^{\prime})=\left[\begin{array}{cccc}k({q}_{1},{q}_{1}^{\prime})& k({q}_{1},{q}_{2}^{\prime})& \dots & k({q}_{1},{q}_{Q}^{\prime})\\ k({q}_{2},{q}_{1}^{\prime})& k({q}_{2},{q}_{2}^{\prime})& \dots & k({q}_{2},{q}_{Q}^{\prime})\\ \vdots & \vdots & \ddots & \vdots \\ k({q}_{Q},{q}_{1}^{\prime})& k({q}_{Q},{q}_{2}^{\prime})& \dots & k({q}_{Q},{q}_{Q}^{\prime})\end{array}\right].$$(3)Likewise,*k*(*q,***q**) and*k*(**q***, q*) denote the corresponding row and column vectors formed by pairing*q*with the elements of**q**. - We indicate elements of a vector
or matrix*v***Σ**with the subscripted, non-bold versions*v*and Σ_{i}. When order of operations is otherwise unclear, we use the notation [_{i,j}]*v*and [_{j}**Σ**]; see immediately below._{i,j} - The use of a Roman T as a superscript denotes the transpose, not the temperature as an exponent. That is, [
**Σ**^{T}]= Σ_{i,j}for matrix_{j,i}**Σ**. - Unless otherwise noted, we distinguish the meaning of indices according to the parent symbol as follows.
*q*refers to the_{k}*k*th state-point in a collection thereof.- Importantly,
*w*,_{n}*f*,_{n}*P*, and variations thereof refer to the_{n}*n*th updates of (or corrections to) the weights, forces, and probabilities. - The double-indexed symbol
*w*(_{m,n}*q*) denotes the*n*th correction of the*m*th mode evaluated at*q*. - Superscripts in braces, such as {
*i*} in*w*^{{}^{i}^{}}, denote the*i*th independent, identically distributed (i.i.d.) realization of a random quantity.

- We sometimes denote a probability over forces using
*P*[*f; q*]. In this expression,*P*depends on*q*but is not normalized with respect to this variable.

Stated imprecisely, our goal is to calibrate *f* on an interval of state-points by ensuring that it (a) is thermodynamically consistent with an atomistic model, and (b) yields the same predictions as some target (*q*), e.g. experimental data or atomistic simulations. Given infinite computational resources, any of the techniques advocated by Refs. ^{5}^{–}^{9}^{, }^{12}^{, }^{15}^{– }^{18}^{, }^{21}^{–}^{23}^{, }^{32}^{–}^{34} determine functions *f*_{0} that satisfy these criteria pointwise in *q*; thus the problem is solved in theory.

In practice, however, several difficulties arise. First, coarse-graining methods generally rely on exhaustive sampling of configurations from fully atomistic simulations. This sampling is expensive, and often one can only realize such *f*_{0} at a small set of anchor points *q _{k}*. Interpolation is sometimes used to overcome this limitation, but such techniques suffer from a second issue: CG simulations may be sensitive to minor changes in

These observations suggest the need to augment a typical calibration procedure with systematic correction and uncertainty quantification (UQ) steps that address the limitations of interpolation. To achieve this, we begin by invoking a Bayesian perspective wherein we model our *belief* that a given *f* is the best choice to satisfy the criteria (a) and (b) above. Operationally, we ask: what is the probability *P*[*f; q*] that *f* is well-calibrated at state-point *q*. We caution that this probability should not be understood in the frequentist sense; i.e. there is no notion of a random draw from a set of optimal force curves. Rather, *P* quantifies belief in a way that allows for analytical and systematic treatment.^{37} This is useful insofar as we can encode information about our initial estimates *f*_{0} (computed, e.g. from coarse-graining or interpolation) and their uncertainties in a prior probability *P*_{0}[*f ; q*], which we update by comparing CG simulations with the target model. In distinctly Bayesian language, we compute a *posterior P*_{1}[*f ; q*] via a relation of the form

$${P}_{n+1}[f;q]\propto P[\mathbb{T}(q)\mid f]{P}_{n}[f;q]=\mathcal{L}(f\mid \delta ,q){P}_{n}[f;q]$$

(4)

where *P*[ (*q*)|*f*] = (*f|δ, q*) is the probability that *f* yields the target predictions (*q*), or equivalently a likelihood function ** expressing our belief that ***f* is well calibrated, given that its CG predictions differ from the target by an amount *δ*(*q*). Importantly, this procedure can be iterated indefinitely by treating the posterior as a new prior and re-updating according to Eq. (4). See Fig. 2 for a flowchart of this correction loop. Within this loop, it is straightforward to extract an average * _{n}* from

In order to generate *transferable* forces, we propose to use this correction loop twice as follows (see also Secs. II B – II D). First, assume that (*q*) has been constructed. Second, use some preferred CG method (such as FM) to estimate *f*_{0} at a small set of anchor points *q _{k}* and correct these estimates according to the algorithm in Fig. 2. Third, interpolate the resulting force curves in state-space between the anchors, and again apply the correction algorithm to update forces pointwise. Note that the only significant difference between the correction loops amounts to how we compute an initial estimate

Process (solid black arrow) and information (dotted red arrow) flow associated with our calibration work flow. Rectangles encompassing many symbols group the corresponding work flow components by the section where they are discussed. Symbol colors have **...**

In Secs. II B–II D, we discuss in more detail the mathematics needed to use Eq. (4). Specifically, we outline procedure by which the ingredients *P*_{0}, , and *P*_{1} can be constructed for one iteration of the Bayesian correction and show results for the work flows in Figs. 2 and and3.3. Before continuing, however, a few comments are in order.

First, we can draw a close analogy between Eq. (4) and other well known optimization techniques. Specifically, in Sec. II C we will adopt the perspective that is the probabilistic counterpart of a Newton's method update equation of the form “
${f}_{n+1}-{f}_{n}=-{\scriptstyle \frac{\mathbb{S}[{f}_{n}]-\mathbb{T}}{{\mathbb{S}}^{\prime}[{f}_{n}]}}$,” where, informally, ′ is a derivative with respect to *f* and the quotes emphasize that this equation is not expressed rigorously. As Newton's method solves *g*(*x*_{0} + Δ*x*) = 0 (for arbitrary *g*) by computing corrections Δ*x* to an initial guess *x*_{0}, similarly ** weights changes Δ***f* according to the likelihood that [*f*_{0} +Δ*f, q*] − (*q*) = 0. In perfect correspondence, *P*_{0} is the probabilistic counterpart of an initial guess, and the product *P*_{1} ∞ *P*_{0} ensures that the posterior *P*1 is (loosely speaking) localized close the corresponding *f*_{0}. In Sec. II D, we make these ideas more precise by showing that, for specific choices of priors and likelihoods, the mean _{1} of *P*_{1} is the solution to a constrained optimization problem whose metric ||*f* − *f*_{0}|| is weighted according to the uncertainty in *f*_{0}.

Second, we assumed that an analytical or surrogate model (*q*) is already known. While the construction of such models is non-trivial, many well-known techniques have been developed to accomplish similar tasks; moreover, the formulation of (*q*) is often problem dependent. For these reasons, we leave the discussion of surrogate building for Appendix C. We note, however, that it is useful to express in the form

$$\mathbb{T}(q)=\stackrel{\sim}{\mathbb{T}}(q)+{\eta}_{\mathbb{T}}(q),$$

(5)

where is the average value of and *η*_{} (*q*) is a random process accounting for uncertainty through its variance
${\sigma}_{\mathbb{T}}^{2}=\text{Var}[{\eta}_{\mathbb{T}}(q)]$. Figure 4 provides an example of what Eq. (5) could look like for density-versus-temperature data. The yellow stripe indicates our level of un-certainty in the black curve, which is fit to atomistic simulation data (
). For the purposes of calibrating *f*, any force curve that yields CG predictions within the yellow stripe is statistically consistent , and thus acceptable. This highlights the usefulness of adopting a probabilistic framework inherent in Eq. (4): we avoid overcorrecting our force curves to match predictions that may be unreliable or have large uncertainties.

Example of a surrogate model (*q*) = *ρ*_{a}(*T* ) for atomistic density-temperature relations *ρ*_{a}(*T* ). The
denote atomistic simulation results, while the black curve is the hyperbola fit to the data. The inset shows the underlying probability **...**

Third, when the total number of modes *M* is finite, it is convenient to express *P*[*f; q*] = *P*[** w**;

Fourth, we have noted that Eq. (4) operates pointwise in *q*; that is, the Bayesian update is agnostic to whether *P*_{0} is constructed from CG forces *f*_{0} or interpolations. For this reason, we combine the discussions of interpolation and building priors in the next section. Gaussian processes, in particular, provide a natural means of accomplishing this task.

From a mathematical perspective, we find it useful to view a non-ideal (or finite-resource) CG algorithm as generating i.i.d. realizations from the statistical model

$${f}_{0}(r,q)={\stackrel{\u2323}{f}}_{0}(r,q)+{\eta}_{f}(r,q),$$

(6)

where _{0}(*r, q*) is a deterministic force curve coming from an ideal (or infinite-resource) GC algorithm, and *η _{f}* (

Assuming orthogonality of the basis functions, we formally decompose *f*_{0}(*r, q*) according to Eq. (1) to yield a set of random weights *w*_{0}(*q*) whose average values **E**[*w*_{0}(*q*)] = _{0}(*q*) reconstruct _{0}(*r, q*). We also assume that for fixed *q*, the functional form of each
${f}_{0}^{\{i\}}(r,q)$ is given on the entire domain 0 ≤ *r* ≤ 1, where *r* is a non-dimensionalized separation, and *r* = 1 is a cutoff distance beyond which particles do not interact;^{39} see Appendix A for details on non-dimensionalization.

In order to generate *P*_{0} for the first correction loop (cf. Fig. 3), we generate *N ≥* 1 realizations of
${f}_{0}^{\{i\}}$ at anchor points *q _{k}* in a vector

$${\stackrel{\sim}{\mathit{w}}}_{0}({q}_{k})=\frac{1}{N}\sum _{i=1}^{N}{\mathit{w}}_{0}^{\{i\}}({q}_{k})$$

(7)

(*which is random* !) and, when we can compute it, the sample covariance,

$$\text{Cov}[{\stackrel{\sim}{w}}_{m,0}({q}_{k}),{\stackrel{\sim}{w}}_{{m}^{\prime},0}({q}_{k})]\approx \frac{1}{N(N-1)}\sum _{i=1}^{N}\left[{w}_{m,0}^{\{i\}}({q}_{k})-{\overline{w}}_{m,0}({q}_{k})\right]\phantom{\rule{0.16667em}{0ex}}\left[{w}_{{m}^{\prime},0}^{\{i\}}({q}_{k})-{\stackrel{\sim}{w}}_{{m}^{\prime},0}({q}_{k})\right],$$

(8)

which quantifies correlations between modes and our uncertainty in _{0} when approximating this quantity with _{0}. Given Eqs. (7) and (8), we then postulate that a reasonable prior at each anchor point is

$${P}_{0}[\mathit{w};{q}_{k}]=\frac{1}{\sqrt{{(2\pi )}^{M}\mid \mathit{C}\mid}}exp\phantom{\rule{0.16667em}{0ex}}\left[\frac{-{[\mathit{w}-{\stackrel{\sim}{\mathit{w}}}_{0}({q}_{k})]}^{\mathrm{T}}{\mathit{C}}^{-1}[\mathit{w}-{\overline{\mathit{w}}}_{0}({q}_{k})]}{2}\right],$$

(9)

where ** C** = Cov[

In order to generate *P*_{0} for input into the second correction loop of Fig. 3, we must interpolate the weights in *q*. In general, we would like for *P*_{0} to account for additional uncertainty that arises when *q* is far from any of the anchor points *q _{k}*. To achieve this, we use GPR, a common tool in the machine learning and uncertainty quantification communities;

$${P}_{0}[\mathit{w};q]=\frac{1}{\sqrt{{(2\pi )}^{M}\mid \mathbf{\sum}\mid}}exp\phantom{\rule{0.16667em}{0ex}}\left\{\frac{-{[\mathit{w}-{\overline{\mathit{w}}}_{0}(q)]}^{T}{\mathbf{\sum}}^{-1}[\mathit{w}-{\overline{\mathit{w}}}_{0}(q)]}{2}\right\}$$

(10)

where **Σ** is a covariance matrix, |**Σ**| is the corresponding determinant, and the _{0}(*q*) are linear combinations of input weights (typically output by the first correction loop) at the anchor points; see also Eqs. (D5) and (D6) in Appendix D. In both Eqs. (9) and (10), the motivation for using Gaussian probabilities arises from the observation that realizations
${f}_{0}^{\{i\}}$ are themselves averages constructed from many i.i.d. samples of a thermodynamic ensemble. As such, it is reasonable to assume that a Central Limit Theorem (and hence a normal distribution) accurately describes the fluctuations of
${f}_{0}^{\{i\}}$ about their mean.

In the special case that the weights *w _{m}*(

$${P}_{0}[\mathit{w};q]=\prod _{m=1}^{M}P[{w}_{m};q],$$

(11)

$$P[{w}_{m};q]=\frac{1}{\sqrt{2\pi {\lambda}_{m}^{2}(q)}}exp\phantom{\rule{0.16667em}{0ex}}\left\{\frac{-{({w}_{m}-{\overline{w}}_{m,0}(q))}^{2}}{2{\lambda}_{m}^{2}(q)}\right\}.$$

(12)

While Eq. (11) neglects dependency between weights, we show in Sec. II D that this example provides an important connection between Eq. (4) and other well known optimization techniques.

In Fig. 5, we illustrate priors calculated according to Eq. (12) for an initial set of force curves *f*^{{}^{i}^{}} that were generated using the FM method. For eight temperatures **q** = (100 K; 200 K, … ,800 K), we generated 100 realizations
${f}_{0}^{\{i\}}(r)$, which were used to construct sample means and variances according to Eqs. (7) and (8); see the top sub-plots, where *ϕ* are the Chebyshev polynomials and the number of modes *M* = 10 (only 7 are shown). The bottom-left sub-plot shows estimates _{0} and their uncertainties according to Eq. (12), where we used the output of the first correction loop (iterated twice) as inputs to the GPR. In the bottom-right sub-plot, we reconstruct the mean force curves from these average weights according to Eq. (1). See Appendix D for more details on how these computations were performed.

We end this section with the caveat that, in general, there are many ways to construct reasonable priors. In essence, the task amounts to one of modeling, and this science is well developed; see, for example, Chapter 3 of Ref. ^{42}. Our analysis is meant to be exploratory, and it is likely that more advanced interpolation methods based on statistical mechanics (e.g. as in Ref. ^{15}) could lead to more suitable priors than we derive here. Nonetheless, we believe the usefulness of our approach here is born out by the results of the following sections and omit more in-depth studies as being beyond the scope of this work.

The key idea underlying the construction of our likelihood function is to invoke a linear correction equation based on Taylor's theorem and Newton's method. To motivate this formalism, fix *q*, denote an initial set of weights as _{0}, and temporarily assume that (*q*) is deterministic.^{43} Provided (*q*) ≠ [_{0}*, q*], we can estimate new weights ** w** via an equation of the form

$$\mathbb{T}(q)-\mathbb{S}[{\mathit{w}}_{0},q]=(\mathit{w}-{\overline{\mathit{w}}}_{0})\xb7{\nabla}_{\mathit{w}}\mathbb{S}[\mathit{w},q]{|}_{\mathit{w}={\overline{\mathit{w}}}_{0}}+\mathcal{O}({|\mathit{w}-{\overline{\mathit{w}}}_{0}\mid}^{2})$$

(13)

where **_{w}** [

In order to make use of Eq. (13), we must estimate the functional derivatives **_{w}** [

$${s}_{j}(f,q)=\underset{\epsilon \to 0}{lim}\frac{\mathbb{S}[f(r)+\epsilon {\varphi}_{j}(r),q]-\mathbb{S}[f(r),q]}{\epsilon}.$$

(14)

In general, it is impossible to evaluate this limit analytically, so we resort to approximating

$${s}_{j}(f,q)\approx \frac{\mathbb{S}[f(r)+\epsilon {\varphi}_{j}(r),q]-\mathbb{S}[f(r),q]}{\epsilon},$$

(15)

for *ε* 1. In practice, evaluating Eq. (15) simply amounts to running CG simulations with the inputs *f*(*r*) + *εϕ _{j}*(

Raw sensitivities and interpolations thereof. *Top plot:* Sensitivities for 8 temperatures as a function of mode number. The dotted lines are meant to aid the eye in seeing trends; they do not indicate values of sensitivities between integer modes. Note **...**

When computing sensitivities for the second correction loop of Fig. 3, it may be necessary to interpolate their values between a grid of anchor points. We omit discussion of this interpolation, as it can be done by the same steps that lead to the probabilistic models of *w*_{0} and ; see Appendix D. The bottom sub-plot of Figure 6 shows an example of interpolated sensitivities for our running example of density-versus-temperature relations. Note that while we use the same anchor points as in Sec. II B, this is done entirely for convenience. In principle any set of anchors is acceptable so long as they are chosen where one has estimates of _{0}. As in the case of (*q*), it is useful to express the sensitivities via an expression of the form

$$\mathit{s}(q)=\overline{\mathit{s}}(q)+\mathbf{D}{\mathit{\eta}}_{s}(q),$$

(16)

where *η** _{s}*(

In order to construct a likelihood function, we first re-write Eq. (13) by replacing deter- ministic quantities with their probabilistic counterparts. Defining Δ*w*_{0} := ** w** −

$$\delta -\mathrm{\Delta}{\mathit{w}}_{0}\xb7\overline{\mathit{s}}(q)=-{\eta}_{\mathbb{T}}(q)+\mathrm{\Delta}{\mathit{w}}_{0}\xb7\mathbf{D}{\mathit{\eta}}_{s}(q),$$

(17)

where *δ* = (*q*) − [_{0}*, q*]; cf. Eq. (5). The above expression states that that we can effect a change *δ* in the CG predictions via Δ*w*_{0}, but subject to the fact that our target and sensitivities are known imprecisely. We identify a likelihood function by noting that the left- hand side of Eq. (17) is normally distributed with mean zero and variance
${\sigma}_{\mathbb{T}}^{2}+{(\mathrm{\Delta}{\mathit{w}}_{0}\xb7\mathbf{D})}^{2}$; viz

$$\mathcal{L}[\mathrm{\Delta}{\mathit{w}}_{0}\mid \delta ]:=P[\delta \mid \mathrm{\Delta}{\mathit{w}}_{0}]=\frac{1}{\sqrt{2\pi [{\sigma}_{\mathbb{T}}^{2}+{(\mathrm{\Delta}{\mathit{w}}_{0}\xb7\mathbf{D})}^{2}]}}exp\phantom{\rule{0.16667em}{0ex}}\left[\frac{-{(\delta -\mathrm{\Delta}{\mathit{w}}_{0}\xb7\overline{\mathit{s}})}^{2}}{2[{\sigma}_{\mathbb{T}}^{2}+{(\mathrm{\Delta}{\mathit{w}}_{0}\xb7\mathbf{D})}^{2}]}\right],$$

(18)

where we equate the likelihood with *P*[*δ |Δ**w*_{0}], the probability of effecting a change δ in the CG predictions conditioned on Δ*w*_{0}.

Given the prior probability in Eq. (10) and likelihood function in Eq. (18), Bayes theorem yields

$${P}_{1}[\mathit{w};q]\propto {P}_{0}[\mathit{w};q]\mathcal{L}[\mathrm{\Delta}{\mathit{w}}_{0}(q)],\propto \frac{1}{\sqrt{2\pi [{\sigma}_{\mathbb{T}}^{2}+{(\mathrm{\Delta}\mathit{w}\xb7\mathbf{D})}^{2}]}}exp\phantom{\rule{0.16667em}{0ex}}\left\{\frac{-1}{2}\left[\frac{{(\delta -\mathrm{\Delta}\mathit{w}\xb7\overline{\mathit{s}})}^{2}}{{\sigma}_{\mathbb{T}}^{2}+{(\mathrm{\Delta}\mathit{w}\xb7\mathbf{D})}^{2}}+\mathrm{\Delta}{\mathit{w}}^{\mathrm{T}}{\mathbf{\sum}}^{-1}\mathrm{\Delta}\mathit{w}\right]\right\}$$

(19)

where a constant of proportionality is required to ensure that *P*_{1} is normalized. Several comments are in order. For one, *δ*, **D**, and *σ*_{} are all functions of *q*. Second, the general form of Eq. (19) makes it difficult to analytically compute the mean *∫*
*w**P*_{1}[** w**] d

In practical applications, the inexpense of CG simulations may allow one to characterize the sensitivities ** s** sufficiently well so that we can approximate

$${P}_{1}[\mathit{w}]\propto exp\phantom{\rule{0.16667em}{0ex}}[{(\mathit{w}-{\overline{\mathit{w}}}_{1})}^{T}\mathbf{\Lambda}(\mathit{w}-{\overline{\mathit{w}}}_{1})]$$

(20)

$${\mathrm{\Lambda}}_{i,j}={[{\mathrm{\sum}}^{-1}]}_{i,j}+\frac{{s}_{i}{s}_{j}}{{\sigma}_{\mathbb{T}}^{2}}$$

(21)

$${\overline{\mathit{w}}}_{1}={\overline{\mathit{w}}}_{0}+\frac{\delta}{{\sigma}_{\mathbb{T}}^{2}}{\mathbf{\Lambda}}^{-1}\mathit{s},\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}{\overline{w}}_{1,i}={\overline{w}}_{0,i}+\delta \frac{{[\mathbf{\sum}\mathit{s}]}_{i}}{{\sigma}_{\mathbb{T}}^{2}+{\mathit{s}}^{\mathrm{T}}\mathbf{\sum}\mathit{s}},$$

(22)

where the _{1,}* _{i}* are computed with the aid of the formula

$${[{\mathbf{\Lambda}}^{-1}]}_{i,j}={\mathrm{\sum}}_{i,j}-\frac{{[\mathbf{\sum}\mathit{s}]}_{i}{[\mathbf{\sum}\mathit{s}]}_{j}}{{\sigma}_{\mathbb{T}}^{2}+{\mathit{s}}^{\mathrm{T}}\mathbf{\sum}\mathit{s}},$$

(23)

due to Miller.^{45} The mean weights _{1} should yield improved predictions when used in place of _{0} in the CG simulations.

We highlight several features of Eq. (22). For one, using the prior Eq. (12) reduces _{1} to the much simpler

$${\overline{w}}_{1,i}-{\overline{w}}_{0,i}=\delta \frac{{s}_{i}{\lambda}_{i}^{2}}{{\sigma}_{\mathbb{T}}^{2}+{(\mathit{s}\xb7\mathit{\lambda})}^{2}}$$

(24)

If we further suppose that a single sensitivity *s _{i}* or uncertainty
${\lambda}_{i}^{2}$ is much greater than all other terms, then the corresponding difference

$${L}_{2,\lambda}(\mathit{w}-{\overline{\mathit{w}}}_{0})=\sum _{i}\frac{{({[\mathit{w}-{\overline{\mathit{w}}}_{0}]}_{i})}^{2}}{{\lambda}_{i}^{2}}$$

(25)

subject to the constraint that Δ*w*_{0} · ** s** =

When ${\sigma}_{\mathbb{T}}^{2}$ is finite, Eq. (24) also optimizes Eq. (25) subject to the modified constraint

$$\mathrm{\Delta}{\mathit{w}}_{0}\xb7\mathit{s}=\delta \frac{{(\mathit{s}\xb7\mathit{\lambda})}^{2}}{{\sigma}_{\mathbb{T}}^{2}+{(\mathit{s}\xb7\mathit{\lambda})}^{2}}\le \delta .$$

(26)

Again, limiting cases bring clarity to this equation. In particular, when
${\sigma}_{\mathbb{T}}^{2}\gg (\mathit{s}\xb7\mathit{\lambda})$, one finds that Δ*w*_{0} ≈ 0; i.e. the uncertainty in the atomistic results is too high to justify forcing the CG simulations to agree with the surrogate. On the other hand, when
${\sigma}_{\mathbb{T}}^{2}\ll (\mathit{s}\xb7\mathit{\lambda})$, the sensitivities and/or uncertainties **λ** are so large that we incur little inconsistency with the prior by requiring exact equality between the CG predictions and their mean atomistic counterparts.

Figure 7 shows the results of correcting force curves for our running example of density-temperature data first discussed in Figs. 5 and and6.6. We used our training data at **q** = (100 K*,* 200 K, … ,800 K) to estimate the mean weights and their uncertainty via Eqs. (7) and (8), which are inputs to the first correction loop in Fig. 3. Priors were computed according to Eq. (9), and sensitivities were recomputed for each iteration. After two passes through this loop, we interpolated the posterior mean forces at 10 K intervals and used these as inputs to the second correction loop in Fig. 3. We assumed independence of the modes to generate an interpolated prior according to Eq. (11) and set *μ*_{0} = 0 and *p* = 2, which are inputs to the GPR [cf. Eqs. (D1) and (D2) in Appendix D]. Approximating **D** = 0, we were able to generate well-calibrated forces after only one iteration through the second correction loop.

As is frequently the case with numerical optimization recipes, the Bayesian correction algorithm requires certain assumptions and inputs for practical implementation, e.g. hyper-parameters (i.e. inputs to GPR) and prior probabilities. As many of these quantities are not typically encountered in the context of coarse-graining, we explore the dependence of our method on its underlying assumptions (Secs. III A and III B). A key theme repeatedly appears: knowledge of the target data is more important than details of the prior. In this vein, we also demonstrate how Eq. (4) can accelerate coarse-graining by correcting initial force curves that are constructed by force matching undersampled atomistic simulations (Sec. III C).

Formally speaking, the posterior probability *P*_{1}[*f ; q*] that we recover from Eq. (4) should be independent of the basis set *ϕ*, provided (i) the latter is complete (in a mathematical sense), (ii) *M* → ∞, and (iii) the linear approximation in Eq. (13) is valid. From a practical standpoint, however, *P*_{1} will depend on the size of *M* and how well a given basis set captures relevant features of the CG forces. While it is beyond the scope of this work to present general statements about convergence, we do demonstrate that calibrating an initial set of force curves with different *ϕ* can lead to approximately the same CG predictions and posterior mean.

Here we consider two additional bases: Bessel functions with *M* = 15 modes and Legendre polynomials with *M* = 10 modes. For each set, we use our training force curves at **q** = (100 K, 200 K, … , 800 K) and update the modes according to the procedure underlying Figs. 3 and and7;7; see Sec. II D. In Fig. 8, we show the results of these computations. In both cases, we only require a few updates of the anchor points, followed by one update of the interpolated curves to yield good agreement with the target data. Thus for the three basis sets considered (including Chebyshev polynomials in Fig. 7), the algorithm converges with roughly the same amount of computational work.

In this section, we verify that corrections according to Eq. (4) are relatively stable to changes in the hyperparameters (i.e. underlying assumptions) of the GPR. We fix {*ϕ _{m}*} to be the Chebyshev polynomials, set

*μ*_{0}= 0 and*p*= 1,*μ*_{0}= 0 and*p*= 1.5,*μ*_{0}=*bT*+*c*and*p*= 2,*μ*_{0}=*aT*^{2}+*bT*+*c*and*p*= 2,

where *p* is the exponent appearing in Eq. (D2), *μ*_{0} is a prior assumption of the functional form of each weight as a function of temperature, and *a, b*, and *c* are free parameters given by a least squares fit to the posterior mean weights output by the first correction loop of Fig. 3.

Given these interpolation methods, we then perform the same steps as in the previous section and Sec. II D, interpolating weights as a function of *T*, computing sensitivities and updating force curves at 10 K increments. As before, we approximate **D** = 0. Figure 9 shows the results of this process. While each interpolation scheme yields somewhat different predictions, the interpolated forces all converge to acceptable targets within a few updates.

In this figure, it is interesting to note that the linear interpolation (*μ*_{0} = *bT* +*c*) behaves far worse than the others. While it is not clear *why* this is so, our work flow nonetheless succeeds in calibrating the underlying forces. This example points to the difficulty of interpolation and need for robust work flows. In such cases, a reasonable approach is to accept the *ad hoc* nature of a given choice and attempt to systematically account for and correct its effects.

In the previous section, we considered how arbitrary changes in the prior affect the correction algorithm. Here we examine more closely the special case in which we change *P*_{0} through the number *N* of force curves realized at each *q _{k}*. Specifically, we demonstrate that it is possible to find a reasonable by constructing priors with one tenth of the atomistic simulation resources used in the FM algorithm.

Figure 10 illustrates the results of this task. To generate this data, we followed the steps leading to Fig. 7, but using one tenth of the i.i.d. force curves
${f}_{0}^{\{i\}}$ available at each of the anchor temperatures when constructing the prior. The target data was kept the same. To arrive at convergence, we require more iterations of the update procedure applied to the interpolated forces. However, given the computational cost associated with running atomistic simulations, significantly less resources are required to calibrate *f* as compared to Fig. 7. In particular, each realization of
${f}_{0}^{\{i\}}$ takes roughly 40 CPU hours to generate, while each CG simulation takes ≈ 10 CPU hours. Thus, Fig. 7 required about 32,000 CPU hours to generate all the
${f}_{0}^{\{i\}}$ and roughly 4,500 CPU hours to run CG simulations, for a total of 36,500 CPU hours.^{46} In contrast, Fig. 10 required 3,200 CPU hours to generate
${f}_{0}^{\{i\}}$ fig and 8,100 CPU hours for CG simulations, for a total of about 11,300 CPU hours, less than a third of the case for Fig. 7.

Given the usefulness of CG simulations, it is not surprising that a vast amount of literature has been devoted to the subject of constructing and calibrating coarse-grained force fields; see, e.g. ^{7}^{, }^{10}^{, }^{11}^{, }^{15}^{–}^{18}^{, }^{20}^{, }^{22}^{, }^{32}^{, }^{33}, and ^{47}. To the best of our knowledge, however, none of these works implement uncertainty quantification as a widespread and integral part of their analysis. Here we wish to emphasize that UQ can play a central role in these analyses by (i) explicitly revealing the subjective aspects of calibration, (ii) recasting the latter as a decision-making problem, and (iii) providing systematic and robust methods for choosing calibrated forces.

To be more specific, the Bayesian framework illustrates that infinitely many plausible force curves yield the target predictions. This degeneracy (in a mathematical sense) is perhaps already foreshadowed by the variety of CG methods that have been successfully developed: empirical Leonard Jones (LJ) potentials,^{11,18,47} iterative Boltzmann inverseion (IBI),^{10,16,17} FM,^{5,6} etc.. But even within the context of a single method, statistical under-sampling, numerical error, and uncertainty in target data suggest that it is difficult, if not impossible to identify a unique, optimal force for a given coarse-graining procedure. Thus, any method for determining force curves involves a subjective choice by the modeler at some level.

Implementation of UQ throughout a CG procedure is useful in that it assesses the extent to which certain choices are more reasonable than others. This is illustrated by Eq. (25), wherein updates Δ** w** are computed by minimizing the norm

$${L}_{2,\lambda}(\mathrm{\Delta}\mathit{w})=\sum _{i}\frac{{({[\mathrm{\Delta}\mathit{w}]}_{i})}^{2}}{{\lambda}_{i}^{2}}.$$

(27)

Within the Bayesian framework, changing modes with small uncertainties would be heavily penalized, and hence implausible. We may interpret this as a compatibility statement: if a given weight output by the CG algorithm is very certain, a force curve with a very different value of this weight is inconsistent with the underlying assumptions or physics of coarse-graining. Thus, the simplest method of calibration (e.g. changing a single mode) could potentially be the least plausible.

Along related lines, our correction algorithm is robust, which is of particular importance when generating transferable potentials; see, e.g. Refs. ^{10}^{, }^{16}^{–}^{18}. These works showed that IBI typically outputs non-transferable forces, and other authors found that individually tailored methods (e.g. empirical potentials or specific interpolation methods) were required to overcome this limitation.^{11,15,18,20,22,47}. Figure 6 illustrates a possible culprit: large sensitivities to small uncertainties in a given mode can lead to wildly incorrect CG predictions. To highlight this issue, consider that for the Chebyshev basis at 800 K, an uncertainty of ±2 GN / mol in the first mode (approximately 2.5% of the weight, corresponding to a vertical translation of the force) leads to a variation of roughly ±0.25 g/cm^{3} (about ±30% !) in the density.^{48} Thus, small differences in reasonable interpolation methods can still lead to vastly different CG predictions. By updating forces according to a paradigm according to Eq. (4), we both (i) acknowledge this problem, while (ii) providing a systematic method for overcoming it; see Figs. 8–10.

Given its similarity to Newton's method, our interpretation and formulation of Eq. (4) lends itself to several extensions. In particular, decades of research have produced a variety of higher-order update equations that are known to have faster rates of convergence than Eq. (13); e.g. see Ref. ^{49}. Generally speaking, it should be straightforward to write expressions equivalent to our likelihood using any of these updates. In many instances, however, the corresponding algorithms may require running additional CG simulations in order to compute changes in sensitivities or higher order corrections to the linear approximation. Moreover, the resulting likelihoods may not be Gaussian, making it harder to analytically compute posterior means. Nonetheless, such extensions could provide promising improvements over the simple linear approximation that we invoke.

A second extension of this work is to enable calibration of multiple physical properties simultaneously. Within the context of the linear approximation, it should be relatively straightforward to solve such a problem. In particular, one can define a sensitivity vector *s** _{k}* for the

$${P}_{1}[\mathit{w},q]\propto {P}_{0}[\mathit{w},q]\prod _{k}{\mathcal{L}}_{k}[{\delta}_{k},q].$$

(28)

Our analysis is also easily generalized to the case in which the target data and sensitivities are both known precisely by reinterpreting likelihoods as Dirac delta functions. In a similar fashion, the above formulas can be generalized to multi-dimensional state points (e.g. temperature - pressure pairs) by recasting Eq. (13) in terms of the appropriate functional derivatives.

A key limitation of our analysis is the need to run a large number of CG simulations. Often times, these simulations are so inexpensive that the practical limitations of the method are negligible. However, in certain cases the cost of a CG simulation may be high enough so that, perhaps only (100) CG runs are available. In this case, the accuracy and/or efficiency of the method is limited by the available resources.

A related issue arises from the fact that, in general, CG simulations are not linear functions of the weights. Figure 7 illustrates the consequences of this; specifically, the correction algorithm sometimes over- or under-corrects the CG predictions, requiring additional iterations of the Bayesian update to fully optimize the forces. While this is a common occurence in numerical root finders that use, e.g. Newton's methods, such algorithms often benefit from the extreme inexpense [often less than (0.01) s per computation] of evaluating a function. Given that our function *is an entire simulation,* the impact of a poor guess *f*_{0} on efficiency cannot be overstated; see Sec. III C. As suggested in the previous section, higher-order methods could also mitigate this issue.

All authors thank Debra Audus, Timothy Burns, Andrew Dienstfrey, and Anthony Kearsley for useful comments during preparation of the manuscript. TR thanks Sergei Izvekov (Army Research Lab) for the use of his FM code implemented on the Garnet cluster located at the ERDC DSRC. Part of this work was completed while PP was a postdoctoral fellow at the Institute the Institute for Mathematics and its Applications (IMA) under a 506 grant from NIST to the IMA. The IMA is a National Science Foundation Math Institute 507 funded under award DMS-0931945.

The reference system used to construct all of the figures consisted of 10 chains of 50 monomers of atomistically resolved atactic polystyrene. To generate an initial configuration we used the amorphous builder within the commercial package MAPS by Scienomics^{50.51} All molecular dynamics simulations were performed with the open source package LAMMPS^{52}; temperature and pressure was maintained with a Nose-Hoover thermostat and barostat with coupling constants of 0.1 ps and 1.0 ps, respectively. The forcefield applied to the atomistic system was PCFF^{53} with a cutoffof 1.2 nm for both Lennard-Jones (van der Waals) and Coulombic nonbonded interactions. Long-range electrostatic interactions were calculated using the particle-particle particle-mesh method with a precision of 1.0 × 10^{−6}. Due to ease of implementation, the tabular force fields generated by the Bayesian calibration found in this work (e.g. Fig. 5) were calculated using dimensionless separations. Because very few interactions for CG sites in the atomistic reference system existed for distances less than 0.48 nm apart, this was defined as r = 0. To determine forces at separations less than 0.48 nm the first 600 values were fit to a linear function and extrapolated to 0 nm. The spacing between force values in real units was 5.0 × 10^{−5} nm. A real space cutoffof 1.4 nm was used in all CG simulations, corresponding to the dimensionless *r* = 1. Bond and bond angle potentials were derived using Inverse Boltzmann inversion at 500 K and applied to all temperatures.

To relax the initial atomistic configuration, we first minimized the energy of the amorphous structure using the Polak-Ribiere conjugate gradient method with a force tolerance of 1.0 × 10^{−8}. Four consecutive NVT simulations at 298 K were then performed for 5.0 × 10^{5} timesteps each, with a timestep of 0.001 fs, 0.01 fs, 0.1 fs, and 1.0 fs, respectively. All following simulations were performed with a timestep of 1.0 fs. A 0.5 ns NVT simulation that continuously increased the temperature from 298 K to 600 K was then performed. Finally, a 2 ns NPT simulation at 600 K was run to ensure an equilibrated structure. To obtain initial configurations for the higher temperatures investigated in this study, we further equilibrated for 0.5 ns at those temperatures. To obtain initial configurations for the lower temperatures we decreased the temperature continuously for 5 ns and then equilibrated for another 0.5 ns. After equilibration, a production run at each temperature was performed for 12 ns with coordinates, velocities, and forces recorded every 1 picosecond. This capture rate determines the number of frames used in the FM averaging of Eq. (B3) in Appendix B.

All CG simulations were performed in the NPT ensemble with a timestep of 1.0 fs. Initial configurations were generated from a direct mapping of the atomistic system at each temperature. Equilibration consisted of 100 ps, followed by a 10 ns production run.

The force matching method applied to coarse-grained systems has been described in detail elsewhere.^{5,6} Here we only outline specifics as they pertain to the initial guess for the CG force curves.

Force matching yields the optimal (in the least-squares sense) pairwise decomposition of the atomistic free-energy surface in CG coordinates. The result is a series of tabular functional-free force curves *f*^{{}^{i}^{}} (*r*, *T*) that we process via the work flow illustrated in Fig. 3.

Given an ensemble of atomic positions, momenta, and forces, ** r_{j}**,

$$\begin{array}{lll}\hfill {\mathit{R}}_{\mathit{J}}& =\hfill & \hfill \sum _{j=1,{a}_{J}}{c}_{jJ}{\mathit{r}}_{\mathit{j}},\\ \hfill {\mathit{P}}_{\mathit{J}}& =\hfill & \sum _{j=1,{a}_{J}}{\mathit{p}}_{\mathit{j}},\hfill \\ \hfill {\mathit{F}}_{\mathit{J}}^{\mathit{atm}}& =\hfill & \sum _{j=1,{a}_{J}}{\mathit{f}}_{\mathit{j}},\hfill \end{array}$$

(B1)

where *a _{J}* is the number of atoms in the

$${\mathit{f}}_{J}^{FM}(\mathrm{\Omega})=\sum _{K\ne J}{f}^{FM}({R}_{JK},\mathrm{\Omega}){\mathit{n}}_{JK},$$

(B2)

where *R _{JK}* = |

$${\chi}^{(\mathrm{\Omega})}=\langle \sum _{J=1}^{N}{\mid {\mathit{F}}_{J}^{\mathit{atm}}-{\mathit{f}}_{J}^{FM}(\mathrm{\Omega})\mid}^{2}\rangle ,$$

(B3)

subject to the constraint that the CG force field also reproduce the atomistic virial, viz.

$$\sum {f}^{FM}R=3V{P}^{\mathit{inst}}$$

(B4)

where *P ^{inst}* is the instantaneous pressure of the reference configuration at a volume

In general, there is no univeral technique for constructing a surrogate or target model; details often depend on the problem at hand. However, certain generic features are useful: (i) (*q*) should provide estimates of the quantity of interest at all state-points of interest; (ii) such estimates should also have a quantified uncertainty; (iii) (*q*) should be quick to evaluate, so that differences δ(*q*) are easy to compute.

In the case of simulated density-vs-temperature curves, we construct a surrogate by fitting a hyperbola to atomistic data.^{54} The motivation for representing the data via hyperbolas stems from the observation that, when a polymer system is deep within the glassy or rubbery regimes, density-temperature relationships *ρ*(*T*) should be approximately linear. The asymptotes of a hyperbola can capture this bilinear behavior, while its transition region can model intermediate temperatures. If we let *γ* denote the parameters of a hyperbola, we determine this quantity for our data by least-squares fit. For more details, see Ref. ^{54}.

Given a hyperbola parameterized by *γ*, we model uncertainty in the atomistic predictions by a simple stochastic process. Specifically, letting *ρ _{}*(

$${\rho}_{a}(T)={\rho}_{\mathscr{H}}(\gamma ,T)+{\eta}_{a}(T),$$

(C1)

where *η _{α}*(

$${\sigma}_{\mathbb{T}}^{2}=\frac{1}{Q}\sum _{k=1}^{Q}R{({T}_{k})}^{2}$$

(C2)

for *Q* anchor points. Given this parameter, we identify our target (*T*) = *ρ _{α}*(

We note in passing that this simple model neglects effects (e.g. correlations) that could be represented by more sophisticated models. Here we do not pursue such models since our focus is on developing tools for updating CG force curves.

In the UQ and machine learning communities, Gaussian-process regression is a well-known set of techniques that can interpolate data while simultaneously assessing confidence in the interpolated results. As there are many references that provide comprehensive reviews of the subject, we only introduce the main ideas and highlight issues that are relevant for constructing force-curves.^{29–31} We confine our discussion to a simple example of interpolating a single weight as a function of temperature, although the formulas we discuss are easily generalized to more complicated scenarios.

To begin, let **q** = (*q*_{1}, …, *q _{Q}*) be a vector of

$$P[w(\mathbf{q})]=\frac{1}{{(2\pi )}^{Q/2}\mid K(\mathbf{q},\mathbf{q})\mid}exp\phantom{\rule{0.16667em}{0ex}}\left\{-\frac{1}{2}{[w(\mathbf{q})-\mu (\mathbf{q})]}^{T}K{(\mathbf{q},\mathbf{q})}^{-1}[w(\mathbf{q})-\mu (\mathbf{q})]\right\},$$

(D1)

where *μ*(*q*) is the average behavior of 𝖜 as a function of *q* and *K*(*q, q′*) is a covariance matrix element indicating the extent to which 𝖜 (*q*) is correlated with 𝖜 (*q′*). Loosely speaking, the goal of GPR is to use the data (**q***, *(**q**)) to determine the the functions *μ* and *K*. In effect, this amounts to finding the stochastic process 𝖜 that is consistent with the available measurements.

Now, it is generally impossible for data to effect a regression by itself; we must invoke some additional or outside information. In the case of GPR, one postulates (i) a prior guess *μ*_{0}(*q*) for the mean, and (ii) a functional form of the covariance *K*(*q, q′;ψ*), parameterized by some unknown constants *ψ* The constants are often referred to as hyperparameters. For all of the regressions in this work, we use an altered form of the squared exponential covariance function

$$K(q,{q}^{\prime};\sigma ,\ell ,p)={\sigma}^{2}\sqrt{q{q}^{\prime}}exp\phantom{\rule{0.16667em}{0ex}}\left[\frac{-{\mid q-{q}^{\prime}\mid}^{p}}{{\ell}^{2}}\right],$$

(D2)

1 ≤ *p* ≤ 2, *σ*^{2}, and ^{2} are free parameters. The exponent *p* alters the smoothness properties of the underlying random functions, while the term
$\sqrt{q{q}^{\prime}}$ is motivated by the observation that MD predictions are noisier at higher temperatures [recall *K* is applied to our running example of *ρ*(*T*)]. The *σ*^{2} determines the scale of the overall noise and determines the temperature scale over which different observations are correlated. In this work, we manually control the values of *p* and *μ*_{0}(*q*), but determine *σ* and by means of a maximum likelihood analysis, which works as follows. First define Λ = diag (Var[(**q**)]) to be a diagonal matrix of the uncertainties in (**q**) [e.g. according to Eq. (8)]. Given this, we insert = *K* + Λ into Eq. (D1) and maximize the probability of finding our data **q***, *(**q**) as a function of these two parameters (assuming that *p* and *μ*_{0} are fixed). Including the term Λ accounts for the uncertainties in the sample estimates of (**q**) when computing the hyperparameters.

Given the covariance function *K*(*q, q′*), it is a straightforward matter to estimate *μ*(*q*) and Var[*w*(*q*)]. Specifically, define = (*q*_{1, …,}
*q _{Q}*

$$P(\stackrel{\u02cb}{\mathit{w}})=\frac{1}{\sqrt{{(2\pi )}^{1+Q}\mid \stackrel{\u02cb}{K}(\stackrel{\u02cb}{\mathbf{q}},\stackrel{\u02cb}{\mathbf{q}})\mid}}exp\phantom{\rule{0.16667em}{0ex}}\left\{-\frac{1}{2}{[\stackrel{\u02cb}{\mathit{w}}(\stackrel{\u02cb}{\mathbf{q}})-{\mu}_{0}(\stackrel{\u02cb}{\mathbf{q}})]}^{\mathrm{T}}\stackrel{\u02cb}{K}{(\stackrel{\u02cb}{\mathbf{q}},\stackrel{\u02cb}{\mathbf{q}})}^{-1}[w(\stackrel{\u02cb}{\mathbf{q}})-{\mu}_{0}(\stackrel{\u02cb}{\mathbf{q}})]\right\},$$

(D3)

where (*q _{k}, q_{k}*) =

$$P[w(q)]\propto exp\phantom{\rule{0.16667em}{0ex}}\left\{-\frac{{(w-{\mu}_{1}(q))}^{2}}{2\text{Var}[w(q)]}\right\},$$

(D4)

where

$${\mu}_{1}(q)={\mu}_{0}(q)+K(q,\mathbf{q})\phantom{\rule{0.16667em}{0ex}}{[K(\mathbf{q},\mathbf{q})+\mathrm{\Lambda}]}^{-1}\phantom{\rule{0.16667em}{0ex}}w(\mathbf{q})$$

(D5)

$$\text{Var}[w(q)]=K(q,q)-K(q,\mathbf{q}){[K(\mathbf{q},\mathbf{q})+\mathrm{\Lambda}]}^{-1}K(\mathbf{q},q).$$

(D6)

Note that this approach updates the prior *μ*_{0}(*q*) by taking into account correlations in the data relative to our initial guess for a functional form. Also, when Λ → ∞, note that *μ*_{1} → *μ*_{0}; that is, data with large uncertainties provide us with no information about how we should update our prior *μ*_{0}.

It is straightforward to generalize these results to the case where the ** w** are correlated across modes and/or additional information about uncertainties in our estimates of are brought to bear. We refer the reader to Refs.

This work is a contribution of the National Institute of Standards and Technology and is not subject to copyright in the United States.

1. Li C, Strachan A. J Polym Sci, Part B: Polym Phys. 2015;53:103.

2. Allison J, Collins P, Spanos G, editors. Proceedings of the 1st World Congress on Integrated Computational Materials Engineering (ICME); John Wiley & Sons, Inc; 2011.

3. Li M, Campbell C, Thornton K, Holm E, Gumbsch P, editors. Proceedings of the 2nd World Congress on Integrated Computational Materials Engineering (ICME); John Wiley & Sons, Inc; 2013.

4. Dienstfrey A, Phelan FR, Jr, Christensen S, Strachan A, Santosa F, Boisvert R. Journal of the Minerals Metals and Materials Society (JOM) 2014;66:1342.

5. Izvekov S, Voth G. J Chem Phys. 2005;123:134105. [PubMed]

6. Izvekov S, Voth G. J Phys Chem B. 2005;109:2469. [PubMed]

7. Moore TC, Iacovella CR, McCabe C. J Chem Phys. 2014;140:224104. [PubMed]

8. Bayramoglu B, Faller R. Macromolecules. 2012;45:9205.

9. Izvekov S, Chung PW, Rice BM. J Chem Phys. 2010;133:064109. [PubMed]

10. Bayramoglu B, Faller R. Macromolecules. 2013;46:7957.

11. Hsu DD, Xia W, Arturo SG, Keten S. J Chem Theory Comput. 2014;10:2514. [PubMed]

12. D'Adamo G, Pelissetto A, Pierleoni C. Soft Matter. 2012;8:5151.

13. Das A, Andersen HC. J Chem Phys. 2012;136:194114. [PubMed]

14. Jacobson LC, Kirby RM, Molinero V. J Phys Chem B. 2014;118:8190. [PubMed]

15. Krishna V, Noid WG, Voth GA. J Chem Phys. 2009;131:024103. [PubMed]

16. Carbone P, Varzaneh HAK, Chen X, Müller-Plathe F. J Chem Phys. 2008;128:064904. [PubMed]

17. Ghosh J, Faller R. Mol Simul. 2007;33:759.

18. Maerzke KA, Siepmann JI. J Phys Chem B. 2011;115:3452. [PubMed]

19. Brini E, Herbers CR, Deichmann G, van der Vegt NFA. Phys Chem Chem Phys. 2012;14:11896. [PubMed]

20. Farah K, Fogarty AC, Boehm MC, Müller-Plathe F. Phys Chem Chem Phys. 2011;13:2894. [PubMed]

21. Martin MG, Siepmann JI. J Phys Chem B. 1998;102:2569.

22. Qian H-J, Carbone P, Chen X, Karimi-Varzaneh HA, Liew CC, Müller-Plathe F. Macromolecules. 2008;41:9919.

23. Ganguly P, van der Vegt NFA. J Chem Theory Comput. 2013;9:5247. [PubMed]

24. Hsu DD, Xia W, Arturo SG, Keten S. J Chem Theory Comput. 2014;10:2514. [PubMed]

25. Hsu DD, Xia W, Arturo SG, Keten S. Macromolecules. 2015;48:3057.

26. The Bayesian interpretation of probabilities is distinct from ones based on ideas of chance and random draws; see, for example, Ref. ^{55}. For a good example of Bayesian calibration techniques applied to computer models, see Ref. ^{56} and the references contained therein. See Sec. II A for more discussion.

27. Das A, Andersen HC. J Chem Phys. 2012;136:194113. [PubMed]

28. Das A, Andersen HC. J Chem Phys. 2009;131:034102. [PubMed]

29. Rasmussen CE, Williams CKI. Gaussian Processes for Machine Learning (Adaptive Computation and Machine Learning) The MIT Press; 2005.

30. Santner TJ, Williams B, Notz W. The Design and Analysis of Computer Experiments. Springer-Verlag; 2003. p. 283.

31. Sacks J, Welch WJ, Mitchell TJ, Wynn HP. Stat Sci. 1989;4:409.

32. Potestio R, Peter C, Kremer K. Entropy. 2014;16:4199.

33. Peter C, Kremer K. Faraday Discuss. 2010;144:9. [PubMed]

34. Louis A. J Phys Condens Matter. 2002;14:9187.

35. Lu J, Qiu Y, Baron R, Molinero V. J Chem Theory Comput. 2014;10:4104. [PubMed]

36. Lobanova O, Avendao C, Lafitte T, Müller EA, Jackson G. Mol Phys. 2015;113:1228.

37. Several additional observations suggest that this perspective is appropriate: (i) there are many reasonable coarse-graining methods for computing *f*_{0}; (ii) it is often difficult to identify a best interpolant; and, as we show, (iii) many plausible force curves can yield the same target predictions. Thus, calibration of CG forces already requires subjective decision-making, which we only seek to acknowledge and systematize.

38. In the unusual case that *η*_{f} → 0, is known to arbitrary accuracy, one can take _{0}, which leads to Dirac delta-function priors following the procedure in this section.

39. Generally speaking, an arbitrary CG algorithm may yield *f*(*r, q*) at a finite, but large number of discretized points on the *r* domain. Here we assume that the spacing between such points is so small that the weights in Eq. (1) can be computed to arbitrary precision using numerical approximations of the appropriate projection integrals.

40. For IBI, it is reasonable that *N* = 1, so that we might need to invoke different methods for determining uncertainties computed in Eq. (8).

41. We use _{0} to denote the mean of a probability density in general so that we do not have to write separate versions of Eq. (4) for the first and second correction loops in Fig. 3.

42. Berger J. Springer Series in Statistics. Springer; 1985. Statistical Decision Theory and Bayesian Analysis.

43. Given that CG simulations are inexpensive, we assume that their outputs are always deterministic.

44. Equation (13) does not specify a unique *w* whenever *M* ≥ 2; in such cases, one could use the Moore-Penrose pseudoinverse to solve for *w*. The Bayesian update that we describe is a variant of this procedure.

45. Miller KS. Math Mag. 1981;54:67.

46. We approximate 900 CPU hours for each iteration of the CG simulations, be it an anchor point update, computation with interpolated forces, or computation with corrected interpolated forces.

47. Mognetti BM, Virnau P, Yelash L, Paul W, Binder K, Müller M, Mac-Dowell LG. J Chem Phys. 2009;130:044101. [PubMed]

48. In this light, it is perhaps no surprise that transferable forces are notoriously difficult to calibrate.

49. Stoer J, Bartels R, Gautschi W, Bulirsch R, Witzgall C. Introduction to Numerical Analysis, Texts in Applied Mathematics. Springer; New York: 2002.

50. “Scienomics, maps platform, version 3.4, 2014, France,”.

51. Certain commercial equipment, instruments, or materials are identified in this paper in order to specify the experimental procedure adequately. Such identification is not intended to imply recommendation or endorsement by the National Institute of Standards and Technology, nor is it intended to imply that the materials or equipment identified are necessarily the best available for the purpose.

52. Plimpton S. J Comput Phys. 1995;117:1.

53. Sun H, Mumby S, Maple J, Hagler A. J A Chem Soc. 1994;116:2978.

54. Patrone PN, Dienstfrey A, Browning AR, Tucker S, Christensen S. Submitted to Polymer. 2015

55. Gelman A, Carlin J, Stern H, Dunson D, Vehtari A, Rubin D. Bayesian Data Analysis. 3. Chapman & Hall/CRC Texts in Statistical Science, Taylor & Francis; 2013.

56. Kennedy MC, O'Hagan A. Journal of the Royal Statistical Society: Series B (Statistical Methodology) 2001;63:425.

PubMed Central Canada is a service of the Canadian Institutes of Health Research (CIHR) working in partnership with the National Research Council's national science library in cooperation with the National Center for Biotechnology Information at the U.S. National Library of Medicine(NCBI/NLM). It includes content provided to the PubMed Central International archive by participating publishers. |