Memory effects are often introduced during coarse-graining of a complex dynamical system. In particular, a generalized Langevin equation (GLE) for the coarse-grained (CG) system arises in the context of Mori–Zwanzig formalism. Upon a pairwise decomposition, GLE can be reformulated into its pairwise version, i.e., non-Markovian dissipative particle dynamics (DPD). GLE models the dynamics of a single coarse particle, while DPD considers the dynamics of many interacting CG particles, with both CG systems governed by non-Markovian interactions. We compare two different methods for the practical implementation of the non-Markovian interactions in GLE and DPD systems. More specifically, a direct evaluation of the non-Markovian (NM) terms is performed in LE-NM and DPD-NM models, which requires the storage of historical information that significantly increases computational complexity. Alternatively, we use a few auxiliary variables in LE-AUX and DPD-AUX models to replace the non-Markovian dynamics with a Markovian dynamics in a higher dimensional space, leading to a much reduced memory footprint and computational cost. In our numerical benchmarks, the GLE and non-Markovian DPD models are constructed from molecular dynamics (MD) simulations of star-polymer melts. Results show that a Markovian dynamics with auxiliary variables successfully generates equivalent non-Markovian dynamics consistent with the reference MD system, while maintaining a tractable computational cost. Also, transient subdiffusion of the star-polymers observed in the MD system can be reproduced by the coarse-grained models. The non-interacting particle models, LE-NM/AUX, are computationally much cheaper than the interacting particle models, DPD-NM/AUX. However, the pairwise models with momentum conservation are more appropriate for correctly reproducing the long-time hydrodynamics characterised by an algebraic decay in the velocity autocorrelation function.

## I. INTRODUCTION

There are many problems in biological systems and soft matter physics that occur on length and time scales beyond the capability of fully atomistic simulations.^{1–4} To this end, coarse-grained (CG) models have been developed to reduce the computational cost by eliminating some degrees of freedom (DOFs) from the all-atom system and creating a coarse system containing only the CG variables.^{5} Commonly, the effect of the unresolved DOFs on the CG variables is approximated by stochastic dynamics.^{6} With much fewer DOFs, the CG model can access much larger spatial and temporal scales, which makes CG modeling a rapidly expanding methodology, especially in the fields of soft matter research,^{7–9} biological cell modeling,^{10,11} macromolecular simulation,^{12–14} and multiscale computation.^{15,16}

In CG approaches, the conservative force defined as the negative gradient of CG potential is responsible for the static properties of the system, e.g., pressure, compressibility, and equilibrium structure often represented by the radial distribution function (RDF). An effective CG potential energy can be obtained numerically using many different coarse-graining methods. Examples include the iterative Boltzmann inversion,^{17} the reverse Monte Carlo,^{18} the force matching method,^{19} the adaptive biasing force method,^{20–23} the multiscale coarse-graining method,^{5,24} and direct ensemble averaging.^{25,26} However, the CG potential alone cannot produce the correct dynamic properties.^{27} More specifically, Izvekov and Voth^{28} found that the 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 representing the influence of the unresolved DOFs on the CG variables are not correctly reproduced.^{5}

Since an effective CG potential can be obtained by numerous methods,^{17–24} in the present work, we focus on reproducing accurately the CG dynamics with stochastic models. Theoretically, the elimination of DOFs from a microscopic system leads to non-conservative forces^{29} that include a random term and a dissipation term. The equation of motion (EOM) of the CG variables is often in the form of a non-Markovian, generalized Langevin equation (GLE), which can be derived from the Mori–Zwanzig (MZ) formalism.^{25,30,31}

By assuming that the non-bonded interactions between neighboring CG clusters in the microscopic system are pairwise decomposable,^{6,32} we can reformulate the GLE into its pairwise version, i.e., non-Markovian dissipative particle dynamics (DPD).^{6,32,33} DPD models then lead to stochastic forces for all pairwise interactions. In the present work, we will compare a GLE model of a single coarse particle with a DPD system consisting of interacting particles.

It is notable that whenever variables are eliminated from a Markovian system of equations, the resultant system is non-Markovian.^{29,34,35} Therefore, both CG systems, GLE and non-Markovian DPD, which are constructed by elimination of certain DOFs from microscopic systems, contain non-Markovian interactions, i.e., a colored random force and a friction force being a convolution of a memory kernel with the particle velocities.

A Markovian approximation is obtained by replacing the correlation of random force $F\u223c(t)$ with the Dirac delta function $\delta (t)$,^{25,26}

where $V$ represents the velocity and *γ* is a constant. This Markovian assumption dramatically simplifies the numerical computation. However, for many stochastic processes in nature, van Kampen^{36} noted that “Non-Markov is the rule, Markov is the exception.” In the present work, we will consider non-Markovian CG systems, where the time scale of unresolved dynamics is comparable with that of CG dynamics. In these cases, the correlation of random force cannot simply be replaced by the Dirac delta function, which would make the Markovian approximation invalid.

To numerically implement the non-Markovian interactions, we can either compute the time-convolution term directly and generate colored noise that satisfies the second fluctuation-dissipation theorem (FDT)^{37} or employ auxiliary variables to replace the non-Markovian dynamics with a Markovian process in a higher dimensional space.^{38} The mapping of a non-Markovian system to a Markovian system using auxiliary variables is approximate since the memory kernel needs to be fitted for this mapping. If the memory kernel of the non-Markovian dynamics can be well approximated by the exponentially decaying oscillatory functions,^{38} the mapped Markovian system would be very close to the original non-Markovian system.

In the following, the direct computation of the time-convolution term is marked by “NM” while the system with auxiliary variables is marked by “AUX.” Therefore, four different models, namely, LE-NM, LE-AUX, DPD-NM, and DPD-AUX will be constructed, where the LE-NM and LE-AUX models simulate the dynamics of non-interacting CG particles governed by GLE, while the DPD-NM and DPD-AUX models simulate the dynamics of interacting particles via stochastic pairwise interactions, as described in Table I. Recently, Davtyan *et al.*^{39} generated non-Markovian dynamics by adding a set of fictitious particles with harmonic interactions among themselves and a special coupling to the CG particles, which is analogous to the idea of the LE-AUX model presented in the present paper.

LE-NM . | LE-AUX . |
---|---|

Direct computation of the | Markovian system with auxiliary |

time-convolution representing | variables coupled to the |

dissipation. The random | momentum, generating |

force needs colored noise. | NM dynamics. |

DPD-NM | DPD-AUX |

Direct computation of the | Markovian system with auxiliary |

time-convolution representing | to the relative momentum |

the pairwise dissipative interaction. | between DPD particles. |

LE-NM . | LE-AUX . |
---|---|

Direct computation of the | Markovian system with auxiliary |

time-convolution representing | variables coupled to the |

dissipation. The random | momentum, generating |

force needs colored noise. | NM dynamics. |

