The Mori-Zwanzig formalism for coarse-graining a complex dynamical system typically introduces memory effects. The Markovian assumption of delta-correlated fluctuating forces is often employed to simplify the formulation of coarse-grained (CG) models and numerical implementations. However, when the time scales of a system are not clearly separated, the memory effects become strong and the Markovian assumption becomes inaccurate. To this end, we incorporate memory effects into CG modeling by preserving non-Markovian interactions between CG variables, and the memory kernel is evaluated directly from microscopic dynamics. For a specific example, molecular dynamics (MD) simulations of star polymer melts are performed while the corresponding CG system is defined by grouping many bonded atoms into single clusters. Then, the effective interactions between CG clusters as well as the memory kernel are obtained from the MD simulations. The constructed CG force field with a memory kernel leads to a non-Markovian dissipative particle dynamics (NM-DPD). Quantitative comparisons between the CG models with Markovian and non-Markovian approximations indicate that including the memory effects using NM-DPD yields similar results as the Markovian-based DPD if the system has clear time scale separation. However, for systems with small separation of time scales, NM-DPD can reproduce correct short-time properties that are related to how the system responds to high-frequency disturbances, which cannot be captured by the Markovian-based DPD model.

## I. INTRODUCTION

Atomistic simulation techniques allow precise representation of molecular structures and have become standard computational tools nowadays for studying molecular systems.^{1} In particular, direct molecular dynamics (MD) simulation can capture all the atomistic details of a molecular system, however, it is computationally prohibitive and impractical to produce large-scale or long-time effects in many applications of biological systems and soft matter physics.^{2–4} The reason is that the size of a MD system is limited by the number of atoms that can be included in the simulation, typically 10^{4}-10^{8} corresponding to a length-scale on the order of tens of nanometers, and the maximum time step in MD simulations is limited by the smallest oscillation period of the fastest atomic motions in a molecule, which is typically several femtoseconds.^{5} However, it may not be necessary to explicitly take into account all the atomistic details of materials when only mesoscopic structures of molecules or their macroscopic properties are of practical interest. To this end, coarse-grained (CG) approaches drastically simplify the atomistic dynamics by eliminating fast degrees of freedom and preserving the behaviors of slow qualities.^{2,3} As a result, the CG modeling provides an economical simulation path to capture observable properties of complex fluids on larger spatial and temporal scales beyond the capability of conventional atomistic simulations. With increasing attention on macromolecules, biomaterials, and soft matter research, CG modeling has become a rapidly expanding methodology in recent years.^{2–4}

Specification of the interactions between CG variables is the basis for construction of CG models. Given a target static property, such as the radial distribution function (RDF) for the center-of-mass (COM) of each molecule in a fully atomistic simulation, an effective CG potential can be obtained numerically by solving an inverse problem via optimization, for example, the iterative Boltzmann inversion method,^{6} the inverse Monte Carlo method,^{7} the force matching method,^{8} or the multiscale coarse-graining methodology.^{9} Moreover, in a particle-based system in a canonical (constant-temperature) ensemble, thermostats are often used with the purpose of maintaining the temperature constant (on average).^{10} In general, the dynamic properties of molecular fluids relate closely to the parameters of applied thermostat. To achieve certain averaged dynamic properties, such as the kinematic viscosity and the diffusion coefficient of the molecule, the parameters of the thermostat are usually fine-tuned in CG simulations.^{11} However, besides the averaged dynamic properties, there is no guarantee that the artificially tuned thermostat reproduces a faithful dynamic property of its corresponding atomistic simulation, such as, the velocity autocorrelation function (VACF) of the molecule. Contrary to the inverse techniques, the Mori-Zwanzig (MZ) projection formalism^{12,13} provides a forward path to constructing CG models directly from microscopic dynamics. By defining CG variables and eliminating irrelevant degrees of freedom, the CG interactions can be directly evaluated from the microscopic dynamics by mapping the microscopic system to a CG system using the MZ projection operator. In particular, the collective dynamics resulting from a projection of an underlying microscopic dynamics is governed by three terms.^{13} The first term represents a conservative interaction due to the change of microscopic configuration. Another term is a random interaction generated by unresolved variables during the coarse-graining process. The last one is a history-dependent term determined by an integral of a memory kernel. The three components of CG interaction are consistent with the framework of dissipative particle dynamics (DPD).^{14} Therefore, in the present work we use DPD as the CG model resulting from a projection of an underlying microscopic dynamics.

In the recent past, several research groups have practiced coarse-graining under the MZ formalism.^{15–21} In particular, a Markovian assumption is often employed so that the implementation of the MZ formalism can be significantly simplified. This simplification is reasonable for many applications, where the time scales are well separated. However, the Markovian assumption is not universally valid and its corresponding coarse-graining implementation may unfaithfully represent the original MZ formalism for certain applications. To this end, Yoshimoto *et al.*^{21} derived the equation of motion for non-Markovian DPD. For construction of non-Markovian CG models being a faithful representation of its underlying microscopic model, an accurate evaluation of the memory kernel is crucial. Furthermore, since the thermodynamic properties of molecular fluids highly depend on the temperature, it is also crucial to implement the second fluctuation-dissipation theorem (FDT)^{22} correctly in non-Markovian CG simulations so that the temperature of the system can be maintained constant. In the present work, we apply two different approaches to evaluate the memory kernel from the trajectories of atomistic simulations. Furthermore, for the CG simulations, we reproduce the non-Markovian property in DPD, and two colored noise generators for satisfying the second FDT are provided and compared. As an example, a MD system of star polymer melts whose molecules can be faithfully coarsened into CG particles is selected as demonstration for the implementation. The VACF of molecules is considered as a crucial dynamic property for validation. At first, we shall show if the non-Markovian DPD (NM-DPD) can recover the results of Markovian DPD and MD at low density (*ρ* = 0.4), where the time scales of VACF and force autocorrelation function (FACF) in the MD system are apparently separated. Subsequently, we shall show if the NM-DPD can represent the MD more faithfully than the Markovian DPD at high density (*ρ* = 0.7), where the two time scales of VACF and FACF are of the same order.

