A comparative study of coarse-graining methods for polymeric fluids: Mori-Zwanzig vs. iterative Boltzmann inversion vs. stochastic parametric optimization

Although the molecular dynamics (MD) simulation has become a standard computational tool for studying molecular systems, an all-atom model is computationally prohibitive to access systems at large spatial-temporal scales.¹ Hence, the atomistic simulation is often unrealistic for many applications of biological systems and soft matter physics even at the mesoscale. To this end, coarse-grained (CG) approaches have been developed to overcome the limitation of scales. The main idea of coarse-graining is to aggregate atoms or molecules into clusters and to average out the irrelevant degrees of freedom for the cluster.^2,3 Consequently, the CG model represents a substantial simplification of the microscopic dynamics, but is still able to capture observable dynamics of complex fluids at larger spatial-temporal scales beyond the capability of conventional MD simulations.⁴ Among these CG methods, dissipative particle dynamics (DPD) conserves the momentum of a system and provides the correct hydrodynamic behavior of fluids at mesoscale,^5,6 which makes DPD one of the currently most popular mesoscopic methods.⁷ 

The DPD method was invented more than two decades ago for simulating complex fluids at mesoscale.⁸ Similar to MD, a DPD system consists of interacting particles and its dynamics is computed by integrating Newton’s equation of motion. But different from MD, it has a softer potential allowing for a much larger time step. Ever since its inception, the DPD method has found a wide spectrum of applications in soft matter including colloidal suspensions,⁹ smart materials,¹⁰ polymer solutions,¹¹ blood rheology,¹² and blood coagulation,¹³ to name but a few. The equations of motion for DPD particles are described as^14,15

where I, J are particle indices and P represents the momentum. The relative displacement and relative velocity of two particles I and J are defined as R_IJ = |R_IJ| = |R_I − R_J| and V_IJ = V_I − V_J, respectively; e_IJ = R_IJ/R_IJ is the unit vector along the radial direction of the two. All three forces are pairwise additive and $FIJC$⁠, $FIJD$⁠, and $FIJR$ are referred to as conservative, dissipative, and random forces, respectively. The total force on particle I is summed over other particles within a cutoff radius R_cut. Also, ω_C, ω_D, and ω_R are the weighting functions while the coefficients a, γ, and δ reflect the individual strength of the three forces. A Gaussian white noise θ_IJ with zero mean and unit variance is generated for the random force. In general, the functional form of $FIJC$ determines the static properties of the system, such as pressure, compressibility, and local structure often represented by the radial distribution function (RDF) of particles. The forms of $FIJD$ and $FIJR$ control the dynamic properties, such as viscosity, diffusivity, and time correlation functions. Moreover, the latter two forces together act as a thermostat and satisfy the fluctuation-dissipation theorem (FDT).¹⁴ This imposes δ² = 2γk_BT and $ωD(R)=ωR2(R)$⁠, where k_B is the Boltzmann constant and T is the temperature of the DPD system.

The essence of a valuable CG model is to obtain the functional forms of $FIJC$⁠, $FIJD$⁠, and $FIJR$ that describe the effective interactions between the DPD particles. To construct the effective DPD force fields from microscopic dynamics, operational strategies generally fall into two classes: a direct forward path by evaluating the microscopic dynamics via the Mori-Zwanzig (MZ) formalism and a reverse iterative path by solving an inverse problem targeting certain properties. In the forward-coarse-graining, the CG force fields between DPD particles are constructed directly from available MD trajectories via the Mori-Zwanzig projection.^16,17 While this procedure may require additional simplifying assumptions for computability, in principle there are no free parameters to be specified. In the reverse-coarse-graining, an initial guess of the CG force fields is posed and a few parameters in the expression of forces are left undetermined. Subsequently, an (iterative) inverse optimization is carried out to correct the free parameters so that target mesoscopic properties are obtained. The majority of available CG methods belong to this category, such as the force matching method,^18–22 inverse Monte Carlo,²³ inverse Boltzmann inversion (IBI),^24,25 stochastic parametric optimization (SPO) using Bayesian inference,^26,27 and minimization of relative entropy.²⁸ The fundamental question of coarse-graining is what parameters in the CG model are free to optimize, or are the parameters interchangeable? In the present work, we shall answer this question by coarse-graining a particular MD simulation as demonstration via both forward and reverse paths.

