Computing the non-Markovian coarse-grained interactions derived from the Mori–Zwanzig formalism in molecular systems: Application to polymer melts

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.⁵ Commonly, the effect of the unresolved DOFs on the CG variables is approximated by stochastic dynamics.⁶ 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,¹⁷ the reverse Monte Carlo,¹⁸ the force matching method,¹⁹ 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.²⁷ More specifically, Izvekov and Voth²⁸ 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.⁵ 

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²⁹ 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∼(t)$ with the Dirac delta function $δ(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³⁶ 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)³⁷ or employ auxiliary variables to replace the non-Markovian dynamics with a Markovian process in a higher dimensional space.³⁸ 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,³⁸ 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.³⁹ 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.

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.

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

where $H(Γ)$ defines the phase space trajectories of the system $Γ≡{𝐫i,𝐩i,i=1,…,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 $𝐑I=MI−1∑mIi𝐫Ii$ and $𝐏I=∑𝐩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 $β=1/kBT$ with $kB$ the Boltzmann constant and T the thermodynamic temperature, $𝐑={𝐑I,I=1,…,K}$ is a phase point in the CG phase space, and $ω(𝐑)$ is defined as a normalized partition function of all the microscopic configurations at a phase point R.³⁰ 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,⁴⁰ 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., $⟨[δ𝐅IQ(t)][δ𝐅JQ(t′)]T⟩=𝜹IJ𝐊(t−t′)$⁠, the generalised Langevin equation describing the evolution of P can be written as

where $δ𝐅Q$ is the random force, and $𝐊(t)=β⟨[δ𝐅Q(t)][δ𝐅Q(0)]T⟩$ is the memory kernel whose definition ensures that the second FDT is satisfied.³⁷ 

In the second CG approach, we first perform a pairwise decomposition,³² i.e., $𝐅I≈∑J≠I𝐅IJ$⁠, and assume negligible many-body correlations between different pairs, i.e., $⟨[δ𝐅IJQ][δ𝐅IKQ]T⟩J≠K≈0$⁠. Then, Eq. (3) can be reformulated into a pairwise form

where $𝐅IJ$ is the instantaneous force between CG particles I and J. Its ensemble average $⟨𝐅IJ⟩$ is taken as the conservative force in the CG model. Also, $𝐕IJ(t)=𝐕I(t)−𝐕J(t)$ denotes the relative velocity between the CG particles I and J. The last term $δ𝐅IJQ$ is the random force resulted from unresolved degrees of freedom. In general, $δ𝐅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 $𝐊IJ(t)=β⟨[δ𝐅IJQ(t)][δ𝐅IJQ(0)]T⟩$⁠.

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∼$ 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Δt$ where $n≥0$ represent discrete time steps,^32,33 also let $Kn=K(tn)$ and $F∼n=F∼(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Δt$⁠, where $Δ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∼n$ satisfying the second FDT given by Eq. (7). Although the Fourier transform⁴³ can be applied to generate colored noise satisfying Eq. (7), the noise represented by Fourier series will introduce an artificial periodicity into the system,⁴⁴ which affects the particle dynamics significantly.³² Here, we adopt an optimization-based scheme proposed in our previous work³² for colored noise generation to avoid the artificial periodicity. More specifically, the random force $F∼n$ at discrete time steps is given by

where $α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∼$ by

where $fn=kBT∑s=0Nαsα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 $αs$ is obtained by performing numerical optimization for minimizing the difference between $fn−m$ and $kBTKn−m$⁠. In practical implementation, the globalized bounded Nelder-Mead algorithm⁴⁵ is applied for the optimization. We note, however, that there are many other optimization methods⁴⁶ that can be used to obtain the coefficients $αs$ by defining the test function $fn=kBT∑s=0Nαsα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.

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

where $𝐀sp=−𝐀psT$⁠, 𝝃 is a vector of uncorrelated Gaussian random numbers with $⟨ξi(t)ξj(0)⟩=δijδ(t)$⁠. Let $𝐱=(𝐏,𝐬)T$⁠, the dynamics of x is the Ornstein-Uhlenbeck (OU) process $𝐱˙=𝐀𝐱+𝐁𝝃$ with $𝐀=((𝟎,𝐀sp)T,(𝐀ps,𝐀ss)T)$ and $𝐁=diag(𝟎,𝐁s)$⁠. Assuming that the non-Markovian dynamics has finite memory, one can safely take $𝐬(−∞)=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³⁸ 

where the matrix $𝐁s$ is determined by $𝐁s𝐁sT=kBT(𝐀ss+𝐀ssT)$ to satisfy the fluctuation-dissipation theorem. In practical implementations, we can always find this $𝐁s$ by designing $𝐀ss$ such that $(𝐀ss+𝐀ssT)$ is positive-definite or positive-semidefinite and using the Cholesky factorization.⁴⁷ When $𝐀ss$ is a real matrix and has complex eigenvalues with positive real part, the memory kernel $𝐊(t)$ is a linear combination of exponentially damped oscillators.³⁸ 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 $𝐊(t)$⁠.

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,⁵⁰ 

where r is the distance between two monomers, and the cut off distance $rc=21/6σ$ 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,⁵¹ 

where k is the spring constant and $r0$ determines the maximum length of the spring. In particular, a parameter set of $k=30ϵ/σ2$ and $r0=1.5σ$ is adopted so that the FENE spring is stiff and short enough to minimize artificial bond crossings.⁵² 

MD simulations of polymer melts are performed in the canonical ensemble (NVT) through the Nośe–Hoover thermostat.⁵³ More specifically, 8000 molecules of star polymer are filled into a periodic cubic box of length $LB=50.095σ$ to model a polymer melt. The number density of monomers is $ρ=0.7σ−3$⁠. 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 $σ=1$⁠, $m=1$⁠, $ϵ=1$⁠, and $τ=σ(m/ϵ)1/2$⁠. In addition, the temperature of the MD system is maintained at $kBT=1.0$ and the time step $δt=0.001τ$ 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 $τv=0.61$ and $τf=0.22$⁠, respectively. Let $κ=τv/τf=2.77$ be the ratio of the time scales, then it is obvious that $τv$ and $τ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,⁵⁴ 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∫0tVACF(τ)dτ$⁠, 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∼t0.57$ before a normal diffusion $MSD∼t$ is reached at long time scales. Such behaviors of the polymer melt shown in VACF and MSD will be reproduced by the CG models.

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 $δ𝐅(t)$ and the memory kernel $𝐊(t)$ should be related by the second FDT. Therefore, if we compute $𝐊(t)$⁠, we can model $δ𝐅(t)$ accordingly using the second FDT.

Multiplying the velocity $𝐕T(0)$ on both sides of Eq. (16) and taking the ensemble average yields

Let us now define the velocity autocorrelation function as $𝐂(t)=⟨[𝐕(t)][𝐕(0)]T⟩$⁠. Since the random force $δ𝐅(t)$ is statistically independent of the velocity $𝐕(0)$⁠, we can rewrite Eq. (17) into a form of Volterra equation

where both $𝐂(t)$ and the force-velocity correlation function (FVCF) $M𝐂˙(t)=⟨[𝐅(t)][𝐕(0)]T⟩$ can be evaluated from MD simulations. Here, $𝐅(t)$ is the total force, and $𝐕(t)$ is the velocity of the COM of star polymers. Then, the memory kernel $𝐊(t)$ can be obtained by solving the Volterra equation (18) with known $𝐂˙(t)$ and $𝐂(t)$⁠.

Fig. 3(a) shows the computed VACF $𝐂(t)$ and FVCF $𝐂˙(t)$ obtained directly from MD simulations, while Fig. 3(b) displays the corresponding memory kernel $𝐊(t)$ and its inset gives a global view of $𝐊(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 $Δ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 $αn$⁠, which are obtained by numerical optimization using the globalized bounded Nelder-Mead algorithm.⁴⁵ 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≤10$⁠. 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.⁵⁵ 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.⁵⁶ 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$∼t0.57$ before a normal diffusion (MSD$∼t$⁠) 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.³⁸ 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)=∑i=13Λiexp(−t/τi)cos(ωit−φi)$ with parameters $(Λi,τi−1,ωi,φi)i=1,2,3$ = ${(355.53$⁠, 23.577, 29.954, $0.666 817 6)$⁠, $(1.935×105$⁠, 9.559, $0.006 42$⁠, $1.570 125)$⁠, $(39 681.12$⁠, $0.248 48$⁠, $7.298×10−5$⁠, $1.570 503)}$⁠. 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,⁵⁵ i.e., $D=kBT/ξ$ where $ξ=∫0∞K(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 $𝐁s𝐁sT=kBT(𝐀ss+𝐀ssT)$ is a sufficient condition to satisfy the FDT for having canonical sampling. The details of construction of the coupling matrices A and $𝐁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 $𝐊(t)=−𝐀psexp(−t𝐀ss)𝐀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,⁵⁴ 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.

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 $δ𝐅IJ(t)$ as the difference between the instantaneous force and the mean force, i.e., $δ𝐅IJ(t)=𝐅IJ(t)−⟨𝐅IJ⟩$⁠. The second equality of Eq. (5) becomes

where $𝐊IJ(t)=β⟨[δ𝐅IJQ(t)][δ𝐅IJQ(0)]T⟩$⁠. By multiplying by the velocity $𝐕IJT(0)$ on both sides of Eq. (20) and taking ensemble average, we have

where the third term $⟨δ𝐅IJQ(t)𝐕IJT(0)⟩$ disappears because the random force $δ𝐅IJQ(t)$ comes from orthogonal dynamics in MZ procedure²⁹ and has no correlation with the CG velocity. Let us now define a function $𝐀IJ(t)=⟨δ𝐅IJ(t)𝐕IJT(0)⟩$ and another function $𝐁IJ(t)=⟨𝐕IJ(t)𝐕IJT(0)⟩$⁠. Then, Eq. (21) can be rewritten as

Considering that the fluctuating force $δ𝐅IJ(t)$ and the relative velocity $𝐕IJ(t)$ have both radial components along vector $𝐞IJ$ and perpendicular components orthogonal to vector $𝐞IJ$⁠, we decompose the fluctuating force as well as the relative velocity into two orthogonal parts,

where the symbol “∥” represents the radial component along $𝐞IJ$ and “⊥” for perpendicular component. Since the two directions of “∥” and “⊥” are orthogonal to each other, the temporal correlations $𝐀IJ(t)$ and $𝐁IJ(t)$ can be approximated by the sum of two orthogonal terms,

where the fluctuating forces $δ𝐅IJ∥(t)$ and $δ𝐅IJ⊥(t)$⁠, the relative velocities $𝐕IJ∥(t)$ and $𝐕IJ⊥(t)$⁠, and their temporal correlations $𝐀IJ(t)$ and $𝐁IJ(t)$ can be computed directly from MD simulations.

Fig. 6(a) illustrates the radial and perpendicular components of the temporal correlation functions $𝐀IJ(t)$ and $𝐁IJ(t)$ at an intermolecular distance $RIJ=2.425σ$⁠. The computed memory kernels $KIJ∥(t)$ and $KIJ⊥(t)$ obtained by solving Eq. (22) are presented in Fig. 6(b). In a pairwise system, the temporal correlation functions, i.e., $𝐀IJ(t)$ and $𝐁IJ(t)$⁠, and the memory kernel $𝐊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 $𝐀IJ(R,t)$ and $𝐁IJ(R,t)$ are obtained by ensemble averaging these correlations over the pairs with a distance $R−Δ/2<RIJ(t)<R+Δ/2$⁠. Meanwhile, we also obtain the conservative force given by $FIJC(R)=946.05(1+4R/3.28)(1−R/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., $𝐊IJ(t)|t>NΔt=0$⁠. Let the friction coefficients $ΓIJ,n∥(R)=KIJ∥(R,nΔt)$ and $ΓIJ,n⊥(R)=0.5KIJ⊥(R,nΔt)$⁠. The factor of 0.5 appears in $ΓIJ,n⊥(R)$ because the modulus of perpendicular dissipation is equally distributed in the plane orthogonal to $𝐞IJ$⁠.^6,32 Then, the EOM of DPD-NM particles is given by^32,33

where $𝛀I$ is the angular velocity of a particle I, $𝐓I$ is the torque, and $𝐋I$ the angular momentum. Here, we include Eq. (28) for the conservation of angular momentum of the system. Also, $ξIJ$ is a Gaussian white noise,⁵⁷ and $d𝐖IJA$ is an antisymmetric noise matrix.⁶ $𝐕IJ∥(t−nΔt)=[𝐕IJ(t−nΔt)⋅𝐞IJ(t−nΔt)]𝐞IJ(t−nΔt)$ is the radial velocity component and $𝐕IJ⊥(t−nΔt)=𝐕IJ(t−nΔt)−𝐕IJ∥(t−nΔt)$ is the perpendicular velocity component. In practical implementation, the temporal dependence and distance dependence of the memory kernel are separated, i.e., $ΓIJ∥(RIJ,t)≈ϕ∥(RIJ)⋅θ∥(t)$ and $ΓIJ⊥(RIJ,t)≈ϕ⊥(RIJ)⋅θ⊥(t)$⁠. Then, the colored noise generator described in Section II B is employed to generate colored noise satisfying the second FDT. The coefficients $αIJ,n∥$ and $αIJ,n⊥$ 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∥=3.10$ is a cut off radius beyond which $ϕIJ∥(R)$ vanishes, and $Rc⊥=3.09$ for $ϕIJ⊥(R)$⁠. The temporal parts shown in Fig. 6(b) are fitted by a linear combination of exponentially damped oscillators, wherein $θIJ∥(t)≈∑i=14Λiexp(−t/τi)cos(ωit−φi)$ and $θIJ⊥(t)≈∑i=12Λi$$exp(−t/τi)cos(ωit−φi)$ with the parameters listed in Table II.

Let us introduce an intermediate variable $d𝐏IJ$ for each pair to assist construction of the DPD-AUX system, where $d𝐏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𝐏I=∑J≠Id𝐏IJ$⁠. Then, $d𝐏IJ$ for each time step is determined by conservative, dissipative, and random forces,