The remainder of this paper is organized as follows: In Section II, we briefly introduce the theoretical background for projecting a microscopic system to a coarse-grained system using the MZ formalism. Section III describes in detail how to practically evaluate the CG interaction and the memory kernel directly from atomistic simulations. Section IV presents the scheme for incorporation of the memory effects into a DPD model and provides quantitative comparisons between Markovian DPD and NM-DPD, and also their accuracy in reproducing the MD system. Finally, we conclude with a brief summary and discussion in Section V.

## II. THEORETICAL BACKGROUND

Consider an atomistic system consisting of numerous atomic particles. The idea of coarse-graining is to divide the atomistic system into *K* clusters where each cluster contains many atomic particles. If the atomistic details are not of practical interest, the dynamics of the system can be described by the evolutions of proper CG variables, i.e., the coordinate **R** and momentum **P** of the COM of the clusters. The equation of motion (EOM) of the CG particles obtained from the Mori-Zwanzig projection is in a form of the generalized Langevin equation,^{18–21,23}

where *β* = 1/*k _{B}T* with

*k*the Boltzmann constant and

_{B}*T*the thermodynamic temperature. Also,

**R**= {

**R**

_{1},

**R**

_{2}, …,

**R**

_{K}} is a phase point in the CG phase space, and

*ω*(

**R**) is defined as a normalized partition function of all the microscopic configurations at phase point

**R**. On the right-hand side of Eq. (1), the first term represents the conservative force due to the change of potential energy. The last term $\delta FIQ(t)$ is the random component, which considers the unresolved degrees of freedom during the CG procedure. The second term of Eq. (1) is the friction force containing an integral of memory kernel $\Psi (t)=\beta [\delta FIQ(t)][\delta FXQ(0)]T$.

In principle, the force on a molecule *I* depends on all the COM coordinates **R** as well as their microscopic configurations. Although the EOM given by Eq. (1) based on the MZ formalism is accurate, a direct computation of the many-body interactions and the memory kernel in Eq. (1) is difficult, even for a simple system derived from a one-dimensional harmonic chain.^{17,24} To this end, we use two approximations to make Eq. (1) resulting from the MZ formalism practical. First, we assume that the non-bonded interactions between neighboring clusters in the microscopic system are pairwise decomposable, and hence the total force consists of pairwise forces, e.g., $FI\u2248\u2211J\u2260IFIJ$ and $\delta FIQ\u2248\u2211J\u2260I\delta FIJQ$. Second, we neglect the many-body correlations between different pairs, and assume that the pairwise force **F**_{IJ} between two clusters *I* and *J* depends only on the relative COM positions **R**_{I} and **R**_{J} and is independent of the positions of the rest of clusters.

With negligible many-body correlation and using the pairwise decomposition, Eq. (1) can be written into pairwise form (the detailed derivation is provided in Appendix A),

where **F**_{IJ} is the instantaneous force and its ensemble average is denoted as 〈**F**_{IJ}〉, and the matrix $KIJ(t)=\beta \u3008[\delta FIJQ(t)][\delta FIJQ(0)]T\u3009$ is the pairwise memory kernel of dissipation. Let the difference of the instantaneous force from the mean force be the fluctuating force *δ***F**_{IJ}(*t*) = **F**_{IJ}(*t*) − 〈**F**_{IJ}〉. By multiplying $VIJT(0)$ on both sides of Eq. (2), we have

where the third term $\u3008\delta FIJQ(t)VIJT(0)\u3009$ disappears because $\delta FIJQ(t)$ has no correlation with velocity. The reason is that the random force $\delta FIJQ(t)$ comes from the orthogonal dynamics in the MZ procedure, and hence it is viewed as noise uncorrelated with resolved components. The temporal correlations $\u3008\delta FIJ(t)VIJT(0)\u3009$ and $\u3008VIJ(t)VIJT(0)\u3009$ can be computed from atomistic simulations, and hence the memory function **K**_{IJ}(*t*) for pairwise systems can be obtained by solving Eq. (3). Moreover, since the random force $\delta FQ(t)=e\u2212QiLt\delta F$ where $Q$ is the orthogonal projection operator in the MZ procedure and *L* the Liouville operator,^{18} using zeroth-order approximation,^{25} e.g., $e\u2212QiLt\u2245e\u2212iLt,$ we can compute the force-force correlation $CIJ(t)=\beta \u3008\delta FIJ(t)\delta FIJT(0)\u3009$ as an alternative approximation of the memory kernel **K**_{IJ}(*t*). We note that **C**_{IJ}(*t*) is the zeroth-order approximation of **K**_{IJ}(*t*), which could be valid for short memory but would become questionable for long memory.^{25,26} In Section IV, the difference between using **K**_{IJ}(*t*) and **C**_{IJ}(*t*) in NM-DPD simulations will be discussed.

The second fluctuation-dissipation theorem^{22} requires that the random force $\delta FIJQ(t)$ in Eq. (2) is temporally correlated with the correlation function **K**_{IJ}(*t*),

which is also the definition of the memory kernel in Eq. (2). Generation of colored noises for the random force satisfying the second FDT is crucial to maintain the NM-DPD system at a desired temperature, which guarantees that the CG simulation is able to reproduce the correct CG dynamics consistent with its underlying atomistic system. In the present paper, two different colored noise generators are tested and their performances in CG simulations are discussed in Appendix B.

## III. COARSE-GRAINING PROCEDURE

To demonstrate the coarse-graining procedure, we construct MD systems of star polymer melts to provide the multiscale dynamics to be coarse-grained. In particular, the entire molecule of a star polymer in the MD system is coarse-grained into a single CG particle during the coarse-graining process, whereby the MD system is projected into a corresponding CG system, as shown in Fig. 1. Then, the effective interactions between CG particles are evaluated through the method described in Section II.

### A. Microscopic model

In practice, MD systems are constructed in a computational box filled with star polymers, which are represented as chains of beads connected by short springs. Each molecule of the star polymer has 10 arms with 1 monomer per arm, hence the total number of monomers per star polymer is *N _{c}* = 11. We note that increasing the arm length of star polymer will introduce many-body correlations between different pairs, which makes the interpretation of results complicated. In the present study, we focus on the failure of the Markovian approximation induced by poorly separated time scales in the multiscale dynamics. Therefore, we simplify the problem and exclude the many-body correlations by choosing the star polymer with short arms.

A purely repulsive Weeks-Chandler-Andersen (WCA) potential^{27} is employed to compute the excluded volume interactions between monomers,

where ε sets the energy scale and *σ* the length scale of the monomers. Also, a finitely extensible nonlinear elastic (FENE) potential^{28} is adopted for the bond interaction between connected monomers,

where *k* = 30ε/*σ*^{2} is the spring constant and *r*_{0} = 1.5*σ* determines the maximum length of the spring.^{29} Such a FENE spring is stiff and short enough so that bond crossing can be prevented in the MD system.^{30}

MD simulations of star polymer melts are performed in a canonical ensemble (*NV T*) with the Nosé-Hoover thermostat.^{31,32} More specifically, 1000 molecules of star polymer are filled into periodic cubic boxes of length *L* = (1000*N _{c}*/

