Estimation of error in observables of coarse-grained models of atomic systems

Abstract

Background

The use of coarse-grained approximations of atomic systems is the most common methods of constructing reduced-order models in computational science. However, the issue of central importance in developing these models is the accuracy with which they approximate key features of the atomistic system. Many methods have been proposed to calibrate coarse-grained models so that they qualitatively mimic the atomic systems, but these are often based on heuristic arguments.

Methods

A general framework for deriving a posteriori estimates of modeling error in coarse–grained models of key observables in atomistic systems is presented. Such estimates provide a new tool for model validation analysis. The connection of error estimates with relative information entropy of observables and model predictions is explained for so-called misspecified models. The relationship between model plausibilities and Kullback-Leibler divergence between the true parameters and model predictions is summed up in several theorems.

Results

Numerical examples are presented in this paper involving a family of coarse-grained models of a polyethylene chain of united atom monomers. Numerical results suggest that the proposed methods of error estimation can be very good indications of the error inherent in coarse-grained models of observables in the atomistic systems. Also, new theorems relating the Kullback-Leibler divergence between model predictions and observations to measures of model plausibility are presented.

Conclusions

A formal structure for estimating errors produced by coarse-graining atomistic models is presented. Numerical examples confirm that the estimates are in agreement with exact errors for a simple class of materials. Errors measured in the D _KL -divergence can be related to computable model plausibilities. The results should provide a powerful framework for assessing the validity and accuracy of coarse-grained models.

Background

Coarse-grained-reduced order models

The most common method of constructing reduced-order models in all of computational science involves the use of coarse-grained models of atomic systems, whereby systems of atoms are aggregated into “beads”, or “super atoms”, or molecules to reduce the number of degrees of freedom and to lengthen the time scales in which the evolution of events are simulated.

The use of coarse-grained (CG) approximations has been prevalent in molecular dynamics (MD) simulations for many decades Comprehensive reviews of a large segment of the literature on CG models was recently given by Noid [1] and Li et al. [2], and an application to semiconductor nano-manufacturing is discussed in Farrell et al. [3]. The issue of central importance in developing CG models is the accuracy with which they approximate key features of the atomistic system. Many methods have been proposed to calibrate CG models so that they qualitatively mimic the all-atom (AA) systems, but these are often based on heuristic arguments.

In this paper, we develop a posteriori estimates of error in CG approximations of observables in the AA system. We focus on standard molecular dynamics models of micro-canonical ensemble (NVE) thermodynamics, and we call upon the theory of model adaptivity and error estimation laid down in [4] and [5]. In this particular setting, new estimates are also obtained when the information entropy of Shannon [6] is used as a quantity of interest. This leads to methods for estimating CG-model parameters that involve the Kullback-Leibler divergence between probability densities of observables in the AA and CG systems.

In the final Results and discussion section of this presentation, we review several statistical properties of parametric models, including asymptotic properties of misspecified models and generalizations of the Bernstein-von Mises theorem advanced by Kleijn and van der Vaart [7]. There, the fundamental role of the Kullback-Leibler distance (the D _KL ) between the true probability distribution and the observations accessible by the model is reviewed. We present results in the form of theorems that relate the D _KL to measures of model plausibility that arise from Bayesian approaches to model selection. The relationships of the a posteriori estimates to the statistical interpretations are summarized in concluding remarks.

Preliminaries, conventions and notations

We generally approach the problem of developing computer models of large atomic systems through the use of any of several hardened molecular dynamics (MD) codes or through equivalent Monte Carlo approximations invoking the ergodic hypothesis. For a system of n atoms, the Hamiltonian is

$$ H\left(\mathbf{r}^{n}, \mathbf{p}^{n}\right) = \sum_{\alpha=1}^{n} \frac{m_{\alpha}}{2} \mathbf{p}_{\alpha} \cdot \mathbf{p}_{\alpha} + u(\mathbf{r}^{n}), $$

((1))

where r ⁿ={r ₁,r ₂,…,r _n } is the set of atomic coordinate vectors, p ⁿ={p ₁,p ₂,…,p _n } is the set of particle momentum vectors, m _α the atomic mass of the α−th atom and u(r ⁿ) is the potential energy or interaction potential. Then (r ⁿ,p ⁿ) defines a point or microstate in the phase space Γ _AA of the all atom (AA) model. In typical MD simulations, the potential is, for example, of the form

$$ u(\mathbf{r}^{n}) = V_{\text{bond}}(\mathbf{r}^{n}) + V_{\text{angle}}(\mathbf{r}^{n}) + V_{\text{dihedral}}(\mathbf{r}^{n}) + V_{\mathrm{non-bonded}}(\mathbf{r}^{n}) + V_{\text{coulomb}}(\mathbf{r}^{n}), $$

((2))

where

$$\begin{array}{@{}rcl@{}}{} V_{\text{bond}}(\mathbf{r}^{n}) = \sum^{N_{b}}_{i=1} \frac{1}{2} k_{r i} |\mathbf{r}_{i} - \mathbf{r}_{0i}|^{2}, \end{array} $$

((3a))

$$\begin{array}{@{}rcl@{}}{} V_{\text{angle}}(\mathbf{r}^{n}) = \sum^{N_{a}}_{i=1} \frac{1}{2} k_{\theta i} (\theta_{i} - \theta_{0i})^{2}, \end{array} $$

((3b))

$$\begin{array}{@{}rcl@{}}{}V_{\mathrm{dihedral}}(\mathbf{r}^{n}) = \sum\limits^{{N_{nb}-1}}_{i=1} \sum\limits_{j=1}^{N_{nb}} \frac{V_{ji}}{2} \left[ 1 + (-1)^{j-1} \cos({j} \phi_{i}) \right], \end{array}$$

((3c))

$$\begin{array}{@{}rcl@{}} V_{\text{non-bonded}}(\mathbf{r}^{n}) = \sum^{N_{nb}}_{i=1} \sum_{j>i} 4 \epsilon_{ij} \left(\left(\frac{\sigma_{ij}}{r_{ij}} \right)^{\alpha} - \left(\frac{\sigma_{ij}}{r_{ij}} \right)^{\beta} \right), \qquad r_{ij} \leq r_{2}, \end{array} $$

((3d))

$$\begin{array}{@{}rcl@{}}{} V_{\text{columb}}(\mathbf{r}^{n}) = \sum^{N_{q}-1}_{i=1} \sum^{N_{q}}_{j>i} 4 \epsilon_{0} \frac{q_{i} q_{j}}{r_{ij}}. \end{array} $$

((3e))

Here covalent bonds are represented by the harmonic potential (1a), changes in bond angles by (1b), torsional potentials by changes in dihedral angles (1c), Lennard-Jones non-bonded potentials by (1d), with r _ij =|r _i −r _j | and r _c the cut-off radius, (α,β) typically = (12,6), and Coulomb potentials between charges q _i at r _i and q _j at r _j (1e). These forms are typical of those implemented in popular MD codes, although several other common potentials could be added. The parameters of the potential model are given by the vector of physical coefficients: {k _i ,k _{θ i} ,V _{b l a c k j i} ,ϕ _i ,ε _ij ,σ _ij ,r _c ,… }. In general, atomic properties and values of parameters for the full all-atom system are supplied by systems calibrated using experimental data or quantum mechanics predictions (see, e.g. the OPLS data in [8,9]).

Given the Hamiltonian (1), Hamilton’s equations of motion are:

$$ \frac{\partial H}{\partial \mathbf{p}_{\alpha}}=\dot{\mathbf{r}}_{\alpha}, \qquad \frac{\partial H}{\partial \mathbf{r}_{\alpha}}=-\dot{\mathbf{p}}_{\alpha}, \qquad 1 \leq \alpha \leq n. $$

((4))

In MD, it is assumed that the atomic system evolves according to the laws of Newtonian mechanics, so we set $\mathbf {p}_{\alpha } = m_{(\alpha)}\dot {\mathbf {r}}_{\alpha }$, and the second Hamiltonian equation in (4) reduces to the system of equations

$$ m_{\alpha \beta}\ddot{r}_{\beta i}(t) + \partial_{\alpha i} u(\mathbf{r}^{n}(t))-f_{\alpha i}(t)=0,\, \qquad 1 \leq \alpha, \beta \leq n \;, \; 1 \leq i \leq 3, $$

((5))

where repeated indices are summed throughout their range, m _{α β} =m _(α) δ _{α β} is the mass of atom α, superimposed dots indicate time derivations, r _{β i} is the component of r _β in direction i, ∂ _{α i} =∂/∂ r _{α i} u(r ⁿ(t)) is the total interatomic potential of the system given, e.g., by (2), and f _{α i} (t) is the ith component of applied force on atom α at time t. We will add initial conditions, $\dot {r}_{\beta i}(0) = v_{\beta i}$, and $r_{\beta i}(0)=r^{0}_{\beta i}$, where v _{β i} and $r^{0}_{\beta i}$, for now, are assumed to be given.

Molecular dynamical equations of the form (5) are typical of those in standard MD codes that are numerically integrated with randomly-sampled initial conditions over time intervals to approximate systems with constant energy and fixed volume and fixed number of particles corresponding to so-called micro-canonical ensembles. Without loss in generality, we confine this development to such thermodynamic scenarios noting that straightforward extensions to, say, constant temperature settings, are covered by replacing (5) with appropriate “thermostat” models, such as the Langevin or Nose–Hoover formulations (see e.g. [10]). The general approach is then applicable to canonical ensembles and more general statistical thermodynamics settings.