DPD-NM | DPD-AUX |

Direct computation of the | Markovian system with auxiliary |

time-convolution representing | to the relative momentum |

the pairwise dissipative interaction. | between DPD particles. |

For our numerical benchmarks, we will perform molecular dynamics (MD) simulations of a star-polymer melt. The corresponding CG system is defined by grouping many bonded atoms into a single cluster. With this definition, the pairwise potential and the memory kernel of friction can be computed from the microscopic trajectories via the MZ formulation.^{6,32} Subsequently, four different CG models, i.e., LE-NM, LE-AUX, DPD-NM, and DPD-AUX, will be developed as the CG representation. We compare the CG dynamics resulting from these models to quantify their accuracy and computational cost. Our numerical benchmarks reveal that the use of auxiliary variables can reduce the computational cost of simulations considerably, compared to the direct NM model (a speedup of 19–34 was observed). Moreover, the DPD model reproduces correctly the algebraic tail of the velocity autocorrelation function (VACF), while the single particle GLE model is accurate at short times, but exhibits a faster decay (exponential with AUX) at long times compared to the MD reference behavior.

The remainder of this paper is organized as follows: Section IIbriefly introduces the theoretical background of coarse-graining and the MZ formalism, as well as the numerical implementation of the non-Markovian dynamics. In Section III, we describe in detail how to practically construct CG models and compute the memory kernel directly from atomistic simulations. Subsequently, we present a quantitative comparison between different CG models. Then, we end up with a brief summary and discussion in Section IV. Moreover, Appendix A provides sensitivity analysis of auxiliary dynamics to the fitting of memory kernels, and Appendix B describes in detail how to construct the coupling matrices for the auxiliary systems.

## II. THEORETICAL BACKGROUND

### A. Coarse-graining via the Mori–Zwanzig formalism

Consider an atomistic *n*-particle system whose microscopic state $\Gamma ={\mathbf{r}n,\mathbf{p}n}$ is identified with the coordinates **r** and momenta **p** of the atomic particles. The Hamiltonian is

where $H(\Gamma )$ defines the phase space trajectories of the system $\Gamma \u2261{\mathbf{r}i,\mathbf{p}i,i=1,\u2026,n}$. When the atomistic details are not of practical interest, by dividing the original *n*-particle system into *K* clusters with each cluster containing many atomic particles, the dynamics of the system can be described by the evolution of proper CG variables, i.e., the coordinate **R** and momentum **P** of the center-of-mass (COM) of the clusters defined by $\mathbf{R}I=MI\u22121\u2211mIi\mathbf{r}Ii$ and $\mathbf{P}I=\u2211\mathbf{p}Ii$, as shown in Fig. 1; $MI$ is the mass of the cluster. Unless otherwise stated, the variables of CG particles are represented with capital symbols, such as *M*, **R** and **P** being mass, position and momentum, respectively, while the corresponding lowercases *m*, **r** and **p** denote the quantities of atomic particles.

Taking those CG quantities as resolved variables, the equation of motion of the CG particles derived from the Mori–Zwanzig formalism is in the form of^{6,25,26,30,32,33}

where $\beta =1/kBT$ with $kB$ the Boltzmann constant and *T* the thermodynamic temperature, $\mathbf{R}={\mathbf{R}I,I=1,\u2026,K}$ is a phase point in the CG phase space, and $\omega (\mathbf{R})$ is defined as a normalized partition function of all the microscopic configurations at a phase point **R**.^{30} The conservative force acting on a CG particle *I* depends on all the COM coordinates **R** shown in the first term, while the dissipative force given by the second term depends on the momenta of all CG particles in the system. We note that the exact memory kernel resulting from the Mori-Zwanzig projection is a complex function of the positions and momenta of CG particles. Although the non-linear terms of memory may appear in theoretical derivations,^{40} it is difficult to compute these non-linear terms for realistic systems. For practical purposes, we adopt the formation of Eq. (3) where the memory kernel is considered independent of the particle’s momentum. Moreover, the terms in Eq. (3) refer to multi-body interactions. A direct computation of the multi-body forces and the memory kernel in Eq. (3) is very challenging, even for a simple system derived from a one-dimensional harmonic chain.^{41,42} In the present work, we will use two different CG strategies to simplify Eq. (3) for computational practicality.

In the first CG approach, we consider an isotropic system free of external fields and hence the total conservative force acting on a CG particle is zero on average. By choosing the momentum **P** of CG particles as the variable of interest and assuming that the random force on different CG particles has no correlation, i.e., $\u27e8[\delta \mathbf{F}IQ(t)][\delta \mathbf{F}JQ(t\u2032)]T\u27e9=\bm{\delta}IJ\mathbf{K}(t\u2212t\u2032)$, the generalised Langevin equation describing the evolution of **P** can be written as

where $\delta \mathbf{F}Q$ is the random force, and $\mathbf{K}(t)=\beta \u27e8[\delta \mathbf{F}Q(t)][\delta \mathbf{F}Q(0)]T\u27e9$ is the memory kernel whose definition ensures that the second FDT is satisfied.^{37}

In the second CG approach, we first perform a pairwise decomposition,^{32} i.e., $\mathbf{F}I\u2248\u2211J\u2260I\mathbf{F}IJ$, and assume negligible many-body correlations between different pairs, i.e., $\u27e8[\delta \mathbf{F}IJQ][\delta \mathbf{F}IKQ]T\u27e9J\u2260K\u22480$. Then, Eq. (3) can be reformulated into a pairwise form

where $\mathbf{F}IJ$ is the instantaneous force between CG particles *I* and *J*. Its ensemble average $\u27e8\mathbf{F}IJ\u27e9$ is taken as the conservative force in the CG model. Also, $\mathbf{V}IJ(t)=\mathbf{V}I(t)\u2212\mathbf{V}J(t)$ denotes the relative velocity between the CG particles *I* and *J*. The last term $\delta \mathbf{F}IJQ$ is the random force resulted from unresolved degrees of freedom. In general, $\delta \mathbf{F}IJQ$ is non-Gaussian and the corresponding dissipative force depends on an integral of history velocities weighted by a memory kernel, which is defined by $\mathbf{K}IJ(t)=\beta \u27e8[\delta \mathbf{F}IJQ(t)][\delta \mathbf{F}IJQ(0)]T\u27e9$.

### B. Direct computation of non-Markovian interactions

In a non-Markovian dynamics described by the GLE of Eq. (4), the dissipative force is history-dependent. To simplify the mathematical notations, we demonstrate the following derivation with the one-dimensional case:

where $K(t)$ is a memory kernel and $F\u223c$ is the random force. In the numerical implementation, both the dissipative force and the random force are imposed at each time step. Let $tn=n\Delta t$ where $n\u22650$ represent discrete time steps,^{32,33} also let $Kn=K(tn)$ and $F\u223cn=F\u223c(tn)$. Then, the second FDT can be rewritten into its discrete form,

