Using the probabilistic language of conditional expectations, we reformulate the force matching method for coarse-graining of molecular systems as a projection onto spaces of coarse observables. A practical outcome of this probabilistic description is the link of the force matching method with thermodynamic integration. This connection provides a way to systematically construct a local mean force and to optimally approximate the potential of mean force through force matching. We introduce a generalized force matching condition for the local mean force in the sense that allows the approximation of the potential of mean force under both linear and non-linear coarse graining mappings (e.g., reaction coordinates, end-to-end length of chains). Furthermore, we study the equivalence of force matching with relative entropy minimization which we derive for general non-linear coarse graining maps. We present in detail the generalized force matching condition through applications to specific examples in molecular systems.
I. INTRODUCTION
Complex molecular systems are materials of amazing diversity ranging from polymers to colloids, hybrid nanocomposites, biomolecular systems, etc., which are directly related with an enormous range of possible applications in nano-technology, bio-technology, food science, drug industry, cosmetics, etc. Due to the above reasons, molecular simulations of complex systems are a very intense research area.1 A main challenge in this field is to predict structure-property relations of such materials and to provide a direct quantitative link between chemical structure at the molecular level and measurable structural and dynamical quantities over a broad range of length and time scales.
On the microscopic (atomistic) level, detailed all-atom molecular dynamics (MD) or Monte Carlo (MC) simulations allow direct quantitative predictions of the properties of molecular systems over a range of length and time scales.1–4 However, due to the broad spectrum of characteristic lengths and times involved in complex molecular systems, it is not feasible to apply them to large realistic systems or molecules of complex structure, such as multi-component biomaterials, polymers of high molecular weight, and colloids. On the mesoscopic level, coarse-grained (CG) models have proven to be very efficient means in order to increase the length and time scales accessible by simulations.1,3,5–22
CG (particle) models can be roughly categorized, based on the way they are developed, into two groups: (a) ad hoc or phenomenological CG models, such as simple bead spring or lattice ones, which are primarily used to study generic behavior (e.g., scaling properties) of complex systems but lack a link to specific systems.1 The interactions between the CG groups in these models are described through semi-empirical functional forms obtained from previous knowledge and with a lot of physical intuition. (b) Systematic CG models, which are usually developed by lumping groups of atoms into “superatoms” and deriving the effective CG interaction potentials directly from more detailed (microscopic) simulations. Such models are capable of predicting quantitatively the properties of specific systems and have been applied with great success to a very broad range of molecular systems (see, for example, Refs. 5–13, 15, 17–19, and 23–25 and references therein).
A main challenge in the later family of CG models is to develop rigorous atomistic to CG methodologies that allow, as accurate as possible, estimation of the CG effective interaction. With such approaches, the hierarchical combination of atomistic and CG models could be in order to study a very broad range of length and time scales of specific molecular complex systems without adjustable parameters, and by that they become truly predictive.11,14,15
The overall procedure of systematic coarse-grained modeling for a given system, based on detailed microscopic data, is shortly described through the following stages: (a) execution of microscopic (e.g., ab initio or atomistic) simulations on small model systems, i.e., usually a relatively small number of molecules with a rather low molecular weight is considered, (b) select a CG map (transformation from the atomistic to the CG description), (c) development of a CG effective interaction potential (or force field), (d) execution of the CG (e.g., MD, Langevin dynamics (LD), or MC) simulations, and (e) re-introduction of the missing atomistic degrees of freedom in the CG structures, in case the properties under study require atomistic detail. From all above stages, the development of the CG force field is the most challenging one. Indeed, an accurate estimation of the way CG “superatoms” interact to each other is a conditio sine qua non in order to understand the behavior and to (quantitatively) predict the properties of the specific complex molecular system under study.
Note that from a mathematical point of view, coarse-graining is a sub-field of the dimensionality reduction;4,26 there are several statistical methods for the reduction of the degrees of freedom under consideration, in a deterministic or stochastic model, such as principal component analysis, polynomial chaos, and diffusion maps.20 Here, we focus our discussion on CG methods based on statistical mechanics, which have been used extensively in the last two-three decades in the theoretical modeling of molecular systems across a very broad range of disciplines, from physics to chemistry and biology as well as in engineering sciences.
There exists a variety of methods that construct a reduced CG model that approximates the properties of molecular systems based on statistical mechanics. These methods usually consider the optimization of proposed parametric models using different minimization principles, that is, considering a pre-selected set of observables {ϕi, i = 1, …, k} and then minimizing (average) values over a parameter set Θ,
where μ(x) and μθ(x) are the atomistic and proposed CG Gibbs measures, respectively. Different methods consider different sets of observables, for example,
in structural based methods, the observable is the pair radial distribution function, g(r), related to the two-body potential of mean force (PMF) (see Section VI C), for the intermolecular interaction potential, and distribution functions of bonded degrees of freedom (e.g., bonds, angles, and dihedrals) for CG systems with intramolecular interaction potential;6,7,9,10,21,22
force matching (FM) methods5,16,27 considered as observable function, the force fj(x) = − ∇xjU(x), j = 1, …, N, for an N-particle system with interaction potential U(x), x ∈ ℝ3N;
- the relative entropy (RE)8,18,19 method employs the minimization of the relative entropy, Kullback-Leibler divergence,These methods, in principle, are employed to approximate a many-body potential describing the equilibrium distribution of CG particles observed in the simulations of atomically detailed models. The many-body potential is defined through the renormalization group map28 that is equivalent to the PMF29 in case the former is differentiable. The force-matching (or multiscale coarse graining (MSCG)) and the relative entropy are minimization methods that construct a best fit of a proposed CG potential for systems in equilibrium. The force-matching method determines a CG potential from atomistic force information through a least-square minimization principle to variationally project the force corresponding to the potential of mean force onto a force that is defined by the form of the approximate potential. The relative entropy approach produces optimal CG potential parameters by minimizing the relative entropy between the atomistic and the CG Gibbs measures sampled by the atomistic model. In addition to the above numerical methods, analytical works for the estimation of the effective CG interaction, based on integral equation theory, have been also developed.30 A brief review and categorization of parametrization methods at equilibrium are given in Ref. 17.
Besides all the above, a classical method for calculating free energy differences using arbitrary reaction coordinates is thermodynamic integration (TI) theory.31–34 Thermodynamic integration is based on writing free energy differences as the integral of free energy derivative and thus computing the derivatives (mean force) instead of directly the free energy.
The purpose of this work is (a) to formulate in the probabilistic language of conditional expectations, the force matching method. In turn, the conditional expectation formulation allows us (b) to reveal the connection of force matching with thermodynamic integration that provides a way to construct a local mean force in order to best approximate the potential of mean force when applying the force matching method and (c) to present, in probabilistic formalism, the equivalence of relative entropy and force matching methods which we derive for general nonlinear coarse graining maps. (d) Furthermore, we discuss structure-based (SB) CG methods in a conditional expectation framework, thus presenting a complete picture of the known numerical (many body) potential of mean force estimation methods for systems at equilibrium and their relation.
The probabilistic formalism gives a geometric representation of the force matching method, i.e., recasts the force matching as a projection procedure onto the space of coarse observables (we refer specifically to Figure 1). The novelty and advantages of our approach is that it allows us to define a generalized force matching minimization problem,
applicable for linear and nonlinear CG mapsξ : ℝ3N → ℝm. The force matching condition introduced
where Jξ(x) = Dξ(x)Dξt(x), Dξ ∈ ℝm×3N (Dξ)ij(x) = ∇xjξi(x), i = 1, …, m, j = 1, …, N ensures the best approximation of the PMF. As in thermodynamic integration, h(x) is called the local mean force. A more general result is available in Theorem 2. This elucidates the direct connection of the above discussed particle CG methodologies with the standard thermodynamic integration approaches.
Geometric representation of the force matching procedure. Projection of the observable h(x) over the set of feasible coarse observables, , that is a subset of all CG observables L2(μ; ξ).
Geometric representation of the force matching procedure. Projection of the observable h(x) over the set of feasible coarse observables, , that is a subset of all CG observables L2(μ; ξ).
The current work is directly related to previous works that concern linear CG mapping schemes.14,16 Here, we recast and extend these works in a probabilistic framework in order to present and compare the relative entropy and force matching methods. This approach allows us to generalize the methodology to nonlinear coarse-graining maps. In the case of a linear CG map, the local mean force is
which reduces to the result in Ref. 14, see the examples in Section V A. Notice that for linear CG maps, the last term in relation (1) vanishes. The proposed formula for the local mean force extends the works14,16 in two aspects: (a) to any non-linear CG map and (b) to the existence of a family of appropriate local mean forces h(x) in closed form.
Finally, we should note that in the above discussion, we focus on molecular systems at equilibrium. The development of CG force fields for systems under non-equilibrium conditions, based on the information theory and path-space tools of relative entropy is the subject of Ref. 35.
The structure of this work is as follows. In Section II, we introduce the atomistic molecular system and its coarse graining through the definition of the CG map and the n-body potential of mean force. The probabilistic formulation of the force matching method and the best approximation of the PMF are presented in Section III, while the force matching condition for approximating the PMF for any, linear or non-linear CG map, is given in Section IV. In Sections V A and V B, we calculate the analytic form of the local mean for examples of linear and non-linear CG maps in molecular systems. A second result of the current work is presented in Section VI where we prove that the relative entropy minimization and the force matching methods are equivalent, producing the same approximation to PMF up to a constant. Furthermore, for completeness, we present the structure based methods in Section VI C. We close with Section VII summarizing and discussing the results of this work.
II. ATOMISTIC AND COARSE-GRAINED SYSTEMS
Assume a prototypical problem of N (classical) molecules in a box of volume V at temperature T. Let x = (x1, …, xN) ∈ ℝ3N describe the position vectors of the N particles in the atomistic (microscopic) description, with potential energy U(x). The probability of a state x at temperature T is given by the Gibbs canonical measure,
where Z = ∫ℝ3Ne−βU(x)dx is the partition function, , and kB the Boltzmann constant. We denote f(x) the force corresponding to the potential U(x), that is, f : ℝ3N → ℝ3N,
For such a system, the n-body, n < N, PMF 29,36 is defined through
where g(n)(x1, …, xn) is the n-body distribution function,
and is the number density.
Coarse-graining is considered as the application of a mapping (CG mapping) ξ : ℝ3N → ℝ3M,
on the microscopic state space, determining the M(<N) CG particles as a function of the atomic configuration x. We denote by z = (z1, …, zM) any point in the CG configuration space ℝ3M and use the bar “” notation for quantities on the CG space. We call “particles” the elements of the microscopic space with positions xj ∈ ℝ3, j = 1, …, N and “CG particles” the elements of the coarse space with positions zi ∈ ℝ3, i = 1, …, M. We should note that a CG mapping ξ(x) does not necessarily map to three-dimensional CG particles, it can be considered in the more general form ξ : ℝ3N → ℝm, for any m ∈ ℕ, with m < 3N. This is the case, for example, when considering some reaction coordinates, like the bending angle and the end-to-end distance (see examples in Section V B).
Assuming a CG mapping, ξ(x), defined by the renormalization group map,28 the expectation (mean) value of any observable ϕ : ℝ3N → ℝ of the form ϕ(x) = ϕ(ξ(x)) can be written as
where 𝔼μ[ϕ] denotes the expectation of ϕ(x) with respect to the probability measure μ(dx),
is defined by
for all z ∈ ℝ3M, being the conditional expectation of the observable quantity ϕ(x) that represents the expectation of ϕ(x) with respect to the Gibbs measure μ(dx) for given z = ξ(x) fixed. For a complete mathematical formulation of conditional expectation, see Appendix A.
Based on the above property, the conditional Helmholtz free energy A(z) is thus defined such that the CG probability density is of Gibbs type, i.e.,
The conditional (M-body) PMF is directly related to the free energy through the reversible work theorem.29,36 In many works, the free energy and potential of mean force are used interchangeably. Here, we use the potential of mean force notation and write . We define the mean force FPMF : ℝ3M → ℝ3M corresponding to the PMF defined by (5), assuming it exists,
The calculation of the potential of the mean force is a task as difficult/costly as is calculating expectations on the microscopic space. Thus, one seeks an effective potential function that is easy to formulate and calculate and approximates as well as possible the PMF. This is the ultimate goal of all the methods (i.e., structural-based methods, force matching, and relative entropy minimization) for systems in equilibrium that we present in the following sections in detail. In all the above mentioned methods, one usually proposes a family of interaction potential functions in a parametrized, , θ ∈ Θ, or a functional form and seeks for the optimal that “best approximates” the PMF. We denote by
the equilibrium probability measure at the coarse grained configurational space for the given CG potential function , where is the CG partition function.
III. CONDITIONAL EXPECTATION AND FORCE MATCHING
A class of the parametrization methods discussed here are based on the observation of a vector field h : ℝ3N → ℝ3M,
from microscopic simulations or experimental observations, and the definition of an optimization problem in order to find an optimal estimator G∗(z) of h(x) as a function of configurations in the coarse space ℝ3M. The optimization problem is to find a G∗ : ℝ3M → ℝ3M such that the mean square error
is minimized, where ‖ ⋅ ‖ denotes the Euclidean norm in ℝ3M. In general, h(x) can be any observable quantity which, eventually, for the force matching method5 that we study, represents the atomistic force on the coarse particles in configuration x ∈ ℝ3N. Hence, the force matching method at equilibrium is seeking for the optimal force G∗(z) as the minimizer of the mean square error over a set of proposed CG forces G(z) at the coarse space.
We recast the force matching optimization problem, as proposed in Refs. 5 and 27, in probabilistic terms using the concept of conditional expectation and its interpretation as a projection on a subspace of observables. First, we present a well known result in probability theory, see Ref. 37. We denote by L2(μ) the space of mean square integrable vector fields with respect to μ(dx), i.e., L2(μ) = {h : ℝ3N → ℝ3M|∫ℝ3N‖h(x)‖2μ(dx) < ∞}. For a given CG map ξ : R3N → R3M, we denote
the space of observables having the properties: (i) are square integrable observables with respect to the Gibbs measure μ(dx) and (ii) are functions of the coarse variable ξ(x). Property (i) ensures the space has a geometry that allows an easy formulation of the concept of projections for functions, e.g., it is a Hilbert space. The later property (ii) is called a “(sub-) σ algebra” in the context of conditional expectations, see Appendix A for further information.
that is the collection of all proposed CG force fields, G(z), see Figure 1. The set may consist of non-parametrized or parametrized elements, i.e., a set of splines, the span of a truncated basis of L2(μ; ξ). When the minimization problem is over the solution is not necessarily the conditional expectation 𝔼μ[h|ξ], defined by relation (4), as the Lemma 1 states, since it is possible that 𝔼μ[h|ξ] ∉ ℰ, rather it is a for which relation (9) holds, see the schematic in Figure 1. In this case, we say that G∗(z), the projection of h(x) on , is a best approximation of the F(z) = 𝔼[h|ξ]. With the following theorem, we state a necessary and sufficient condition that the observed quantity h(x) should satisfy so that mean force FPMF(z) (6) is best approximated with the force matching method.
- if ,
- if , the minimizer of satisfies
The statement is a direct consequence of Lemma 1.
The h(x) that satisfies (10) is called local mean force, as in the thermodynamic integration theory. Theorem 1 suggests that h(x) should have a specific form which is not at all obvious. In thermodynamic integration theory, approaches such as blue moon sampling and generalizations33,34 provide expressions for the local mean force for which the derivative of the free energy (the mean force) is the conditional expectation of the local mean force with respect to the reaction coordinates (coarse variables), that is, FPMF(z) = 𝔼μ[h|z]. This observation in combination with the statement of Theorem 1 that the optimal solution of the force matching problem is the conditional expectation of the local mean force, suggest that we can use the expressions for the local mean force given in the thermodynamic integration literature, see Theorem 2, to establish well defined force matching problems for general linear and non-linear CG mappings ξ : ℝ3N → ℝ3M. Note that this extends the validity of the force matching method to non-linear coarse graining mappings as the local mean force is chosen appropriately. The purpose of Sec. IV is to provide closed form representations for the local mean force h(x).
IV. CONSTRUCTION OF THE LOCAL MEAN FORCE AND SYSTEMATIC FORCE MATCHING
In this section, we give a closed form of the local mean force h(x), appearing in the force matching problem (8) for which the mean force is best approximated, based on the statement of Theorem 1 and results from TI theory.31–34 We introduce the derived form of h(x) as the appropriate observable to be used in a force matching method implementation in order to best approximate the mean force. In thermodynamic integration, the goal is to calculate free energy differences for a given reaction coordinate using the derivative of the free energy, see Chapter 3 in Ref. 38. We think of the coarse grained variable ξ(x) as a reaction coordinate, even though in the later case, one does not necessarily consider coarse graining of the system. Then, we use the result that the derivative of free energy (the mean force) is given as the conditional expectation on ξ of a local mean force that has a specific form that we state and prove here for completeness.
Before we state the result, we introduce some notations and assumptions. We denote Dξ(x) the 3M × 3N matrix with block elements (Dξ)ij(x) = ∇xjξi(x), i = 1, …, M, j = 1, …, N and Jξ(x) = Dξ(x)Dξt(x) the Jacobian matrix of the transformation. For a matrix A, At denotes its transpose, detA the determinant, and A−1 its inverse. We assume that the map ξ is smooth and such that This assumption ensures that the Jacobian matrix of the transformation Jξ(x) is non-degenerate, i.e., detJξ(x)≠0 and its inverse Jξ−1(x) exists.
The assumption in Theorem 1 is that h(x) must satisfy
which, as the following theorem states, is not unique rather it is parametrized by a family of vector valued functions W : ℝ3N → ℝ3M×3N related to the coarse graining map.
From the practical point of view, through the above theorem, someone can numerically estimate the mean force and then, by integrating obtain the PMF, which is the one used in the CG simulations. However, the above theorem also states that the choice of the local mean force h(x), that is how we construct the total force for each CG particle that corresponds to the PMF, is not unique; nevertheless, the PMF is well defined. For different choices of W(x), we can consider various force matching minimization problems; however, the corresponding PMF is the same. Furthermore, some of the problems may be better than others, simpler, and cheaper to implement. At first glance, formula (12) seems complicated though a suitable choice of W(x) can introduce major simplifications. Note also that in the low temperature regime, where 1/β ≪ 1, the term is not contributing significantly and can be neglected.
A W(x) always exists, at least in the case of smooth coarse graining map that we consider, since choosing W(x) = Dξ(x), we have that GW(x) = Jξ(x) is invertible. In thermodynamic integration, a well studied choice is W(x) = Dξ(x)31,33 that we present in the sequel as a corollary of Theorem 2.
Note that the second term in (13) depends on the curvature ∇x ⋅ Jξ−1(x)Dξ(x) of the sub-manifold Ω(z) = {x : ξ(x) = z}.
The coarse graining maps that are mainly considered in the equilibrium parametrization methods, the force matching and the relative entropy minimization, are linear mappings ξ : ℝ3N → ℝ3M,
for which the corresponding curvature ∇x ⋅ Jξ−1(x)Dξ(x) = 0, since Dξ(x) = T, where T = [ζijI3]i=1,…,M,j=1,…,N is independent of x and I3 denotes the 3 × 3 identity matrix. The form of the local mean force is thus simplified as described in the following corollary.
Therefore, we are looking for W such that
where O3M×3N is the 3M × 3N matrix with zero entries. The above system of equations has non-trivial solution, i.e., a non-zero matrix W, if B is such that I3N − TtB is a singular matrix.
In Sec. V, we study representative examples of molecular systems and coarse graining mappings and show in detail the application of the results of the current section.
V. FORCE MATCHING FORMULATION FOR LINEAR AND NON-LINEAR CG MAPS
The subject of this section is to present analytically the form of the local mean force h(x) for specific examples of molecular systems and for linear and nonlinear CG mappings. Specifically, we study the coarse graining of an N particle system with linear CG mappings, and the bending angle and end-to-end distance in a three atom system that serve as toy examples for presenting the nonlinear coarse graining case. Moreover, through the examples, we verify that even for linear but complex CG, the typical choice of the local mean force in force matching will not lead to approximations of the PMF, see the example in Section V A 1 b. Note that additional examples are given in Appendix C.
A. N-particle system under linear coarse graining maps
Let us consider a microscopic system of N particles with masses mj, j = 1, …, N, and position vectors x = (x1, …, xN) ∈ ℝ3N. In Sec. V A 1 and Appendix C, we consider different linear coarse graining maps ξ(x) for which we derive explicit forms of the local mean force h(x). Define the linear mapping ξ : ℝ3N → ℝ3M by
for ζij ∈ ℝ such that for all i = 1, …, N. The corresponding matrix is the 3M × 3N matrix , where Tij are the 3 × 3 blocks,
and I3 denotes the 3 × 3 identity matrix.
1. Two CG particles coarse space
We consider an example where the coarse space consists of M = 2 CG particles. The corresponding coarse graining map is defined by
with the matrix
a. Case 1.
For the CG, where each particle is contributing only to one CG particle, we verify that a W such that h∗(x) = hW(x),
typically used in force matching, does exist based on condition (14).
We consider now the (complex) mapping where there is one particle that contributes to both CG particles. We prove that for a W such that h∗(x) = hW(x) does not exist and thus, one must construct a local mean force based on Corollary 2, an example is given in case V A 1 b. We provide the computation details in Appendix C.
b. Case 2.
Let W = T, then
Furthermore, if each particle is contributing only to one CG particle, the form of the local mean force simplifies to
Note that when , and Ci = {j : ζij≠0}, or equivalently, Ci = {j : particlej contributes to CG particle i}, then from relation (17), we have
denoting ; thus, only when all particles have equal mass mj = m, j = 1, …, N,
Note that in all the above examples with linear CG maps, the form of h(x) can also be written in the form
B. Force matching formulation and non-linear CG maps
In this section, we examine the application of the force matching method with examples where we consider that the coarse graining mapping corresponds to a reaction coordinate, that is, in principle, a nonlinear mapping ξ : ℝ3N → ℝm.
We borrow the example from Ref. 32, where the corresponding free energy differences and PMF were calculated explicitly using generalized coordinates. Here, we only consider the mapping to the reaction coordinate and a proper selection of W(x) appearing in (12), as is also remarked in Ref. 39 Section 4.4. In this example, the microscopic model is a single molecule consisting of three atoms. Let xj ∈ ℝ3, j = 1, 2, 3 denote the position vectors of the atoms, see Figure 2.
1. Bending angle
The coarse variable is the bending angle θ = < x1x2x3, see Figure 2,
where and Jξ(x) = Dξ(x)Dξt(x), with
for j = 1, 2, 3 and Jξ(x) = ‖Dξ(x)‖2 ∈ ℝ. Thus,
where ξ(x) is given by (18).
2. End to end distance
Let us now choose the end to end distance ‖ℓ13‖ as a coarse variable,
for which and Jξ = DξDξt = 2. Applying Corollary 1, we have
since ∇x ⋅ Dξ(x) = 6, following definition (19).
A coarse variable that is of interest in molecular systems is the end-to-end vector ℓ13 = x3 − x1 ∈ ℝ3 for which we present the associated local mean force in Appendix C 4.
VI. FORCE MATCHING AND INFORMATION-BASED PROJECTIONS
In this section, we show a strong link between coarse-graining viewed as minimization of relative entropy and CG derived from the force matching optimization principle in L2(μ) presented in Sections III and IV. We first start the discussion with a brief outline of the relative entropy minimization and continue with its relation to force matching. Finally, we include a brief description of structural based methods in order to provide a complete view of the methods for potential of mean force approximations in coarse graining.
A. Relative entropy
The relative entropy approach8,40 considers the minimization of the relative entropy functional
over a space of interaction potentials. The minimization problem is based on the properties of the relative entropy (a) for all probability measures μ, π and (b) if and only μ ≡ π. The relative entropy is a pseudo-distance between the microscopic Gibbs measure μ(x) ∝ e−βU(x)dx and a back-mapping of the proposed Gibbs measure at the CG space ,
associated with the proposed interaction potential , where
with γ any non-negative L1(ℝ3N) function, and
Recall that Ω(z) = {x ∈ ℝ3N : ξ(x) = z}. The measure ν(x|z) is a normalized conditional probability of sampling an atomistic configuration x given a CG configuration z (microscopic reconstruction). A mathematical formulation of microscopic reconstruction is presented in our work,41 while probabilistic reconstruction methodologies are proposed and tested in Refs. 7, 10, and 41–44.
The difference in relative entropy between μ(x) and is given as
where is the exact coarse grained measure (3), and μ(x|z) is the unique measure such that . Relation (22) shows that the difference is composed of two parts: (a) the error in the approximation of the exact Gibbs measure corresponding to the by , , and (b) the error in reconstruction, , that is the error in approximating μ(x|z) by ν(x|z).
In the relative entropy minimization method, as defined by Shell et al.8,18 γ(x) = 1 is assigning the same probability to all atomistic configurations x that map to the same z. The reconstruction measure is the uniform distribution ν(x|z) = 1/|Ω(z)|, where |Ω(z)| is the volume of the set Ω(z), and the error introduced is
Note that for this choice of reconstruction, the error does not depend on the proposed approximating potential , the error is constant for any . In the ideal case where γ(x) = μ(x), the reconstruction is considered exact, there is no reconstruction error since ν(x|z) = μ(x|z) and , and the minimization problem is equivalent to
In view of the last two observations, it is verified that the relative entropy minimization method, with uniform or exact reconstruction, is indeed approximating the potential of mean force since the minimization problem is equivalent to the
B. Relative entropy and force matching
The goal of the last part of this section is to compare the force matching method with the relative entropy minimization method through their relation to the PMF. The relative entropy is directly related with the PMF through relation (22), while the force matching method approximates the PMF if the local mean force h(x) is such that , as stated in Theorem 2.
As discussed in the previous section, a reasonable choice for the reconstruction is γ(x) = μ(x), the equilibrium Gibbs measure;16 thus,
In practice, this choice of γ(x) means that we sample from the Gibbs measure using constraints on z. One can easily check that the relative entropy for γ(x) = μ(x) is rewritten as
Based on the above equality and the properties of the relative entropy, we can see that the minimum value of is given when corresponding to the PMF , under the assumption that the reconstruction probability ν(x|z) is exact, i.e., .
With the following theorem, we compare the relative entropy minimization and the force matching methods under the assumptions that both approximate the PMF , in the sense discussed in Secs. III and IV, i.e., in force matching and γ(x) = μ(x) in relative entropy minimization.
Observing relations (23) and (24), we notice that the leading term at the relative entropy approach minimizes the average of the square of potential difference , i.e., it is an error. On the other hand, the force matching minimizes the average of , an error, where . Thus, assuming the minimization problems have unique optimal solutions, and for the relative entropy and force matching methods, respectively, these solutions differ by a constant.
C. Structural based parametrization methods
For completeness, in the last part of this section, we discuss structural-based or correlation-based methods for obtaining CG interaction potential, such as the inverse Boltzmann, direct,6 and iterative7,21 and the inverse Monte Carlo methods,22 using a conditional expectation framework. The relation between these methods and the potential of mean force is clear. Indeed, theoretically, if one can compute, from the microscopic simulations, the n-body correlation function in the CG state space, under the specific CG mapping scheme, , n < M,
then according to the relation,29
where z(n) = (z1, …, zn), the computation of is straightforward. Hence, the correlation-based methods, in principle, can provide exactly the (M-body) potential of the mean force, as is also the case of exact convergence for the relative entropy and force matching methods.
However, in practice, the computation of is not feasible for large n, and what is used at inverse Boltzmann and inverse Monte Carlo methods is the pair correlation
In homogeneous systems, depends on the relative position between two particles r = ‖z1 − z2‖, , called the radial distribution function,
that is, the average density of finding the CG particle 1 at a distance r from the particle 2.
Moreover, all structure based methods rely on Henderson’s uniqueness theorem,45 which states that for a given radial distribution function, provided the PMF can be written as a sum of pair potentials, v(r), such that
for z = (z1, …, zM) then, this potential is unique, up to a constant.
In practice, structure based methods look for an optimal pair potential v(r) as the solution of the minimization problem is
where
with . The error introduced with a structure based method is due to the assumption of pair potentials, that is,
which is analogous to the error introduced in a force matching method presented in relation (11).
Overall, the methods discussed in this section (force matching, relative entropy, and correlation-based) in principle can estimate the unique potential of mean force for a specific CG mapping. The structure based methods with the use of the pair radial distribution function are comparable to the force matching and relative entropy when the later ones consider the family of proposed potentials, in Theorem 3, to consist of pair interaction potentials. The numerical comparison of all methods for molecular systems under equilibrium and non-equilibrium conditions is the subject of the future work.46
VII. CONCLUSIONS
The main goal of all systematic CG approaches, based on statistical mechanics, is in principle, to derive effective CG interactions as a numerical approximation of the many-body potential of the mean force, which for realistic molecular complex systems cannot be calculated exactly.
In this work, we have presented a general formalism for the development of CG methodologies for molecular systems. Below, we summarize the main outcomes of the detailed analysis presented in Secs. I–VI.
The probabilistic formalism discussed allows us to define a generalized force matching, as a CG minimization problem both for linear and nonlinear CG maps. This probabilistic formulation gives a geometric representation of the force matching method, as is schematically depicted in Figure 1.
A practical outcome of (a) is the connection of force matching with thermodynamic integration that provides a way on how to construct a local mean force in order to best approximate the potential of mean force with force matching. Specifically, this connection introduces a family of corresponding (to the CG map) coarsening transformations of the microscopic forces (local mean force) that ensure the best approximation of the PMF. This approach extends the work in Refs. 14 and 16, for any nonlinear CG map.
CG methods based on relative entropy and force matching are in principle asymptotically equivalent, in the sense of Theorem 3, both for the case of linear and nonlinear coarse-graining maps.
Furthermore, we prove, for linear CG maps in a specific example of a system with N molecules, that the (un-weighted) total force exerted at each CG particle satisfies the force matching condition when each particle is contributing to a single CG particle, see the example in V A 1 a. This fact, along with the example of the nonlinear CG map studied in Section V B, suggests that for complicated linear and nonlinear CG mappings, one can use appropriately formula (12) and achieve the best approximation of PFM with the force matching method.
Acknowledgments
The research of E.K. and V.H. was supported by the European Union (ESF) and Greek national funds through the Operational Program Education and Lifelong Learning of the NSRF-Research Funding Program: THALIS. The research of V.H. was also partially supported by ARISTEIA II. The research of M.K. was supported in part by the Office of Advanced Scientific Computing Research, U.S. Department of Energy under Contract No. DE-SC0010723. The research of P.P. was partially supported by the U.S. Department of Energy Office of Science, Office of Advanced Scientific Computing Research, Applied Mathematics program under Award No. DE-SC-0007046.
APPENDIX A: CONDITIONAL EXPECTATION AND COARSE GRAINING
Let (ℝ3N, 𝒢, μ) be the probability space induced by the random variable X of atomic configuration. is the σ-algebra generated by the random variable X, i.e., it is the collection 𝒢 = {A ∈ ℝ3N : ∃B Borel in ℝ3N s.t. X−1(B) = A}. Consider the coarse-grained random variable ξ = ξ(X) and define the sub σ-algebra of , induced by ξ,
i.e., any function ϕ : ℝ3N → ℝ, that is, -measurable is of the form
Denote Ω(z) = {x ∈ ℝ3N : ξ(x) = z}, the sub-manifold of ℝ3N corresponding to configurations x at a fixed value of the coarse grained variable z ∈ ℝ3M. The conditional expectation with respect to is the random variable , defined by
for any z and for any -measurable ϕ, with
that is, the average of ϕ keeping z fixed.
APPENDIX B: PROOFS
1. Proof of Theorem 2
Let the sub-manifold Ω(z) = {x ∈ ℝ3N : ξ(x) = z} of ℝ3N have the co-dimension 3M, i.e., dim(ℝ3N) − dim(Ω(z)) = 3M. The δ measure is defined as follows, for any smooth test function ϕ : ℝ3N → ℝ:
where ⋅t denotes the matrix transpose, det(⋅) the matrix determinant, and ΣΩ(z) denotes the surface measure on Ω(z). We define a mollifier on ℝ3N,
We have
recalling the notation (Dξ)ij = ∇xjξi, i = 1, …, M, j = 1, …, N, then we can write
and, since we assume that WDξt(x) is invertible, we can write
from which we have
Taking the limit as ϵ → 0, in view of Lemma 2 in Appendix B, we have
We recall the definition of potential of mean force (5),
which we rewrite as
Therefore, in view of relation (B1),
Thus, .□
APPENDIX C: EXAMPLE APPLICATIONS
1. Two CG particles coarse space
a. Proof of (16) and (17)
Applying Corollary 2, we have that
where , and the local mean force is given by
To prove (17), note that if we consider that each particle is contributing only to one CG particle, i.e., ζ1jζ2j = 0, the form of h(x) is simplified, since
and relation (17) is proved.
b. Proof of Sec. V A 1 a
Indeed, in this case, ζ1jζ2j = 0 for all j = 1, …, N, and following relation (14), we write
where δζij = 1 if ζij≠0 and δζij = 0 if ζij = 0 and calculate the det(I3N − TtB), where
Using the assumption that and properties of matrix determinants, we have that
Thus, we see that det(I3N − TtB) = 0 if 1 − δζ1j − δζ2j = 0 for all j = 1, …, N, i.e., when ζ1jζ2j = 0, or in other words when particle jcontributes only to one CG particle. Therefore, for this case, the typical choice of force matching methods to use the total force acting on each CG particle as the local mean force is, as expected, indeed a correct choice.
Assume now that there exists only one particle k that contributes to both CG particles. Let for simplicity choose k = 1, and ζ21 = 1, ζ2j = 0, j≠1, and ζ1j≠0 for all j = 1, …, N s.t. , then
which is nonzero. Thus, we found an example of CG map for which there does not exist any W s.t. This suggests that for this CG map, one should choose a W and then construct the h(x) in order to achieve the PMF approximation with the force matching, as is calculated, for example, in case (16).
2. Center of mass of N particles
In this example, the coarse grained variable is the center of mass of the N particles, see Figure 3. That is, M = 1 and the elements of the 3 × 3N coarse graining mapping matrix T are T1j = ζ1jI3, j = 1, …, N, as in (15), with
such that
Coarsening a many particle system to one CG particle, the center of mass of the N particles.
Coarsening a many particle system to one CG particle, the center of mass of the N particles.
We distinguish two cases: the first, choosing a specific W and looking for the form of h(x) and the second, by choosing a local mean force h(x) and looking for the existence of W such that the force matching indeed approximates the PMF.
a. Case 1
If we choose W(x) = T, then the local mean force h(x) is given by
This is a result of the application of Corollary 2, indeed we have that
where
Thus,
b. Case 2
We look for W such that the local mean force is the total force exerted at the center of mass . Note that hW(x) is nonzero if external forces are present. In view of relation (14), such a 3 × 3N W exists if I3N − TtB is singular, where B is the 3 × 3N matrix for which holds hW(x) = Bf(x). We have that
Since we assume that , the sum of all column elements is zero and det(I3N − TtB) = 0. Thus, there are infinitely many nontrivial solutions of W(I − TtB) = O3M×3N, W = [w11I3, …, w1NI3], w1j ∈ ℝ, j = 1, …, N. For example, one solution is given for w1j = 1, j = 1, …, N, that is, matrix W in Corollary 2 is for which , the typical local mean force for the force matching problem.
3. Two particles contributing to each CG particle
With this example, we examine the coarse graining where each CG particle position vector is the average of two particles position vectors, which contribute only to that CG particle, Figure 4. That is, assuming that the number of particles N is even, the number of CG particles is M = N/2 and the mapping is defined by
for ζij ∈ ℝ, and ζij = 0 if j≠2i − 1, 2i, such that ζi,2i−1 + ζi,2i = 1 for all i = 1, …, M. The 3M × 3N matrix of the linear mapping ξ is
where Tij = ζijI3, i = 1, …, M, j = 1, …, 2M.
Coarsening a many particle system with two particles per CG particle.
a. Example
Let W = T, applying Corollary 2, we have that
for i = 1, …, M. Let mj denote the mass of the jth particle and set
for i = 1, …, M, then if m2i−1 = m2i, i = 1, …, M, we can have
We show that there exists a family of 3M × 3N matrices W appearing in Theorem 2, such that
Let W with block entries wijI3, i = 1, …, M j = 1, …, N. The later equality holds if, in view of Corollary 2,
which we rewrite as ,where
Thus det(I3N − TtB) = 0 since ζi,2i−1 + ζi,2i = 1 for all i = 1, …, M. Therefore, according to (14), there exist infinitely many M × (2M) matrices W such that W(I3N − TtB) = 0 that gives
and .
4. Three atom molecule: End-to-end vector
A coarse variable that is of interest in molecular systems is the end-to-end vector ℓ13 = x3 − x1 ∈ ℝ3, the corresponding map is linear with
The mapping has the 3 × 9 corresponding matrix
Since the mapping is linear, we can apply Corollary 2, which states that the average of the microscopic forces is given by