In the remainder of this paper, at first we will describe a reference MD system in Section II. Subsequently, in Section III we will construct CG models via the Mori-Zwanzig formalism, iterative Boltzmann inversion, and stochastic parametric optimization methods. Section IV presents the results and discussion. In particular, we examine how close the values for the pressure, compressibility, radial distribution function, viscosity, and diffusivity of the CG system match the corresponding values of the reference MD system. Finally, we conclude with a brief summary in Section V.

Molecular dynamics (MD) simulations of star polymer melts are performed to provide multiscale dynamics for coarse-graining. More specifically, star polymer molecules in the MD system are represented as chains of beads connected by short springs.^29,30 Each molecule has ten arms inter-connected with a center bead with two identical monomers per arm, and hence the total number of beads per molecule is N_c = 21 (see Fig. 1). Excluded volume interaction between monomers is described by Lennard-Jones (LJ) potential truncated for only repulsion, also known as the Weeks-Chandler-Andersen (WCA) potential,³¹ 

where r is the distance between two monomers and ϵ and σ set the energy and length scales, respectively. Moreover, neighboring monomers are connected by a finite extensible nonlinear elastic (FENE) spring³² 

where k = 30ϵ/σ² is the spring constant and r₀ = 1.5σ determines its maximum length.²⁹ Such a FENE spring is short and stiff enough to prevent neighboring bonds from crossing each other.³⁰ 

In the MD simulation, the polymer melt is modeled by 1000 molecules filled in a periodic cubic box of length L_B = (1000N_c/ρ)^1/3 and the number density of monomers is ρ = 0.4σ⁻³. Unless otherwise indicated, quantities of both MD and DPD simulations are reported in reduced LJ units, that is, the length, mass, energy, and time units are set as σ = 1, m = 1, ϵ = 1, and τ = σ(m/ϵ)^1/2. The system is coupled with an external heat bath via the Nośe-Hoover thermostat to sample a canonical ensemble.³³ In addition, the temperature of the system is maintained at k_BT = 1.0 and the velocity-Verlet integrator with a time step δt = 0.001τ is adopted.³³ 

In this section, we will describe how to construct mesoscopic models using different CG procedures, namely the Mori-Zwanzig formalism, iterative Boltzmann inversion, and stochastic parametric optimization methods.

A coarse-grained model can be constructed directly from the microscopic dynamics by applying the Mori-Zwanzig (MZ) projection operator.¹⁷ Starting from a Hamiltonian of the microscopic system and defining proper CG variables, the evolution of the CG variables can be described by a generalized Langevin equation (GLE) after projecting the microscopic system to the coarse-grained system.^{16,17,34–36} Upon a pairwise decomposition, this GLE can be cast into the DPD equation as^37,38

where F_IJ is the instantaneous force between CG particles I and J. The ensemble average 〈F_IJ〉 is taken as the conservative force in the CG model. The last term $δFIJQ$ is the random force due to unresolved degrees of freedom. In general, $δFIJQ$ is non-Gaussian and the corresponding dissipative force depends on an integral of the past history of motion weighted by a memory kernel given as $KIJ(t)=(kBT)−1〈[δFIJQ(t)][δFIJQ(0)]T〉$⁠, which ensures that the CG system satisfies the second FDT.³⁹ 

The momentum of center-of-mass (COM) is often a slow variable due to its inertia, while the random force is a fast variable due to frequent atomic collisions. When the relaxation time scales of the COM’s momentum and the random force are clearly separable, the convolution integral in the expression of dissipative force can be further simplified using a Markovian approximation by replacing the memory kernel with the Dirac delta function, that is, $∫0tKIJ(t−s)VIJ(s)ds=γIJVIJ(t)$⁠. Furthermore, if only the radial components of CG interactions are considered, Eq. (4) reduces to the classical DPD formulation in Eq. (1).