To avoid prohibitive computational cost in practical implementation, we assume that the memory is finite, i.e., a history length $N\Delta t$, where $\Delta t$ is the time step in the simulations. Therefore, the time convolution representing the dissipative force can be computed by a summation over *N* discrete time steps,

Next, we need to generate the random force $F\u223cn$ satisfying the second FDT given by Eq. (7). Although the Fourier transform^{43} can be applied to generate colored noise satisfying Eq. (7), the noise represented by Fourier series will introduce an artificial periodicity into the system,^{44} which affects the particle dynamics significantly.^{32} Here, we adopt an optimization-based scheme proposed in our previous work^{32} for colored noise generation to avoid the artificial periodicity. More specifically, the random force $F\u223cn$ at discrete time steps is given by

where $\alpha s$ is undetermined coefficient and $Wn$ is a sequence of independent identically distributed Gaussian random variables with zero mean and unit variance. Then, we can compute the correlation of the random force $F\u223c$ by

where $fn=kBT\u2211s=0N\alpha s\alpha n+s$. The second FDT requires the correlation of random force related to the memory kernel by Eq. (7). To achieve this, a set of coefficients $\alpha s$ is obtained by performing numerical optimization for minimizing the difference between $fn\u2212m$ and $kBTKn\u2212m$. In practical implementation, the globalized bounded Nelder-Mead algorithm^{45} is applied for the optimization. We note, however, that there are many other optimization methods^{46} that can be used to obtain the coefficients $\alpha s$ by defining the test function $fn=kBT\u2211s=0N\alpha s\alpha n+s$ and the target function $kBTKn$.

With the discretization of the time convolution of dissipative force in form of Eq. (8) and the colored noise generator given by Eq. (9), the simulation of the generalized Langevin dynamics can be performed directly. It is worth noting that the direct computation of non-Markovian forces by Eqs. (8) and (9) requires the storage of historical information of particles, which could dramatically increase the computational cost and complexity.

### C. Auxiliary variables for non-Markovian dynamics

In order to bypass the complexity of computing the non-Markovian interactions directly, we supplement the system with auxiliary variables $\mathbf{s}={si}$, which are linearly coupled to the momentum and among themselves.^{38} The corresponding stochastic differential equation is given as

where $\mathbf{A}sp=\u2212\mathbf{A}psT$, *𝝃* is a vector of uncorrelated Gaussian random numbers with $\u27e8\xi i(t)\xi j(0)\u27e9=\delta ij\delta (t)$. Let $\mathbf{x}=(\mathbf{P},\mathbf{s})T$, the dynamics of **x** is the Ornstein-Uhlenbeck (OU) process $\mathbf{x}\u02d9=\mathbf{A}\mathbf{x}+\mathbf{B}\bm{\xi}$ with $\mathbf{A}=((\U0001d7ce,\mathbf{A}sp)T,(\mathbf{A}ps,\mathbf{A}ss)T)$ and $\mathbf{B}=diag(\U0001d7ce,\mathbf{B}s)$. Assuming that the non-Markovian dynamics has finite memory, one can safely take $\mathbf{s}(\u2212\u221e)=0$. By solving the equation of **s** in Eq. (11) we arrive at

Substituting Eq. (12) into Eq. (11), we eliminate the auxiliary variables **s** from the dynamics of momentum **P** and have^{38}

where the matrix $\mathbf{B}s$ is determined by $\mathbf{B}s\mathbf{B}sT=kBT(\mathbf{A}ss+\mathbf{A}ssT)$ to satisfy the fluctuation-dissipation theorem. In practical implementations, we can always find this $\mathbf{B}s$ by designing $\mathbf{A}ss$ such that $(\mathbf{A}ss+\mathbf{A}ssT)$ is positive-definite or positive-semidefinite and using the Cholesky factorization.^{47} When $\mathbf{A}ss$ is a real matrix and has complex eigenvalues with positive real part, the memory kernel $\mathbf{K}(t)$ is a linear combination of exponentially damped oscillators.^{38} As a result, we have Markovian dynamics in the form of Eq. (11), which is equivalent to the non-Markovian dynamics given by Eq. (13) for a certain form of $\mathbf{K}(t)$.

## III. COARSE-GRAINING PROCEDURE

### A. Reference microscopic system

A reference MD system of star-polymer melts is constructed to generate multiscale dynamics for coarse-graining. In practice, the star-polymers are represented as chains of monomers connected by short springs. We note that increasing the length of polymer chains enhances the many-body correlations between different pairs.^{6,48,49} In the present work, we focus on the implementation of non-Markovian interactions in CG simulations. The investigation of many-body correlations induced by long arms of star polymers is beyond the scope of this paper. Therefore, a polymer melt consisting of star polymers with short arms is employed, in which each star polymer has 10 arms with one monomer per arm, and hence the total number of monomer per molecule is $Nc=11$.

The excluded volume interactions between monomers are represented through a Lennard-Jones (LJ) potential truncated for repulsion only, which is also known as the Weeks–Chandler–Andersen (WCA) potential,^{50}

where *r* is the distance between two monomers, and the cut off distance $rc=21/6\sigma $ is chosen so that only the repulsive part of the LJ potential is considered. Also, *ϵ* and *σ* set the energy and length scales, respectively. The bond interactions between neighboring monomers are described by a finitely extensible nonlinear elastic (FENE) potential,^{51}

where *k* is the spring constant and $r0$ determines the maximum length of the spring. In particular, a parameter set of $k=30\u03f5/\sigma 2$ and $r0=1.5\sigma $ is adopted so that the FENE spring is stiff and short enough to minimize artificial bond crossings.^{52}

MD simulations of polymer melts are performed in the canonical ensemble (NVT) through the Nośe–Hoover thermostat.^{53} More specifically, 8000 molecules of star polymer are filled into a periodic cubic box of length $LB=50.095\sigma $ to model a polymer melt. The number density of monomers is $\rho =0.7\sigma \u22123$. Unless otherwise indicated, all the results of both MD and CG simulations are reported in the reduced LJ units, i.e., the length, mass, energy, and time units are set as $\sigma =1$, $m=1$, $\u03f5=1$, and $\tau =\sigma (m/\u03f5)1/2$. In addition, the temperature of the MD system is maintained at $kBT=1.0$ and the time step $\delta t=0.001\tau $ is used in the MD simulations.

To obtain non-Markovian dynamics for coarse-graining, the force and momentum in the MD system should have comparable correlation time scales so that the memory effect is important for the dynamics. Fig. 2(a) plots the velocity autocorrelation function (VACF) and force autocorrelation function (FACF) of the COM of star polymers. The correlation time scales of VACF and FACF are $\tau v=0.61$ and $\tau f=0.22$, respectively. Let $\kappa =\tau v/\tau f=2.77$ be the ratio of the time scales, then it is obvious that $\tau v$ and $\tau f$ are of the same order and there is no apparent separation of time scales between the force and the velocity. Thus, a non-Markovian dynamics of the polymer melt is expected for coarse-graining.

