Relative Resolution: A Multipole Approximation at Appropriate Distances

Recently, we introduced Relative Resolution as a hybrid formalism for fluid mixtures [1]. The essence of this approach is that it switches molecular resolution in terms or relative separation: While nearest neighbors are characterized by a detailed fine-grained model, other neighbors are characterized by a simplified coarse-grained model. Once the two models are analytically connected with each other via energy conservation, Relative Resolution can capture the structural and thermal behavior of (nonpolar) multi-component and multi-phase systems across state space. The current work is a natural continuation of our original communication [1]. Most importantly, we present the comprehensive mathematics of Relative Resolution, basically casting it as a multipole approximation at appropriate distances; the current set of equations importantly applies for all systems (e.g, polar molecules). Besides, we continue examining the capability of our multiscale approach in molecular simulations, importantly showing that we can successfully retrieve not just the statics but also the dynamics of liquid systems. We finally conclude by discussing how Relative Resolution is the inherent variant of the famous"cell-multipole"approach for molecular simulations.

1 Introduction most significant challenge with all of these multiscale strategies is the fact that a very undesirable computational step is always required: This step involves performing a molecular simulation of a particular FG model before even using a respective CG model; realize that in general, if we can construct a molecular simulation of the former, we do not really need the latter. Finally, we note that because of the reduction in the degrees of freedom, time-scales naturally change during multiscale optimization, and thus, kinetic features will inherently deviate from their true behavior.
The in-parallel multiscale methods generally contain FG and CG information simultaneously in a unified system: A single molecular simulation contains both FG and CG models, with vital aspects described in detail by the latter, and trivial aspects described with simplicity by the former. Diverse strategies have been developed that follow this hybrid computational path, and we categorize them here in four main classes. The most straightforward approach is the "group-based" multiscale method: Here, essential molecular groups (e.g., solute molecules) are described via FG interactions, and auxiliary molecular groups (e.g., solvent molecules) are described via CG interactions [24,25,26,27]. The next class of strategies can be thought of as the "time-based" multiscale method: A given molecular simulation spends some of its time as the FG system and some of its time as the CG system, with the results being recorded for the former and discarded for the latter [28,29]. A powerful variation of this approach actually exchanges among FG and CG replicas of the molecular simulation, which enhances the overall efficiency of the computational procedure [30,31].
In either case, this hybrid approach is obviously convenient for overcoming various time-scale barriers in a system of interest. Perhaps the most successful class of the in-parallel multiscale methods has been the one which is commonly called Adaptive Resolution. With much research done on it for over a decade, this approach is a "one-body distance-based" multiscale method: In a single molecular simulation, adaptive molecules switch their resolution in terms of absolute position [32,33,34,35]. For example, as molecules move around a system, they adopt a FG model if near to the origin and a CG model if far from the origin. In practice, Adaptive Resolution essentially constructs hybrid molecules that embody both FG and CG models, and in fact, the main difference between all variations of Adaptive Resolution has been the mixing rule for the hybrid interaction as a function of absolute position. In the original publication for Adaptive Resolution, Praprotnik et al. introduced the hybrid interaction as a linear combination of forces, and this has been the ideal choice for evolving Newtonian trajectories [32]; in that same journal issue, Abrams presented an alternative route that switches between the models in a stochastic manner [33]. Conversely, Ensing et al. formulated the hybrid interaction as a linear combination of potentials [34]; in turn, Potestio et al. showed that such a Hamiltonian version of Adaptive Resolution naturally corresponds with many aspects in statistical mechanics [35]. On a practical level, Adaptive Resolution is especially useful if a given system has a specific region of special interest (e.g., a protein at the origin, with everything else being water). As such, Adaptive Resolution has been already applied in some biological scenarios, having quite successful results [36,37].
Very recently, we developed another type of multiscale methods, termed Relative Resolution (RelRes) [1]. Inspired by Adaptive Resolution, this approach is a "two-body distance-based" multiscale method: In a single molecular simulation, hybrid molecules switch their resolution in terms of relative separation; in particular, molecules that are near neighbors (i.e., their pairwise distance is small) interact via their FG models, and molecules that are far neighbors (i.e., their pairwise distance is large) interact via their CG models [1]. Importantly, RelRes is the sole class of multiscale simulations which describes all molecules, for all times, with both FG and CG models; the resolution of a given molecule is always relative of its observer. While other strategies, reminiscent of RelRes, have been also developed [38,39], our formalism is unique in that it mathematically finds a natural connection between the FG and CG potentials. While Ref. [38] numerically parametrizes between the two models via "force-matching", Ref. [39] does not make a clear connection between the two models. Conversely, we developed an analytical parameterization between the FG and CG models, which is just based on energy conservation [1]. Unlike those other strategies [38,39], we were consequently able of correctly retrieving across state space, the structural and thermal behavior of several nonpolar mixtures [1]. Besides, we showed that our hybrid approach can be considered as a generalized extension of established theories for uniform liquids, which assume a mean field for interactions beyond a certain distance [40,41,42,43].
In this current work, we build on our original communication of RelRes [1]. Foremost, we mathematically present the comprehensive framework for RelRes. Above all by performing a multipole expansion for an arbitrary potential of inverse powers, we show that its zero-order term, which is sufficient for nonpolar systems, is identical with the succinct expression we presented in our initial publication; in consideration of polar systems, the current work also presents the first-order and second-order terms of the Taylor series, both of which involve molecular orientation. By use of this multipole approximation, we introduce RelRes as a Hamiltonian that switches between the FG and CG models at an appropriate distance. Building on the computational results of our original work, we consequently continue examining RelRes with molecular simulations of nonpolar systems. Interestingly, we find that the ideal switching distance between the FG and CG models is roughly equivalent with the location at which a particular orientational function essentially decorrelates. Besides, we also demonstrate that RelRes can correctly describe not just static behavior but also dynamic behavior. Together with our previous publication [1], we suggest that RelRes can be a powerful multiscale tool for efficiently studying soft matter via molecular simulations.

Theoretical Foundation
Contrary to our original publication [1], we present here our hybrid formalism in a reversed order. In the previous work, we started by defining RelRes, and we ended by parameterizing between the FG and CG potentials. In the current work, we start by parameterizing between the FG and CG potentials, and we end by defining RelRes. Regardless, the mathematical notation is identical between our two publications. Note also that the current framework is the comprehensive one, applying for all molecular systems. The original publication is the terse version, which is only applicable for nonpolar scenarios.