However, a molecule consists of many discrete monomers in the MD system, therefore, the total force F_IJ contains not only the radial components but also the perpendicular components. Including both radial and perpendicular interactions in the CG description leads to a full DPD (FDPD) model.³⁷ It is worth noting that the CG interactions are evaluated directly from MD simulations and there are no iteratively optimized parameters in the MZ-guided DPD models. For further technical details we refer to Refs. 35–38.

Iterative Boltzmann inversion is primarily used to obtain non-bonded interactions for reproducing the RDF of a reference system.²⁴ Given a RDF g_ref(R), which is often obtained from all-atom MD simulations or from experiments, the potential of mean force (PMF) U_PMF(R) between pairs of CG particles as a function of their distance R can be obtained through the Boltzmann inversion²⁴ 

which is a free energy and cannot be directly used as a pairwise potential in a CG model.⁴ However, U_PMF(R) is usually sufficient to serve as an initial guess, $U0CG(R)$⁠, for the pairwise CG potential during an iterative procedure. The first CG simulation performed with $U0CG(R)$ yields a corresponding RDF g₀(R), which is different from the target g_ref(R). Thereafter, a correction term $αkBTlng0(R)/gref(R)$ is added to improve the potential and this procedure is performed iteratively as⁴⁰ 

where the subscript i denotes the iteration number, and α is a scaling factor that improves the convergence rate and stability of the IBI process.⁴⁰ Due to a fundamental theorem,⁴¹ for a given RDF the pairwise potential is unique up to a constant. Therefore, the iterative procedure reaches a converged potential that generates the target RDF well. However, in practical implementations the choice of cutoff radius R_cut is arbitrary and as a consequence, there can be as many CG potentials as the choices of R_cut, each of which reproduces the given target RDF well.

If the pressure is also of concern, only one of those potentials with the correct R_cut can generate the correct pressure as that of the reference MD. If a different R_cut is selected initially, a typical strategy is to add a linear term to the potential^24,25,42 to correct the pressure at the expense of the accuracy of the RDF. However, this correction leads to a significant deviation of the compressibility from the target value.^25,42 From the discussion above, it is reasonable to hypothesize that for a given MD system the R_cut in the CG potential is not an arbitrary parameter. Once R_cut is selected falsely, the desired thermodynamic properties may not be recovered well by other means.

To evaluate this hypothesis explicitly, we will select three different choices of R_cut = 5.23, 7.35, and 10.0, which are denoted by IBI-1, IBI-2, and IBI-3, respectively. The first value, R_cut = 5.23, is chosen based on the MZ-guided interaction range, the second one R_cut = 7.35 takes the value of the second peak of the given RDF g_ref(R), and the last one, R_cut = 10.0, is the longest range of RDF that we consider in the present study. The reference RDF, g_ref(R), of COM of star polymers is obtained by running MD simulations. Then, the IBI procedure is implemented in the VOTCA⁴³ package while coarse-grained runs are performed using GROMACS.⁴⁴ 

Stochastic parametric optimization (SPO) is another reverse method to obtain an effective CG potential. The SPO method differentiates itself from IBI by assuming empirical functions with undetermined parameters for the interaction potential. Then, a stochastic parametric optimization is performed to determine these parameters so that designed target properties can be obtained by the optimized CG potential.^19,27,45 For static properties we set the RDF of COM of star polymers and the pressure of the reference MD system as our target properties. Here, we choose two different empirical functions for the CG potential: SPO-1 has a form of U₁(R) = 5.23a₁(s₁ + 1)⁻¹(1 − R/5.23)^s₁+1 with two free parameters a₁ and s₁, while SPO-2 has U₂(R) = 0.5a₂R_cut(1 − R/R_cut)² with two free parameters a₂ and R_cut. The potential U₂(R) corresponds to a conservative force being a linear function of distance, which is widely used in classic DPD simulations.¹⁵ 