As shown in Fig. 2(a), at short time scales, the VACF decays exponentially, and then turns into negative values because of the backscattering effect,^{54} which refers to the phenomenon of a particle bouncing back and forth in a temporary cage formed by its neighboring particles. Since the diffusion coefficient is an integral of VACF, i.e., $D(t)=13\u222b0tVACF(\tau )d\tau $, the negative VACF leads to a decrease in the diffusion coefficient. Correspondingly, the mean squared displacement (MSD) of the COM of star polymers in Fig. 2(b) shows that the star polymers undergo a subdiffusion process distinguished by $MSD\u223ct0.57$ before a normal diffusion $MSD\u223ct$ is reached at long time scales. Such behaviors of the polymer melt shown in VACF and MSD will be reproduced by the CG models.

### B. Generalized Langevin dynamics

In the first CG approach introduced in Section II A, the EOM of a CG particle is described by the generalized Langevin equation, which is a linear integro-differential equation

where the conservative force disappears because an isotropic system free of external fields is considered and hence the total conservative force acting on a CG particle is zero on average. If the system is in thermal equilibrium, the random force $\delta \mathbf{F}(t)$ and the memory kernel $\mathbf{K}(t)$ should be related by the second FDT. Therefore, if we compute $\mathbf{K}(t)$, we can model $\delta \mathbf{F}(t)$ accordingly using the second FDT.

Multiplying the velocity $\mathbf{V}T(0)$ on both sides of Eq. (16) and taking the ensemble average yields

Let us now define the velocity autocorrelation function as $\mathbf{C}(t)=\u27e8[\mathbf{V}(t)][\mathbf{V}(0)]T\u27e9$. Since the random force $\delta \mathbf{F}(t)$ is statistically independent of the velocity $\mathbf{V}(0)$, we can rewrite Eq. (17) into a form of Volterra equation

where both $\mathbf{C}(t)$ and the force-velocity correlation function (FVCF) $M\mathbf{C}\u02d9(t)=\u27e8[\mathbf{F}(t)][\mathbf{V}(0)]T\u27e9$ can be evaluated from MD simulations. Here, $\mathbf{F}(t)$ is the total force, and $\mathbf{V}(t)$ is the velocity of the COM of star polymers. Then, the memory kernel $\mathbf{K}(t)$ can be obtained by solving the Volterra equation (18) with known $\mathbf{C}\u02d9(t)$ and $\mathbf{C}(t)$.

Fig. 3(a) shows the computed VACF $\mathbf{C}(t)$ and FVCF $\mathbf{C}\u02d9(t)$ obtained directly from MD simulations, while Fig. 3(b) displays the corresponding memory kernel $\mathbf{K}(t)$ and its inset gives a global view of $\mathbf{K}(t)$. In the model of LE-NM, we compute the non-Markovian interactions directly by approximating the convolution integral by a summation over many discrete time steps, as described in Section II B. The curve of VACF contains a long-time hydrodynamic tail characterized by an algebraic decay, which yields the corresponding algebraic decay in the memory kernel. Theoretically, the temporal length of such a memory kernel is infinite, but, in practical implementation of LE-NM simulations, we need to assume that the memory is finite. Otherwise, the computers will run out of physical memory quickly for long simulations.

In the present LE-NM simulation, we truncate the memory kernel at $t=10$. In this case, a CG particle carries its historical information for $N=201$ time steps at a time step of $\Delta t=0.05$. The dissipative force of the generalized Langevin dynamics is computed by Eq. (8), while the corresponding random force is generated by Eq. (9) with a set of 201 coefficients $\alpha n$, which are obtained by numerical optimization using the globalized bounded Nelder-Mead algorithm.^{45} The resultant VACF of the LE-NM system and the reference VACF of MD system are plotted in Fig. 4(a). Fig. 4(b) zooms in the initial part of VACF to examine the short time dynamics, while Fig. 4(c) compares the long-time tail of VACF in double-logarithmic coordinates.

We observe in Fig. 4(a) that the LE-NM model is able to generate correct non-Markovian dynamics consistent with the reference MD system for $t\u226410$. In particular, the LE-NM model successfully captures the zero initial slope of VACF, which cannot be observed in Langevin dynamics with the Markovian assumption that ignores the temporal correlations of random forces. The reason is that the VACF of a Langevin dynamics decays exponentially with a non-zero initial slope.^{55} Also, since the memory kernel in LE-NM/AUX is not the Dirac delta function, the LE-NM and LE-AUX models can have more complex dynamics than normal diffusion.^{56} Fig. 5 shows the MSD of CG particle in the LE-NM system with comparison to that of the COM of star polymer in the MD system. At short time scales, the molecules in MD system first experience a ballistic motion where the MSD approaches $(3kBT/M)t2$. Then, the star polymers undergo a subdiffusion process distinguished by MSD$\u223ct0.57$ before a normal diffusion (MSD$\u223ct$) is reached at long time scales. It is observed in Fig. 5 that all the three stages of diffusion process can be reproduced by the LE-NM model. However, since the memory kernel has been truncated at $t=10$, the curve of temporal correlation contains an artificial peak at $t=10$ as shown in Fig. 4(c), which also raises $D(t)$ at $t=10$ displayed in the inset of Fig. 4(a). Hence, we conclude that although the LE-NM model with truncated memory kernel can generate correct non-Markovian dynamics for short time ($t<10$), the long-time dynamics (longer than the truncate length of memory kernel) cannot be correctly reproduced because of the truncation.

Alternatively, we construct a LE-AUX model based on the description of Section II C to avoid memory truncation. Instead of computing the non-Markovian interaction directly, the LE-AUX model employs auxiliary variables to model the non-Markovian dynamics by a set of Ornstein–Uhlenbeck processes.^{38} In particular, the computed memory kernel shown in Fig. 3(b) is fitted by a linear combination of three exponentially damped oscillators, i.e., $f(t)=\u2211i=13\Lambda iexp(\u2212t/\tau i)cos(\omega it\u2212\phi i)$ with parameters $(\Lambda i,\tau i\u22121,\omega i,\phi i)i=1,2,3$ = ${(355.53$, 23.577, 29.954, $0.666\u2009817\u20096)$, $(1.935\xd7105$, 9.559, $0.006\u200942$, $1.570\u2009125)$, $(39\u2009681.12$, $0.248\u200948$, $7.298\xd710\u22125$, $1.570\u2009503)}$. The validity and accuracy of the extended dynamics depend on the fitting of the computed memory kernels. In Appendix A, we present the sensitivity analysis of the LE-AUX system to the fitting of memory kernels. Since the self-diffusion coefficient of a Brownian particle is determined by the integral of the memory kernel,^{55} i.e., $D=kBT/\xi $ where $\xi =\u222b0\u221eK(t)dt$, matching the integration value of the fitting function $f(t)$ to the integration value of $K(t)$ yields the correct diffusion coefficient. Based on the fitting function, the corresponding matrices in Eq. (11) are given by