A fundamental concept in the molecular theory of matter is that macroscopic properties actually observed in experiments, the so-called observables, must be the result of averages over a time interval τ of some phase function q(r ⁿ,p ⁿ) that depends on the phase-point positions (r ⁿ,p ⁿ) in phase space Γ _AA , as all measurements require a finite duration. Moreover, for thermodynamic systems in equilibrium, this average, denoted 〈q〉, must be independent of the starting time t ₀ and it must attain a value from an essentially infinite time duration. Thus, our goal in constructing the all-atom model is to compute observables of the form (cf [10,11]),

$$ \langle q \rangle = {\lim}_{\tau \rightarrow \infty} \tau^{-1} \int_{t_{0}}^{t_{0}+\tau} q(\mathbf{r}^{n}(t), \mathbf{p}^{n}(t)) \text{d} t. $$

((6))

Here we shall confine our attention to phase functions that depend only on the configurations of systems in thermodynamic equilibrium. Our quantities of interest are then written,

$$ Q_{\mathbf{r}} = {\lim}_{\tau \rightarrow \infty} \tau^{-1} \int_{t_{0}}^{t_{0}+\tau} q(\mathbf{r}^{n}(t)) \text{d} t. $$

((7))

In all but the simplest applications, it is impossible to solve the dynamic system (5) owing to its enormous size. Therefore, reduced-order models must be developed. The process involves aggregating groups of atoms into beads or molecules or “super atoms” so as to create a coarse-grained (CG) molecular model. The CG model has N coordinate vectors R ^N(t)={R ₁(t),R ₂(t),…,R _N (t)}, N<n; and the corresponding equations of motion are

$$ M_{AB}\ddot{R}_{B i}(t) + \partial_{A i} U(\mathbf{R}^{N}(t),\boldsymbol{\theta})-F_{A i}(t)=0,\, \qquad 1 \leq A, B \leq N \;, \; 1 \leq i \leq 3. $$

((8))

M _AB defining the CG mass matrix, ∂ _Ai =∂/∂ R _Ai ,U(·,·) the interaction potential energy of the CG system, θ a vector of parameters defining the CG model, and F _Ai (t) the ith component of applied force at bead A at time t. Initial conditions are $\dot {R}_{A i}(0) = V_{A i}$, and $R_{A i}(0)=R^{0}_{A i}$. The unknown parameters with a potential of the form (2) are denoted, for example,

$$ \boldsymbol{\theta} = (K_{Ri}, R_{0i}, K_{\theta i}, \theta_{0 i}, V_{0 i}, \varepsilon_{ii}, \sigma_{ii}, \cdots), $$

((9))

the notation, in analogy with (1), being chosen to indicate parameters of the CG model.

It is important to establish a kinematic and algebraic relation between coordinates of particles in the AA system and those in the CG system. A very large literature exists on various coarse-graining mapping schemes, and choices of the appropriate map from the AA to the CG system or vice versa are often based on heuristic methods (see e.g. [2]). Our general approach can be adapted to any such well-defined AA-to-CG or CG-to-AA map, but for definiteness, we describe one such family of mappings.

Let $\mathcal {J}_{A}$ be the index set of AA-coordinate labels of atoms aggregated into a single bead A with CG-coordinate vector R _A emanating from the origin to a reference point labeled A within the bead (e.g. R _A could be chosen to define the center of mass, $\mathbf {R}_{A} = \sum _{\alpha \in \mathcal {J}_{A}} m_{\alpha }\mathbf {r}_{\alpha }/M_{A}, M_{A} = \sum _{\alpha \in \mathcal {J}_{A}} m_{\alpha }$). Let a _{A α} (t) be a vector from the reference center of bead A to the end point of AA-coordinate vector $\mathbf {r}_{\alpha }(t), \alpha \in \mathcal {J}_{A}$. Let $G^{\cdot A}_{\alpha }$ be component of the n×N array,