Technically, we employ generalized polynomial chaos (gPC)⁴⁶ to construct a surrogate model for DPD systems using a linear combination of a set of special basis functions defined in the parameter space.²⁷ The two parameters of DPD model to be optimized are given by

where $(ā1,s̄1)=(150,3.5)$ with (δ_a1, δ_s1) = (30, 1.0) for SPO-1, and $(ā2,R̄cut)=(55,4.5)$ with (δ_a2, δ_Rcut) = (50, 1.0) for SPO-2. Here, ζ₁ = (ξ₁, ξ₂) and ζ₂ = (ξ₃, ξ₄) are i.i.d uniform random variables distributed on [ − 1, 1]. To obtain the optimal parameter set, we first generate 65 samples of random variables based on sparse grid method.⁴⁷ Then, we construct the surrogate model using the 65 samples of DPD simulations to infer two parameters in the DPD model for achieving the pressure and the RDF of the reference MD system.

In this section, a quantitative evaluation of the three different CG models by comparing the CG systems with the reference MD system will be presented and discussed, including both static properties (i.e., the radial distribution function, pressure, and compressibility) and dynamic properties (i.e., diffusivity, viscosity, and velocity autocorrelation function) of the polymeric fluid.

The conservative force is responsible for the static properties of molecular fluids. Figure 2 shows the distance-dependent pairwise potential from the MZ formalism, which is defined as the spatial integration of the mean force 〈F_IJ(R)〉, that is, $UMZ(R)=∫R∞〈FIJ(r)〉dr$⁠. Here, the cutoff radius R_cut = 5.23 is determined as the distance beyond which the pairwise interactions between CG particles vanish.^37,38 〈F_IJ(R)〉 has no data available at short distances due to the fact that pairs become improbable at small distance,³⁷ as shown in Fig. 3. In practice, we use a bell-shaped function to extrapolate 〈F_IJ(R)〉 at short distances, which is presented in a global view of the MZ-guided pairwise potential in the inset of Fig. 2.

The resultant RDF of the MZ-DPD system as well as the RDF of COM obtained from the reference MD system are plotted in Fig. 3. To quantify the deviation of a RDF from the reference RDF, we define a ℓ₂-norm

where L is a length chosen for comparison of the two RDFs. Here, we take the value of the cutoff distance L = R_cut = 5.23 to penalize deviations at small distances. The MZ-guided pairwise potential results in $ℓ2RDF=1.20%$⁠, which reveals that the reference RDF of the MD system has been well reproduced by the MZ-guided DPD model, as shown in Fig. 3. The values of $ℓ2RDF$ and other comparisons for all cases are summarized in Table I.

In general, the deviation of RDF obtained by IBI from the target RDF decreases with the iteration step, but eventually the deviation saturates after some number of iterations.⁴³ Specifically, we perform 300 iterations for updating the IBI potential, and the resultant potentials for different cutoff radii are shown in Fig. 2. The inset of Fig. 2 gives a global view of these potentials, where the negative parts of the potentials of IBI-2 and IBI-3 are displayed by lines without symbols, while the potential of IBI-1 with a cutoff radius of R_cut = 5.23 has no negative part. The potentials obtained by IBI are assumed exponential functions for R_cut < 2.2 because the RDF becomes zero at short distances.

It is shown in Fig. 3 that all the IBI potentials (IBI-1, IBI-2, IBI-3) can reproduce well the RDF of the reference MD system with $ℓ2RDF$ within 0.3% as listed in Table I. However, the pressures from the three CG systems are apparently different. The IBI-1 has a pressure 4.0% higher than that of the reference MD system, while pressures of the IBI-2 and IBI-3 deviate from the reference by +10.1% and −38.9%, respectively. The errors on the static properties of different CG systems are summarized in Table I. These results reveal that the pressure of the reference system can be reproduced in the CG system with a correct cutoff radius without deteriorating the quality of RDF or the compressibility.