where $\mathbf{B}s\mathbf{B}sT=kBT(\mathbf{A}ss+\mathbf{A}ssT)$ is a sufficient condition to satisfy the FDT for having canonical sampling. The details of construction of the coupling matrices **A** and $\mathbf{B}s$ from the fitting functions are described in Appendix B. As shown in Fig. 4(b), the initial value $VACF(0)=0.273=3kBT/M$ confirms that the FDT is satisfied in the LE-AUX simulations and the temperature of the LE-AUX system is maintained at $kBT=1.0$.

The memory kernel of a multivariate OU process is in the form of $\mathbf{K}(t)=\u2212\mathbf{A}psexp(\u2212t\mathbf{A}ss)\mathbf{A}sp$, which is in principle an arbitrary combination of complex exponentials. As shown in Fig. 3(b) and its inset, although the leading part of the computed memory kernel $K(t)$ can be approximated well by the fitting function $f(t)$, it is difficult to fit the long-time algebraic decay with a combination of exponential functions. As a result, the short time dynamics of the reference MD system can be reproduced by the LE-AUX model. In particular, the zero initial slope of VACF shown in Fig. 4(b) and the three stages of diffusion process shown in Fig. 5 are correctly captured by the LE-AUX model. However, since the memory kernel of the LE-AUX model is a combination of complex exponentials, the long-time tail of VACF shows an exponentially decay rather than an algebraic decay of the MD system. Therefore, it can be seen in Fig. 4(c) that the VACF of the LE-AUX system diverges from the algebraically decaying VACF of the MD system for $t>5$.

For problems where only short time dynamics are concerned, both LE-NM and LE-AUX can be applied to produce the underlying multiscale dynamics at a CG level. However, if the long-time hydrodynamics are also of importance and have to be correctly preserved, the non-interacting particle model, i.e., LE-NM and LE-AUX, is not accurate enough. The reason is that the momentum conservation is known to be essential for the appearance of long-time tail of VACF,^{54} but neither LE-NM nor LE-AUX conserves the momentum and hence the correct long-time tail of VACF cannot be expected. In the following, we will construct CG models based on pairwise interactions so that the momentum conservation is rigorously satisfied.

### C. Non-Markovian dissipative particle dynamics

In the second CG approach introduced in Section II A, with the assumption of negligible many-body correlation and using the pairwise decomposition, the EOM of a CG particle can be written as Eq. (5) in a pairwise form. Let us define a fluctuating force $\delta \mathbf{F}IJ(t)$ as the difference between the instantaneous force and the mean force, i.e., $\delta \mathbf{F}IJ(t)=\mathbf{F}IJ(t)\u2212\u27e8\mathbf{F}IJ\u27e9$. The second equality of Eq. (5) becomes

where $\mathbf{K}IJ(t)=\beta \u27e8[\delta \mathbf{F}IJQ(t)][\delta \mathbf{F}IJQ(0)]T\u27e9$. By multiplying by the velocity $\mathbf{V}IJT(0)$ on both sides of Eq. (20) and taking ensemble average, we have

where the third term $\u27e8\delta \mathbf{F}IJQ(t)\mathbf{V}IJT(0)\u27e9$ disappears because the random force $\delta \mathbf{F}IJQ(t)$ comes from orthogonal dynamics in MZ procedure^{29} and has no correlation with the CG velocity. Let us now define a function $\mathbf{A}IJ(t)=\u27e8\delta \mathbf{F}IJ(t)\mathbf{V}IJT(0)\u27e9$ and another function $\mathbf{B}IJ(t)=\u27e8\mathbf{V}IJ(t)\mathbf{V}IJT(0)\u27e9$. Then, Eq. (21) can be rewritten as

Considering that the fluctuating force $\delta \mathbf{F}IJ(t)$ and the relative velocity $\mathbf{V}IJ(t)$ have both radial components along vector $\mathbf{e}IJ$ and perpendicular components orthogonal to vector $\mathbf{e}IJ$, we decompose the fluctuating force as well as the relative velocity into two orthogonal parts,

where the symbol “∥” represents the radial component along $\mathbf{e}IJ$ and “⊥” for perpendicular component. Since the two directions of “∥” and “⊥” are orthogonal to each other, the temporal correlations $\mathbf{A}IJ(t)$ and $\mathbf{B}IJ(t)$ can be approximated by the sum of two orthogonal terms,

where the fluctuating forces $\delta \mathbf{F}IJ\u2225(t)$ and $\delta \mathbf{F}IJ\u22a5(t)$, the relative velocities $\mathbf{V}IJ\u2225(t)$ and $\mathbf{V}IJ\u22a5(t)$, and their temporal correlations $\mathbf{A}IJ(t)$ and $\mathbf{B}IJ(t)$ can be computed directly from MD simulations.

Fig. 6(a) illustrates the radial and perpendicular components of the temporal correlation functions $\mathbf{A}IJ(t)$ and $\mathbf{B}IJ(t)$ at an intermolecular distance $RIJ=2.425\sigma $. The computed memory kernels $KIJ\u2225(t)$ and $KIJ\u22a5(t)$ obtained by solving Eq. (22) are presented in Fig. 6(b). In a pairwise system, the temporal correlation functions, i.e., $\mathbf{A}IJ(t)$ and $\mathbf{B}IJ(t)$, and the memory kernel $\mathbf{K}IJ(t)$ are functions of both time and intermolecular distance. To compute the pairwise correlations and memory kernel at various distances from MD simulations, the intermolecular distance is divided into many bins with width of Δ. Then, the values of $\mathbf{A}IJ(R,t)$ and $\mathbf{B}IJ(R,t)$ are obtained by ensemble averaging these correlations over the pairs with a distance $R\u2212\Delta /2<RIJ(t)<R+\Delta /2$. Meanwhile, we also obtain the conservative force given by $FIJC(R)=946.05(1+4R/3.28)(1\u2212R/3.28)4$ with a cut off radius $Rc=3.28$ beyond which $FIJC(R)$ vanishes. For more technical details, we refer to our previous works.^{6,32}

The DPD-NM model computes the non-Markovian interactions directly as described in Section II B. To avoid prohibitive computational cost in DPD-NM simulations, we assume that the pairwise memory kernel is finite, i.e., $\mathbf{K}IJ(t)|t>N\Delta t=0$. Let the friction coefficients $\Gamma IJ,n\u2225(R)=KIJ\u2225(R,n\Delta t)$ and $\Gamma IJ,n\u22a5(R)=0.5KIJ\u22a5(R,n\Delta t)$. The factor of 0.5 appears in $\Gamma IJ,n\u22a5(R)$ because the modulus of perpendicular dissipation is equally distributed in the plane orthogonal to $\mathbf{e}IJ$.^{6,32} Then, the EOM of DPD-NM particles is given by^{32,33}