$$ G^{\cdot A}_{\alpha}=\left\{ \begin{array}{ccc} 1 & \text{if} & \alpha \in \mathcal{J}_{A} \\ 0 & \text{if} & \text{otherwise}. \end{array}\right. $$

((10))

Then we have,

$$ \mathbf{r}_{\alpha}= \sum_{A=1}^{N} G^{\cdot A}_{\alpha} \left(\mathbf{R}_{A} + \mathbf{a}_{A \alpha} \right), \qquad 1 \leq \alpha \leq n. $$

((11))

Here we assume that each AA coordinate vector r _α belongs to only one bead identified with CG vector R _A , but this is not a necessary assumption. In the case of bonded systems in which r _α is associated with, say, two index sets $\mathcal {J}_{A}$ and $\mathcal {J}_{B}$, we simply choose either $\mathcal {J}_{A}$ or $\mathcal {J}_{B}$ as the representative of r _α and associate r _α with only one bead to avoid double counting.

These kinematical conventions are illustrated in Figure 1. The AA coordinates r _α of atoms assigned to molecular bead A remain with that bead throughout all possible motions of the CG system. The Boolean array $G^{\cdot A}_{\alpha }$ merely adjusts labels of the AA system to agree with labels assigned beads in the CG system.

Figure 1

A Lagrangian description of the motion of a molecular bead A and its associated AA atoms $(i,j,k,l) \in \mathcal {J}_{A}$.

Returning to (7), it is clear that the CG approximation of the quantities of interest (QoI) are of the form,

$$ Q_{\mathbf{R}}(\boldsymbol{\theta}) = {\lim}_{\tau \rightarrow \infty} \tau^{-1} \int_{t_{0}}^{t_{0}+\tau} q\left(G\left(\mathbf{R}^{N}(t),\boldsymbol{\theta}\right)\right) \text{d} t $$

((12))

where G is the AA-to-CG map defined in (11), where we denote r ⁿ=G(R ^N(t),θ), and where we specifically present the dependence of the QoI on the CG potential parameters θ.

Now it is obvious that the evolution of the CG system defined by the coordinate vectors R _A (t) satisfying (8) do not satisfy the “true” equations of motion (5). Indeed, if R ^N(t)={R ₁(t),R ₂(t),…,R _N (t)} is the set of N-vectors satisfying (8), they constrain the motion of the AA system via (11), so upon introducing the constrained motions into (5) we obtain the residual,

$$\begin{array}{@{}rcl@{}} \rho_{\alpha i} (\boldsymbol{\theta},t) &=& m_{\alpha \beta} \sum^{N}_{A=1} G^{\cdot A}_{\beta} \left(\ddot{R}_{Ai} + \ddot{a}_{A \beta i} \right) (t) \\ & & + \partial_{\alpha i} u \left(G(\mathbf{R}^{N}(t), \mathbf{a}^{N}(t)) \right) - \sum_{A=1}^{N} G^{\cdot A}_{\beta} \omega^{\beta}_{\alpha} F_{Ai}(t), \\ & & 1 \leq \alpha, \beta \leq n \;, \; 1 \leq i \leq 3, \end{array} $$

((13))

where $\sum _{\beta } \omega ^{\beta }_{\alpha }=1, \sum _{\alpha \in \mathcal {J}_{A}} \mathbf {f}_{\alpha } = \mathbf {F}_{A}$, and

$$\begin{array}{@{}rcl@{}} G(\mathbf{R}^{N}(t), \mathbf{a}^{N}(t)) &=& \left\{ \sum_{A} G^{\cdot A}_{1} (\mathbf{R}_{A} + \mathbf{a}_{A1}), G^{\cdot A}_{2} (\mathbf{R}_{A} + \mathbf{a}_{A2})\right. \\ & & \left.\dots, G^{\cdot A}_{n} (\mathbf{R}_{A} + \mathbf{a}_{An}) \vphantom{\sum_{A} G^{\cdot A}_{1} (\mathbf{R}_{A} + \mathbf{a}_{A1}), G^{\cdot A}_{2} (\mathbf{R}_{A} + \mathbf{a}_{A2})}\right\}, \\ &=& \mathbf{r}^{n}(t). \end{array} $$

((14))

On the left-hand side of (13), we have taken note that the residual depends on the CG-model parameters θ, not explicitly presented on the right for simplicity. The AA force f _α at coordinate α is a fraction $\omega ^{\beta }_{\alpha }$ of the CG force F _{A i} . In general, as a first-order approximation, one can take the vectors a _{A α} as time-independent constant vectors equal to their value in a given reference configuration $(\mathbf {a}_{A \alpha }(t) = \mathbf {a}_{A \alpha }(t_{0}) ; \dot {\mathbf {a}}_{A \alpha } = \ddot {\mathbf {a}}_{A \alpha }=\mathbf {0})$. Then the approximate residual results of the form

$$ \hat{\rho}_{\alpha i} (\boldsymbol{\theta},t) = G^{\cdot A}_{\alpha} m_{\alpha \beta} \ddot{\mathbf{R}}_{Ai} + \partial_{\alpha i} u \left(G(\mathbf{R}^{N}(t)), \mathbf{a}^{N}(t_{0}) \right) - f_{\alpha i}, $$

((15))

with $f_{\alpha i} = \omega ^{\beta }_{\alpha } G^{\cdot A}_{\beta } F_{Ai}$, and with repeated indices summed, 1≤α,β≤n,1≤A≤N.

It is possible to define a reverse or “push back” relationship that assigns to every AA coordinate vector r _α the CG coordinate vector R _A of the bead to which r _α belongs. Given r _α , set

$$ \mathbf{R}_{A} = G^{\cdot \alpha}_{A} (\mathbf{r}_{\alpha} - \mathbf{a}_{A \alpha}), \qquad 1 \leq A \leq N, $$

((16))

where $G^{\cdot \alpha }_{A}$ is the transpose of $G^{\cdot A}_{\alpha }$ (no sum on α).

Thus, for each time t∈[t ₀,τ], one can select a sample ω(t) of AA coordinates $\{ \mathbf {r}^{*}_{1}(t), \mathbf {r}^{*}_{2}(t), \dots, \mathbf {r}^{*}_{n}(t) \}$ and employ (16) to generate the corresponding CG coordinates $\{ \mathbf {r}^{*}_{1}(t), \mathbf {r}^{*}_{2}(t), \dots, \mathbf {r}^{*}_{n}(t) \}$. We use the simplified notation,

$$ \mathbf{R}^{N}(t) = G(\omega(t)) $$

((17))

to define the image of this sample in the CG system.

Methods

Weak forms of the dynamical problem

It is convenient to recast the molecular dynamics problem into a weak or variational form. Toward this end, we introduce the Banach spaces

$$ \begin{aligned} \mathcal{V} &= \left\{\vphantom{\int_{0}^{\tau} \left[ \ddot{\mathbf{r}}_{\alpha} \cdot \ddot{\mathbf{r}}_{\alpha} + \dot{\mathbf{r}}_{\alpha} \cdot \dot{\mathbf{r}}_{\alpha} + \mathbf{r}_{\alpha} \cdot \mathbf{r}_{\alpha} \right] \text{d} \tau < \infty} \mathbf{r}^{n}(t); \mathbf{r}_{\alpha} \in \Omega \subset \mathbb{R}^{3} \times H^{2}(0, \tau); 1 \leq \alpha \leq n; \right.\\ \left\| \mathbf{r}^{n}(\cdot) \right\|^{2} &=\left. \int_{0}^{\tau} \left[ \ddot{\mathbf{r}}_{\alpha} \cdot \ddot{\mathbf{r}}_{\alpha} + \dot{\mathbf{r}}_{\alpha} \cdot \dot{\mathbf{r}}_{\alpha} + \mathbf{r}_{\alpha} \cdot \mathbf{r}_{\alpha} \right] \text{d} \tau < \infty \right\}, \end{aligned} $$

((18))

and the semilinear and linear forms $\mathcal {B}: \mathcal {V} \times \mathcal {V} \rightarrow \mathbb {R}, \mathcal {F}: \mathcal {V} \rightarrow \mathbb {R}$, given by

$$ \mathcal{B} \left(\mathbf{r}^{n} ; \mathbf{v}^{n} \right) = \int_{0}^{\tau} \left[ m_{\alpha \beta} \ddot{\mathbf{r}}_{\alpha}(t) \cdot \mathbf{v}_{\beta}(t) \right. + \left. \partial_{\alpha i} u \left(\mathbf{r}^{n} (t)\right) v_{\alpha i} (t) \right] \text{d}t + I_{0}(\mathbf{v}^{n}), $$

((19))

$$\begin{array}{@{}rcl@{}}{} \mathcal{F} \left(\mathbf{v}^{n} \right) = \int_{0}^{\tau} f_{\alpha i}(t) v_{\alpha i}(t) \text{d}t + v_{\alpha i}(0) m_{\alpha \beta} V_{\beta} - \dot{v}_{\alpha i}(0) m_{\alpha \beta} r^{0}_{\beta}, \end{array} $$

((20))

where

$$ I_{0}(\mathbf{v}^{n}) = v_{\alpha i}(0) m_{\alpha \beta} \dot{r}_{\beta}(0) - \dot{v}_{\alpha i} (0) m_{\alpha \beta} r_{\beta}(0). $$

((21))

The notation $\mathcal {B(\cdot ; \cdot)}$ is intended to mean that $\mathcal {B(\cdot ; \cdot)}$ is possibly nonlinear in entries to the left of the semi-colon and linear in the entries to the right of it.

The problem of finding $\mathbf {r}^{n} \in \mathcal {V}$ such that

$$ \mathcal{B} \left(\mathbf{r}^{n} ; \mathbf{v}^{n} \right) = \mathcal{F} \left(\mathbf{v}^{n} \right) \qquad \forall \mathbf{v}^{n} \in \mathcal{V}, $$

((22))

is equivalent to (8) in the sense that every solution of (8) with appropriate initial conditions, satisfies (22), and any sufficiantly smooth solution of (22) satisfies (8).

The adjoint problem

Let

$$ \mathcal{B}' (\mathbf{r}^{n}; \mathbf{z}^{n}, \mathbf{v}^{n}) = {\lim}_{\theta \rightarrow 0} \theta^{-1} \left[ \mathcal{B} (\mathbf{r}^{n} + \theta \mathbf{z}^{n}; \mathbf{v}^{n}) - \mathcal{B}(\mathbf{r}^{n}; \mathbf{v}^{n}) \right] $$

((23))

and

$$ Q' (\mathbf{r}^{n}; \mathbf{v}^{n}) = {\lim}_{\theta \rightarrow 0} \theta^{-1} \left[ Q(\mathbf{r}^{n} + \theta \mathbf{v}^{n}) - Q(\mathbf{r}^{n}) \right], $$

((24))

where Q is a functional on , and both $\mathcal {B}'(\cdot ; \cdot)$ and Q ^′(·;·) are assumed to exist and be finite (i.e. $\mathcal {B}(\cdot ; \cdot)$ and Q(·) are Gateaux differentiable). Then the adjoint or dual problem associated with (22) is

Find $\mathbf {z}^{n} = \{ \mathbf {z}_{1}, \mathbf {z}_{2}, \dots, \mathbf {z}_{n} \} \in \mathcal {V}$ such that

$$ \mathcal{B}'(\mathbf{r}^{n}; \mathbf{z}^{n}, \mathbf{v}^{n}) = Q'(\mathbf{r}^{n}; \mathbf{v}^{n}) \qquad \forall \mathbf{v}^{n} \in \mathcal{V}. $$

((25))

Introducing (19) into (23) gives, after some manipulations,

$$\begin{array}{@{}rcl@{}} \mathcal{B}'(\mathbf{r}^{n}; \mathbf{z}^{n}, \mathbf{v}^{n}) &=& \int_{0}^{\tau} \left(m_{\alpha \beta} \ddot{z}_{\beta i} - H_{\alpha i \beta j} (\mathbf{r}^{n}(t)) z_{\beta j} \right) v_{\alpha i} \text{d}t \\ & & + \dot{m}_{\alpha \beta} z_{\beta i} ({\tau}) v_{\alpha i}(\tau) \\ & & - m_{\alpha \beta} \dot{z}_{\beta i} ({\tau}) v_{\alpha i}(\tau), \end{array} $$

((26))

where H _{α i β j} is the Hessian,

$$ H_{\alpha i \beta j}(\mathbf{r}^{n}(t)) = \frac{\partial^{2} u (\mathbf{r}^{n}(t)) }{\partial r_{\alpha i} \partial r_{\beta j}}. \qquad 1 \leq \alpha, \beta \leq n \;, \; 1 \leq i, j \leq 3. $$

((27))

Likewise, if $Q(\mathbf {r}^{n}) = \int _{0}^{\tau } q(\mathbf {r}^{n}(t)) \text {d}t$, then

$$ Q'(\mathbf{r}^{n}(t) ; \mathbf{v}^{n}) = \int_{0}^{\tau} \partial_{\alpha i} q(\mathbf{r}^{n}(t)) v_{\alpha i} \text{d}t. $$

((28))

Note that (25) is solved “backward in time;” the forward problem (22) is solved for r ⁿ(t), which determines the coefficients in $\mathcal {B}'(\mathbf {r}^{n}; \mathbf {z}^{n}, \mathbf {v}^{n})$ which marches the adjoint solution from t=τ to t=0. The dynamical system corresponding to (25) is:

$$ m_{\alpha \beta} \ddot{z}_{\beta i} + H_{\alpha i \beta j} \left(\mathbf{r}^{n}(t) \right) z_{\beta j}(t) = \partial_{\alpha i} q \left(\mathbf{r}^{n}(t) \right). $$

((29))

Theory of a posteriori estimation of modeling error

Let us now review the theory of a posteriori estimation of modeling error advanced in [4] and expanded in [5]. We consider an abstract variational problem of finding an element u in a topological vector space such that,

$$ \mathcal{B} (u ; v) = \mathcal{F}(v), \qquad \forall v \in \mathcal{V}, $$

((30))

where $\mathcal {B(\cdot ; \cdot)}$ is a semilinear form from $\mathcal {V} \times \mathcal {V}$ into and is a linear functional on . Problem (30) is equivalent to the problem of finding a solution u of the problem A(u)=F in the dual space $\mathcal {V}'$, where A is the map induced by $\mathcal {B(\cdot ; \cdot)} : \langle A(u),v \rangle = \mathcal {B} (u ; v) = \mathcal {F}(v) = \langle \mathcal {F}, v \rangle$, 〈·;·〉 denoting duality pairing in $\mathcal {V}' \times \mathcal {V}$. Assuming (30) is solvable for u, we wish to compute the value Q(u) of a functional $Q: \mathcal {V} \rightarrow \mathbb {R}$ representing a quantity of interest, or an observable of interest.

We assume that the semilinear form $\mathcal {B} (\cdot ; \cdot)$ and the functional Q(·;·) are three times Gateaux differentiable on with respect to u. In particular, the following limits exist (recall (23) and (24)),

$$ \left. \begin{array}{ccccc} \mathcal{B}'(u; w, v) & = & {\lim}_{\theta \rightarrow 0} & \theta^{-1} \left[ \mathcal{B}(u+\theta w, v) - \mathcal{B}(u, v) \right] \\ Q' (u;v) & = & {\lim}_{\theta \rightarrow 0} & \theta^{-1} \left[ Q(u+\theta v) - Q(u) \right ] \end{array} \right\}. $$

((31))

with similar definitions of higher-order derivatives, e.g. $\mathcal {B}^{\prime \prime }(u; w_{1}, w_{2},v)$, $\mathcal {B}^{\prime \prime \prime }(u; w_{1}, w_{2}, w_{3}, v)$, Q ^′′(u;w,v), Q ^′′(u;v ₁,v ₂,v ₃), etc. See [5] for details.