*ρ*)

^{1/3}to model the polymer melts, in which

*ρ*is the number density of monomers. Unless otherwise specified, all the results in the present paper are expressed in the reduced LJ units, i.e., the length, mass, energy, and time units are set as

*σ*= 1,

*m*= 1,

*ϵ*= 1, and

*τ*=

*σ*(

*m*/ε)

^{1/2}, respectively. Also, the temperature of the MD systems is maintained at

*k*= 1.0 and the time step used in the velocity-Verlet integrator

_{B}T^{10}is

*δt*= 1.0 × 10

^{−3}

*τ*.

To examine the importance of memory effects in the coarse-graining procedure, two cases are considered: one has apparent time scale separation while the other has small separation of time scales. In particular, the MD simulations of polymer melts at two different monomer densities *ρ* = 0.4 and 0.7 are performed. Figure 2 plots the VACF and FACF of COM at two different monomer densities. It is clearly shown that, at lower monomer density *ρ* = 0.4, the correlation time scales of VACF and FACF are *τ _{v}* = 10.30 and

*τ*= 0.44, respectively. Here, the magnitude of

_{f}*τ*and

_{v}*τ*takes the time when normalized VACF and FACF decrease to 10

_{f}^{−2}. In this case, the ratio of the time scales defined as

*κ*=

*τ*/

_{v}*τ*= 23.41 indicates a clearly separated time scale of momentum and force, in which the memory effects are less important and a Markovian behavior can be expected. However, at higher monomer density

_{f}*ρ*= 0.7, the time scales of

*τ*= 0.61 and

_{v}*τ*= 0.22 result in

_{f}*κ*=

*τ*/

_{v}*τ*= 2.77, which means that the correlation time scales of VACF and FACF become comparable and the Markovian approximation becomes inaccurate.

_{f}Moreover, the molecules in the two MD systems follow different diffusion processes. Figure 3 shows the mean squared displacement (MSD) of the COM of star polymer at monomer densities *ρ* = 0.4 and *ρ* = 0.7. At short time scales, the molecules in both systems experience the ballistic regime where the MSD approaches (3*k _{B}T*/

*M*)

*t*

^{2}. Then, the molecules in the MD system of

*ρ*= 0.4 have normal diffusion in which the MSD is proportional to the length of time, i.e.,

*MSD*∼

*t*. However, in the MD system of

*ρ*= 0.7, the star polymers undergo a sub-diffusion process distinguished by

*MSD*∼

*t*

^{0.57}before a normal diffusion is reached at long time scales. With these two significantly different cases, we will practically coarsen the MD systems of polymer melts into corresponding CG systems to demonstrate how to evaluate the memory kernel from atomistic simulations and how to incorporate the memory effects in CG simulations correctly.

### B. Coarse-grained force field

Consider a pair of molecules *I* and *J* in MD systems and let **e**_{IJ} be the radial unit vector directed along center-to-center from *J* to *I*. In general, the total force **F**_{IJ} between the two molecules *I* and *J* is not parallel to the radial vector **e**_{IJ} because the molecule of star polymer consists of many discrete monomers. Here, we define three directions as shown in Fig. 4 for description of pairwise interactions between clusters *I* and *J*, where the symbol “∥” denotes the direction parallel to **e**_{IJ}, “⊥_{1}” for the direction along the perpendicular velocity component $VIJ\u22a51=VIJ\u2212(VIJ\u22c5eIJ)eIJ$ and “⊥_{2}” for the direction orthogonal to both **e**_{IJ} and **V**_{IJ}.

As shown in Eq. (2), the evolution of a CG particle is governed by three types of interaction. First of all, we need to sample the intermolecular mean force from atomistic simulations by averaging the total forces between different molecules of star polymer. Since the instantaneous pairwise force **F**_{IJ} has no preference in the plane orthogonal to **e**_{IJ}, the average pairwise force 〈**F**_{IJ}〉, which is taken as the conservative force, has only radial component^{23}