For the CG procedure using SPO, we define the relative error of static properties as $Δ=|P−Pref|/Pref+ℓ2RDF$⁠. Figure 4(a) shows the response surface of the relative error function Δ(a₁, s₁) of SPO-1 in the parameter space, while Δ(a₂, R_cut) of SPO-2 is shown in Fig. 4(b). The optimal parameter set for achieving the target static properties is (a₁, s₁) = (144.72, 3.27) for SPO-1 and (a₂, R_cut) = (21.00, 4.38) for SPO-2, where the relative error function reaches the minimum in the parameter space, as shown in Fig. 4.

We find that the best result of SPO-2 is Δ_min = 16%, which implies significant deviations of pressure and RDF from the targets, i.e., P = 0.208 is 5.0% higher than the pressure P = 0.198 of MD system and $ℓ2RDF=11.3%$ indicates an obvious deviation of RDF. However, when the cutoff radius is properly selected, SPO-1 with R_cut = 5.23 yields its best result with Δ_min = 2% indicating a good reproduction of both pressure and RDF. The only difference between SPO-1 and SPO-2 is that for SPO-1 we fixed R_cut while for SPO-2 we took R_cut as a tunable parameter. Results show that the best CG model of SPO-1 reproduces well the reference MD system on static properties but the best SPO-2 contains significant deviations of pressure and RDF from targets, which implies that the cutoff radius should not be taken as a free parameter in CG modeling of polymeric fluids, similarly to the findings in the construction of CG model using IBI.

The dimensionless compressibility is defined by κ⁻¹ = [L]³/ρk_BTκ_T, where κ_T is the usual isothermal compressibility of the fluid, [L] the length scale, and ρ the number density of particles in a volume of [L]³. In particle-based systems, the value of κ⁻¹ can be computed by¹⁵ 

where ρ^∗ = ρ/N_c is the number density of DPD particles or the number density of molecules in the MD system. To obtain the equation of state (EOS) of these systems the monomer density is varied from ρ = 0.3 to ρ = 0.6 in a step of δρ = 0.025. Figure 5 shows the dependence of pressure P on the monomer density ρ for the MD system and the CG systems. We note that all the CG potentials are obtained based on the MD system at ρ = 0.4. According to Eq. (9), the compressibility κ⁻¹ of these systems can be obtained based on the gradient ∂P/∂ρ of the EOS at ρ = 0.4, as listed in Table I.

Figure 5 shows that the CG potentials of MZ, IBI-1, and SPO-1 have similar performance that is consistent with the equation of state of the reference MD system for ρ < 0.45. Although the coarse-graining strategies of MZ, IBI, and SPO are significantly different, the pairwise potentials obtained by MZ, IBI-1, and SPO-1 converge to the optimal one when the cutoff radius of CG interactions is defined properly. Here, we note that both IBI-1 and SPO-1 use the cutoff radius R_cut = 5.23, which is obtained from the MZ formulation as a guidance for construction of effective CG models. Therefore, the CG models of IBI-1 and SPO-1 yield good results. However, when IBI-2 and IBI-3 use an arbitrarily selected R_cut while SPO-2 takes R_cut as a free parameter to be optimized, none of them can yield a good CG model. To this end, the MZ-guided CG model has no free parameters and is able to provide accurate information for CG representations, which can be taken as necessary guidance for performing reverse CG procedures more effectively.

It is observed in Fig. 5 that all the EOSs of CG models using the pairwise potentials obtained at ρ = 0.4 diverge from the MD system as the monomer density increases. The reason is that if the CG clusters are soft in the MD system and can change their configurational morphology as the polymer melt becomes dense, the effective interactions between molecules depend on all the neighboring molecules of the star polymer, and hence the many-body interactions should be considered in a CG model. However, when we construct the CG models and compute the pairwise potentials at ρ = 0.4, we do not include the many-body effect in the pairwise potential. To quantify the deformation of CG clusters in the MD system as the monomer density changes, we compute the probability density functions (PDFs) of their gyration radius R_g defined by

