Incorporation of memory effects in coarse-grained modeling via the Mori-Zwanzig formalism

Atomistic simulation techniques allow precise representation of molecular structures and have become standard computational tools nowadays for studying molecular systems.¹ 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⁴-10⁸ 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.⁵ 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,⁶ the inverse Monte Carlo method,⁷ the force matching method,⁸ or the multiscale coarse-graining methodology.⁹ 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).¹⁰ 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.¹¹ 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.¹³ 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).¹⁴ 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.²¹ 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)²² 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.

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_BT with k_B the Boltzmann constant and T the thermodynamic temperature. Also, R = {R₁, R₂, …, 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 $δ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 $Ψ(t)=β[δFIQ(t)][δ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≈∑J≠IFIJ$ and $δFIQ≈∑J≠Iδ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)=β〈[δFIJQ(t)][δFIJQ(0)]T〉$ 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 $〈δFIJQ(t)VIJT(0)〉$ disappears because $δFIJQ(t)$ has no correlation with velocity. The reason is that the random force $δ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 $〈δFIJ(t)VIJT(0)〉$ and $〈VIJ(t)VIJT(0)〉$ 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 $δFQ(t)=e−QiLtδF$ where $Q$ is the orthogonal projection operator in the MZ procedure and L the Liouville operator,¹⁸ using zeroth-order approximation,²⁵ e.g., $e−QiLt≅e−iLt,$ we can compute the force-force correlation $CIJ(t)=β〈δFIJ(t)δFIJT(0)〉$ 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²² requires that the random force $δ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.

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.

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²⁷ 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²⁸ is adopted for the bond interaction between connected monomers,

where k = 30ε/σ² is the spring constant and r₀ = 1.5σ determines the maximum length of the spring.²⁹ Such a FENE spring is stiff and short enough so that bond crossing can be prevented in the MD system.³⁰ 

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 = (1000N_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_BT = 1.0 and the time step used in the velocity-Verlet integrator¹⁰ is δt = 1.0 × 10⁻³τ.

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 τ_f = 0.44, respectively. Here, the magnitude of τ_v and τ_f takes the time when normalized VACF and FACF decrease to 10⁻². In this case, the ratio of the time scales defined as κ = τ_v/τ_f = 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 ρ = 0.7, the time scales of τ_v = 0.61 and τ_f = 0.22 result in κ = τ_v/τ_f = 2.77, which means that the correlation time scales of VACF and FACF become comparable and the Markovian approximation becomes inaccurate.

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 (3k_BT/M) t². 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.

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, “⊥₁” for the direction along the perpendicular velocity component $VIJ⊥1=VIJ−(VIJ⋅eIJ)eIJ$ and “⊥₂” 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²³ 

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⋅eIJ$ over all those pairs with a center-to-center distance between R_IJ − Δ/2 and R_IJ + Δ/2. Figure 5 shows the distance dependence of the conservative force $FIJC(RIJ)$ at the monomer densities ρ = 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−RIJ/3.32)4$ for ρ = 0.4 and $FIJC(RIJ)=946.05(1+4RIJ/3.28)(1−RIJ/3.28)4$ for ρ = 0.7. Here, the cutoff radius beyond which the pairwise interactions vanish is set to be R_cut = 3.32 for ρ = 0.4 and R_cut = 3.28 for ρ = 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=∫0tKIJ(t−s)VIJ(s)ds$ and a random force $δFIJQ(t)$⁠, as shown in Eq. (2). For a particle-based system in thermal equilibrium, the random force $δ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)=β〈δFIJ(t)δFIJT(0)〉$⁠, or obtained by solving Eq. (3) with computed correlations $〈δFIJ(t)VIJT(0)〉$ and $〈VIJ(t)VIJT(0)〉$⁠.

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 ⊥₁ and ⊥₂, and decompose the fluctuating force δF_IJ into two orthogonal parts,

where $δFIJ∥$ is the radial component along vector e_IJ and $δFIJ⊥$ is the perpendicular component whose modulus is equally distributed on directions ⊥₁ and ⊥₂. Since the two components $δFIJ∥$ and $δFIJ⊥$ are orthogonal to each other, the temporal correlations are approximated by the sum of two orthogonal terms,

where the fluctuating forces $δFIJ∥(t)$ and $δFIJ⊥(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)≈KIJ∥(t)eIJeIJT+KIJ⊥(t)(1−eIJeIJT)$⁠. Since the directions ⊥₁ and ⊥₂ have equivalent contribution to $CIJ⊥(t)$ and $KIJ⊥(t)$⁠, we have $CIJ⊥1(t)=CIJ⊥(t)/2$ and $KIJ⊥1(t)=KIJ⊥(t)/2$⁠.