where **e**_{IJ} = (**R**_{I} − **R**_{J})/*R _{IJ}* with

*R*= |

_{IJ}**R**

_{I}−

**R**

_{J}|. The distance-dependent function $FIJC(RIJ)$ represents the magnitude of the conservative force $FIJC$.

To compute the value of $FIJC$ at various distances from MD trajectories, the distance between two molecules is divided into many bins with width of Δ. Then, the magnitude of $FIJC(RIJ)$ is obtained by averaging $FIJ\u22c5eIJ$ over all those pairs with a center-to-center distance between *R _{IJ}* − Δ/2 and

*R*+ Δ/2. Figure 5 shows the distance dependence of the conservative force $FIJC(RIJ)$ at the monomer densities

_{IJ}*ρ*= 0.4 and 0.7. It indicates that the conservative force changes slightly when the monomer density increases from

*ρ*= 0.4 to

*ρ*= 0.7. For easier numerical implementation of CG simulations, the raw data of $FIJC(RIJ)$ obtained directly from MD simulations are fitted by a bell-shaped function, i.e., $FIJC(RIJ)=759.69(1+4RIJ/3.32)(1\u2212RIJ/3.32)4$ for

*ρ*= 0.4 and $FIJC(RIJ)=946.05(1+4RIJ/3.28)(1\u2212RIJ/3.28)4$ for

*ρ*= 0.7. Here, the cutoff radius beyond which the pairwise interactions vanish is set to be

*R*= 3.32 for

_{cut}*ρ*= 0.4 and

*R*= 3.28 for

_{cut}*ρ*= 0.7. The inset of Fig. 5 provides a global view of these fitting functions, which will be used as the conservative force for DPD simulations in Section IV.

The non-conservative forces between two CG particles consist of a dissipative force $FIJD=\u222b0tKIJ(t\u2212s)VIJ(s)ds$ and a random force $\delta FIJQ(t)$, as shown in Eq. (2). For a particle-based system in thermal equilibrium, the random force $\delta FIJQ(t)$ is related to the dissipative force by the second FDT given by Eq. (4). To obtain CG force field that governs the motion of CG particles, the memory kernel **K**(*t*) should be evaluated from MD trajectories. In practice, the memory kernel **K**(*t*) can be either approximated by the temporal autocorrelation function of fluctuating force $CIJ(t)=\beta \u3008\delta FIJ(t)\delta FIJT(0)\u3009$, or obtained by solving Eq. (3) with computed correlations $\u3008\delta FIJ(t)VIJT(0)\u3009$ and $\u3008VIJ(t)VIJT(0)\u3009$.

In general, the fluctuating force *δ***F**_{IJ} between clusters *I* and *J* is three-dimensional as shown in Fig. 4. However, *δ***F**_{IJ} has no preference in the plane orthogonal to **e**_{IJ} because of the spherical symmetry of star polymers. Therefore, we ignore the difference between the perpendicular directions ⊥_{1} and ⊥_{2}, and decompose the fluctuating force *δ***F**_{IJ} into two orthogonal parts,

where $\delta FIJ\u2225$ is the radial component along vector **e**_{IJ} and $\delta FIJ\u22a5$ is the perpendicular component whose modulus is equally distributed on directions ⊥_{1} and ⊥_{2}. Since the two components $\delta FIJ\u2225$ and $\delta FIJ\u22a5$ are orthogonal to each other, the temporal correlations are approximated by the sum of two orthogonal terms,

where the fluctuating forces $\delta FIJ\u2225(t)$ and $\delta FIJ\u22a5(t)$ as well as their temporal correlations can be evaluated directly from the MD simulations. Similarly, the memory kernel **K**_{IJ}(*t*) is also approximated by $KIJ(t)\u2248KIJ\u2225(t)eIJeIJT+KIJ\u22a5(t)(1\u2212eIJeIJT)$. Since the directions ⊥_{1} and ⊥_{2} have equivalent contribution to $CIJ\u22a5(t)$ and $KIJ\u22a5(t)$, we have $CIJ\u22a51(t)=CIJ\u22a5(t)/2$ and $KIJ\u22a51(t)=KIJ\u22a5(t)/2$.

Figure 6 illustrates the temporal correlations of fluctuating force $CIJ\u2225(t)$ and $CIJ\u22a5(t)$, and also the computed memory kernels $KIJ\u2225(t)$ and $KIJ\u22a5(t)$ for polymer melts at a monomer density of *ρ* = 0.4. These correlations and memory kernels for the MD system of *ρ* = 0.7 are presented in Fig. 7. It should be noted that the force-force autocorrelation function **C**_{IJ}(*t*) and the memory kernel **K**_{IJ}(*t*) are functions of both time and intermolecular distance. To compute the pairwise correlations and memory kernels at different distances, the intermolecular distance *R _{IJ}* is divided into many bins with a width of Δ, which is similar to the computation of the conservative force $FIJC(RIJ)$. Then, the values of

**C**

_{IJ}(

*R*,

*t*) and

**K**

_{IJ}(

*R*,

*t*) are obtained by computing the force-force correlations and cross correlations appear in Eq. (3) over the pairs with a distance

*R*− Δ/2 <

*R*(

_{IJ}*t*) <

*R*+ Δ/2. For instance, $CIJ\u2225(R,t)$ is computed through

We note that the bin size Δ should be selected properly. Because of the diffusion of molecules, the intermolecular distance *R _{IJ}* changes continuously. On the one hand, using small bin size means that only short time correlations can be obtained, and on the other hand, large bin sizes will reduce the spatial resolution of

**C**

_{IJ}(

*R*,

*t*) and

**K**

_{IJ}(

*R*,

*t*). To balance the time length of correlations and their spacial resolution, a bin size Δ = 0.05 is used in the present study.

## IV. INCORPORATING MEMORY EFFECTS IN DPD SIMULATIONS

Once the conservative force and the memory kernel of dissipation are obtained from atomistic simulations, a CG system can be constructed that takes into account such a force field. In the present work, DPD is employed as such a CG system. More specifically, DPD is a particle-based mesoscopic method in which the particles interact through conservative, dissipative, and random forces.^{14} Since the interactions between DPD particles are pairwise, the total momentum of the system is strictly conserved. Moreover, Español^{33} and Marsh *et al.*^{34} reported that the hydrodynamic equations of a DPD system recover the continuity and Navier-Stokes equations, which means that DPD can provide the correct hydrodynamic behavior of fluids at mesoscale. Therefore, DPD and its extensions^{35–37} have been successfully applied to simulations of biological systems^{37,38} and complex fluids.^{39} However, these DPD models cannot consider the correlations of fluctuating force. In this section, we first apply the Markovian approximation to simplify the dissipative term in Eq. (2), which results in an EOM consistent with the formulations of the conventional DPD model. Subsequently, we discuss how to incorporate the memory effects in a NM-DPD model when the Markovian assumption becomes inaccurate.

### A. Markovian DPD model

Usually, if the typical time scales of random force and momentum are well separated, the Markovian assumption can be applied, where the time correlation of random force is approximated by the Dirac delta function. Then, we have

where the *γ*_{IJ} is the friction tensor. In practice, *γ*_{IJ} is decomposed into two orthogonal terms, which are computed by

*τ*

_{0}and

*τ*

_{1}are the cutoff times for the integral of

*K*(

*t*). More specifically,

*τ*

_{0}= 0.45 and

*τ*

_{1}= 0.25 are used for the cases

*ρ*= 0.4 and 0.7. We note here that the time step

*δt*= 0.001 is used for MD simulations and Δ

*t*= 0.005 for DPD and NM-DPD simulations.

The radial and perpendicular components of the computed friction coefficients $\gamma IJ\u2225(RIJ)$ and $\gamma IJ\u22a51(RIJ)$ versus the distance *R _{IJ}* are presented in Fig. 8. Their values obtained from the integral of the memory kernel

**K**(

*t*) and the force correlation function

**C**(

*t*) are quite similar. It is worth noting that there is no symbol for

*R*< 2.2 because the radial distribution function (RDF) is zero when

_{IJ}*R*< 2.2, as shown in Fig. 9. To simplify the numerical implementation of DPD simulations, the symbols in Fig. 8 are fitted by a function

_{IJ}*f*(

*R*) = Λ(1 +

*χR*/