Defining our System
Consider that we are performing a molecular simulation of a pair of molecules in vacuum; a particular configuration of this trivial system is depicted in Fig. 1. We label the gravitational centers by the Latin indices i and j, and we label the atomistic coordinates by the Greek indices µ and ν. In our notation throughout, we denote n i or n j as the number of sites on a molecule i or j, respectively. The relative separation between the centers is − → r ij , and the relative separation between the coordinates is − → r µiνj . Furthermore, the distance between a certain atomistic coordinate and a respective gravitational center is − → ∆ µi or − → ∆ νj . By Fig. 1, it is clear that these distances are related as follows, and we defined in this expression the following variable: Note that we introduce here all distances in vector form, and their scalar definitions naturally follow. We also define dimensionless variables that will compact all of our ensuing mathematics: With these definitions, together with some rearrangement, we attain this useful relation which basically stems in the dot product of Eq. 1 with itself: Finally for clarity, we must make a subtle distinction with a variable, ∆r µiνj , that appears in our original communication, but that we do not use here: ∆r µiνj = r µiνj − r ij = ∆ µiνj .
We now focus on the energetics of this rudimentary system. The governing function here is defined as follows: Here, u µiνj is the intrinsic potential between atomistic coordinates µ i and ν j , and u ij is the resultant potential between gravitational centers i and j; the latter is obtained by the respective summation of the former. Notice that we are allowing for absolute nonuniformity in our mixtures considering the indices on the potentials, u µiνj and u ij . Importantly, we assume that the former is isotropic, being exclusively a function of the scalar r µiνj . Conversely, we expect that the latter is geometric, being chiefly a function of the vector − → r ij , as well as the two sets of intramolecular distances, − → ∆ µi and − → ∆ νj . Importantly, we ignore intramolecular energetics during most of our derivation: These are accounted for once we present the entire Hamiltonian of RelRes.
One of the focal assumptions in our analysis is that the basis function of the potential can be cast as an inverse power of the isotropic distance: Obviously, m is the power associated with this inverse law. In most scenarios, like with the Coulomb potential (m = 1), as well as with the Lennard-Jones (LJ) potential (m = 12 and m = 6), the inverse law is a natural choice. The proportionality coefficient c µiνj corresponds with the unique parameters involved in the interaction between sites µ i and ν j , and it is usually known empirically for use in molecular simulations. We will later make an important assumption of geometric mixing for c µiνj (i.e., we will express it as a product of two separate parameters).
In the remainder of the derivation, our goal is to reformulate u ij of Eq. 6 in terms of r ij , together with variations of − → ξ µiνj and cos θ µiνj , as defined by Eqs. 3 and 4. Such a reformulation holds a promise for the efficient calculation of the energy of the pair in Fig. 1; obviously, computation of a single distance, r ij , is preferable over the computation of n i n j distances, r µiνj , which the summation of Eq. 6 requires. Of course, such a reformulation may involve approximations, and thus, we must make valid assumptions for − → ξ µiνj and cos θ µiνj so that we maintain the correct energy for this pair.

Introducing the Multipole Expansion
We substitute Eq. 7 in Eq. 6, which yields the following: For clarity throughout most of the ensuing analysis, − → ξ µiνj and cos θ µiνj are omitted yet implied in the functionalities involved in Eq. 8; the same also goes for the bar on − → r ij . Consider now the radial distance which appears here together with the power of m; we correspondingly exponentiate Eq. 5, and we present the ensuing reciprocal: We can use this expression to substitute the computationally inexpensive r ij for the computationally expensive r µiνj in Eq. 8. We consequently attain the following expression: Analogous with the usual multipole expansion, we now reformulate, via a Taylor series, the cumbersome expression which appears above, with Λ m,ℵ basically defined by the appropriate derivative in terms of the index ℵ: Clearly, Λ m,ℵ is only a function of cos θ µiνj , with the functionality of ξ µiνj appearing as a power in Eq. 11. Realize that such a Taylor expansion is justified by the fact that usually, ξ µiνj < 1 (i.e., ∆ µiνj < r ij ), since in most cases, bonds are negligible in magnitude compared with intermolecular distances. Notice also that the exponent m of the inverse law is just a parameter in Λ m,ℵ . Besides, if m = 1 (i.e., the Coulomb potential), Eq. 11 simply becomes the common definition of the Legendre polynomials. Substituting Eq. 11 in Eq. 10 yields the following power series for our energy function: and it is particularly convenient to cast the above as follows, in which we obviously defined the following: One may think of c ij,ℵ as an effective coefficient for a given configuration of Fig. 1; analogous with u ij , ξ µiνj and cos θ µiνj are consequently omitted yet implied in the functionality of c ij,ℵ throughout most of our derivation. In any case, notice that these expressions are especially useful for computational purposes: By summing over all sites µ i and ν j , one can separately calculate the energetic coefficient of any ℵ term by Eq.
15, and by including as many ℵ terms in Eq. 14 as desired, one can have a valid approximation for the overall interaction between the molecular pair of Fig. 1. Of particular importance is the leading (nonvanishing) ℵ term: For making a connection with the terminology of the original work, we term it as the Molecular Pair in the Infinite Limit (MPIL).

Evaluating the Multipole Terms
By not even performing any differentiation, the ℵ = 0 term is readily determined by Eq. 12 as the following: Conveniently, this is always unity, regardless of the value of m. Plugging this in Eq. 15, we attain this expression: Notice here that this is not dependent at all on the configuration of the molecular pair: This means c ij,0 is a constant that can be calculated before any molecular simulation is performed. If this is nonzero, we can approximate our energy in Eq. 14: In the later part of this expression, we obviously invoked Eq. 17. Notice the similarity between Eqs. 8 and 18: We start with an m inverse law in the former, and we end with an m inverse law in the latter. In Eq. 8, we sum over the many distances between the atomistic coordinates, with their intrinsic interaction coefficients, and in Eq. 18, we take in the single distance between the gravitational centers, with its resultant interaction coefficient. This is in effect the monopole term; once we introduce another assumption, we will show that this is specifically equivalent with the monopole-monopole interaction.
This is in essence the sole term that we derived in our previous work. We simply attained it by vanishing all bond lengths (i.e., ∆ µiνj = 0) in Fig. 1; this is basically equivalent to setting all ξ µiνj = 0 in Eq. 10, and in turn, it also eliminates the corresponding functionality of all cos θ µiνj in our energy expression. The only discrepancy of Eq. 18 with the corresponding one in our original work is that in this one, we have the assumption of the inverse power for the interaction. We can remove this assumption temporarily by invoking Eq. 7 for Eq. 18, and we in turn get the exact expression which we termed as the MPIL in our original publication: Realize that that this is a special case of our current definition of the MPIL: As mentioned before, the general definition of our MPIL is the leading (nonvanishing) ℵ term in Eq. 13. Importantly in our original work, we showed that if Eq. 18 is implemented with RelRes, it can successfully describe nonpolar mixtures that are based on the LJ potential. Theoretically, this approximation shall also work with the Coulomb potential for ionic molecules that have a net charge. In fact for RelRes, we expect that Eq. 18 shall be adequate for any m, given that this monopole term does not vanish. Nevertheless, once the above summation for a certain molecular pair contains both positive and negative values for c µiνj , there may be some issues; specifically, if all c µiνj in the summation cancel each other, Eq. 18 is strictly zero. The most obvious such case is the zwitterionic scenario, in which the molecules have Coulombic charges, yet they are overall neutral (e.g., water). It is obvious that the interaction between such polar pairs is finite, and consequently, one must go beyond the monopole term of Eq. 18 in order to describe the corresponding energy, at least approximately. This is in fact the purpose of the ensuing mathematics. Note that the equations below may be also useful as correction terms for the interaction between the molecular pair, even if their monopole term is nonzero.
Obviously, all the terms of the Taylor series can be evaluated by successive differentiation as presented in Eq. 12. Here is the first-order term: Consequently, the ℵ = 1 term of Eq. 15 is the following, Here is the second-order term: Consequently, the ℵ = 2 term of Eq. 15 is the following, Here, we introduced this definition: We also emphasize here that if m = 1 (e.g., the Coulombic scenario), Λ m,1 and Λ m,2 become the familiar first and second Legendre polynomials, respectively. Of course, these c ij,ℵ can be readily used for evaluating the energy function of Eq. 14: They may serve as correction terms for Eq. 18, and if c ij,1 or c ij,2 is actually the leading term of this multipole expansion, it can even become the sole approximation for the energy (i.e., the MPIL). Notice that unlike c ij,0 , these current expressions are functions of cos θ µiνj , and well as ξ µiνj . Thus, while Eq. 18 is isotropic, involving just the scalar r ij in its functionality, the energy which corresponds with other c ij,ℵ is geometric, involving also the vector − → ∆ µiνj in its functionality. Thus for ℵ = 0, the coefficient must be computed at each step of a molecular simulation. As such, going beyond the monopole term obviously requires significant computational cost, and it is only recommended if it is necessary (i.e., for polar systems in which these are the leading terms in the Taylor series).