Figure 6 illustrates the temporal correlations of fluctuating force $CIJ∥(t)$ and $CIJ⊥(t)$⁠, and also the computed memory kernels $KIJ∥(t)$ and $KIJ⊥(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∥(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.

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.¹⁴ Since the interactions between DPD particles are pairwise, the total momentum of the system is strictly conserved. Moreover, Español³³ and Marsh et al.³⁴ 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.³⁹ 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.

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

The radial and perpendicular components of the computed friction coefficients $γIJ∥(RIJ)$ and $γIJ⊥1(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_IJ < 2.2 because the radial distribution function (RDF) is zero when R_IJ < 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 f(R) = Λ(1 + χR/R_cut)(1 − R/R_cut)^χ with parameters listed in Table I.

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.²³ Therefore, we consider the pairwise interactions in all the three directions e^∥, e^⊥₁, and e^⊥₂ 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=〈∑mi(yˆi2+zˆi2)〉=23〈MRg2〉$⁠. Here, R_g is the gyration radius of star polymers.²³ Also, ξ_IJ is a Gaussian white noise that is independent for different pairs and at different times.⁴² $dWIJA=12(dWIJμν−dWIJνμ)$ is an antisymmetric noise matrix. The relationships [σ^∥(R_IJ)]² = 2γ^∥(R_IJ) k_BT and [σ^⊥₁(R_IJ)]² = 2γ^⊥₁(R_IJ) k_BT are required so that the fluctuation-dissipation theorem is satisfied and the DPD system can be maintained at constant temperature.⁴² 

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 $ΓIJ,n∥(R)=KIJ∥(R,nΔt)$ and $ΓIJ,n⊥1(R)=0.5KIJ⊥(R,nΔt)$⁠, where the memory kernels $KIJ∥(R,nΔt)$ and $KIJ⊥(R,nΔt)$ should be replaced by the temporal correlations of fluctuating force $CIJ∥(R,nΔt)$ and $CIJ⊥(R,nΔt)$ when C(t) is taken as an approximation of K(t). Then, the EOM of NM-DPD particles is given by

where $VIJ∥(t−nΔt)=[VIJ(t−nΔt)⋅eIJ(t−nΔt)]eIJ(t−nΔt)$ is the parallel velocity component and $VIJ⊥(t−nΔt)=VIJ(t−nΔt)−VIJ∥(t−nΔt)$ the perpendicular velocity component. The coefficients $αIJ,n∥$ and $αIJ,n⊥1$ are used to generate colored noise satisfying the second FDT given by Eq. (4). These coefficients are related to the friction kernels by

With the constructed CG force field, DPD and NM-DPD simulations are performed with 1000 DPD particles in periodic cubic boxes of length L = (1000N_c/ρ)^1/3, where N_c = 11 and ρ = 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_BT = 1.0, which indicates correct implementations of FDT and second FDT.

The 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 τ_f = 0.44, as shown in Fig. 2. The ratio of the time scales κ = τ_v/τ_f = 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 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/τ_f = 2.77 implies a small separation of time scales. It is found that using 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.

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.

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.

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=〈FIJ〉$ 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) $δFIJQ(t)$ is also decomposed into pairwise random forces,

According to the Mori-Zwanzig formalism, the random force $δ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 $δ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)=β〈δFIJ(t)δFIJT(0)〉$ can be computed from microscopic dynamics as an approximation of the memory kernel $KIJ(t)=β〈[δFIJQ(t)][δFIJQ(0)]T〉$⁠.

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)=β〈[δFIJQ(t)][δFIJQ(0)]T〉$ is the pairwise memory kernel involved in Eq. (2).

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.²² If the dissipative force is history-dependent, the equation of motion of the particle is given in the form of generalized Langevin equation,¹⁶ 

where θ(t) is a memory kernel and R is the random force. Then, the second FDT²² 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_n) and θ_n = θ(t_n). As a result, the second FDT given by Eq. (B2) can be rewritten into its discrete form,

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 $θˆk=|αˆk|2$ is greater or equal to zero, by defining $αˆk=θˆk1/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_nR_m〉 = k_BTθ_n−m.

A memory kernel θ(t) = exp(−19.30t)cos(28.25t) 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 α_s are shown in Fig. 12(b). As expected, the correlation of the colored noises and the memory kernel coincide in Fig. 12(a).

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∥(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_BT = 1.0 and the mass of each DPD particle is m = 11. Also, $FijC(r)=795.69(1+4r/rc)(1−r/rc)4$ and ϕ(r) = 3.66 × 10⁴(1 − r/r_c)^3.84 with r_c = 3.32, and the time step δ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).