where $\mathbf{\Omega}I$ is the angular velocity of a particle *I*, $\mathbf{T}I$ is the torque, and $\mathbf{L}I$ the angular momentum. Here, we include Eq. (28) for the conservation of angular momentum of the system. Also, $\xi IJ$ is a Gaussian white noise,^{57} and $d\mathbf{W}IJA$ is an antisymmetric noise matrix.^{6} $\mathbf{V}IJ\u2225(t\u2212n\Delta t)=[\mathbf{V}IJ(t\u2212n\Delta t)\u22c5\mathbf{e}IJ(t\u2212n\Delta t)]\mathbf{e}IJ(t\u2212n\Delta t)$ is the radial velocity component and $\mathbf{V}IJ\u22a5(t\u2212n\Delta t)=\mathbf{V}IJ(t\u2212n\Delta t)\u2212\mathbf{V}IJ\u2225(t\u2212n\Delta t)$ is the perpendicular velocity component. In practical implementation, the temporal dependence and distance dependence of the memory kernel are separated, i.e., $\Gamma IJ\u2225(RIJ,t)\u2248\varphi \u2225(RIJ)\u22c5\theta \u2225(t)$ and $\Gamma IJ\u22a5(RIJ,t)\u2248\varphi \u22a5(RIJ)\u22c5\theta \u22a5(t)$. Then, the colored noise generator described in Section II B is employed to generate colored noise satisfying the second FDT. The coefficients $\alpha IJ,n\u2225$ and $\alpha IJ,n\u22a5$ involved in Eq. (9) are related to the memory kernel by

More specifically, the computed data of memory kernels can be fitted by

where $Rc\u2225=3.10$ is a cut off radius beyond which $\varphi IJ\u2225(R)$ vanishes, and $Rc\u22a5=3.09$ for $\varphi IJ\u22a5(R)$. The temporal parts shown in Fig. 6(b) are fitted by a linear combination of exponentially damped oscillators, wherein $\theta IJ\u2225(t)\u2248\u2211i=14\Lambda iexp(\u2212t/\tau i)cos(\omega it\u2212\phi i)$ and $\theta IJ\u22a5(t)\u2248\u2211i=12\Lambda i$$exp(\u2212t/\tau i)cos(\omega it\u2212\phi i)$ with the parameters listed in Table II.

. | Index i
. | $\Lambda i$ . | $\tau i\u22121$ . | $\omega i$ . | $\phi i$ . |
---|---|---|---|---|---|

$\theta IJ\u2225(t)$ | 1 | 0.07932 | 9.8295 | 59.7436 | 0.163 067 |

2 | 0.44711 | 13.1517 | 19.9223 | 0.583 477 | |

3 | 0.45455 | 11.7759 | 37.6419 | 0.303 194 | |

4 | 2.63894 | 7.1461 | 0.31109 | 1.527 280 | |

$\theta IJ\u22a5(t)$ | 1 | 1.10320 | 25.8467 | 28.9060 | 0.729 580 |

2 | 0.21980 | 7.5003 | 10.2789 | 0.630 370 |

. | Index i
. | $\Lambda i$ . | $\tau i\u22121$ . | $\omega i$ . | $\phi i$ . |
---|---|---|---|---|---|

$\theta IJ\u2225(t)$ | 1 | 0.07932 | 9.8295 | 59.7436 | 0.163 067 |

2 | 0.44711 | 13.1517 | 19.9223 | 0.583 477 | |

3 | 0.45455 | 11.7759 | 37.6419 | 0.303 194 | |

4 | 2.63894 | 7.1461 | 0.31109 | 1.527 280 | |

$\theta IJ\u22a5(t)$ | 1 | 1.10320 | 25.8467 | 28.9060 | 0.729 580 |

2 | 0.21980 | 7.5003 | 10.2789 | 0.630 370 |

Let us introduce an intermediate variable $d\mathbf{P}IJ$ for each pair to assist construction of the DPD-AUX system, where $d\mathbf{P}IJ$ represents the increment of momentum contributed from an individual pair $IJ$. By ignoring the many-body correlation between different pairs, the differential of total momentum comes from the pairwise contributions, i.e., $d\mathbf{P}I=\u2211J\u2260Id\mathbf{P}IJ$. Then, $d\mathbf{P}IJ$ for each time step is determined by conservative, dissipative, and random forces,

Based on the description in Section II C, Eq. (31) can be transformed into a Markovian equation in which the momentum is coupled to a set of auxiliary variables,

As a result, the computation of non-Markovian interactions in Eq. (31) is transferred to the evaluation of $d\mathbf{P}IJ$ for each time step, wherein $d\mathbf{P}IJ$ has both radial and perpendicular components. Since the directions $\mathbf{e}IJ\u2225$ and $\mathbf{e}IJ\u22a5$ are orthogonal to each other, we can compute the radial and perpendicular components independently. The EOM of a DPD-AUX is given by

Based on the fitting function, the corresponding matrices in Eq. (32) are given by

where $\lambda =(\varphi IJ\u2225(R))1/2$ and $\mu =(\varphi IJ\u22a5(R))1/2$. Matrix $\mathbf{B}s$ is determined by $[\mathbf{B}s\u2225][\mathbf{B}s\u2225]T=kBT(\mathbf{A}ss\u2225+[\mathbf{A}ss\u2225]T)$ and $[\mathbf{B}s\u22a5][\mathbf{B}s\u22a5]T=2kBT(\mathbf{A}ss\u22a5+[\mathbf{A}ss\u22a5]T)$ to satisfy the FDT.

With the constructed CG force field and corresponding matrices given by Eq. (35), DPD-NM and DPD-AUX simulations are performed with 8000 DPD particles in periodic cubic boxes of length $LB=50.095$. The velocity-Verlet algorithm is employed for the numerical integration in CG simulations with a time step of $\Delta t=0.005$. Since the DPD-NM model computes the non-Markovian interactions directly based on the algorithm introduced in Section II B, it has to store historical information for each pair. To avoid prohibitive computational cost, the memory kernel $KIJ\u2225(t)$ shown in Fig. 6(b) is truncated at $t=0.5$ and $KIJ\u22a5(t)$ is truncated at $t=0.25$ in practical simulations, which correspond to storing the historical information of each pair for 101 and 51 time steps in radial and perpendicular directions, respectively. By monitoring the temperature of DPD-NM and DPD-AUX systems, we observed that both of their temperatures are maintained at $kBT=1.0$ and the initial value of VACF(t = 0) shown in Fig. 7 is $0.273=3kBT/M$, which indicates correct implementations of the second FDT in both DPD-NM and DPD-AUX simulations.