*R*)(1 −

_{cut}*R*/

*R*)

_{cut}^{χ}with parameters listed in Table I.

ρ
. | Based on . | γ
. | Λ . | χ
. | R
. _{cut} |
---|---|---|---|---|---|

0.4 | C(t) | γ^{∥}(R) | 146.18 | 3.00 | 3.32 |

γ^{⊥1}(R) | 110.76 | 3.95 | 3.32 | ||

K(t) | γ^{∥}(R) | 84.13 | 2.02 | 3.13 | |

γ^{⊥1}(R) | 112.85 | 3.51 | 3.13 | ||

0.7 | C(t) | γ^{∥}(R) | 112.21 | 2.17 | 3.11 |

γ^{⊥1}(R) | 102.30 | 3.27 | 3.11 | ||

K(t) | γ^{∥}(R) | 82.07 | 1.78 | 3.10 | |

γ^{⊥1}(R) | 42.28 | 2.62 | 3.09 |

ρ
. | Based on . | γ
. | Λ . | χ
. | R
. _{cut} |
---|---|---|---|---|---|

0.4 | C(t) | γ^{∥}(R) | 146.18 | 3.00 | 3.32 |

γ^{⊥1}(R) | 110.76 | 3.95 | 3.32 | ||

K(t) | γ^{∥}(R) | 84.13 | 2.02 | 3.13 | |

γ^{⊥1}(R) | 112.85 | 3.51 | 3.13 | ||

0.7 | C(t) | γ^{∥}(R) | 112.21 | 2.17 | 3.11 |

γ^{⊥1}(R) | 102.30 | 3.27 | 3.11 | ||

K(t) | γ^{∥}(R) | 82.07 | 1.78 | 3.10 | |

γ^{⊥1}(R) | 42.28 | 2.62 | 3.09 |

Since the pairwise interactions between DPD particles have perpendicular components, the rotational motions of DPD particles should be included to conserve the angular momentum of the DPD system.^{23} Therefore, we consider the pairwise interactions in all the three directions **e**^{∥}, **e**^{⊥1}, and **e**^{⊥2} shown in Fig. 4, and also include the rotational motion of DPD particles. As a result, the time evolution of a DPD particle is given by^{23,40,41}

where **Ω**_{I} is the angular velocity of particle *I*, **T**_{I} is the torque and **L**_{I} = *I*_{RI}**Ω**_{I} the angular momentum. The magnitude of the rotational inertia of a star polymer is 6.55, which is obtained by a computation of $IRx=\u3008\u2211mi(y\u02c6i2+z\u02c6i2)\u3009=23\u3008MRg2\u3009$. Here, *R _{g}* is the gyration radius of star polymers.

^{23}Also,

*ξ*is a Gaussian white noise that is independent for different pairs and at different times.

_{IJ}^{42}

^{,}$dWIJA=12(dWIJ\mu \nu \u2212dWIJ\nu \mu )$ is an antisymmetric noise matrix. The relationships [

*σ*

^{∥}(

*R*)]

_{IJ}^{2}= 2

*γ*

^{∥}(

*R*)

_{IJ}*k*and [

_{B}T*σ*

^{⊥1}(

*R*)]

_{IJ}^{2}= 2

*γ*

^{⊥1}(

*R*)

_{IJ}*k*are required so that the fluctuation-dissipation theorem is satisfied and the DPD system can be maintained at constant temperature.

_{B}T^{42}

### B. Non-Markovian DPD model

Without the Markovian approximation, we have to preserve the temporal memory of each pair in the NM-DPD model. To avoid prohibitively computational cost, we assume that the memory on time is finite, e.g., history length *N*Δ*t* where Δ*t* is the time step of NM-DPD simulations. Therefore, the time correlations of the fluctuating and random forces are zero when the time interval is larger than *N*Δ*t*,