where M is the mass of a CG cluster, m is the mass of a monomer, and $rˆi=ri−R$ is the relative displacement of a monomer with respect to the COM of the cluster.

Figure 6 illustrates the PDF of R_g for two coarse-graining levels, N_c = 11 and 21, as the monomer density ρ increases from 0.3 to 0.6. For the star polymer with short arms N_c = 11, the cluster has less deformability in the MD system. Therefore, the PDF of R_g does not change as the monomer density increases, which indicates that the many-body effect is not important for the system of N_c = 11. Consequently, the pairwise potential of N_c = 11 obtained at ρ = 0.4 can be safely applied to the system at other monomer densities and reproduces well the reference MD system, as shown in the inset of Fig. 5. However, the star polymer with long arms, N_c = 21, is soft and deforms easily as the monomer density changes, which corresponds to an obvious shift on PDF of R_g shown in Fig. 6. The deformability of cluster makes the many-body effect non-negligible, and hence all the CG models of N_c = 21 using pairwise potentials without many-body corrections diverge from the MD system shown in Fig. 5. This result suggests that many-body corrections should be included for achieving a transferable CG potential when CG clusters are soft and can deform significantly. Otherwise, the constructed CG potential is only valid near the training point and cannot be applied to other systems.

In general, the CG potential alone cannot produce the correct dynamic properties, which are irrelevant for the reproducibility of structural correlations.⁴⁸ Izvekov and Voth²⁰ reported that Hamiltonian mechanics on a CG potential surface yields faster diffusion dynamics than its underlying all-atom system. This faster CG dynamics arises from the fact that the effective frictional forces, which are induced by the effects of unresolved degrees of freedom, are not considered in the CG representation. Having obtained the interaction potential, the dynamic properties of a CG system are determined by the dissipative and random forces. In this respect, the IBI method starts from RDF to obtain effective CG potentials that produce correct static properties but does not work for dynamic properties, while the SPO method can work for targeting both static and dynamic properties. Similarly to the optimization of two parameters for desired pressure and RDF, we can optimize the dissipative force in the form of $FIJD(R)=−γ(1−R/Rcut)s2 (eIJ⋅VIJ)eIJ$ with two tunable parameters γ and s₂ to target the diffusivity and viscosity of the MD system.

We define the relative error of dynamic properties as Δ = |ν − ν_ref|/ν_ref + |D − D_ref|/D_ref. Figure 7 shows the response surface of Δ(γ, s₂) in the parameter space, which is obtained by performing the stochastic parametric optimization process. Let φ(R) = γ(1 − R/R_cut)^s₂, we have the optimal parameter set giving φ(R) = 144.0(1 − R/5.23)^2.42 for SPO-1 and φ(R) = 56.0(1 − R/4.38)^1.01 for SPO-2, which result in (D, ν) = (0.065, 1.433) of SPO-1 and (D, ν) = (0.064, 1.456) of SPO-2 with respect to (D, ν) = (0.061, 1.413) of the reference MD system.

According to the Green-Kubo relations,⁴⁹ the transport properties of fluids can be expressed in terms of integrals of time correlation functions. Although the integral properties such as the diffusivity and the viscosity of SPO-1 and SPO-2 approximate the values of the MD system, as shown in Table I, Figure 8 reveals that the time correlations in SPO-1 and SPO-2 systems are significantly different from that of the MD system. When the inverse methods using optimization are employed to construct CG models, even if the target properties are achieved, one should not expect that other behavior besides the targeted ones can be correct automatically. In contrast, the MZ-based CG procedure does not set target properties before constructing a CG model as it uses a forward path to compute the effective CG interactions from MD simulations. Since the CG interactions are correctly considered in the MZ-guided model, the time correlations and also the integral properties of the CG system agree well with its underlying MD system. For CG models derived from the reverse CG strategies, only if the time correction functions were defined as target properties, the dynamics of these CG systems would have correct time evolution.