Assuming a Proportionality Coefficient
Let us now discuss the proportionality coefficient c µiνj of the inverse law of Eq. 7. In most cases, such a parameter is empirically known for use in molecular simulations. In fact, it is usually expressed in terms of some mixing rules between two separate coefficients c µi and c νj of sites µ i and ν j , respectively. For example with the Coulomb potential, c µiνj is fundamentally the product of partial charges; an analogous product is also sometimes used for the LJ potential. For the remainder of our work, we assume that the following equation holds: For completeness, this means that Eq. 8 (i.e., the initial energy function in our derivation that introduced the inverse power), becomes the following: In any case, Eq. 26 substantially facilitates our computational algorithm for RelRes, especially since a variant of c ij,ℵ , for any ℵ, becomes a constant during a molecular simulation.
Let us begin with the zero-order term. Substituting Eq. 26 in Eq. 17, we attain the following equation in which we have introduced the definition of the monopole c i,0 or c j,0 , for each molecule i or j, respectively: If both of these monopoles are nonzero, the MPIL is the following according to Eq. 14: Analogous with our comparison of Eqs. 8 and 18 above, notice the similarity between Eqs. 27 and 30: Performing the approximation, we start with several inverse coefficients c µi c νj , and we end with one inverse We now consider the benefit of the product assumption of Eq. 26. You may notice that Eq. 30 may be slightly more computationally efficient than Eq. 18. The former involves two summations over n i and n j parameters, while the latter involves one summation over n i n j parameters. Regardless, all of these summations can be performed before a molecular simulation, and thus, they have negligible computational cost. The computational superiority introduced via Eq. 26 becomes very apparent once we move to the first-order and second-order terms of the multipole expansion.
We may now apply the assumption of Eq. 26 for evaluating any c ij,ℵ , as presented in Eq. 15, after the appropriate differentiation of Eq. 12. As an important part of the ensuing mathematics, we introduce the familiar definition for the dipole − → p i,1 or − → p j,1 , of each molecule i or j, respectively: Moreover, we introduce components of the familiar definition for the quadrupole − → − → q i,2 , q i,0 or − → − → q j,2 , q j,0 , of each molecule i or j, respectively: By casting a linear combination of these, we retrieve the conventional definition of the quadruple for each molecule i or j, respectively; this is most apparent for m = 1 (e.g, the Coulombic scenario). As in the usual case, we are permitting here for generalized tensor algebra (e.g., − → − → I is the identity tensor); the amount of bars on these parameters, together with their last index for extra emphasis, corresponds with the order of the tensor. For compactness, we introduce here a compact notation for the set of dipoles (i.e., , as well as for the set of quadrupoles (i.e., As we all know, the elegance of each of these parameters is that it is only a function of the topology of each molecule, and it does not at all depend on the overall configuration of the molecular pair in Fig. 1.
By invoking these, we specifically show in the Appendix for ℵ = 1 and ℵ = 2 how the following expressions for c ij,ℵ can be derived. Here is the first-order term: Here is the second-order term: Of course, keep in mind the linear combination of Eq. 25, which we cast as follows here: Besides, note that the "hats" above r denote a distance tensor of a unity magnitude; their amount is equivalent with the order of the tensor. While usually omitting yet implying the functionality of c ij,ℵ in the remainder of our mathematics, we emphasize here that the coefficients above are actually dependent on the set of dipoles − → p ij,1 , together with the set of quadrupoles − → − → q ij,2 and q ij,0 ; of course, these are also functions of r ij , as well as of the monopoles. Besides, note that for m = 1, the correspondence of these expressions with the familiar Coulombic case becomes apparent, considering our definitions of the quadruples in Eqs. 32 and 33.
Importantly in our general formalism, these coefficients can be employed in Eq. 14 for evaluating the interaction between the molecular pair of Fig. 1: These can be used as correction terms for Eq. 30, and for cases in which the zero-order term vanishes, c ij,1 or c ij,2 may even serve as the MPIL. Let us now consider the computational implementation of Eqs. 34 and 37, especially in relation with our corresponding discussion of Eqs. 21 and 25, respectively (i.e., the discussion before we assumed the product rule of Eq. 26). Above all, Eqs. 34 and 37 are dependent on parameters of gravitational centers i and j (i.e., − → p ij,1 , as well as − → − → q ij,2 and q ij,0 ), while Eqs. 21 and 25 are dependent on parameters of atomistic coordinates µ i and ν j (i.e., ξ µiνj , as well as cos θ µiνj ). Of course, the latter requires a summation over n i n j sites, while the former does not; besides, − → p ij,1 , together with − → − → q ij,2 and q ij,0 , can be evaluated analytically at no computational cost irrespective of − → r ij , while ξ µiνj together cos θ µiνj must be calculated at each step of a molecular simulation with consideration of − → r ij . Still, just as with Eqs. 21 and 25, the energy corresponding with Eqs. 34 and 37 involves vectors, which is unlike the very convenient energy function of the monopole-monopole term of Eq. 30, which solely depends on scalars. In turn on a practical level, we do not recommend going beyond the zero-order term in the multipole expansion unless it is necessary (e.g., water).

Defining Relative Resolution
Up until now, we have been essentially looking only at the infinite limit associated with Fig. 1. How do we deal with an arbitrary distance between the molecular pair? RelRes can resolve this issue, yet we must begin by clearly defining the FG and CG potentials, which mostly amounts to introducing the appropriate labels on the potentials we have been working with above. Considering Fig. 1, the FG potential u F G µiνj , only a function of r µiνj , is the fundamental interaction between atomistic coordinates µ i and ν j , and the CG potential u CG ij , notably a function of r ij , is the apparent interaction between gravitational centers i and j; the latter is an approximation of the appropriate summation of the former: Compare this expression with Eq. 6; importantly, the right sides but not the left sides are identical between the two expressions (i.e., u F G µiνj = u µiνj but u CG ij = u ij ), which is in fact the reason for the approximation in Eq. 38. With an inverse assumption analogous with Eq. 7, we can clearly define the approximation here via the multipole expansion we have performed above (i.e., from Eq. 8 to Eq. 13). The following then ensues: Remember of course the relation between the coefficients which appear here as given by Eq. 15, and ℵ * is the leading (nonvanishing) term in the Taylor series of Eq. 12; this is in fact equivalent with the MPIL.
Besides, if we invoke the product assumption of Eq. 26, we can cast these potentials in terms of the ensuing multipoles (e.g., Eq. 28).
Let us now cast the interaction between the molecular pair of Fig. 1 as a single function by using Eqs.
39 and 40. For this purpose, we must again remember the MPIL that we thoroughly discussed above. If the molecules of Fig. 1 are near to each other, the approximation of MPIL is unreasonable, and the FG potential of Eq. 39 is the relevant one. If the molecules of Fig. 1 are far from each other, the approximation of MPIL is legitimate, and the CG potential of Eq. 40 is the relevant one. We illustrate these ideas for nonpolar molecules in the top panel of Fig. 2. This is in fact the main idea of RelRes, and we define its potential as follows: In essence, RelRes is a linear combination ofũ F G µiνj andũ CG ij , which are slight modifications of the FG and CG potentials of Eqs. 39 and 40, respectively; we thoroughly discuss these functions below. Besides, note that we emphasize in Eq. 41 that the RelRes potential is a function of various pairwise distances (i.e., atomistic coordinates, r µi νj , and gravitational centers, r ij ); in such a way, RelRes maintains all degrees of freedom. Note also that the dependence on ξ µiνj and cos θ µiνj has been omitted here for clarity, especially since these variables do not have any effect on nonpolar molecules.
The modifications of the FG and CG potentials must ensure that for small relative separations,ũ F G µiνj = u F G µiνj while the derivative ofũ CG ij is zero, and for large relative separations,ũ CG µiνj = u CG µiνj while the derivative ofũ F G ij is zero. Their exact functionalities are not unique, yet we choose the following piecewise functions for our purposes:ũ The switching distance r s is a constant that represents the relative separation at which the MPIL supposedly switches from being unreasonable to being legitimate. Note that the vertical shifts u F G µiνj (r s ) and u CG ij (r s ) merely ensure that the corresponding forces are continuous throughout the domains; through the linear combination of RelRes, these vertical shifts approximately cancel each other out, with the extent of the cancellation depending on how many nonzero terms the summation of Eq. 41 includes.
Eq. 41 is actually the pairwise version of RelRes, which only applies for a molecular pair in vacuum. Extending RelRes for a multiscale simulation with many molecules, we define the complete version of RelRes through a summation over all different pairs of molecules in the system: We omit the functionality of this complete version of RelRes throughout our work for clarity. This Hamiltonian introduces a notable conceptual complexity, as compared with Eq. 41 of the pairwise case. For a molecular pair in vacuum, the entire system is basically either a pure FG or CG scenario: If the molecules are near to each other, they are both described by the FG model, and if the molecules are far from each other, they are both described by the CG model (of course, there are subtleties for moderate relative separations). Nevertheless, in the case of a molecular simulation of a realistic liquid, the situation is completely different, with the system being neither purely FG nor CG, but instead having a hybrid nature: Each molecule simultaneously embodies both models, and a given molecule interacts with its near neighbors via the FG potential and with its far neighbors via the CG potential. The bottom panel of Fig. 2 illustrates a RelRes system for a nonpolar scenario, which emphasizes that all molecules embody both molecular resolutions.
Let us also discuss RelRes in terms of its r s parameter, labeling such a dependence asŨ rs . The two limits of r s are of particular importance for emphasis:Ũ ∞ is identical with the corresponding Hamiltonian of the pure FG system, andŨ 0 is identical with the corresponding Hamiltonian of the pure CG system. The former is the basis for reference mixtures. Such a system has r s → ∞ (i.e., a system with just FG and no CG interactions), and it is depicted in the bottom panel of Fig. 2. An inherent goal for RelRes with an arbitrary r s is of approximately capturing all behavior of its corresponding reference mixture (e.g.,Ũ rs ≈Ũ ∞ ). For completeness, let us also define the entire system Hamiltonian: K is the kinetic energy, while U 0 accounts for all intramolecular interactions. Note that these do not have a tilde above them. This is because RelRes maintains all degrees of freedom, and in turn, its K, as well as its U 0 , is strictly unaltered with the corresponding functions of the reference system (e.g., K rs = K ∞ ).

Computational Validation
We now turn our attention for testing the efficacy of our multiscale approach in describing reference systems . We do this via molecular simulations, complementing in turn the initial examination of our previous publication [1]. Our original work was a proof of concept that RelRes is very successful in describing (nonpolar) multi-component and multi-phase systems across state space; by invoking a tuning parameter in our Hamiltonian, we also showed that MPIL is in fact the ideal choice for the potential beyond nearest neighbors. Consequently for simplicity in this current work, we only examine uniform liquids, taking MPIL for granted in all of our systems. Here, we instead focus on systematically varying the distance at which we switch between the FG and CG potentials. Besides, we also complement our initial work by showing that RelRes with MPIL works well not only for static behavior but also for dynamic behavior.

Constructing the Molecular Simulations
Let us begin by thoroughly discussing the general implementation of our molecular simulations. It is almost identical with the implementation we had in our original publication [1]; we emphasize important aspects which are different. Above all, we again perform our molecular simulations with the GROMACS package, specifically its 4.6 version [44]. If we do not specify here certain features of the molecular simulations, that basically means that the default values are used.
In the current publication, all of our systems contain a total of 2000 (identical) molecules: All of them are dumbbell-like molecules (e.g., ethenes), with two FG sites and one CG site, meaning that n = 2. All FG atomistic coordinates have mass m, and all CG gravitational centers have mass zero; the latter is constructed via equal weights of the former. In all cases, the LJ potential is the sole component for the interactions between molecules; this in turn sets the length and energy parameters, σ and ǫ, respectively, for our work. For a practical implementation of RelRes so that there are no singularities in Newtonian trajectories, we modify the step functionality of Eqs. 42 and 43: In each of these in the interval [r s − δr s , r s + δr s ], we introduce the following four-term polynomial in terms of the inverse distance: The four coefficients γ are determined by ensuring that the potential, as well as its derivative, is continuous at its two boundaries, r s ±δr s . We vary r s between 1.1σ and 1.9σ at intervals of 0.2σ, and we set δr s = 0.0625σ.
On an analogous note, rather than sharply truncating all potentials at 2.5σ with a step function, we invoke the sigmoid function above between 2.25σ and 2.75σ. Note that all interactions in our reference system, denoted by r s → ∞, practically vanish beyond this distance.
The intramolecular energetics are mostly governed by elastic bonds, which is unlike our previous publication that had all bonds being rigid. Here, the distance of each bond is 0.5σ, while its spring is 2000ǫ/σ 2 .
Comparing with a conventional bond, our distance is shorter by a significant fraction, while our spring is about an order of magnitude weaker. This is purposefully done as a stringent test: If RelRes can work on our large flexible molecules, it can also work on small stiff molecules.
In all cases, the protocol that we use involves a sequence of molecular simulations at a temperature of 1.0ǫ/k, with k being Boltzmann's constant. The initial molecular simulation is of the reference liquid, and it is distinct in that it is coupled to a barostat, whose pressure is 1.0ǫ/σ 3 ; the purpose of this molecular simulation is to fix the system size for the entire sequence, while all of its results are disregarded. The rest of the molecular simulations are subsequently in the canonical ensemble: One of them is of the reference liquid, while the rest are RelRes systems of different switching distances. Besides, all molecular simulations evolve via Newtonian equations of motion. Each molecular simulation starts with an equilibration phase of 5, 000 steps of size 0.002τ and ends with a production phase of 1, 000, 000 steps of size 0.001τ ; note that τ = mσ 2 /ǫ is our unit of time. Realize that all molecular simulations that we perform in this work are of single-component and single-phase systems.

Computing Static and Dynamic Features
Above all, we thoroughly examine in our work several structural correlations, usually in terms of relative separation r ij ; we frequently omit the ij indices throughout most of our work. For this purpose, we record the positions of all molecules every 100 steps of the molecular simulation. The domain of our r goes completely from the origin to the edge of the system box; we discretize it by 1000 bins.
Foremost, we look at radial distributions, g, between the midpoints of bonds. For comparison between radial distributions, it is convenient to invoke a functional for them. We particularly choose the Jensen-Shannon entropy: with the conventional entropy taking on this definition: κ is just the normalization constant (i.e.,´κgr 2 = 1), which is inversely proportional with the system volume. In Eq. 47, the term on the left (i.e., S 1 2 (g ∞ + g rs ) ) is one entropy, altogether for an average of two distributions, and the term on the right (i.e., 1 2 (S [g ∞ ] + S [g rs ])) is an average of two entropies, each for one distribution. As such, the Jensen-Shannon entropy is a functional of two radial distributions, which just measures the disparity between them: Obviously if g rs ≈ g ∞ , S JS ≈ 0, and S JS increases as the discrepancy between g rs and g ∞ intensifies.
Note that while in the previous work, we presented g between atomistic coordinates, in the current work, we basically focus on g between gravitational centers. In turn, it is important to complement these radial distributions with orientational correlations. By letting a bond be a vector, we define s i and s j as the respective directions of molecules i and j. We particularly focus on computing the moments of their corresponding dot product, s i · s j , as a function of the relative separation between the midpoints of the bonds. Most of our case studies have symmetry in their molecules, and thus for s i · s j , the average vanishes but the variance persists. We consequently present in our work only the latter, ( s i · s j ) 2 ; note that ( s i · s j ) 2 = 1 3 for decorrelated cases. Moreover, the bond vectors can be cast as a linear decomposition in terms of two components s i = s i + s ⊥ i and s j = s j + s ⊥ j : In terms of the relative separation between molecules i and j, s i , s j are the components parallel to r ij , and s ⊥ i , s ⊥ j are the components perpendicular to r ij . In turn, we have the following: The approximation is introduced here because in most of our case studies, we observe that the cross term in the above expression is rather negligible, and thus, we only present the two main components, s i s j 2 and s ⊥ i s ⊥ j 2 , of this orientational function; realize that s i s j 2 = 1 9 and s ⊥ i s ⊥ j 2 = 2 9 for an ideal gas.
Besides these static functions, we also examine dynamic functions, something that we did not at all do in our original publication. Specifically, we look at the squared displacement ( r t − r 0 ) 2 , as well as its derivative, ∂ t ( r t − r 0 ) 2 . r t is the position of molecule i or j at a given time t, and the index 0 is of course for zero time. Besides, we also examine an orientational correlation s t · s 0 . s t is analogous with the dumbbell directions introduced earlier; just as with the squared displacement, we have here the time index, while the molecule indices, i or j, are omitted for clarity. Importantly in calculating these transport functions, we employ an algorithm of multiple time origins, in fact using each recorded step for such purposes.
We furthermore examine thermal properties in our work. We foremost look at the average U and ; as we frequently do, we omit the index of r s for these thermal properties. Besides, we also look at a transport function of the total energy, δ E t δ E 0 / δ E 2 , computing it again using multiple time origins. All of our thermal properties are based on probing our molecular simulations every 10 steps.

Examining Elementary Dumbbells
We begin by examining a system of elementary dumbbells. These can be considered as ethenes, and one such molecule is depicted in Fig. 3. Such a dumbbell has two identical FG sites, with σ and ǫ for its LJ parameters. By the MPIL, this means that the respective CG site has σ and 4ǫ for its LJ parameters since n = 2. The inherent bond of these elementary dumbbells is set at at a distance of 0.5σ. The density for this system is 1.00σ 3 /m.
The radial distribution between the dumbbells is given in the top panel of Fig. 3. This structural correlation for the reference liquid is given here as the solid black curve; the remaining dashed curves are for the RelRes systems, with each color representing a different switching distance. All of these switching distances give a sufficient description of the radial distribution of the reference system, and it is clear that as r s increases, the capability of RelRes improves, in an apparently asymptotic manner: While for r s = 1.3σ (i.e., the blue curve), the replication is fairly decent, for r s = 1.7σ (i.e., the red curve), the replication is essentially perfect; in fact, r s = 1.5σ (i.e., the violet curve) already captures the radial distribution of the reference system as well as one may desire. The systematic observation here is reminiscent with the one we had in the original publication [1]: RelRes captured structural correlations with r s ≈ 1.1σ equitably yet with r s ≈ 1.6σ splendidly.
For the purpose of better clarifying this asymptotic behavior, we invoke the Jensen-Shannon entropy, as defined by Eq. 47. For each of the relevant curves in Fig. 3, we calculate S JS , and we plot it (i.e., an indigo circle) in terms of the switching distance in Fig. 4. Clearly, this entropy is monotonic throughout most of the domain of r s , and considering the logarithmic scale here, S JS is roughly an exponential decay of the distance; this is especially true around r s = 1.5σ, which we mentioned in the context of Fig. 3 as the critical switching distance for adequately capturing the radial distribution of the reference system. Thus, we confirm here our earlier observation: RelRes improves asymptotically as r s approaches infinity. Importantly, it appears that the asymptote is practically reached at r s ≈ 1.5σ: At this distance, the discrepancy between the radial distributions almost vanishes, and in turn, S JS becomes rather negligible.
So why is 1.5σ the critical distance for switching between the FG and CG potentials in RelRes? For this purpose, we now look at orientational correlations in the bottom panel of Fig. 3. Let us begin by discussing the black solid curve, which is ( s i · s j ) 2 of the reference liquid. Notice that while it clearly has a maximum and a minimum early on, it quickly flattens out after the ensuing inflection at 1.42σ. This consequently explains why 1.5σ is an excellent choice for r s in RelRes: At this distance, the dumbbell directions apparently become decorrelated, and the influence of these orientational degrees of freedom becomes negligible on the various structural correlations of the system (e.g., radial distributions). Interestingly, the flattening of ( s i · s j ) 2 in the bottom panel happens just around the middle between the respective maximum and minimum in g of the corresponding top panel (i.e., the relative separation of 1.54σ); this is also around the respective inflection of g (i.e., the relative separation of 1.46σ), which we mark by a vertical brown line in both panels. This means that as soon as nearest neighbors depart from each other, their directions quickly become decorrelated.
Let us continue examining the two main components of ( s i · s j ) 2 , as given in Eq. 49. For the reference system, these are plotted as gray curves in the bottom panel of Fig. 3, the lower one is for s i s j 2 and the higher one is for s ⊥ i s ⊥ j 2 ; importantly, realize that their summation essentially yields the black curve.
Surprisingly, the behavior of ( s i · s j ) 2 is very distinct from its components: While ( s i · s j ) 2 becomes mostly decorrelated as the molecules leave the primary coordination shell, its components are still correlated even as the molecules enter the secondary coordination shell. In fact, s i s j Importantly, RelRes captures the orientational correlation with its components perfectly well, regardless of which switching distances we use. This is in contrast with the radial distribution, which is noticeably influenced by r s . We now move on to examine thermal properties of these elementary dumbbells. We specifically look at functionals of the configurational Hamiltonian of RelRes, presenting, as a function r s , its normalized average U and variance δ U 2 in the bottom and top panels of Fig. 5, respectively. The coloring here is analogous with Fig. 4. Focus again on the indigo circles: Because of normalization, they fluctuate about the horizontal brown line which is of course set at unity. In fact, for all switching distances, the average and variance are essentially both within a fraction of 0.05 from their reference values. Importantly, we notice that a dampening effect of the fluctuation, with both the average and variance eventually reaching unity.
Besides, the fluctuating characteristic of the symbols in this figure is in striking contrast to the observation we made for the functional of the structural correlations in Fig. 4, which showed an asymptotic behavior with r s . Bearing also in mind Fig. 3, this means that if one is fine with a decent yet rough description of radial distributions in a liquid, one can still adequately capture thermal properties, together with orientational correlations, while using just a modest switching distance in RelRes that does not entirely account for all nearest neighbors. Finally, these systematic findings are analogous with the observations we made in our previous work for the pressure, together with its corresponding response function. For both of our switching distances there, r s ≈ 1.1σ and r s ≈ 1.6σ, RelRes was very successful once we used it together with MPIL; [1].
In general, the excellent replication of thermal properties reiterates the validity of the MPIL approximation, which is essentially based on energy conservation.
So far, considering our original publication as well, we have only examined static behavior; we now turn our attention to dynamic behavior. We consequently plot various transport functions in terms of time, t, in Fig. 6. The coloring here is again analogous with Fig. 3: While all solid curves are for the reference liquid, the dashed and dotted lines are for the RelRes systems with varying r s . In the top panel for the reference liquid, we give ( r t − r 0 ) 2 as the black line, with its negative derivative −∂ t ( r t − r 0 ) 2 as the gray one.
In the bottom panel for the reference liquid, the orientational correlation s t · s 0 is given as the black line.
Besides these structural functions, we also present here a thermal function, particularly for the total energy, δ E t δ E 0 / δ E 2 , as the gray line. Of course, the functionality of the colored lines is the same as that of their respective neighboring curves.
In an all cases, irrespective of the value of r s , RelRes satisfactorily captures the transport functions of the reference system. Of course, there are some nuances between the various transport functions in terms of r s , and these dynamic characteristics are reminiscent of our observations for the static features. For example, r s = 1.7σ is required for a flawless replication of the translational ( r t − r 0 ) 2 , yet for the orientational s t · s 0 , only r s = 1.3σ can not excellently capture the behavior of the reference liquid. Besides, for the thermal function δ E t δ E 0 / δ E 2 , all switching distances yield a perfect description of this transport function. In summary, in terms of r s , we notice the same trends for dynamic behavior as we do with static behavior: Thermal properties can be superbly captured even with a surprisingly small r s , but structural correlations require a relatively large r s ; for the latter, we notice that its orientational correlations are rather feasible of capturing as compared with their translational counterparts. Specifically for these elementary dumbbells, it appears that r s = 1.5σ can capture the entire behavior of the reference liquid very well, and thus, this is our recommended switching distance for these molecules. This in turn reiterates one of our main findings of the original publication: RelRes works best if molecules interact with each via a FG potential between nearest neighbors and a CG potential between other neighbors.

Varying the Bond Length
We continue by examining dumbbells in which we vary the length of the bonds. In particular, we construct two systems, one with short dumbbells (i.e., ℓ = 0.3σ) and one with long dumbbells (i.e., ℓ = 0.7σ); representative sketches are given in Figs. 7 and 8, respectively. The modification of the bond length does not alter the LJ parameters of the FG and CG sites: They are again σ and ǫ for the former, and σ and 4ǫ for the latter. The system density is 1.34 for the short molecules and 0.78 for the long molecules, in units of m/σ 3 . The structural correlations for the short and long dumbbells are given in Figs. 7 and 8, respectively, which are basically analogous with those of Fig. 3. All orientational functions of the reference liquids are perfectly captured by RelRes, regardless of the switching distance that we use; we consequently give in the respective bottom panels the r s = 1.5σ curve. Conversely, the radial distributions are given in the respective top panels; notice importantly the cyan curve (i.e., r s = 1.1σ) for the ℓ = 0.3σ case and the magenta curve (i.e., r s = 1.9σ) for the ℓ = 0.7σ. In essence for each case, we have shifted our systematic r s examination by 0.2σ. This is because we observe that RelRes captures the radial distributions at a switching distance that differs by roughly 0.2σ as compared with the ℓ = 0.5σ of Fig. 3: The apparently asymptotic r s is 1.3σ in Fig. 7 and 1.7σ in Fig. 8. It is of course natural that the ideal r s for the short dumbbells is less than the ideal r s for the long dumbbells, since the former are more "sphere-like" than the latter. This further means that RelRes works better as the molecular size decreases. As a further analysis, we compute S JS for these radial distributions, plotting them in Fig. 4 as orange (downward) triangles for ℓ = 0.3σ and as green (upward) triangles for ℓ = 0.7σ. For these two cases, we again observe an essentially exponential decay. We also reaffirm that the efficacy of RelRes is more deficient for the long molecules than for the short molecules, considering that S JS for ℓ = 0.3σ is less than that for ℓ = 0.7σ by about an order of magnitude; remember that S JS → 0 means that we are approaching perfect replication of the radial distribution.
In the context of Fig. 3 for ℓ = 0.5σ, we have thoroughly discussed that the recommended value of r s is signaled by structural correlations. We find that this is still the case here once we vary the bond length of the dumbbells. Most importantly, in the bottom panels of Figs. 7 and 8, we again notice the flattening of ( s i · s j ) 2 beyond a certain distance, which stems in the mirroring behavior of the maxima and minima of its components s i s j 2 and s ⊥ i s ⊥ j 2 that cancel each other. We summarize all the critical distances associated with these orientational correlations, as well as with the translational ones in the respective top panels, in Tab. 1. Among these different bond lengths, we generally notice the same trends: The inflection in ( s i · s j ) 2 comes just before the inflection in g, and this is followed consecutively by the the minimum in s ⊥ i s ⊥ j 2 and the maximum in s i s j 2 , with the midpoint between the extrema in g being around there as well. For a given bond length, the most important aspect of these distances is that they almost all occur within 0.1σ of each other, and they are roughly the same as the estimate for the asymptotic r s via our systematic RelRes examination. As such, it appears that if one is interested in determining the ideal switching distance, there may be no need of constructing many RelRes systems with different r s : One can just look at the structural correlations of the reference liquid, with the radial distribution being perhaps the most feasible choice. As such, we consequently again draw vertical lines in Figs. 7 and 8 that represent their inflections in g. Nevertheless, we can even go further: For such dumbbells of an arbitrary bond length, we may not even need to do any molecular simulations of the reference liquid at all for determining r s , since it is quite obvious here that the following linear approximation holds: This further suggests that for an arbitrary system, one may be able to do a rough estimate for the ideal r s in RelRes just by considering the size of the respective molecules. Finally, we also calculate thermal properties for these two systems. Analogous with the elementary dumbbells, we plot their average U and variance δ U 2 in the bottom and top panels, respectively, of Fig. 5, with the orange (downward) triangles for ℓ = 0.3σ and the green (upward) triangles for ℓ = 0.7σ.
We again notice the fluctuating trend in all sets here. Importantly, notice that the magnitude of these fluctuations decreases with bond length. This reiterates the fact the RelRes works better for "sphere-like" molecules.

Examining Nonuniform Dumbbells
We now construct another system, whose molecules again have a bond length of 0.5σ . In this case however, we are dealing with nonuniform dumbbells in terms of their LJ parameters, although they do still have equal mass across their sites; a sketch of such a molecule is given in Fig. 9. Elaborating on this nonuniformity, one FG site is small yet strong, and one FG site is large yet weak: Their respective length parameters are 6/ 1 + √ 37 ≈ 0.85 and 1 + √ 37 /6 ≈ 1.18, in units of σ, and their respective energy parameters are 1 + √ 5 /2 ≈ 1.61 and 2/ 1 + √ 5 ≈ 0.62, in units of ǫ; notice that the mixed interaction between these two still has the standard σ and ǫ parameters. Using MPIL, we obtain that the ideal CG site has LJ parameters that are roughly 1.08σ and 2.76ǫ. Importantly, realize that even though these dumbbells are nonuniform, the system itself is still uniform, being a single-component and single-phase liquid. Besides, the density for this system is 0.88σ 3 /m. For this dumbbell system, we again perform a systematic examination for the switching distance. We present the corresponding structural correlations in Fig. 9, whose format is identical with that of Fig. 3. We again notice the same trends that we have observed for our elementary dumbbells. For the orientational functions of the bottom panel, the replication is perfect irrespective of the value of r s , just as we found in all of our other scenarios. Most importantly, for the radial distributions in the top panel, we find that the asymptotic r s is at 1.5σ, just as we have for our elementary dumbbells (i.e., ℓ = 0.5σ). Once we determine all the critical distances in g, as well as in ( s i · s j ) 2 with its components, we again find that they are all roughly the same as the ideal r s for RelRes; these are summarized in Tab. 1. Realize now that Eq. 50 still applies for this case, even though we established it while only considering uniform dumbbells. Finally, once we invoke our functionals for the radial distributions, we show in Fig. 4 that for all r s , S JS for this nonuniform scenario (i.e., the almost hidden yellow diamonds) is basically identical with its uniform counterpart (i.e., the indigo circles). All of these observations for the structural correlations of the nonuniform dumbbells, particularly the fact that their relationship with the switching distance of RelRes is identical with that of our elementary dumbbells (i.e., those with the same bond length of 0.5σ, yet which are uniform across their sites), has significant ramifications. For an arbitrary molecule which has much uniformity within it, we do not need to perform any molecular simulations for it for determining its ideal r s ; all one must do is to construct a system of similar molecules that are actually uniform across their sites (perhaps using the overall mean values for σ and ǫ of the nonuniform molecule). According to our observations, these uniform molecules, whose structural correlations are fairly feasible for computation, would roughly have the same r s as the nonuniform molecules. This fact becomes especially useful once we are dealing with a set of molecules that roughly have the same topology, yet whose sites notably vary in their LJ parameters: If the molecules do not have the linear estimate for r s as given in Eq. 50, we suggest that all one does is measure, via a single molecular simulation, the radial distribution of the reference liquid of the uniform version of these molecules, and we recommend using its inflection in g as the switching distance in RelRes for the entire family of molecules. We now proceed with the thermal properties of the nonuniform scenario in Fig. 5, which is presented as yellow diamonds. In this case, we actually do not observe the fluctuating trend which we have for all other dumbbell scenario, and instead we have a slight asymptotic behavior; of course, it is very modest as compared with the entropic functionals of Fig. 4. Still, this is not necessarily a drawback for RelRes: If we look at our recommended r s of 1.5σ, the variance in the top panel is off by a negligible fraction, yet the average of the bottom panel is off by a ∼ 0.1 fraction. Of course, the latter is still satisfactory, especially if one is not too interested in exact values of thermal properties. The discrepancy between the ability of RelRes in retrieving the average and variance may stem in the fact that the latter are related with response functions, and it is well known that these can be perfectly captured if the structural correlations are also perfectly captured, which is our case in this scenario. In fact, the main message here is that while structural correlations of nonuniform molecules can be as feasible of retrieving as those of uniform molecules, nonuniform energetics involve much intricacies with them, and thus a perfect replication is rather difficult to achieve.
Finally, we examine the transport functions of this nonuniform scenarios; we present them in Fig. 10, which is identical in format with Fig. 6. We again reaffirm our earlier observations. Foremost, RelRes perfectly captures the time behavior of the thermal property of the reference system, regardless of the switching distance. While the transport functions of the structural correlations are also adequately retrieved, there are some subtleties: For perfect replication, ( r t − r 0 )

Conclusion
In this current work, we continue our initial effort [1], presenting in turn a comprehensive picture of RelRes with MPIL. As shown in our original communication [1], our hybrid formalism can be thought of as an extension of the common approach in statistical mechanics which describes all interactions beyond a certain distance via a mean field [40,41,42,43]. The main aspect of RelRes is that molecular resolution switches in terms of relative separation: Each molecule embodies both FG and CG models, and a given molecule then interacts with near neighbors via its FG potential and with far neighbors via its CG potential.
Via energy conservation, we relate the FG and CG potentials by an arithmetic parameterization, which we call MPIL.
In our current publication, we specifically cast MPIL as a multipole approximation. Reminiscent of the familiar scenario, we do this derivation by assuming that the potential can be cast as a linear combination of inverse powers for the relative separation, and we then perform the common Taylor expansion. We call the leading nonzero term the MPIL, and we show that it naturally fits with RelRes. Note that in our original publication, we naively assumed infinite separation, in turn readily attaining MPIL [1]; this current work is in fact the comprehensive derivation, in which we specifically present the first and second terms of the Taylor series. In our original publication, we already showed that RelRes with MPIL works incredibly well for nonpolar molecules, as complex as a generic tetrachloromethane-thiophene fluid mixture which contains a vacuum cavity [1]. In the current work, we take a better examination at the role of the switching distance, finding an ideal value for it, which can be usually obtained by running a single molecular simulation of the reference mixture. Besides, we also show that RelRes can describe not only statics but also dynamics.
You may notice that RelRes with MPIL bears much resemblance with the famous "cell-multipole" formalism, the algorithm which is frequently mentioned as one of the most computationally powerful methods [45]. For an arbitrary potential of an inverse power, the "cell-multipole" approach also invokes a multipole approximation at appropriate distances: Inside a "cell" at small separations, one evaluates interactions between discrete sites, and outside a "cell" at large separations, one evaluates interactions between continuum patches. A crucial distinction with our multiscale scheme is that the "cell-multipole" strategy considers interaction sites that move freely of each other, in turn lumping them together in arbitrary "cells" in space.
In molecular simulations, interaction sites do not move freely of each other as they are constrained by bonds; RelRes, as well as many other multiscale methods [32,33,34,35], uses this signature of molecules, invoking a natural mapping from the FG model to the CG model (i.e., from atomistic coordinates to gravitational centers, respectively). As such, one can think of RelRes with MPIL as an natural modification of the "cellmultipole" approach for molecular simulations. This can be extremely useful, since, despite the fact that the "cell-multipole" method is frequently used in a myriad of systems (e.g., galaxies), it is rarely employed in molecular systems.
As a next step, we are interested in implementing RelRes for polar molecules, specifically for water.
Importantly, all the necessary mathematical ingredients have been already derived here. Unlike the nonpolar case, our model for water will have an orientational component, as manifested by its dipole. As we move to biology, we also expect to model various "chains". Of course, RelRes is only one route for enhancing the efficiency of molecular simulations, and we recommend that one employs it concurrently with other computational approaches for optimal results. For biological systems (e.g., a protein in water), we expect that combining RelRes with explicit-implicit strategies may be especially powerful: Reminiscent of Adaptive Resolution [46], these algorithms switch from an explicit solute near to the origin to an implicit solvent far from the origin [47]; the former is the one that is ideal for being modeled by RelRes. In summary, we emphasize that RelRes with MPIL can be of much use in aiding the study of soft matter via molecular simulations.

Acknowledgments
Foremost, we are thankful for partial funding by the Alexander von Humboldt Foundation, as well as the financial support of the Melvin J. Berman Hebrew Academy. We also appreciate the insightful conversations with Denis Andrienko regarding polar systems. Finally, we are grateful for the special services of Mark Wilson.
This is compactly presented in the main text as Eq. 34. For ℵ = 2 + , substituting Eqs. 66 and 68 in Eqs. 60 and 61, respectively, we obtain the following: For ℵ = 2 − , substituting Eqs. 67 and 69 in Eqs. 64 and 65, respectively, we obtain the following: Invoking Eq. 57, these are compactly presented in the main text as Eqs. 35 and 36. Besides, in an analogous manner which we performed in this Appendix for the first and second terms of c ij,ℵ , other terms of this coefficient can be evaluated as well. Figure 1: A pair of molecules in vacuum. The gravitational centers are represented by hollow rings, and the atomistic coordinates are represented by replete disks. The left side is for molecule i with its µ i sites, and the right side is for molecule j with its ν j sites. The various colors of these mean that they can all be inherently different in terms of their interaction parameters. The gray shading does not have much of a physical meaning: It just helps us delineate the approximate boundary of a molecule. The various arrows are distance vectors. − → r ij is the distance between centers i and j, and r µiνj is the distance between coordinates µ i and ν j . Within each molecule, − → ∆ µi or − → ∆ νj is a vector between an atomistic coordinate, µ i or ν j , and a gravitational center, i or j, respectively.

Reference System
RelRes System Figure 2: A schematic representation of our multiscale approach for nonpolar systems. The red and blue colors mean that we formally have two different molecules here. The top panel characterizes the geometrical FG and isotropic CG models on its left, and on its right, it illustrates that the former applies between atomistic coordinates, if their relative separation is small, and the latter applies between gravitational centers, if their relative separation is large. The bottom panel basically shows two molecular simulations, for the reference system on the left and for the RelRes system on the right. The arrow represents the MPIL parameterization, which makes the two systems approximately equal to each other. These molecules have uniform LJ parameters for their two sites, σ and ǫ, and the bond length between them is 0.5σ. Here, all the black curves, as well as the gray ones, correspond with the reference system. The colored curves, correspond with the RelRes strategy, with each color representing a different switching distance r s . Note particularly the violet curve, which is of the ideal switching distance in capturing the behavior of the reference liquid; the blue or red curve has a shorter or longer switching distance, respectively. Everything is plotted as a function of the distance between the midpoints of the dumbbells, r. The top panel plots the corresponding radial distribution g, with the solid curve being for the reference system, while the dashed curves being for the RelRes scenarios. The bottom panel, for the reference liquid, plots the orientational correlation ( s i · s j ) 2 as the black curve, together with its two main components s i s j 2 and s ⊥ i s ⊥ j 2 as the lower and higher gray curves, respectively; the functionality of the dotted curves corresponds with the functionality of their neighboring solid curves; note that all of these curves are almost identical. In both panels, the vertical brown line represents the inflection of g, which is located here at 1.46. Elementary Short Long Nonuniform Figure 4: The Jensen-Shannon entropy, S JS , for the radial distribution, plotted as a function of the switching distance r s in RelRes. Each curve is for a different system. The indigo color is for the molecules which have identical LJ parameters between their two sites, ǫ and σ, and the corresponding bond length is 0.5σ. The orange and green curves refer to the systems with bond lengths of 0.3σ and 0.7σ, while their LJ parameters are also uniform. The almost hidden yellow curve corresponds with {1.61ǫ, 0.62ǫ} and {0.85σ, 1.18σ} for its LJ parameters, while its bond is identical with that of the base case. Importantly, note that the entropy is in logarithmic scale. Besides, the lines here just serve as guides.  Figure 5: Energetic functionals plotted as a function of the switching distance r s in RelRes. Particularly, the average U and variance δ U 2 of the configurational energy is given in the bottom and top panels respectively; the asterisk denotes that these functionals are normalized by the respective value of the reference system. Besides, the color coding here is identical with that of Fig. 4. Everything is plotted in terms of time t. The top panel plots the squared displacement ( r t − r 0 ) 2 of a given molecule as the black curve for the reference system, together with its neighboring dashed lines for the RelRes systems; the corresponding negative derivative of this function −∂ t ( r t − r 0 ) 2 is given as the gray curve, together with its dotted lines. In an analogous manner, the bottom panel plots the orientational function s t · s 0 as the black curve, together with its neighboring dashed lines; the gray curve, together with its dotted lines, represents the energetic function δ E t δ E 0 / δ E 2 .    Table 1: Signature distances in the structural correlations of the various dumbbell scenarios. Rows 1-2 give us the molecule type, as well as its respective bond length, which each column corresponds with. Column 2 goes together with Fig. 3, columns 3 and 4 refer to Figs. 7 and 8, respectively, and column 5 goes with Fig.  9. Column 1 then tells us which distance value we are dealing with in each case. Row 3 is the switching distance that we recommend for use in RelRes, based on our systematic examination of varying r s . Rows 4-5 deal with signatures of the radial distributions, while rows 6-8 deal with signatures of the orientational functions. Realize that all numbers here are given in units of σ.