Let the friction coefficients $\Gamma IJ,n\u2225(R)=KIJ\u2225(R,n\Delta t)$ and $\Gamma IJ,n\u22a51(R)=0.5KIJ\u22a5(R,n\Delta t)$, where the memory kernels $KIJ\u2225(R,n\Delta t)$ and $KIJ\u22a5(R,n\Delta t)$ should be replaced by the temporal correlations of fluctuating force $CIJ\u2225(R,n\Delta t)$ and $CIJ\u22a5(R,n\Delta t)$ when **C**(*t*) is taken as an approximation of **K**(*t*). Then, the EOM of NM-DPD particles is given by

where $VIJ\u2225(t\u2212n\Delta t)=[VIJ(t\u2212n\Delta t)\u22c5eIJ(t\u2212n\Delta t)]eIJ(t\u2212n\Delta t)$ is the parallel velocity component and $VIJ\u22a5(t\u2212n\Delta t)=VIJ(t\u2212n\Delta t)\u2212VIJ\u2225(t\u2212n\Delta t)$ the perpendicular velocity component. The coefficients $\alpha IJ,n\u2225$ and $\alpha IJ,n\u22a51$ are used to generate colored noise satisfying the second FDT given by Eq. (4). These coefficients are related to the friction kernels by

*η*

^{∥}be the integral of

*θ*

^{∥}(

*t*), that is, $\eta \u2225=\u222b0\tau 0\theta \u2225(t)dt$. Then, the distance-dependent function

*ϕ*

^{∥}(

*R*) =

_{IJ}*γ*

^{∥}(

*R*)/

*η*

^{∥}where

*γ*

^{∥}(

*R*) is the function listed in Table I. The details of colored noise generator are described in Appendix B, in which two different schemes for obtaining the coefficients $\alpha IJ,n\u2225$ and $\alpha IJ,n\u22a51$ are tested and their performances are compared.

With the constructed CG force field, DPD and NM-DPD simulations are performed with 1000 DPD particles in periodic cubic boxes of length *L* = (1000*N _{c}*/

*ρ*)

^{1/3}, where

*N*= 11 and

_{c}*ρ*= 0.4 or 0.7. Here, each DPD particle corresponds to a molecule of star polymer in MD systems. The velocity-Verlet algorithm

^{10,14}is employed for the numerical integration with a time step of Δ

*t*= 0.005. By monitoring the temperature in the simulations, it has been verified that the temperatures of DPD and NM-DPD systems are maintained at

*k*= 1.0, which indicates correct implementations of FDT and second FDT.

_{B}TThe comparisons on the RDF between Markovian-based DPD systems, NM-DPD systems, and the MD systems at different monomer densities are made in Fig. 9. In particle-based systems, RDF is the primary linkage between macroscopic thermodynamic properties, such as pressure and compressibility, and microscopic interactions. Figure 9 shows that the CG force fields obtained from MD trajectories reproduce excellent RDFs consistent with the original MD systems at both *ρ* = 0.4 and *ρ* = 0.7. Including non-Markovian memory in the DPD system does not change its performance on reproducing RDF. The reason is that the RDF of a DPD system is only determined by the conservative force, and the changes of non-conservative force do not affect RDF, even if the non-Markovian memory is introduced in the simulations. Moreover, the pressure is not affected by the inclusion of non-Markovian memory either, for example, the MD system of *ρ* = 0.4 has a pressure of *P* = 0.191 while both DPD and NM-DPD give *P* = 0.193.

Now, we examine the dynamics properties of MD, DPD, and NM-DPD systems by comparing their VACFs. The behavior of VACF implicitly represents the dynamical properties of the system. For example, a typical VACF shows an exponential decay determined by the viscosity, while the diffusion constant can be obtained by the integral of VACF via Green-Kubo relations. To examine the VACF of a CG system, comparing with its underlying MD system will quantify how well the dynamics of the MD system is reproduced by the CG model.

For the MD system of star polymer melts at a monomer density *ρ* = 0.4, the typical time scale of momentum is *τ _{v}* = 10.30 while the time scale of fluctuating force is

*τ*= 0.44, as shown in Fig. 2. The ratio of the time scales

_{f}*κ*=

*τ*/

_{v}*τ*= 23.41 indicates an apparent separation of time scales, and hence the Markovian assumption should be valid for this MD system. Figure 10 shows the performance of DPD and NM-DPD models in reproducing the MD system on VACF, in which Fig. 10(a) uses the force-force autocorrelation function

_{f}**C**

_{IJ}(

*t*) for computation of the dissipative coefficients while Fig. 10(b) uses the memory kernel

**K**

_{IJ}(

*t*). Each figure has a zoom-in view near

*t*= 0 as an inset. Results indicate that the VACF of MD system is well reproduced by the Markovian DPD model. Incorporation of the memory effects into DPD simulations yields small improvement on VACF at short time, but the improvement is not obvious because the deviation of Markovian DPD from MD is small.

It is worth noting that the VACFs of MD system start with a zero initial slope. However, applying the Markovian assumption completely ignores the temporal correlations of random forces, which results in an exponential decay of VACF at short time scales. Therefore, the VACFs of Markovian DPD model in all the insets of Figs. 10 and 11 have non-zero slopes at *t* = 0. The failure to capture short time VACF indicates that the dynamics of the system involved in high-frequency processes cannot be correctly produced. However, by including the non-Markovian memory in the pairwise interactions, NM-DPD generates much better short time VACF consistent with the VACF of MD system. Figure 11 shows that this improvement on short time VACF is obvious in the case of *ρ* = 0.7, where the ratio of time scale *κ* = *τ _{v}*/

*τ*= 2.77 implies a small separation of time scales. It is found that using

_{f}**C**

_{IJ}(

*t*) or

**K**

_{IJ}(

*t*) results in little difference on VACF, which indicates that the force-force autocorrelation function

**C**

_{IJ}(

*t*) can be a good approximation of the memory kernel

**K**

_{IJ}(

*t*) in pairwise systems.

## V. SUMMARY AND DISCUSSION