The performance of DPD-NM and DPD-AUX models in reproducing the reference MD system on VACF is shown in Fig. 7(a), in which the inset shows the time integral of VACF defined by $D(t)=13\u222b0tVACF(\tau )d\tau $. Based on the Green-Kubo relations, the plateau of $D(t)$ indicates the value of diffusion constant. Comparison of time-correlation function, i.e., VACF, is more stringent requirement than macroscopic properties (diffusion and viscosity) for reproducing the microscopic dynamics. Diffusion and viscosity are integrals of corresponding time-correlation functions, which can be matched even if the time-correlation functions are not correctly reproduced.^{49} Since the entire VACF curve of the MD system is well reproduced by the DPD-NM and DPD-AUX models, as shown in Fig. 7, the macroscopic properties should be correctly produced as well. In particular, let *ν* br rate conditions, the MD systee the kinematic viscosity under low sheam of polymer melt has ${D,\nu}={6.18\xd710\u22123,13.39}$, while the DPD-NM and DPD-AUX systems have ${D,\nu}={6.82\xd710\u22123,12.96}$ with relative errors ${+10%,\u22123.2%}$ and ${5.98\xd710\u22123,13.73}$ with relative errors ${\u22123.2%,+2.5%}$, respectively.

Fig. 7(b) zooms in the initial part of VACF to examine short time dynamics, while Fig. 7(c) plots the VACF in double-logarithmic coordinates to compare the long-time hydrodynamics. These results show that the VACF of MD system is well reproduced by both DPD-NM and DPD-AUX models. Since the memory kernels used in the DPD-NM model have been well approximated by linear combinations of exponentially damped oscillators, as shown in Fig. 6(b), the DPD-AUX model has the same (better in D(t)) performance as that of the DPD-NM model. In particular, both DPD-NM and DPD-AUX models can correctly reproduce the zero initial slope of VACF consistent with the reference MD system, while the backscattering effect yielding negative VACF is also correctly reproduced.

In Fig. 8, we compare the MSD of CG particles in DPD-NM, DPD-AUX systems, and the MSD of the COM of star polymers in reference MD system. We found that all three stages of diffusion process, i.e., ballistic motion, subdiffusion, and normal diffusion, are captured by DPD-NM and DPD-AUX models. Moreover, the diffusion constants of DPD-NM and DPD-AUX system are in good agreement with that of MD system, which can be also seen in the curve of $D(t)$ plotted in the inset of Fig. 7(a). It is obvious that the DPD-NM and DPD-AUX models have better performance than the LE-NM and LE-AUX models in reproducing the reference MD system. The reason is that, unlike the memory kernel $K(t)$ shown in Fig. 3(b), the pairwise memory kernels $KIJ\u2225(t)$ and $KIJ\u22a5(t)$ shown in Fig. 6(b) decay faster and have no long-time term, indicating that they can be safely truncated and well approximated by linear combinations of exponentially damped oscillators.

For long-time hydrodynamics characterized by an algebraic decay in VACF, the momentum conservation is known to be essential for the appearance of long-time tail in VACF. Fig. 7(c) shows that both DPD-NM and DPD-AUX models can reproduce well the algebraically decaying long-time tail of VACF consistent with the reference MD system, which benefits from the rigorous conservation of momentum guaranteed by the pairwise interactions. Hence, the DPD-AUX model has both the advantages of auxiliary variable for convolution-free and pairwise interaction for momentum conservation. Using this coarse-graining strategy, the short-time dynamics can be correctly generated by coupling momentum to these auxiliary variables, while the correct long-time hydrodynamics can be captured through the pairwise interactions.

## IV. SUMMARY AND DISCUSSION

We constructed coarse-grained (CG) models by applying the Mori–Zwanzig (MZ) projection operator to microscopic dynamics generated by molecular dynamics (MD) simulations. In general, elimination of degrees of freedom from complex dynamics introduces non-negligible memory effects, resulting in a generalized Langevin equation (GLE) for the CG system in the context of the MZ formalism. Upon a pairwise decomposition, we reformulated the GLE into its pairwise version, i.e., non-Markovian dissipative particle dynamics (DPD). Both CG systems contain non-Markovian interactions.

To generate microscopic dynamics to be coarsened, we performed MD simulations of star-polymer melts. Subsequently, GLE and non-Markovian DPD models were constructed by grouping many bonded atoms into a single cluster, wherein the CG potential and the memory kernel of dissipation are computed from the MD trajectories. In particular, the GLE is the original formulation resulted from applying the MZ projection operator, while non-Markovian DPD is the pairwise form of GLE under a pairwise decomposition. For numerical implementation of the non-Markovian interactions, we can either perform the non-Markovian formulation directly with colored noise, or introduce auxiliary variables to model the non-Markovian interactions. Therefore, we tested four different strategies to implement the CG simulation, i.e., LE-NM, LE-AUX, DPD-NM, and DPD-AUX. Here, LE-NM and LE-AUX simulate non-interacting particle dynamics governed by GLE, while DPD-NM and DPD-AUX are interacting systems with pairwise interactions. All these CG systems were compared with each other as well as with the original MD system of polymer melt to quantify their difference and performance on reproducing the reference MD system.

Our simulation results demonstrate that a Markovian dynamics with auxiliary variables successfully generates accurate non-Markovian dynamics consistent with the reference MD system at short time scales. All the three stages of diffusion process of the MD system, i.e., ballistic motion, subdiffusion, and normal diffusion can be captured by the CG models that generate non-Markovian dynamics. For problems where only short time dynamics are concerned, both LE-NM and LE-AUX can be applied to reproduce an underlying multiscale dynamics at a CG level. However, if the long-time hydrodynamics are also important and have to be correctly preserved, the LE-NM and LE-AUX models will be inaccurate. The reason is that the momentum conservation is known to be essential for the appearance of long-time tail of VACF, but neither LE-NM nor LE-AUX conserves the momentum and hence the correct long-time hydrodynamics cannot be expected. The CG systems with pairwise interactions, i.e., DPD-NM and DPD-AUX, conserve the momentum rigorously and can reproduce the correct long-time hydrodynamics characterized by an algebraic decay in VACF. Quantitative comparison between different CG models reveals that introducing auxiliary variables into a CG system with pairwise interactions combines both advantages of auxiliary system and pairwise formulation, where the former yields accurate short-time dynamics while the latter guarantees momentum conservation for capturing the correct long-time hydrodynamics.