We have constructed coarse-grained (CG) models for polymeric fluids using both forward and reverse coarse-graining strategies. The forward-coarse-graining procedure eliminates irrelevant atomistic variables by using the Mori-Zwanzig (MZ) projection operators and computes the CG interactions directly from a microscopic dynamics. In contrast, the reverse-coarse-graining procedure obtains the effective CG interactions by solving inverse problems, such as the iterative Boltzmann inversion (IBI), starting from a given radial distribution function (RDF) and iteratively ending with an effective pairwise potential, and the stochastic parametric optimization (SPO) optimizing undetermined parameters for achieving targeted properties.

In particular, we performed molecular dynamics (MD) simulations of star polymer melts to provide the fields to be coarse-grained. In the CG representation, the internal degrees of freedom of each star polymer are averaged out and the entire molecule of the star polymer is coarsened into a single CG particle. Then, different CG methods, i.e., MZ, IBI, and SPO, were employed to construct an effective CG model for representing the MD system at a cheaper coarse-grained level. Quantitative comparison between these CG models indicates that both the forward-coarse-graining and reverse-coarse-graining methods are able to yield an effective CG model that recovers well the reference MD system if the free parameters in a CG model are properly selected. Some quantities, however, such as the cut-off radius, have a strong influence and cannot be considered as arbitrarily tunable parameters.

To construct the effective pairwise potential, the three CG techniques have similar computational efficiency. In particular, the MZ method requires several MD simulations for ensemble average of the mean force, while the IBI method needs many iterations with several DPD simulations in each iteration, and the SPO method needs many samples of DPD simulations to construct the surrogate model. For the effective dissipative and random forces resulting in correct dynamic properties, only the MZ and the SPO techniques can be applied. In general, the MZ method is more expensive than the SPO method. The reason is that we have to run a lot of MD simulations for ensemble average in the MZ procedure so that an accurate memory kernel can be obtained. However, despite the inaccurate time correlation functions, the SPO procedure only requires many samples of DPD simulations to construct the surrogate model. Moreover, recently developed methods can substantially reduce the number of DPD simulations required for constructing the surrogate model in SPO.^50–52 

We also explained how the coarse-graining procedure introduces the many-body effect that makes pairwise potentials highly specific to conditions under which the potentials are obtained. We have demonstrated that coarse-graining rigid CG clusters with less deformability does not introduce the obvious many-body effect, and consequently the CG pairwise potential can be safely applied to other conditions beyond the training point. However, in aggressive coarse-graining, the CG clusters are soft and can deform significantly. As a result, the configurational morphology of soft CG clusters depends on all neighboring clusters and their configurations, which introduces the many-body effect that cannot be ignored in CG representations. In order to obtain CG models being transferable to a different density, temperature, or composition, the many-body corrections should be incorporated into the CG models at aggressive coarse-graining levels. To this end, it would be interesting to consider the density-dependence of the effective CG potentials by using the framework of many-body DPD,^53,54 where the pairwise interaction depends not only on the distance but also on the density of neighboring particles.

We acknowledge support from the DOE Center on Mathematics for Mesoscopic Modeling of Materials (CM4). This work was also sponsored by the U.S. Army Research Laboratory and was accomplished under Cooperative Agreement No. W911NF-12-2-0023. An award of computer time was provided by the Innovative and Novel Computational Impact on Theory and Experiment (INCITE) program. This research used resources of the Argonne Leadership Computing Facility, which is a DOE Office of Science User Facility supported under Contract No. DE-AC02-06CH11357. This research also used resources of the Oak Ridge Leadership Computing Facility, which is a DOE Office of Science User Facility supported under Contract No. DE-AC05-00OR22725. Z. Li would like to thank Professor Bruce Caswell and Dr. Xuejin Li for helpful discussions.