We applied the MZ projection operator to a microscopic system to construct corresponding CG systems. In particular, the CG dynamics resulted from a MZ projection of an underlying microscopic dynamics is governed by three terms. The conservative term results from the change of potential energy. The random term is generated by unresolved degrees of freedom. The dissipative term depends on history variables and an integral of the memory kernel of the random term. In general, including the memory effects is complicated and the Markovian assumption is often employed to simplify the formulation of CG models and corresponding implementation. However, when the time scales of a system are not clearly separated, the Markovian assumption becomes inaccurate and the memory effects should be taken into account. To this end, we studied two significantly different systems to demonstrate how to evaluate the memory kernel directly from atomistic simulations and how to incorporate the memory effects in CG simulations when the Markovian assumption fails.

For practical implementation, MD simulations of star polymer melts were performed to provide multiscale dynamics to be coarse-grained. By grouping many bonded atoms of the MD system into single CG particles, both the conservative and non-conservative interactions were evaluated from MD trajectories. For the MD system having clear time scale separation, we used the Markovian approximation to obtain a DPD model. Otherwise, we preserved the temporal memory of pairwise interactions in a NM-DPD model, where colored noises are generated by two different schemes to satisfy the second FDT. Results show that including memory effects changes the dynamic properties of DPD systems but does not affect static properties, i.e., radial distribution function and pressure. Also, it was demonstrated that the Markovian assumption works well for the system with clear time scale separation, where NM-DPD model has little improvement on the VACF compared with Markovian DPD model. However, when the time scales of a system are not fully separated, the NM-DPD can reproduce correct short-time properties that are related to how the system responds to high-frequency disturbances, which cannot be captured by the Markovian-based DPD model.

The present work provides a clear path to constructing a CG force field for DPD simulations from atomistic simulations. On the one hand, we demonstrated that the Markovian assumption can be safely applied to simplify the computation of dissipative force when the time scales of microscopic dynamics are apparently separated. On the other hand, for the systems without clear separation of time scales, we suggested a methodology to compute the memory kernel directly from microscopic dynamics. Moreover, we presented the procedure for incorporation of the memory effects in the CG simulations with the correct implementation of the second FDT. We note that the present results are only a first step in the direction of including finite memory effects in coarse-grained modeling. In future work, we plan to apply the non-Markovian DPD model to macromolecules at a more aggressive coarse-graining level so that larger clusters may introduce long memory effects in their effective dynamics.

## Acknowledgments

This work was supported by the DOE Center on Mathematics for Mesoscopic Modeling of Materials (CM4) and the National Institutes of Health Grant No. U01HL114476. Computational resources were provided by the Innovative and Novel Computational Impact on Theory and Experiment (INCITE) program. Z. Li would like to acknowledge Professor Eric Darve and Dr. Heesun Lee for helpful discussions and Dr. Yu-Hang Tang for computational support.

### APPENDIX A: PAIRWISE DECOMPOSITION

In this appendix, we present the detailed derivations of the pairwise decomposition for obtaining Eq. (2) from Eq. (1). By assuming that the non-bonded interactions between clusters are pairwise decomposable and ignoring the many-body correlations between different pairs, the first term on the right-hand side (RHS) of Eq. (1) is approximated by

where $FIJC(RIJ)eIJ=\u3008FIJ\u3009$ is the mean force between two clusters *I* and *J* at a separation distance of *R _{IJ}*. Here,

**e**

_{IJ}is the unit vector from cluster

*J*to

*I*,

**F**

_{IJ}is the instantaneous force exerted by cluster

*J*on cluster

*I*, and 〈

**F**

_{IJ}〉 is the ensemble average of

**F**

_{IJ}. It is worth noting that the instantaneous pairwise force

**F**

_{IJ}has no preference in the plane orthogonal to

**e**

_{IJ}, therefore, the average pairwise force 〈

**F**

_{IJ}〉 has only non-zero component along the radial direction

**e**

_{IJ}.

Similarly, the third term on the RHS of Eq. (1) $\delta FIJQ(t)$ is also decomposed into pairwise random forces,

According to the Mori-Zwanzig formalism, the random force $\delta FIQ(t)$ is generated from the orthogonal dynamics, which depends on the full initial conditions of the microscopic system being projected. We note that although the random force $\delta FIJQ$ is difficult to be evaluated in a microscopic dynamics, a direct computation of the fluctuating force *δ***F**_{IJ} is relatively straightforward. Here, the fluctuating force *δ***F**_{IJ} is defined as the difference of the instantaneous force **F**_{IJ} from the mean force 〈**F**_{IJ}〉, i.e., *δ***F**_{IJ} = **F**_{IJ} − 〈**F**_{IJ}〉. Then, the force autocorrelation function $CIJ(t)=\beta \u3008\delta FIJ(t)\delta FIJT(0)\u3009$ can be computed from microscopic dynamics as an approximation of the memory kernel $KIJ(t)=\beta \u3008[\delta FIJQ(t)][\delta FIJQ(0)]T\u3009$.

Now, we consider the second term on the RHS of Eq. (1), which has memory and contains an integral over past values of variables. If the many-body correlations between different pairs are ignored, only the terms that contain the correlations of same pair can survive in the derivation. Therefore, we have

where **V**_{IJ} = **V**_{I} − **V**_{J} is the relative velocity of CG particle *I* to *J*, and $KIJ(t)=\beta \u3008[\delta FIJQ(t)][\delta FIJQ(0)]T\u3009$ is the pairwise memory kernel involved in Eq. (2).

### APPENDIX B: COLORED NOISE GENERATOR

The thermodynamic properties of a particle-based system highly depend on its temperature. For a particle in a thermally equilibrium environment, the dissipative and random forces exerted on the particle should satisfy the FDT.^{22} If the dissipative force is history-dependent, the equation of motion of the particle is given in the form of generalized Langevin equation,^{16}