Although the LE-NM and LE-AUX models have difficulty in capturing the long-time algebraic decay in VACF because of violation of momentum conservation, they reproduced the short-time dynamics accurately at a much lower computational complexity than the pairwise DPD-NM and DPD-AUX models. We note that keeping longer memory in LE-NM or using more auxiliary variables in LE-AUX can gradually improve their long-time dynamics. Also the exact memory kernel resulting from the Mori-Zwanzig projection contains non-linear terms, which are not included in models for practical purposes. Incorporating the non-linear terms of memory and non-Gaussian noise in GLE will be a topic for future work. To test the computational cost of different models, we performed simulations with a system containing 8000 CG particles on a single core of Intel Xeon E5-2637 CPU. We found that the LE-NM simulation needs 0.68 s per time step that is reduced to 0.02 s per time step using the LE-AUX model, while the DPD-NM simulation needs 8.3 s per time step that is reduced to 0.42 s per time step using the DPD-AUX model. As shown in Table III, the non-interacting particle models (LE-NM and LE-AUX) are computationally much cheaper than the pairwise models (DPD-NM and DPD-AUX). Therefore, for the problems where the long-time tail in VACF is important, the DPD-AUX model can be applied as a replacement of DPD-NM to achieve the same accuracy at a lower computational cost. However, if the algebraic decay of the long-time tail in VACF is not important, the non-interacting particle models can be applied with much less computational cost.

. | Model . | Cost . | Speedup . |
---|---|---|---|

Non-interacting particle model | LE-NM | 0.68 | $\u2026$ |

LE-AUX | 0.02 | 34.0 | |

Pairwise interaction model | DPD-NM | 8.3 | $\u2026$ |

DPD-AUX | 0.42 | 19.8 |

. | Model . | Cost . | Speedup . |
---|---|---|---|

Non-interacting particle model | LE-NM | 0.68 | $\u2026$ |

LE-AUX | 0.02 | 34.0 | |

Pairwise interaction model | DPD-NM | 8.3 | $\u2026$ |

DPD-AUX | 0.42 | 19.8 |

Because the dynamics of complex fluids often involves different time scales,^{58} mesoscopic models bridging the atomistic and continuum scales can be a promising candidate for tackling challenges in multiscale modeling. In this work, we have demonstrated that a multiscale dynamics, which spans multiple time scales, can be generated by a combination of different algorithms. More specifically, the effects of atomistic collisions at small time scales on collective dynamics of CG clusters can be captured by coupling CG momentum to auxiliary variables, while the hydrodynamics at large time scales can be generated through implementation of pairwise interactions.

Moreover, the present work provides quantitative comparisons between different CG models derived from the MZ formalism, and also their accuracy in reproducing the reference non-Markovian dynamics. With a better understanding of the non-interacting particle model and the pairwise interaction model, their benefits and limits for modeling non-Markovian dynamics have been clarified, which would help choose an appropriate model for practical applications as a compromise between computational complexity and performance. We note that the present results are based on the non-Markovian dynamics with negligible many-body correlations. However, in an aggressive coarse-graining, the correlations between different pairs of CG clusters may be strong and cannot be ignored during the coarse-graining procedures.^{33,49} In future work, we plan to explore how to incorporate the many-body effects into effective CG dynamics at more aggressive coarse-graining levels, where we can take advantage of the framework of many-body DPD^{59,60} whose pairwise potential depends not only on the distance but also on the density of neighboring particles.

## ACKNOWLEDGMENTS

We would like to thank the anonymous referees whose insightful comments helped us to improve our paper. This work was primarily supported by 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 with the resources of the Argonne Leadership Computing Facility (ALCF) and the resources of the Oak Ridge Leadership Computing Facility (OLCF). Z. Li would like to acknowledge Dr. Xin Bian, Dr. Changho Kim, and Dr. Yuta Yoshimoto for helpful discussions.

### APPENDIX A: FITTING THE MEMORY KERNEL

The extended dynamics using auxiliary variables is an approximation of the original non-Markovian dynamics as described in Section II C. The validity and accuracy of the extended dynamics depend on the fitting of the computed memory kernels. To analyze the sensitivity of the auxiliary system to the fitting of memory kernels, we compare the results of the LE-AUX model with different fitting functions. Figure 9 shows the normalized memory kernel obtained from the data of MD simulations by solving (18).

The computed memory kernel $K(t)$ is fitted by a linear combination of 2, 3, and 5 exponentially damped oscillators in the form of $f(t)=\u2211i=1n\Lambda iexp(\u2212t/\tau i)cos(\omega it\u2212\phi i)$, which means that the momentum of each particle in the LE-AUX system is coupled with 4, 6, and 10 auxiliary variables, respectively. The inset of Fig. 9 gives a global view of $K(t)$ and these fitting functions. We observe that the fitting function with 2 terms diverges significantly from the memory kernel $K(t)$. Therefore, its corresponding LE-AUX system does not reproduce well the VACF of the reference MD system. However, it can be seen in Fig. 9 that the convergence of the fitting is fast. More specifically, the fitting function with 3 terms already represents the memory kernel $K(t)$ well, and the fitting function of 5 terms improves slightly further. As a result, the VACFs of their extended dynamics show good agreement with the VACF of the reference MD system.

As shown in Fig. 10, although the long-time tail of VACF in extended dynamics generated by the LE-AUX model will eventually show exponential decay, using more terms in the fitting function $f(t)=\u2211i=1n\Lambda iexp(\u2212t/\tau i)cos(\omega it\u2212\phi i)$, which corresponding to employing more auxiliary variables in the extended dynamics, the long-time tail of VACF can be captured more and more accurately.

### APPENDIX B: CONSTRUCTION OF COUPLING MATRICES

Suppose we have the coupling matrix **A** in the form of

where $\mathbf{A}ss=((a,\u2212b)T,(b,0)T)$ with a positive *a* ensures that the matrix $(\mathbf{A}ss+\mathbf{A}ssT)$ is positive semidefinite. The matrix $\mathbf{B}s$ in Eq. (11) is determined by $\mathbf{B}s\mathbf{B}sT=kBT(\mathbf{A}ss+\mathbf{A}ssT)$ to satisfy the fluctuation-dissipation theorem. We do eigen-decomposition of $\mathbf{A}ss$ and obtain

where $\chi 1=12(a\u2212a2\u22124b2)$ and $\chi 2=12(a+a2\u22124b2)$ are the two eigenvalues of the matrix $\mathbf{A}ss$. The memory kernel in the form of Eq. (13b) can be computed by

where $\lambda 1=c12+c22$, $\lambda 2=(c22\u2212c12)a/4b2\u2212a2$ and $\omega =4b2\u2212a2/2$. When a memory kernel $K(t)$ is obtained from MD simulations by solving Eq. (18) or (22), we can fit the computed memory kernel by $K(t)=pexp(\u2212qt)cos(rt\u2212s)$. Then, we have all the elements of the coupling matrix **A** of Eq. (B1) given by

In practical implementations, we take $c1=0$ so that the derivative of the memory kernel is zero at $t=0$, i.e., $dK(t)/dt|t=0=0$. When there are more than one term in fitting of the computed memory kernel $K(t)$, we just extend the block of **A** in the form of Eq. (B1), as shown in Eqs. (19) and (35).

## REFERENCES

The spurious periodicity can be mitigated by using a noise with longer memory but at an extra computational and memory cost.