The adjoint problem associated with (30) and the quantity of interest Q consists of finding $z \in \mathcal {V}$ such that

$$ \mathcal{B}' (u; z, v) = Q' (u; v), \qquad \forall v \in \mathcal{V}. $$

((32))

Now let u ₀ be an arbitrary element selected in . The residual functional (or “residuum”) associated with u ₀ is defined as the semilinear functional $\mathcal {R}: \mathcal {V} \times \mathcal {V} \rightarrow \mathbb {R}$,

$$ \mathcal{R} (u_{0}; v) = \mathcal{F}(v) - \mathcal{B}(u_{0};v), $$

((33))

which, for each $u_{0} \in \mathcal {V}$, is a linear functional on .

Obviously, if u ₀=u, the solution of (30), $\mathcal {R}(u;v)=0 \;\; \forall v \in \mathcal {V}$. Thus, $\mathcal {R}(u_{0};v)$ describes the degree to which the vector u ₀ fails to satisfy the central problem (30).

We now recall the basic theorem in [5]:

Theorem 1.

Let the semilinear form $\mathcal {B}(\cdot ; \cdot)$ in (30) and the quantity of interest Q be three-times continuously Gateaux differentiable on . Let u ₀ be an arbitrary element of . Then the error in Q(u)produced by replacing u by u ₀ is given by:

$$ Q(u) - Q(u_{0}) = \mathcal{R}(u_{0};z) + \Delta $$

((34))

where Δis a remainder involving higher-order terms in e ₀=u−u ₀ and ε ₀=z−z ₀, z ₀ being an approximation of z.

An explicit form of Δ is given in the appendix.

If u ₀ is not an arbitrary vector taken from but is a solution of a surrogate problem approximating (30) (such as a coarse-grained model approximating an AA model), then it often happens that Δ is negligible compared to the residual. Then (34) reduces to the approximation,

$$ Q(u) - Q(u_{0}) \approx \mathcal{R}(u_{0};z). $$

((35))

This relation is the basis for many successful methods of a posteriori error estimation of both modeling error and numerical error. Whenever $\mathcal {B}(\cdot ;\cdot)$ is a bilinear form and Q(·) is linear, Δ≡0.

A–posteriori estimation of error in CG approximations

The CG approximations of the “ground truth” AA system are characterized by a parametric class $\mathcal {P}(\boldsymbol {\theta })$ of molecular dynamics models, one model corresponding to each choice of the vector θ in a space Θ of parameters defining the CG intermolecular potential U(R ^N(t),θ). For a given value of θ, observables of interest in states of thermodynamic equilibrium of the CG system are typically generated as averages of samples of the observables taken over subintervals [t _k ,t _k+1]⊂[0,τ], for a distribution of initial conditions (see, e.g. [10]).

If we employ the approximation (35) to the AA and CG models, then an estimate of the error in CG approximations of the observable is immediate. Let q(r ⁿ) be a phase function whose ensemble average 〈q〉 is an observable of interest, denoted Q _r as in (7). The CG approximation is Q _R (θ) and the error, given by (35), is

$$ \varepsilon(\boldsymbol{\theta}) = Q_{\mathbf{r}} - Q_{\mathbf{R}}(\boldsymbol{\theta}) \approx \mathcal{R}(\mathbf{R}^{N}(\boldsymbol{\theta}); \mathbf{z}^{n}), $$

((36))

where $\mathcal {R}(\mathbf {R}^{N}(\boldsymbol {\theta }); \mathbf {z}^{n})= \int _{0}^{\tau } \rho _{\alpha i} (\boldsymbol {\theta }, t) z_{\alpha i} (t) \text {d}t$, ρ _{α i} (θ,t) is the residual in (13) (or (15)), and z ⁿ is the solution to the corresponding adjoint problem, and is generally unknown. If Z ⁿ is an approximation of z ⁿ, then

$$ | \varepsilon(\boldsymbol{\theta}) | \leq | \mathcal{R}(\mathbf{R}^{N}(\boldsymbol{\theta}); \mathbf{Z}^{n}) | + C(\boldsymbol{\theta}) \| \mathbf{z}^{n} - \mathbf{Z}^{n} \| $$

((37))

with