where *θ*(*t*) is a memory kernel and **R** is the random force. Then, the second FDT^{22} should be satisfied so that a constant temperature is guaranteed,

where *k _{B}* is the Boltzmann constant and

*T*the temperature. Let

*K*≥ 0 be the dimension of the random force

**R**(

*t*), and so the memory kernel

**(**

*θ**t*) would be a

*K*×

*K*matrix. To simplify the notations, the following derivations are demonstrated for the one-dimensional case.

In the numerical implementation, the random forces are imposed at each time step. Considering discrete time steps *t _{n}* =

*n*⋅

*δt*where

*n*≥ 0, we define

*R*=

_{n}*R*(

*t*) and

_{n}*θ*=

_{n}*θ*(

*t*). As a result, the second FDT given by Eq. (B2) can be rewritten into its discrete form,

_{n} which should be satisfied for all the integers *m* and *n*. For finite number of steps, i.e., 0 ≤ *n* ≤ *N*, let us write the memory kernel *θ _{n}* as follows:

where *α _{s}* is a set of undetermined coefficients. Then, we apply the discrete Fourier transform (DFT) on the memory kernel,

Since $\theta \u02c6k=|\alpha \u02c6k|2$ is greater or equal to zero, by defining $\alpha \u02c6k=\theta \u02c6k1/2$ the coefficients *α _{s}* can be obtained by the inverse Fourier transform

Let *W _{n}* be a sequence of independently and identically distributed normal random variables, we generate the random forces by

Then, the resultant random forces are correlated and satisfy the second FDT, i.e., 〈*R _{n}R_{m}*〉 =

*k*

_{B}Tθ_{n−m}.

A memory kernel *θ*(*t*) = exp(−19.30*t*)cos(28.25*t*) for 0 ≤ *t* ≤ 0.39, as plotted in Fig. 12(a), is used to test the DFT-based colored noise generator. With a periodic extension, the period Π = 0.78 is used in DFT. Let the time step *δt* = 0.005, and so a set of 156 coefficients *α _{s}* can be obtained directly by an inverse Fourier transform described in Eq. (B6). These coefficients

*α*are shown in Fig. 12(b). As expected, the correlation of the colored noises and the memory kernel coincide in Fig. 12(a).

_{s}To further test the DFT-based colored noise generator in practical simulations, we use a simpler version of non-Markovian DPD model, in which only radial interactions between DPD particles are considered,

where *R*(*t*) is generated by Eq. (B7), and $Vij\u2225(t)=[Vij(t)eij(t)]eij(t)$ is radial component of the relative velocity **V**_{ij}(*t*).

In our test, a NM-DPD simulation is performed in a periodic cubic box of length *L* = 26.37 filled with 1000 DPD particles. The temperature is set to *k _{B}T* = 1.0 and the mass of each DPD particle is

*m*= 11. Also, $FijC(r)=795.69(1+4r/rc)(1\u2212r/rc)4$ and

*ϕ*(

*r*) = 3.66 × 10

^{4}(1 −

*r*/

*r*)

_{c}^{3.84}with

*r*= 3.32, and the time step

_{c}*δt*= 0.005. Figure 13(a) shows time evolution of temperature in the NM-DPD simulation and its inset presents a FFT analysis of the temperature. It is found that the temperature is maintained at a constant temperature of

*T*= 1.0. Therefore, the second FDT is correctly implemented in the NM-DPD simulation. However, the evolution of temperature contains a dominant frequency of

*f*= 1.2821 corresponding to an obvious periodicity with a period of Π = 0.78, which is the period introduced in DFT. Moreover, such a periodicity is also observed in the VACF as shown in Fig. 13(b).

We note the imposed periodicity by DFT-based scheme is known and hence the extra effects in temperature and VACF can be removed easily in postprocessing, but the analysis of artificial effects in particle trajectories would be complicated. Thus, although the DFT-based scheme can generate colored noises with correct correlations, the resultant random forces introduce artificial periodicity that significantly changes the dynamics of the NM-DPD system. To avoid this problem, we consider another approach without using DFT.

When the random force *R _{n}* is generated by

we can compute the correlation of the random forces by

where $fn=kBT\u2211s=0N\alpha s\alpha n+s$. It is known that the second FDT requires a correlation of random forces given by 〈*R _{n}R_{m}*〉 =

*k*

_{B}Tθ_{n−m}. To achieve this, numerical optimization techniques can be used to obtain a set of coefficients

*α*by minimization of difference between

_{s}*f*

_{n−m}and

*k*

_{B}Tθ_{n−m}. In practice, the globalized bounded Nelder-Mead algorithm

^{43}is used for the optimization. We note, however, that many other optimization methods can be also used to obtain

*α*by defining the test function $fn=kBT\u2211s=0N\alpha s\alpha n+s$ and the target function

_{s}*k*.

_{B}Tθ_{n}Figure 14(a) shows the target function *k _{B}Tθ*(

*t*) together with the realistic correlations of colored noise

*f*(

*t*), while the coefficients

*α*obtained by optimization-based scheme (OPT) are shown in Fig. 14(b). It can be found that the target function

_{s}*k*(

_{B}Tθ*t*) is obtained by the optimization scheme. Actually, the iterative strategies that take previous results as an initial guess can be applied to increase the accuracy of optimization further. Since the coefficients

*α*are chosen from random variables, there are no apparently dominant frequencies involved.

_{s}Next, the OPT-based colored noise generator is tested in a practical NM-DPD simulation. Here, we use the same DPD system of previous test for the DFT-based scheme. The evolution of temperature in the NM-DPD simulation is presented in Fig. 15(a). It can be observed that the temperature fluctuates slightly around *T* = 1.0, which indicates a correct implementation of the second FDT. Moreover, a smooth VACF is shown in Fig. 15(b). No detectable periodicity is introduced in the NM-DPD system when OPT-based scheme is used. Therefore, all the NM-DPD simulations in the present work adopt the OPT-based colored noise generator.