$$ C(\boldsymbol{\theta}) = \underset{\mathbf{v}^{n} \in \mathcal{V}}{sup} \frac{\| \mathcal{R}(\mathbf{R}^{N}(\boldsymbol{\theta},t)) \|_{\mathcal{V}'}}{\| \mathbf{v}^{n} \|_{\mathcal{V}}} $$

((38))

$\| \cdot \|_{\mathcal {V}'}$ being the norm on the dual space $\mathcal {V}'$. The problem of error estimation thus reduces to one of developing efficient procedures to compute the residual (ρ) and to compute reasonable approximations of z ⁿ.

It is clear that a quantitative estimate such as (36) (or an approximation with z ⁿ replaced by Z ⁿ) could be a powerful tool for determining validity of the CG model or in designing validation experiments for CG models. In theory, it also provides a basis for selecting optimal parameters for a given model so as to manage ε(θ). We elaborate on this notion in the final part of the Results and discussion section.

Results and discussion

Numerical example: estimation of error in CG-approximation of a polyethelyne chain

We describe in this section an application of a poasterior error estimation described in the previous section, involving CG approximations of a well-known model of polyethylene. For the base “AA” model, we consider a united atom model of a polyethylene chain containing 200 CH2 (methyl) monomers, meaning we have aggregated hydrogen and carbon atoms into an “AA” bead for simplicity. The united-atom coordinates r _i define the locations of each particle on an r-axis, and the displacement is denoted u _i (t). As an additional simplification, we assume that the interatomic potential is characterized by harmonic bonds of the form (1a), with parameter k _l =k= 350 k C a l/m o l, bond length of l= 1.5 Å, and atomic mass m= 14.026 g r/m o l. Initially, each united atom is separated by bond length l, the initial velocities are zero, and the system is assigned an initial displacement field u _i (0)=f(r _i ), where $f(r_{i}) = 1.2 \text {e}^{-0.1r_{i}(0)}$. Under these conditions, the AA system (5) reduces to

$$ \boldsymbol{m}\ddot{\boldsymbol{u}} - \boldsymbol{k}\boldsymbol{u} = 0, \qquad \boldsymbol{u}(0)=\mathbf{u}_{0},\; \dot{\boldsymbol{u}}(0) = \mathbf{v}_{0}, $$

((39))

where m and k are the mass matrix and the stiffness matrix of the united atom system and are of the form

$$ \boldsymbol{m} = m\left[ \begin{array}{lll} 1 & & \\ & \ddots & \\ & & 1 \\ \end{array}\right], \quad \boldsymbol{k} = k\left[ \begin{array}{cccr} 2 & -1 & & \\ -1 & 2 & \ddots & \\ & \ddots & \ddots & -1 \\ & \ddots & -1 & 2 \end{array}\right]. $$

((40))

A family of CG approximations of this model is obtained by aggregating the atoms into beads, with CG models in distinguished by the number P of atoms per CG bead. The resulting CG system (8) is of the form

$$ \boldsymbol{M}\ddot{\boldsymbol{U}} - \boldsymbol{K}\boldsymbol{U} = \mathbf{0}, \qquad \boldsymbol{U}(0)=\mathbf{U}_{0},\; \dot{\boldsymbol{U}}(0) = \mathbf{V}_{0}, $$

((41))

with the mass of each bead set to M=P m and the bond stiffness K=α k/P; $\alpha \in \mathbb {R}^{+}$.

Upon solving (41) for the CG displacement trajectory U(t), we compute the residual trajectory

$$ \boldsymbol{\rho} = \boldsymbol{m}\ddot{\boldsymbol{U}}_{\text{CG}} - \boldsymbol{k}\boldsymbol{U}_{\text{CG}}, $$

((42))

where U _CG(t) is the projection Π U(t) onto AA atom locations.

The bilinear and linear forms described in (19)-(21) reduce, in this case, to

$$\begin{array}{@{}rcl@{}} \mathcal{B} (\boldsymbol{u} ; \boldsymbol{v}) &=& \int_{0}^{\tau} \boldsymbol{v}^{T} \left(\boldsymbol{m}\ddot{\boldsymbol{u}} - \boldsymbol{k}\boldsymbol{u} \right) \text{d}t -\boldsymbol{v}^{T}(0) \boldsymbol{m} \dot{\boldsymbol{u}}(0) -\dot{\boldsymbol{v}}^{T}(0) \boldsymbol{m} \boldsymbol{u}(0), \end{array} $$

((43))

$$\begin{array}{@{}rcl@{}} \mathcal{F} (\boldsymbol{v}) &=& -\boldsymbol{v}^{T}(0) \boldsymbol{m} \mathbf{v}_{0} -\dot{\boldsymbol{v}}^{T}(0) \boldsymbol{m} \mathbf{u}_{0}. \end{array} $$

((44))

and (25) yielding

$$\begin{array}{@{}rcl@{}} \mathcal{B}' (\boldsymbol{u} ; \boldsymbol{v}, \boldsymbol{z}) &=& \int_{0}^{\tau} \left(\boldsymbol{m}\ddot{\boldsymbol{z}} - \boldsymbol{k}\boldsymbol{z} \right)^{T} \boldsymbol{v} \text{d}t\\ &+& \left(\boldsymbol{m} \boldsymbol{z}(\tau)\right)^{T} \dot{\boldsymbol{v}}(\tau) - \left(\boldsymbol{m} \dot{\boldsymbol{z}}(\tau)\right)^{T} \boldsymbol{v}(\tau) \\ &=& \mathcal{Q}'(\boldsymbol{v}). \end{array} $$

((45))

As an example of a QoI, we take Q _r to be the locally-averaged displacement,

$$ \mathcal{Q}_{\bold{r}} = \int_{0}^{\tau} \zeta(t) \text{d}t; \; \; \zeta (t) = \frac{1}{N_{0}} \sum_{i \in \mathcal{N}} u_{i}(t); \; \; \mathcal{N} = \{i: x_{i} \leq \beta l ; \beta \in \mathbb{R}^{+}\}, \; \; N_{0} = card\; \mathcal{N}, $$

((46))

for which the strong form of the dual problem is

$$ \boldsymbol{m}\ddot{\boldsymbol{z}} - \boldsymbol{k}\boldsymbol{z} = \boldsymbol{q}, \qquad m \dot{\boldsymbol{z}}(\tau) = \boldsymbol{0}, \quad \boldsymbol{z}(\tau)=\boldsymbol{0},\, $$

((47))

where N ₀ is the number of united atoms considered in set and q is the vector defined such that $\mathcal {Q}(\boldsymbol {u}) = \boldsymbol {q}^{T}\boldsymbol {u}$. Given the QoI (46), q will be as a n×1 vector,

$$ q_{i} = \begin{cases} 1 & \text{if } x_{i} \leq \beta l \\ 0 & \text{otherwise } \end{cases} \qquad i=1, \cdots, n. $$

((48))

The residual function $\mathcal {R} (\cdot, \cdot)$ of (36) in this example is of the form

$$ \mathcal{R} (\boldsymbol{U}_{\text{CG}},\boldsymbol{z}^{n})= \int_{0}^{\tau} \eta_{t}(t) \text{d}t; \qquad \eta_{t}(t) = \sum^{n}_{i=1} {z}_{i}(t) \cdot {\rho}_{i}(t) \; \text{d}t. $$

((49))

The estimated error in the QoI is then

$$ \mathcal{E}_{\text{est.}} = \mathcal{R} (\boldsymbol{U}_{\text{CG}},\boldsymbol{z}^{n}) \approx \mathcal{Q}_{\mathbf{r}} - \mathcal{Q}_{\mathbf{R}}, $$

((50))

where $Q_{\mathbf {R}} = \int _{0}^{\tau } \eta _{t}(t)\text {d}t$ and then the exact error is

$$ \mathcal{E}_{\text{exact}} = \mathcal{E}_{\text{est.}} + \Delta, $$

((51))

Δ being the blackremainder in (34). Since the forms in (43)-(46) are linear in their respective arguments, the exact blackremainder Δ should be zero, but the error introduced by the numerical integration schemes employed generally leads to an additional numerical error Δ _{Δ t} ≠0. We employ a converted Runga-Kutta algorithm here to integrate (39), (41), and (47).

The results of the coarse-grained model for the case of P=4 are presented in Figure 2. Figure 2a shows the coarse-scale displacement U=U(t) at different times, obtained from solution of (41) over the time domain t∈[0,τ]. The local residual of Figure 2b is then computed from (42). Figure 2c shows the solution of z(t) at different times. It observed that the adjoint solution propagates in time in the opposite direction to the primal solution, (47) being integrated backward in time.

Figure 2

Numerical solutions of the coarse-grained model. Coarse-grained model of a polyethylene chain containing 200 monomers with 4 united atoms per coarse-grained bead, P=4, α=1, and β=4: (a) coarse-scale displacement, U(t); (b) local residual ρ(t); (c) adjoint solution, z(t); (d) time evolution in estimated and exact error in QoI, η _t .

It is known that in general, the solution of the base model is not available. However, in order to show the effectiveness of the method presented here, the equations of motion for the united atom system is also solved in this example. Having the solution of the united atom model, u(t), the evolution of the exact ζ and estimated η _t over time is shown in Figure 2d.

Numerical approximations to the exact error are compared with the estimated error for various CG approximation of the united atom model in Figure 3a. The estimated error $\mathcal {R} = \mathcal {R}(\boldsymbol {U}_{\textit {CG}}(\boldsymbol {\theta }); \boldsymbol {z}^{n})$ for various values of P, with θ=α k/P, are indicated in Figure 3b for α=1. The computed estimated error ($\mathcal {E}_{\text {est}}=\mathcal {R}$) versus the parameter α are indicated in Figure 3c.

Figure 3

Numerical approximations of estimated error. Different CG-mapping G defined by number P of united atoms per bead and different CG parameters corresponding to θ=K=α k/P: (a) evolution of η _t over time (α=1, β=4), (b) values of versus P (α=1, β=4), and (c) values of versus α (P=4, β=4).

In general, the solution of the base model is not available, but the effectiveness of the method presented here is determined by comparing the CG solutions with the exact united atom model. The exact ζ and estimated η _t over time are shown in Figure (2d).

Maximum entropy principle for atomic systems

Among features of the AA system that could qualify as quantities of interest, we consider a special measure of uncertainty content embodied in the so-called information entropy. In 1948, Shannon [6] introduced the concept of information entropy as a real-valued function H(p) of probability distributions (densities) p as a logical measure of uncertainty content in p that satisfied four rather straight forward “common-sense” desiderata (see also [12] for full details). For a discrete pdf p={p ₁,p ₂,…,p _n }, the entropy is defined by

$$ H(p) = - \sum_{i=1}^{n} p_{i} \log p_{i}, $$

((52))

and for a continuous density, $p \in L^{2} (\mathbb {R})$, we write

$$ H(p) = - \int_{\mathbb{R}} p(y) \log p(y) \text{d}y. $$

((53))

Given two probability densities p and q, with non-empty support of domains, the relative entropy between p and q is given by the Kullback-Leibler divergence,

$$\begin{array}{@{}rcl@{}} D_{KL}(p \| q) &=& \int_{\mathbb{R}} p(y) \log \frac{p(y)}{q(y)} \text{d}y \\ &=& H(p,q) - H(p), \end{array} $$

((54))

where H(p,q) ($=-\int _{\mathbb {R}} p \log q \text {d}y$) is the cross entropy and it is understood that $0 \log \frac {0}{0}=0$ and $0 \log \frac {0}{q}=0$.

Shannon’s principle of maximum entropy asserts that in the set of all possible probability distributions relevant to a random field, the correct probability p corresponds to the maximum entropy:

$$ H(p) = \underset{q \in \mathcal{P}}{max} \; H(q). $$

((55))

Errors in information entropy

The connection with the statistical mechanics characterization of the AA and CG models can be established by choosing as a quantity of interest the infinite–time average of the phase function q(r ⁿ(t)) over [0,∞]. For this we invoke the ergodic hypothesis,

$$\begin{array}{@{}rcl@{}} Q_{\mathbf{r}} &=& {\lim}_{\tau \rightarrow \infty} \tau^{-1} \int_{t_{0}}^{t_{0}+\tau} q(\mathbf{r}(t)) \text{d}t \end{array} $$

((56))

$$\begin{array}{@{}rcl@{}} &=& \int_{\Gamma} \rho(\mathbf{r}^{n}) q(\mathbf{r}^{n}) \text{d} \mathbf{r}^{n} \end{array} $$

((57))

$$\begin{array}{@{}rcl@{}} &=& \langle q \rangle \end{array} $$

((58))

ρ(r ⁿ) being the distribution function for the ensemble under study and Γ the corresponding phase space subdomain. The corresponding CG approximation is

$$\begin{array}{@{}rcl@{}} Q_{\mathbf{R}}(\boldsymbol{\theta}) &=& {\lim}_{\tau \rightarrow \infty} \tau^{-1} \int_{t_{0}}^{t_{0}+\tau} q(G(\mathbf{R}^{N}(t);\boldsymbol{\theta})) \text{d}t \end{array} $$

((59))

$$\begin{array}{@{}rcl@{}} &=& \int_{\Gamma} \rho(\mathbf{r}^{n}) q(G(\mathbf{R}^{N};\boldsymbol{\theta})) \text{d} \mathbf{r}^{n} \end{array} $$

((60))

where the notation G(R ^N;θ) represents the relation (11). Setting

$$ q(\mathbf{r}^{n}) = \log \rho (\mathbf{r}^{n}), $$

((61))

gives immediately

$$ Q_{\mathbf{r}} - Q_{\mathbf{R}}(\boldsymbol{\theta}) = D_{KL} \left(\rho(\mathbf{r}^{n}) \| \rho (G(\mathbf{R}^{N};\boldsymbol{\theta})\right), $$

((62))

where D _KL (·∥·) is the Kullback-Leibler divergence defined in (54). Thus, if z ⁿ is the equilibrium solution of (25) with $Q'(\bold {r}^{n} ; \bold {v}^{n}) = \int _{0}^{\tau } \partial _{\alpha _{i}} q(\mathbf {r}^{n}(t)) v_{\alpha i} \text {d}t$, then

$$ D_{KL}(\rho(\mathbf{r}^{n}) \| \rho(G (\mathbf{R}^{N};\boldsymbol{\theta}))) \cong \mathcal{R}(\mathbf{R}^{N} (\boldsymbol{\theta}); \mathbf{z}^{n}). $$

((63))

The specification (60) of the CG approximation (with ρ(r ⁿ) as opposed to ρ(R ^N(θ))) requires some explanation. In interpreting (60), one assumes the role of an observer who resides in the AA system and, instead of the true phase function q(r ⁿ), observes a corrupted version for each choice of θ constrained to reside only in microstates accessible by the CG-model. This is also the interpretation of the residual described in (13) and (15). It is also noted that the estimate (63) is reminiscent of the minimum relative entropy method suggested by Shell [13].

A fundamental question arises at this point: given estimates (36) or (63), is it possible to find a special parameter vector θ ^∗ that makes the error ε(θ ^∗)=0? This question is related to the so-called well-specification or missspecification of the CG model. We believe the answer to this question is generally “no.”

Model misspecification and statistical analysis

A fundamental concept in the mathematical statistics literature on parametric models is the notion of a well-specified model, one that has the property that a special parameter vector θ ^∗ exists that the model maps into the truth; i.e. the true observational data. If no such parameter exists, the model is said to be misspecified.

More generally, we consider a space of physical observables (in our case, the values of appropriate observables sampled from the AA model) and a set $\mathbb {M(\mathcal {Y})}$ of probability measures μ on . As always, a target quantity of interest $Q:\mathbb {M} \rightarrow \mathbb {R}$ is identified (e.g. Q(μ)=μ[X≥a], X being a random variable and a a threshold value). We seek a particular measure μ ^∗ which yields the “true” value of the quantity of interest Q(μ ^∗). We wish to predict Q(μ ^∗) using a parametric model $\mathcal {P}: \Theta \rightarrow \mathbb {M}(\mathcal {Y})$, Θ being the space of parameters. Again, if a θ ^∗∈Θ exists such that $\mathcal {P}(\boldsymbol {\theta }^{*})=\mu ^{*}$, the model is said to be well-specified; otherwise, if $\mu ^{*} \notin \mathcal {P}(\Theta)$, the model is misspecified. See, e.g., Geyer [14], Kleijn and van der Vaart [7], Freedman [15], Nickl [16]. In our model discussed in Section ‘Preliminaries, conventions and notations’, we seek a parameter θ ^∗ of the CG model such that ε(θ ^∗) of (36) is zero, an unlikely possibility for most choices of Q.

To recast the issue of error estimation into a statistical setting, we presume that our goal is to determine (predict) a probability distribution of a random variable, an observable q in the AA system, using a CG model , given a set y ₁,y ₂,⋯,y _n of iid (independent, identically-distributed) random variables representing samples y _i =q(ω _i ) ($\omega _{i} = \mathbf {r}_{i}^{n}$ is meant to denote a particular point in phase space). We denote by π(y _i |θ) the conditional probability density p of the distance between the random data y _i and the parameter-to-observation map d _i (θ),

$$ p(y_{i} - d_{i}(\boldsymbol{\theta})) = \pi (y_{i}|\boldsymbol{\theta}) ; \quad i = 1, 2, \cdots, n, $$

((64))

where π(y _i |θ) is the ith component of the likelihood function. The joint density of the data vector y ⁿ=y ₁,y ₂,⋯,y _n is then,

$$ \pi_{n}(y_{1}, y_{2}, \cdots, y_{n} | \boldsymbol{\theta}) = \pi (\mathbf{y} | \boldsymbol{\theta}) = \prod_{i=1}^{n} \pi (y_{i}, \boldsymbol{\theta}). $$

((65))

The log-likelihood function is

$$ L_{n}(\boldsymbol{\theta}) = \log \pi(\mathbf{y} | \boldsymbol{\theta}) = \sum_{i=1}^{n} \log \pi(y_{i} | \boldsymbol{\theta}). $$

((66))

Let π(θ) be any prior probability density on the parameters θ (computed, for instance, using the maximum entropy method of Jaynes [12], as described for CG models in [3]); then the posterior density satisfies,

$$ \pi_{n}(\boldsymbol{\theta}|\mathbf{y}) = \pi (y_{1}, y_{2}, \cdots, y_{n} | \boldsymbol{\theta}) \pi (\boldsymbol{\theta}) / Z(\boldsymbol{\theta}), $$

((67))

where $Z(\boldsymbol {\theta }) = \int _{\Theta } \pi (\mathbf {y} | \boldsymbol {\theta }) \pi (\boldsymbol {\theta }) \text {d}\boldsymbol {\theta }$ is the model evidence.

The following definitions and theorems follow from these relations:

• The Maximum Likelihood Estimate (MLE) is the parameter $\hat {\boldsymbol {\theta }}^{n}$ that maximizes L _n (θ):

  $$ \hat{\boldsymbol{\theta}}^{n} = \underset{\boldsymbol{\theta} \in \Theta}{\text{argmax}} \; L_{n}(\boldsymbol{\theta}). $$

  ((68))

• The Maximum A Posterior Estimate (MAP) is the parameter $\tilde {\boldsymbol {\theta }}^{n}$ that maximizes the posterior pdf:

  $$ \tilde{\boldsymbol{\theta}}^{n} = \underset{\boldsymbol{\theta} \in \Theta}{\text{argmax}} \; \pi_{n}(\boldsymbol{\theta} | \mathbf{y}). $$

  ((69))

• The Bayesian Central Limit Theorem for well-specified models under commonly satisfied smoothness assumptions (also called the Bernstein-von Mises Theorem [7,16,17]) asserts that

  $$ \pi_{n}(\boldsymbol{\theta} | \mathbf{y}) \overset{\mathcal{P}}{\rightarrow} \mathcal{N}(\boldsymbol{\theta}^{*}; \mathbf{I}^{-1}(\boldsymbol{\theta}^{*})), $$

  ((70))

  where convergence is convergence in probability, $\mathcal {N}(\boldsymbol {\mu }, \boldsymbol {\Sigma })$ denotes a normal distribution with mean μ and covariance matrix Σ, $\hat {\boldsymbol {\theta }}$ is the generalized MLE, and I(θ) is the Fisher information matrix,

  $$ I_{ij}(\boldsymbol{\theta}) = -\sum^{n}_{k=1} \left[ \frac{\partial^{2}}{\partial \theta_{i} \partial \theta_{j}} \log \pi({y}_{k} | \boldsymbol{\theta}) \right]_{\boldsymbol{\theta} = \boldsymbol{\theta}^{*}} $$

  ((71))

• Given a set of parametric models, $\mathcal {M} = \{ \mathcal {P}_{1} (\boldsymbol {\theta }_{1}), \mathcal {P}_{2}(\boldsymbol {\theta }_{2}), \cdots, \mathcal {P}_{m} (\boldsymbol {\theta }_{m}) \}$, the posterior plausibility of model j is defined through the applications of Bayesian arguments by (see [3,18])

  $$ \rho_{j} = \pi(\mathcal{P}_{j} | \mathbf{y}, \mathcal{M}) = \frac{\int_{\Theta_{i}} \pi (\mathbf{y} | \boldsymbol{\theta}_{j}, \mathcal{P}_{j}, \mathcal{M}) \pi (\boldsymbol{\theta}_{j} | \mathcal{P}_{j}, \mathcal{M}) \text{d} \boldsymbol{\theta}_{j} \pi (\mathcal{P}_{j} | \mathcal{M})} {\pi(\mathbf{y} | \mathcal{M})} $$

  ((72))

  with $\sum _{j=1}^{m} \rho _{j} = 1$, and the largest ρ _j ∈[0,1] corresponds to the most plausible model for data $\mathbf {y} \in \mathcal {Y}$.

Finally, we come to the case of misspecified parametric models in which $\mu ^{*} \notin \mathcal {P}(\Theta)$; i.e. no parameter θ ^∗ exists such that the truth $\mu ^{*} = \mathcal {P}(\boldsymbol {\theta }^{*})$. This situation, we believe, is by far the most common encountered in the use of CG models.

We remark that in the (rare?) case of a well-specified CG model, for any continuous functional $Q: \Theta \rightarrow \mathbb {R}$ and if Θ is compact, if θ ^∗ is the unique minimizer of Q and if

$$ \underset{\boldsymbol{\theta} \in \Theta}{\text{sup}} \; | Q(\boldsymbol{\theta}; y_{1}, y_{2}, \cdots, y_{n}) - Q(\boldsymbol{\theta}) | \overset{\mathcal{P}}{\rightarrow} 0, $$

((73))

as n→∞, then the sequence

$$ \hat{\boldsymbol{\theta}}^{n} = \underset{\boldsymbol{\theta} \in \Theta}{\text{argmin}} \; Q_{n} (\boldsymbol{\theta}; y_{1}, y_{2}, \cdots, y_{n}) $$

((74))

converges to θ ^∗ in probability as n→∞. This is proved in Nickl [16]. In particular, under mild assumptions on the smoothness of the log-likelihood L _n (θ),

$$ Q(\boldsymbol{\theta}^{*}) - Q(\boldsymbol{\theta}) = -D_{KL} (\pi(\cdot | \boldsymbol{\theta}) \; \| \; \pi(\cdot | \boldsymbol{\theta}^{*})) $$

((75))

D _KL (·∥·) being the Kullback-Leibler distance defined in (54). By Jensen’s inequality (see, e.g. [16]), Q(θ ^∗)≤Q(θ)∀θ∈Θ; i.e. θ ^∗ is the minimizer of Q.

The asymptotic results for the finite misspecified case is summed up in the powerful result of Kleijn and van der Vaart [7,19]: let g(y) denote the probability density associated with the true distribution μ ^∗. Then the posterior density π _n (θ|y) converges in probability to the normal distribution,

$$ \pi_{n}(\boldsymbol{\theta}| \mathbf{y}) \overset{\mathcal{P}}{\rightarrow} \mathcal{N} (\boldsymbol{\theta}^{\dagger}, \mathbf{V}^{-1}(\boldsymbol{\theta}^{\dagger})), $$

((76))

where

$$ \mathbf{V}_{ij}(\boldsymbol{\theta}) = -\mathbb{E}_{g} \left[ \frac{\partial^{2}}{\partial \theta_{i} \partial \theta_{j}} D_{KL} \left(\cdot \; | \; \pi(\cdot|\boldsymbol{\theta}) \right) \right]_{\boldsymbol{\theta} = \boldsymbol{\theta}^{\dagger}}. $$

((77))

Thus, the best approximation to g in $\mathcal {P}(\Theta)$ is the model with the parameter

$$ \boldsymbol{\theta}^{\dagger} = \underset{\boldsymbol{\theta} \in \Theta}{\text{argmin}} \; D_{KL} \left(g \| \pi (\cdot | \boldsymbol{\theta}, \mathcal{P}, \mathcal{M}) \right) $$

((78))

being a class of parametric models to which belongs.

It is easily shown that θ ^† is a maximum likelihood estimate, i.e. it maximizes the expected value of the log-likelihood relative to the true density g:

$$\begin{array}{@{}rcl@{}} \boldsymbol{\theta}^{\dagger} & = & \underset{\Theta}{\text{argmin}} \left[ \int_{\mathcal{Y}^{n}} g(\textbf{y}) \log g(\textbf{y}) \; d\textbf{y} - \int_{\mathcal{Y}^{n}} g(\textbf{y}) \log \pi(\textbf{y} | \boldsymbol{\theta}) \; d\textbf{y} \right] \\ & = & \underset{\Theta}{\text{argmin}} \left[ - \int_{\mathcal{Y}^{n}} g(\textbf{y}) \log \pi(\textbf{y} | \boldsymbol{\theta}) \; d\textbf{y} \right] \\ & = & \underset{\Theta}{\text{argmax}} \int_{\mathcal{Y}^{n}} g(\textbf{y}) \log \pi(\textbf{y} | \boldsymbol{\theta}) \; d \textbf{y} \\ & = & \underset{\Theta}{\text{argmax}} \mathbb{E}_{g} \left[ \log \pi(\textbf{y} | \boldsymbol{\theta}) \right], \end{array} $$

((79))

where the negative self-entropy $\int g \log g \; d\textbf {y}$ was eliminated since it does not depend on θ and therefore does not affect the optimization.

Plausibility- D _KL theory

Let us now suppose that we have two misspecified models, $\mathcal {P}_{1}$ and $\mathcal {P}_{2}$. We may compare these models in the Bayesian setting through the concept of model plausibility: if $\mathcal {P}_{1}$ is more plausible than $\mathcal {P}_{2}$, ρ ₁>ρ ₂. In the maximum likelihood setting, the model that yields a probability measure closer to μ ^∗ is considered the “better” model. That is, if

$$ D_{KL}(g \| \pi(\textbf{y} | \boldsymbol{\theta}^{\dagger}_{1}, \mathcal{P}_{1}, \mathcal{M})) < D_{KL}(g \| \pi(\textbf{y} |\boldsymbol{\theta}^{\dagger}_{2}, \mathcal{P}_{2}, \mathcal{M})), $$

((80))

it can be said that model $\mathcal {P}_{1}$ is better than model $\mathcal {P}_{2}$. The theorems presented here define the relationship between these two notions of model comparison.

However, Bayesian and frequentist methods fundamentally differ in the way they view the model parameters. Bayesian methods consider parameters to be stochastic, characterized by probability density functions, while frequentist approaches seek a single, deterministic parameter value. To bridge this gap in methodology, we note that considering parameters as deterministic vectors, for example θ ₀, is akin to assigning them delta functions as their posterior probability distributions, which result from delta prior distributions. In this case, the model evidence is given by

$$ \pi(\textbf{y} | \mathcal{P}_{i}, \mathcal{M}) = \int_{\Theta} \pi(\textbf{y} | \boldsymbol{\theta}, \mathcal{P}_{i}, \mathcal{M}) \delta(\boldsymbol{\theta} - \boldsymbol{\theta}_{0}) \; d\boldsymbol{\theta} = \pi(\textbf{y} | \boldsymbol{\theta}_{0}, \mathcal{P}_{i}, \mathcal{M}). $$

((81))

In particular, if we consider the optimal parameter $\boldsymbol {\theta }^{\dagger }_{i}$ for model $\mathcal {P}_{i}$, $\pi (\textbf {y} | \mathcal {P}_{i}, \mathcal {M}) = \pi (\textbf {y} | \boldsymbol {\theta }_{i}^{\dagger }, \mathcal {P}_{i}, \mathcal {M})$. We can take the ratio of posterior model plausibilities,

$$ \frac{\rho_{1}}{\rho_{2}} = \frac{\pi(\textbf{y} | \mathcal{P}_{1}, \mathcal{M}) \pi(\mathcal{P}_{1} | \mathcal{M})}{\pi(\textbf{y} | \mathcal{P}_{2}, \mathcal{M}) \pi(\mathcal{P}_{2} | \mathcal{M})} = \frac{\pi(\textbf{y} | \boldsymbol{\theta}_{1}^{\dagger}, \mathcal{P}_{1}, \mathcal{M}) \pi(\mathcal{P}_{1} | \mathcal{M})}{\pi(\textbf{y} | \boldsymbol{\theta}_{2}^{\dagger}, \mathcal{P}_{2}, \mathcal{M}) \pi(\mathcal{P}_{2} | \mathcal{M})} = \frac{\pi(\textbf{y} | \boldsymbol{\theta}_{1}^{\dagger}, \mathcal{P}_{1}, \mathcal{M})}{\pi(\textbf{y} | \boldsymbol{\theta}_{2}^{\dagger}, \mathcal{P}_{2}, \mathcal{M}) } O_{12}, $$

((82))

where $O_{12} = \pi (\mathcal {P}_{1}|\mathcal {M})/\pi (\mathcal {P}_{2}|\mathcal {M})$ is the prior odds and is often assumed to be one. With these assumptions in force, we present the following theorems.

Theorem 2.

Let (82) hold. If $\mathcal {P}_{1}$ is more plausible than $\mathcal {P}_{2}$ and O ₁₂≤1, then (80) holds.

Proof.

If $\mathcal {P}_{1}$ is more plausible than $\mathcal {P}_{2}$,

$$ 1 < \frac{\rho_{1}}{\rho_{2}} = \frac{\pi(\textbf{y} | \boldsymbol{\theta}_{1}^{\dagger}, \mathcal{P}_{1}, \mathcal{M})}{\pi(\textbf{y} | \boldsymbol{\theta}_{2}^{\dagger}, \mathcal{P}_{2}, \mathcal{M})} O_{12} \leq \frac{\pi(\textbf{y} | \boldsymbol{\theta}_{1}^{\dagger}, \mathcal{P}_{1}, \mathcal{M})}{\pi(\textbf{y} | \boldsymbol{\theta}_{2}^{\dagger}, \mathcal{P}_{2}, \mathcal{M})} $$

((83))

Equivalently, the reciprocal of the far right-hand side is less than one, so

$$ \log \frac{\pi(\textbf{y} | \boldsymbol{\theta}_{2}^{\dagger}, \mathcal{P}_{2}, \mathcal{M})}{\pi(\textbf{y} | \boldsymbol{\theta}_{1}^{\dagger}, \mathcal{P}_{1}, \mathcal{M})} < 0. $$

((84))

Since g(y) is a probability measure, it is always non-negative. Thus

$$ g(\textbf{y}) \log \frac{\pi(\textbf{y} | \boldsymbol{\theta}_{2}^{\dagger}, \mathcal{P}_{2}, \mathcal{M})}{\pi(\textbf{y} | \boldsymbol{\theta}_{1}^{\dagger}, \mathcal{P}_{1}, \mathcal{M})} < 0 \Rightarrow \int_{\mathcal{Y}^{n}}g(\textbf{y}) \log \frac{\pi(\textbf{y} | \boldsymbol{\theta}_{2}^{\dagger}, \mathcal{P}_{2}, \mathcal{M})}{\pi(\textbf{y} | \boldsymbol{\theta}_{1}^{\dagger}, \mathcal{P}_{1}, \mathcal{M})}\; d\textbf{y} < 0. $$

((85))

This can be expanded into

$$ \int_{\mathcal{Y}^{n}}g(\textbf{y}) \log \pi(\textbf{y} | \boldsymbol{\theta}_{2}^{\dagger}, \mathcal{P}_{2}, \mathcal{M}) \; d\textbf{y}- \int_{\mathcal{Y}^{n}}g(\textbf{y}) \log \pi(\textbf{y} | \boldsymbol{\theta}_{1}^{\dagger}, \mathcal{P}_{1}, \mathcal{M}) \; d\textbf{y}< 0, $$

((86))

which means

$$ - \int_{\mathcal{Y}^{n}}g(\textbf{y}) \log \pi(\textbf{y} | \boldsymbol{\theta}_{1}^{\dagger}, \mathcal{P}_{1}, \mathcal{M}) \; d\textbf{y} < - \int_{\mathcal{Y}^{n}}g(\textbf{y}) \log \pi(\textbf{y} | \boldsymbol{\theta}_{2}^{\dagger}, \mathcal{P}_{2}, \mathcal{M}) \; d\textbf{y}. $$

((87))

By adding the quantity $\int _{\mathcal {Y}^{n}} g \log g \; d\textbf {y}$ to both sides, the desired result (80) immediately follows.

This theorem demonstrates that if model $\mathcal {P}_{1}$ is “better” than model $\mathcal {P}_{2}$ in the Bayesian sense, it is also a “better” deterministic model in the sense of (80). However, the reverse implication requires much stronger conditions. The assertion (80) can be equivalently written as

$$ \int_{\mathcal{Y}^{n}} g(\textbf{y}) \log \frac{\pi(\textbf{y} | \boldsymbol{\theta}^{\dagger}_{2}, \mathcal{P}_{2}, \mathcal{M})}{ \pi(\textbf{y} | \boldsymbol{\theta}^{\dagger}_{1}, \mathcal{P}_{1}, \mathcal{M})} \; d \textbf{y} < 0. $$

((88))

For this inequality to hold, the relationship

$$ \frac{\pi(\textbf{y} | \boldsymbol{\theta}^{\dagger}_{2}, \mathcal{P}_{2}, \mathcal{M})}{ \pi(\textbf{y} | \boldsymbol{\theta}^{\dagger}_{1}, \mathcal{P}_{1}, \mathcal{M})} < 1 $$

((89))

does not necessarily need to be true for every point $\textbf {y} \in \mathcal {Y}^{n}$.

One perhaps naive way to proceed is to invoke the Mean Value Theorem: if $\left | \mathcal {Y}^{n} \right | < \infty$ and under suitable smoothness conditions, there exists some $\bar {\textbf {y}} \in \mathcal {Y}^{n}$ such that

$$ \int_{\mathcal{Y}^{n}} g(\textbf{y}) \log \frac{\pi(\textbf{y} | \boldsymbol{\theta}^{\dagger}_{2}, \mathcal{P}_{2}, \mathcal{M})}{ \pi(\textbf{y} | \boldsymbol{\theta}^{\dagger}_{1}, \mathcal{P}_{1}, \mathcal{M})} \; d \textbf{y} = \left| \mathcal{Y}^{n} \right| g(\bar{\textbf{y}}) \log \frac{\pi(\bar{\textbf{y}} | \boldsymbol{\theta}^{\dagger}_{2}, \mathcal{P}_{2}, \mathcal{M})}{ \pi(\bar{\textbf{y}} | \boldsymbol{\theta}^{\dagger}_{1}, \mathcal{P}_{1}, \mathcal{M})}. $$

((90))

Then, combining (88) and (90) yields,

$$ \left| \mathcal{Y}^{n} \right| g(\bar{\textbf{y}}) \log \frac{\pi(\bar{\textbf{y}} | \boldsymbol{\theta}^{\dagger}_{2}, \mathcal{P}_{2}, \mathcal{M})}{ \pi(\bar{\textbf{y}} | \boldsymbol{\theta}^{\dagger}_{1}, \mathcal{P}_{1}, \mathcal{M})} < 0. $$

((91))

Since $\left | \mathcal {Y}^{n} \right | > 0$ and g(y)>0,

$$ \log \frac{\pi(\bar{\textbf{y}} | \boldsymbol{\theta}^{\dagger}_{2}, \mathcal{P}_{2}, \mathcal{M})}{ \pi(\bar{\textbf{y}} | \boldsymbol{\theta}^{\dagger}_{1}, \mathcal{P}_{1}, \mathcal{M})} < 0 \Rightarrow \frac{\pi(\bar{\textbf{y}} | \boldsymbol{\theta}^{\dagger}_{2}, \mathcal{P}_{2}, \mathcal{M})}{ \pi(\bar{\textbf{y}} | \boldsymbol{\theta}^{\dagger}_{1}, \mathcal{P}_{1}, \mathcal{M})} < 1 \Rightarrow \frac{\pi(\bar{\textbf{y}} | \boldsymbol{\theta}^{\dagger}_{1}, \mathcal{P}_{1}, \mathcal{M})}{ \pi(\bar{\textbf{y}} | \boldsymbol{\theta}^{\dagger}_{2}, \mathcal{P}_{2}, \mathcal{M})} > 1. $$

((92))

If O ₁₂≥1,

$$ \frac{\pi(\bar{\textbf{y}} | \boldsymbol{\theta}^{\dagger}_{1}, \mathcal{P}_{1}, \mathcal{M})}{ \pi(\bar{\textbf{y}} | \boldsymbol{\theta}^{\dagger}_{2}, \mathcal{P}_{2}, \mathcal{M})} O_{12}> 1 \Rightarrow \frac{\rho_{1} }{\rho_{2}} > 1. $$

((93))

Thus $\mathcal {P}_{1}$ is more plausible than $\mathcal {P}_{2}$ for given data $\bar {\textbf {y}}$.

In summary, we have:

Theorem 3.

If $D_{\textit {KL}}(g \| \pi (\textbf {y} | \boldsymbol {\theta }^{\dagger }_{1}, \mathcal {P}_{1}, \mathcal {M})) < D_{\textit {KL}}(g \| \pi (\textbf {y} | \boldsymbol {\theta }^{\dagger }_{2}, \mathcal {P}_{2},\mathcal {M}))$ and if $\left | \mathcal {Y}^{n} \right | < \infty$ and if (90) holds, then there exists a $\bar {\textbf {y}} \in \mathcal {Y}^{n}$ such that $\mathcal {P}_{1}$ is more plausible than $\mathcal {P}_{2}$, given that O ₁₂≥1.

Conclusions

The formal structure of a posterior estimates for errors in quantities of interest in CG approximations of atomistic systems is given by (36) if the CG model is sufficiently close to the AA model in some sense, and this error depends upon the CG model parameter θ. Numerical experiments presented in Section ‘Numerical example: estimation of error in CG-approximation of a polyethelyne chain’ involving a family of CG models of a polyethylene chain of united atom monomers suggest that these estimates can be very good indications of the error inherent in CG models of observables in the AA system.

In section ‘Errors in information entropy’ an example of special interest arises in the comparison of the information entropy of AA and CG models. This leads to estimates (62) and (63) involving the Kullback-Leibler divergence, D _KL .

When the CG model is misspecified in a statistical sense, which is generally the case, the “ D _KL -distance” between the AA truth and the best possible approximation of any CG model is defined by choosing θ=θ ^†, the minimizer of the D _KL as indicated in (78). Under special assumptions, one can relate the D _KL distance to Bayesian posterior model plausibility, as stated in our Theorem 2, which provides sufficient conditions for the most plausible model among a class of models to be in fact closest to the AA model in the D _KL sense. The possible role of estimates such as (36), (62), and (63) in model validation should be noted.

For each map G:AA→CG of the type defined by (10), a set of parametric model classes $\{ \mathcal {P}_{1} (\boldsymbol {\theta }_{1}), \mathcal {P}_{2} (\boldsymbol {\theta }_{2}), \cdots, \mathcal {P}_{m} (\boldsymbol {\theta }_{m}) \}$ is defined, each with undetermined and possibly random parameter vectors θ _i . For a calibration scenario, AA calibration data y _c ={y ₁,y ₂,⋯,y _n } are sampled, and a series of Bayesian updates is performed using an expanded form of Bayes’s rule that recognizes prior choices of the set and the class $\mathcal {P}_{j}$ within :

$$ \pi(\boldsymbol{\theta}_{j} | \mathbf{y}_{c}, \mathcal{P}_{j}, \mathcal{M}) \propto \pi (\mathbf{y}_{c} | \boldsymbol{\theta}_{j}, \mathcal{P}_{j}, \mathcal{M}) \pi (\boldsymbol{\theta}_{j} | \mathcal{P}_{j}, \mathcal{M}), \qquad 1 \leq j \leq m $$

((94))

The marginalization of the right-hand side of this relation is the model evidence, which serves as a likelihood function for a higher level of Bayes’s rule. The corresponding posteriors are the model plausibilities of (72). We remark that the notion of model plausibilities is an extension of the idea of Bayes factors prevalent in statistic literature (see e.g. [12] for discussion of the ideas) and was introduced to the best of our knowledge in [18]. The development of algorithms involving Bayesian plausibilities to study model selection in CG models of complex atomic system is discussed in [3,20].

It has been demonstrated, the most plausible model in a set will, under stated assumptions, involve parameters that minimize the D _KL −distance between the model and the so-called truth parameters. Whether that “best” model is valid for the intended purpose depends on tolerances set of error in key observables, the QoIs of the validation scenario (see [3]).

Appendix

A surrogate pair of equations approximating (30) and (32) may be embodied in the problem of finding the pair $(u_{0}, z_{0}) \in \mathcal {V} \times \mathcal {V}$ such that

$$ \left. \begin{array}{llll} \mathcal{B}_{0}(u_{0};v) & = & \mathcal{F}_{0}(v) & \forall v \in \mathcal{V} \\ \mathcal{B}'_{0}(u_{0};z_{0},v) & = & Q'_{0}(u_{0};v) & \forall v \in \mathcal{V} \end{array} \right\}. $$

((A-1))

The remainder Δ in (34) can, in this case, be shown (see [4]) to be given by:

$$\begin{array}{@{}rcl@{}} \Delta &=& \frac{1}{2} {\int_{0}^{1}} \left\lbrace \mathcal{B}^{\prime\prime} (u_{0}+\theta e_{0}; e_{0}, z_{0}+\theta \varepsilon_{0}) \right. \\ & & - \left. Q^{\prime\prime} (u_{0}+\theta e_{0}; e_{0}, e_{0}, \varepsilon_{0}) \right\rbrace \text{d} \theta \\ & & + \frac{1}{2} {\int_{0}^{1}} \left\lbrace Q^{\prime\prime\prime} (u_{0}+\theta e_{0}; e_{0}, e_{0}, \varepsilon_{0}) -3 \mathcal{B}^{\prime\prime} (u_{0}+\theta e_{0}; e_{0}, \varepsilon_{0}) \right. \\ & & - \mathcal{B}^{\prime\prime\prime} (u_{0}+\theta e_{0}; e_{0}, e_{0}, e_{0}, z_{0}+\theta) \\ & & - \left. \mathcal{B}^{\prime\prime\prime} (u_{0}+\theta e_{0}; e_{0}, e_{0}, e_{0}, z_{0}+\theta \epsilon_{0}) \theta (1-\theta) \text{d} \theta \right\rbrace, \end{array} $$

((A-2))

where

$$ e_{0} = u - u_{0} \qquad \text{and} \qquad \varepsilon_{0} = z - z_{0}. $$

((A-3))

The theory and estimates reduce to finite element a posterior error estimates in the special case in which u ₀=u _h and z ₀=z _h are finite element approximation of solutions (u,z) to partial differential equations (see e.g. [21]).