We propose a new coarse-grained (CG) molecular simulation technique based on the Mori–Zwanzig (MZ) formalism along with the iterative Boltzmann inversion (IBI). Non-Markovian dissipative particle dynamics (NMDPD) taking into account memory effects is derived in a pairwise interaction form from the MZ-guided generalized Langevin equation. It is based on the introduction of auxiliary variables that allow for the replacement of a non-Markovian equation with a Markovian one in a higher dimensional space. We demonstrate that the NMDPD model exploiting MZ-guided memory kernels can successfully reproduce the dynamic properties such as the mean square displacement and velocity autocorrelation function of a Lennard–Jones system, as long as the memory kernels are appropriately evaluated based on the Volterra integral equation using the force-velocity and velocity-velocity correlations. Furthermore, we find that the IBI correction of a pair CG potential significantly improves the representation of static properties characterized by a radial distribution function and pressure, while it has little influence on the dynamic processes. Our findings suggest that combining the advantages of both the MZ formalism and IBI leads to an accurate representation of both the static and dynamic properties of microscopic systems that exhibit non-Markovian behavior.

Soft matter and biological systems intrinsically involve hierarchical phenomena that are governed by interactions and processes on a wide range of spatiotemporal scales.1,2 Therefore, it is of great importance to elucidate static and dynamic properties at each scale and seamlessly link them towards a better understanding of such multiscale phenomena. Molecular dynamics (MD)3,4 is a powerful tool for analyzing these properties in a microscopic system where the constitutive equation for a continuum description is unknown. However, all-atom MD simulations of complex systems are computationally expensive because they explicitly deal with fine-grained atomistic motions according to Newton’s law, and hence they are typically impractical for the application to mesoscale phenomena. Coarse-grained (CG) molecular simulations have been often employed to bridge the scales from a microscopic level to a mesoscopic one.5,6 In CG simulations, a single CG particle represents multiple atoms or molecules, allowing for tracing mesoscopic trajectories at significantly larger spatiotemporal scales compared to all-atom MD. Obviously, the qualities of static and dynamic properties of a CG system are determined by the description of the CG force field and its fluctuations.

In general, there are two methods for the construction of CG models. The first strategy is the so-called reverse coarse-graining, where a CG force field is determined by tuning functional forms or parameters so that target properties such as radial distribution functions (RDFs), pressure, isothermal compressibility, and mean force acting on centers-of-mass (COMs) of molecules are reproduced as accurately as possible. Indeed, most CG models that have been proposed thus far fall into this category. Among such reverse coarse-graining techniques are, for instance, the force-matching method,7,8 iterative Boltzmann inversion (IBI),9,10 reverse Monte Carlo,11,12 stochastic parameter optimization using Bayesian inference,13,14 and minimization of relative entropy.15,16 Certainly, these CG models can accurately represent the target properties, but the qualities of other properties are not generally guaranteed. In particular, an effective CG potential alone cannot yield correct dynamic properties such as diffusivity, viscosity, or, more generally, time correlation functions, unless the dynamic properties themselves are chosen as target properties. More specifically, a coarse-grained molecular dynamics (CGMD) simulation generally gives rise to faster diffusion dynamics compared to MD17,18 since the effective CG potential is not accompanied by friction and fluctuating forces responsible for dynamic properties, which should be taken into account to compensate for the lost degrees-of-freedom (DOFs) as a consequence of coarse-graining. By contrast, in the second strategy—referred to as forward coarse-graining—a CG force field is directly constructed from microscopic trajectories of an underlying fine-grained system based on the Mori–Zwanzig (MZ) formalism,19–30 inherently involving friction and fluctuating terms. In the forward coarse-graining, once relevant variables are chosen, the generalized Langevin equation (GLE) to which they are subject can be straightforwardly derived using the projection operator method while retaining an explicit relation to the microscopic dynamical variables. Therefore, there are no free parameters to be tuned in principle, although a few approximations are required for practical implementation.

Recently, we have proposed non-Markovian dissipative particle dynamics (NMDPD) models,25,29 which are derived from the MZ-guided GLE on the assumption that forces acting on CG particles are pairwise decomposable and that many-body correlations between different pairs of CG particles are negligible. The GLE-based NMDPD with pairwise interactions allows for an accurate representation of hydrodynamics because it rigorously satisfies mass and momentum conservation. Moreover, in contrast to the conventional DPD based on the Markovian approximation,31–33 the NMDPD models explicitly incorporate memory effects, corresponding to the extension of the conventional DPD. In fact, we found that a typical time scale of CG particle motion is comparable to that of memory kernels in dense liquid systems, thus necessitating the incorporation of the memory effects into CG modeling.25,29 Subsequently, we demonstrated that incorporating the memory effects based on the MZ formalism improves the dynamic properties remarkably. More recently, Li et al.30 have proposed MZ-guided NMDPD employing auxiliary variables, which enable the replacement of a non-Markovian equation with a Markovian one in a higher dimensional space. The technique of introducing auxiliary variables has also been investigated by other researchers to address the complexity in computing the memory term in the GLE.34–36 Unlike the original NMDPD which needs to store historical variables to explicitly compute memory kernels, the NMDPD employing auxiliary variables yields a convolution-free Markovian equation, leading to orders of magnitude less memory usage and faster computation in terms of numerical simulation.

The MZ-guided modeling of memory kernels is indeed a powerful choice for the faithful representations of dynamic properties including time correlation functions. However, static properties such as a RDF and pressure, which are predominantly determined by an effective CG potential, may not be generated as accurately compared to the dynamic ones since our NMDPD model utilizes a pair CG potential obtained from a direct ensemble average of instantaneous forces between two clusters composed of microscopic particles.24,25 Indeed, the pair CG potentials obtained in that manner fail to yield RDFs of corresponding MD systems where Lennard–Jones (LJ) clusters are greatly overlapping, implying non-negligible many-body effects.25 To tackle this issue, we propose a new coarse-graining strategy that combines forward and reverse coarse-graining methods. More specifically, we demonstrate that incorporating memory effects based on the MZ-guided GLE, followed by the IBI correction of a CG potential, is an effective technique for constructing a CG model that can successfully represent both static and dynamic properties of an underlying microscopic system. In addition, we compare the memory kernels constructed from the Volterra integral equation using force-velocity and velocity-velocity correlations with their zeroth-order approximate counterparts, which are also referred to as a Q-approximation.26,27,37

The remainder of this paper is organized as follows. In Sec. II, we detail the MZ-guided NMDPD employing auxiliary variables followed by a brief review on the IBI method. Section III describes a reference MD system to be coarse-grained, where the COM positions and momenta of each LJ cluster are regarded as relevant variables. Here, we present in detail the procedure for extracting a pair CG potential and memory kernels together with their zeroth-order approximate forms from the MD simulation. In Sec. IV, we describe the method for implementation of the NMDPD models exploiting auxiliary variables based on the CG potential and memory kernels. Subsequently, we evaluate the NMDPD models by comparing static and dynamic properties with those of the MD system. Finally, we discuss the effect of the IBI correction of the CG potential on the static and dynamic properties. Concluding remarks are presented in Sec. V.

The MZ projection operator method19,20,38 provides us a systematic way for extracting relevant dynamical variables from a very high DOF system. Let us consider a microscopic system composed of N interacting particles. The microscopic particles are grouped into K clusters, and the COM positions and momenta of each cluster {RI, PI; I = 1, …, K} are regarded as the relevant variables. By introducing the projection operator P onto the subspace spanned by the relevant variables and the orthogonal counterpart Q=(1P), the equation of motion (EOM) for the COMs equivalent to CG particles can be derived in the form of the GLE21,23–25,28–30

dPIdt=1βRIlnωRβJ=1K0tdtδFIQtt×δFJQ0TPJtMJ+δFIQt,
(1)

where β = (kBT)−1 with kB as the Boltzmann constant and T as the thermodynamic temperature, R = {RI; I = 1, …, K} is a phase point of the CG configuration space, and MI is the mass of CG particle I; also, ω(R) is a normalized partition function associated with all the microscopic configurations at phase point R. For a detailed definition of the projector operator P, see Eq. (17) in Ref. 21. Each term of the right-hand side of Eq. (1) corresponds to mean force, friction force, and time-correlated fluctuating force, respectively. The friction and fluctuating forces intrinsically appear in Eq. (1) to compensate for the lost DOFs as a consequence of coarse-graining. The friction force involves the time convolutions of the memory kernels βδFIQtδFJQ0T and the momenta of CG particles, satisfying the second fluctuation-dissipation theorem (FDT).21,39 More generally, nonlinear terms appear in the memory when the stochastic force is not fully decoupled from the past CG coordinates, resulting in position- and momentum-dependent memory functions.26 However, it is very demanding to compute the nonlinear terms for realistic systems. In a practical sense, the linear approximation of the memory in Eq. (1) is considered to be reasonable unless the system is far from equilibrium.

Although GLE (1) accurately represents the trajectories in the CG phase space, direct computation of Eq. (1) is very demanding even for one-dimensional harmonic chain systems40,41 since the mean force depends on all the positions of CG particles, while the friction force depends on all the momenta of CG particles, thus requiring the evaluation of multi-body interactions. For practical implementation, we employ two approximations. The first approximation is a pairwise decomposition, where the force acting on CG particle I is given in a pairwise additive form FItJIFIJt with FIJt as the instantaneous force acting from CG particle J to I. The second approximation is to neglect many-body correlations between different pairs δFIJQtδFIXQ0TXJ0. These two approximations render Eq. (1) a GLE in a pairwise additive form corresponding to NMDPD25,29

dPIdt=JIFIJt=JIFIJ0tdtKIJttVIJt+δFIJQt,
(2)

where FIJ is an ensemble average of FIJ, δFIJQt is a pairwise fluctuating force acting from CG particle J to I, KIJt=βδFIJQtδFIJQ0T is a pairwise memory kernel satisfying the second FDT, and VIJ=VIVJ is a relative velocity with VI as the velocity of CG particle I. If KIJ(t) reduces to 2γIJδt in Eq. (2), where γIJ0KIJtdt, then the EOM for Markovian DPD is recovered.

To implement the EOM for NMDPD (2), we introduce the intermediate variables defined for each pair of CG particles dPIJ, which represent the increments of momenta contributed from CG particle J to I.30 Accordingly, the time variation of PI is given by

dPIdt=JIdPIJdt,
(3a)
dPIJdt=FIJ0tdtKIJttVIJt+δFIJQt.
(3b)

We note that the pairwise intermediate variables dPIJ must be calculated for each pair I and J, namely, dPIJ/dt ≠ dPI/dt − dPJ/dt. Then, instead of explicitly evaluating the non-Markovian equation (3b),25,29,42 we introduce auxiliary variables sIJ for each pair I and J, which are linearly coupled to PIJ and themselves. Consider the following stochastic differential equations with respect to PIJ and sIJ:30,34

ṖIJṡIJ=FIJ00ApsAspAssPIJsIJ+000Bs0ξIJ,
(4)

where Asp=ApsT and ξIJ is a vector of uncorrelated Gaussian random numbers with each element ξIJi satisfying ξIJit=0 and ξIJitξKLj0=δijδIKδJL+δILδJK×δt. The matrices A0,AspT,Aps,AssT and B diag0,Bs are termed drift and diffusion matrices. Assuming that non-Markovian dynamics has finite memory, we can safely take sIJ=0 and solve Eq. (4) with respect to sIJ,

sIJt=tdtettAssAspPIJt+BsξIJt.
(5)

Substituting Eq. (5) into Eq. (4) allows for the elimination of sIJ in the resultant equation for PIJ,30,34

ṖIJt=FIJtdtK̃IJttPIJt+ζIJt,
(6a)
K̃IJt=ApsetAssAspt0,
(6b)
ζIJt=tdtApsettAssBsξIJt.
(6c)

The lower limit of the integration in Eq. (3b) can also be taken as −∞ since the pairwise memory kernel is finite, i.e., KIJtt>t=0. Therefore, Eq. (6) corresponds to Eq. (3b), meaning that solving Eq. (3b) in a non-Markovian form is equivalent to solving Eq. (4) in a Markovian form. The second FDT reads

BsBsT=kBTAss+AssT
(7)

for the Markovian dynamics given by Eq. (4).30,34 Equation (7) can be seen as an extension of the standard second FDT for free-particle Langevin dynamics, i.e., b2 = 2akBT with a being the friction coefficient and b being the intensity of random force. Practically, Bs can be obtained through the Cholesky factorization43 by designing Ass+AssT to be positive-definite or positive-semidefinite. This constraint also ensures that the stochastic process is consistent with the second law of thermodynamics since the Fourier transform of the memory kernel is positive for all real frequencies.44 It has been reported that this modest loss of generality does not significantly affect the accuracy or flexibility of the representation of the memory kernel.34 Furthermore, when Ass is a real matrix and has complex eigenvalues with positive real parts, K̃IJt can be expressed as a linear combination of exponentially damped oscillators.30,34

In contrast to the MZ formalism, the IBI is a structure-based inverse coarse-graining method for obtaining a CG pair potential that reproduces the RDF of a reference system.9 The simple Boltzmann inversion yields

UPMFR=kBTlngrefR,
(8)

where R is the separation between CG particles, grefR is the target reference RDF, and UPMF(R) is the potential of mean force (PMF). Unfortunately, this PMF cannot be directly adopted as a pair potential for a CG model since it inherently involves multibody contribution from all the particles in the system.45 Therefore, the iterative update of the pair potential is performed in the standard IBI as

Ui+1CGR=UiCGRkBTlngiRgrefR,
(9)

where the subscript i denotes the iteration number. The RDF after the ith update gi(R) is usually obtained through a CGMD simulation exploiting the pair potential UiCGR. This iterative update is performed until the converged pair potential that reproduces gref(R) is attained. Typically, UPMF(R) given by Eq. (8) is adopted as an initial guess of the pair potential U0CGR, although in the present study, we utilize a MZ-guided pair potential as an initial guess as described in Sec. IV C. In addition, the practical implementation of the IBI requires smoothing and extrapolation of the potential to provide continuous forces.

In principle, the Henderson theorem46 guarantees a one-to-one correspondence between the target RDF gref(R) and CG pair potential. However, in practical molecular simulations, the cutoff distance for pairwise interactions breaks down the uniqueness of the CG pair potential. More specifically, different cutoff distances lead to different CG pair potentials, all of which yield identical RDFs but different pressures. In the present study, we choose the cutoff distance obtained from the MZ-guided pair potential because it ensures the accurate representations of both the RDF and pressure.47 

Based on the MZ formalism, we extract the pair potential and memory kernels from a LJ system, where 10 LJ particles are grouped into a single cluster and the COM positions and momenta of each cluster {RI, PI; I = 1, …, K} are regarded as CG variables (Fig. 1). The system size is 158 × 158 × 158 Å3 with periodic boundary conditions in all the directions. The system consists of 8000 clusters and thus 80 000 LJ particles. The interatomic potential is given by the LJ (12-6) form

ϕr=4εσr12σr6,
(10)

where r is the interatomic distance, σ = 3.4 Å, and ε = 1.67 × 10−21 J corresponding to argon atoms with the mass of 6.63 × 10−26 kg. The cutoff radius is 5σ. Each LJ cluster is subject to the following constraint condition:24,25,48

1Nmi=1NmrIiRI2=Rg2,
(11)

where Nm = 10 is the number of LJ particles per cluster, rIi is the position of LJ particle i in cluster I, and Rg = 4.985 Å is the constant radius of gyration. The constraint condition (11) prevents extreme expansion of constituent LJ particles, thus allowing for the evaluation of intercluster forces. Here, Rg is determined such that the inner number density, defined by 3Nm/(4πRg3), matches the number density of the system. The MD simulation is performed in a microcanonical ensemble with a 2-fs time step. The average temperature of the system is T = 121 K.

FIG. 1.

The reference LJ system with 8000 clusters, each of which contains 10 LJ particles. The system size is 158 × 158 × 158 Å3 with the periodic boundary conditions in all the directions. The clusters are tinted in twelve colors for easy identification.

FIG. 1.

The reference LJ system with 8000 clusters, each of which contains 10 LJ particles. The system size is 158 × 158 × 158 Å3 with the periodic boundary conditions in all the directions. The clusters are tinted in twelve colors for easy identification.

Close modal

The instantaneous force FIJ acting from cluster J to I is evaluated by summing the interatomic forces between the LJ particles in both clusters, and it is decomposed into two components that are parallel and orthogonal to RIJ = RIRJ. Then, the distance-dependent conservative force is obtained by averaging the parallel component of FIJ as

FIJ=FIJCRIJeIJ,
(12a)
VCGRIJ=RIJFIJCRdR,
(12b)

where eIJ=RIJ/RIJ with RIJ=RIJ, FIJCRIJ is the distance-dependent magnitude of the conservative force, and VCG(RIJ) is the pair CG potential. The pair potential shown in Fig. 2 has a significantly soft repulsion core and yields a finite value even around RIJ ≈ 0, although there are larger statistical errors for smaller RIJ, where it is more improbable to find cluster pairs. Accordingly, the pair potential is fitted by the linear combination of cubic spline functions as

          VfitRIJ=a1wRIJ;Rc1,n1a2wRIJ;Rc2,n2,
(13a)
wRIJ;Rc,n=1+34n12RIJRc36n12RIJRc2+7n12RIJRc30RIJ<Rc2n22RIJRc3Rc2RIJ<Rc0RcRIJ,
(13b)

where Rc is the cutoff distance of the smoothing function and n is a dimensionless parameter. Equation (13b) with n = 1/4 reduces to the smoothing function commonly used in smoothed particle hydrodynamics (SPH).49,50 Unlike the cubic spline function used in SPH, where the first derivative of the function is necessarily zero, Eq. (13) allows for the extrapolation of the pair potential, whose first derivative is a finite negative value at RIJ = 0, as shown in Fig. 2.

FIG. 2.

The pair potentials obtained from Eq. (12) and fitted by Eq. (13) with (a1, Rc1, n1) = (5.229 × 10−20 J, 13.31 Å, 1.802 × 10−1) and (a2, Rc2, n2) = (1.670 × 10−20 J, 22.32 Å, 8.082 × 10−2).

FIG. 2.

The pair potentials obtained from Eq. (12) and fitted by Eq. (13) with (a1, Rc1, n1) = (5.229 × 10−20 J, 13.31 Å, 1.802 × 10−1) and (a2, Rc2, n2) = (1.670 × 10−20 J, 22.32 Å, 8.082 × 10−2).

Close modal

Let us define a deviation of the instantaneous force from the mean value as δFIJtFIJtFIJ. Then, Eq. (2) yields29,30

δFIJt=0tdtKIJttVIJt+δFIJQt.
(14)

By multiplying both sides with VIJT0 and taking an ensemble average, Eq. (14) can be further reformulated into a form of the Volterra integral equation of the first kind,

δFIJtVIJT0=0tdtKIJttVIJtVIJT0.
(15)

Here, we note that δFIJQtVIJT0 vanishes because δFIJQt originates from the orthogonal dynamics,19,21 and thus the correlation with the relevant variable VIJT0 is zero. For simplicity, we define AIJtδFIJtVIJT0 and BIJtVIJtVIJT0, leading to

AIJt=0tdtKIJttBIJt.
(16)

For practical implementation, it is useful to decompose each function in Eq. (16) into two components that are parallel and perpendicular to eIJ. δFIJ, δFIJQ, and VIJ can be rewritten as

δFIJ=eIJeIJTδFIJ+1eIJeIJTδFIJ=δFIJ||+δFIJ,
(17)
δFIJQ=eIJeIJTδFIJQ+1eIJeIJTδFIJQ=δFIJQ,||+δFIJQ,,
(18)
VIJ=eIJeIJTVIJ+1eIJeIJTVIJ=VIJ||+VIJ,
(19)

where the superscripts “||” and “” denote parallel and perpendicular components. Therefore, AIJ, BIJ, and KIJ can be approximated by

AIJt=δFIJ||tVIJ||0T+δFIJtVIJ0T=AIJ||teIJeIJT+AIJt1eIJeIJT,
(20)
BIJt=VIJ||tVIJ||0T+VIJtVIJ0T=BIJ||teIJeIJT+BIJt1eIJeIJT,
(21)
KIJt=βδFIJQ,||tδFIJQ,||0T+βδFIJQ,tδFIJQ,0T=KIJ||teIJeIJT+KIJt1eIJeIJT,
(22)

where the scalar quantities are given by

AIJ||t=δFIJ||tVIJ||0,
(23a)
AIJt=12δFIJtVIJ0,
(23b)
BIJ||t=VIJ||tVIJ||0,
(24a)
BIJt=12VIJtVIJ0,
(24b)
KIJ||t=βδFIJQ,||tδFIJQ,||0,
(25a)
KIJt=12βδFIJQ,tδFIJQ,0.
(25b)

Note that the orthogonal scalar quantities (23b), (24b), and (25b) involve factors of 1/2; this can be easily confirmed by taking the traces of both sides of, for instance, AIJt1eIJeIJT=δFIJtVIJ0T. Substitution of Eqs. (20)–(22) into Eq. (16) yields

AIJ||t=0tdtKIJ||ttBIJ||t,
(26a)
AIJt=0tdtKIJttBIJt.
(26b)

Both the force-velocity correlations AIJ||,t and velocity-velocity correlations BIJ||,t can be directly evaluated from the MD simulation, as shown in Fig. 3. We note that these correlations are dependent on the inter-cluster separation as well as the time, and representative correlations at RIJ = 9.5 Å are displayed in Fig. 3. Subsequently, the memory kernels KIJ||,t can be obtained through the deconvolutions of Eq. (26), as shown in Fig. 4. The velocity autocorrelation function (VACF) of the cluster COMs is also displayed, indicating that its time scale is comparable with that of the memory kernels and hence the Markovian approximation is not applicable to the present system. In addition, we compute an approximate form of KIJ(t) given by

CIJt=βδFIJ||tδFIJ||0T+βδFIJtδFIJ0T=CIJ||teIJeIJT+CIJt1eIJeIJT,
(27)

where

CIJ||t=βδFIJ||tδFIJ||0,
(28a)
CIJt=12βδFIJtδFIJ0.
(28b)

Since the force deviation is related to the fluctuating force as δFIQt=expQiLtδFI with L as the Liouville operator and Q as the orthogonal projection operator,21,CIJ(t) can be seen as a zeroth-order approximation of KIJ(t),26,29 i.e., expQiLtexpiLt, and it has been adopted as a memory kernel in some studies.24,25,28,51 However, Fig. 4 indicates the clear difference between KIJ(t) and CIJ(t) in contrast to a MD system of star polymer melts,29 where there is no significant difference between them and their time integrals. The validity of the approximation KIJtCIJt will be confirmed in Sec. IV B.

FIG. 3.

Force-velocity correlations (a) and velocity-velocity correlations (b) at RIJ = 9.5 Å.

FIG. 3.

Force-velocity correlations (a) and velocity-velocity correlations (b) at RIJ = 9.5 Å.

Close modal
FIG. 4.

Memory kernels in the radial direction (a) and perpendicular direction (b) at RIJ = 9.5 Å. KIJ||,t are obtained through the deconvolutions of Eq. (26), while CIJ||,t are computed via Eq. (28). The velocity autocorrelation function (VACF) of the cluster COMs is also displayed as a measure of time scale separation.

FIG. 4.

Memory kernels in the radial direction (a) and perpendicular direction (b) at RIJ = 9.5 Å. KIJ||,t are obtained through the deconvolutions of Eq. (26), while CIJ||,t are computed via Eq. (28). The velocity autocorrelation function (VACF) of the cluster COMs is also displayed as a measure of time scale separation.

Close modal

The distance-dependent friction coefficients in the radial and perpendicular directions are defined by

γ||,RIJ=0tcKIJ||,RIJ;tdt,
(29)

where tc is the upper limit of the time integration. Basically, a self-diffusion coefficient which characterizes long-time dynamics at the hydrodynamic spatiotemporal scale is mainly determined by the friction coefficients, while short-time dynamics originating from molecular interaction significantly depends on the functional forms of the memory kernels as well as their integrals. Figure 5 shows the friction coefficients in the radial and perpendicular directions computed from KIJ||,RIJ;t and CIJ||,RIJ;t with tc = 1 ps. Again, there are remarkable differences between the friction coefficients obtained from KIJ||, and CIJ||,. The friction coefficients are fitted by

γ||,RIJ=1αγmainRIJ+αγtailRIJ,
(30a)
γmainRIJ=Λ1RIJRf1s10RIJ<Rf10Rf1RIJ,
(30b)
γtailRIJ=ηRIJRps20RIJ<Rf20Rf2RIJ,
(30c)
α=expχ1RIJRf2s3,
(30d)

where Rp = 10 Å, s3 = 5, and the other parameters (Λ, s1, Rf1, η, s2, Rf2, and χ) for each case are listed in Table I. γRIJ based on KIJt can be well fitted by Eq. (30b) alone. The function γtail(RIJ) is responsible for the algebraically decaying tails of the friction coefficients (see the insets of Fig. 5), which would slightly influence the diffusivity of CG particles.

FIG. 5.

Friction coefficients obtained through Eq. (29) with tc = 1 ps for the radial direction (a) and the perpendicular direction (b). The memory kernels are taken as KIJ||,t and CIJ||,t. Each friction coefficient is fitted by Eq. (30) with the fitting parameters listed in Table I. The insets show global views of the fitting functions.

FIG. 5.

Friction coefficients obtained through Eq. (29) with tc = 1 ps for the radial direction (a) and the perpendicular direction (b). The memory kernels are taken as KIJ||,t and CIJ||,t. Each friction coefficient is fitted by Eq. (30) with the fitting parameters listed in Table I. The insets show global views of the fitting functions.

Close modal
TABLE I.

Fitting parameters used for fitting Eq. (30) to the friction coefficients shown in Fig. 5. γRIJ based on KIJt is fitted by Eq. (30b) alone.

Memory kernelγ (kg/s)Λ (10−12 kg/s)s1Rf1 (Å)η (10−11 kg/s)s2Rf2 (Å)χ
KIJ||,t γ||RIJ 3.168 2.8 16.7 20 19 2 × 103 
γRIJ 3.271 4.8 18.8 … … … … 
CIJ||,t γ||RIJ 1.690 2.2 16.6 20 19 3 × 103 
γRIJ 1.699 2.4 15.6 21 18 3 × 103 
Memory kernelγ (kg/s)Λ (10−12 kg/s)s1Rf1 (Å)η (10−11 kg/s)s2Rf2 (Å)χ
KIJ||,t γ||RIJ 3.168 2.8 16.7 20 19 2 × 103 
γRIJ 3.271 4.8 18.8 … … … … 
CIJ||,t γ||RIJ 1.690 2.2 16.6 20 19 3 × 103 
γRIJ 1.699 2.4 15.6 21 18 3 × 103 

To take into account the friction and fluctuating forces orthogonal to RIJ, we define the local coordinates (eIJ, eIJ1, and eIJ2) assigned to each pair of CG particles as depicted in Fig. 6. Here, the unit vectors eIJ1 and eIJ2 are orthogonal to each other and eIJ. Then, the EOM (3a) can be rewritten as

dPIdt=JIdPIJdt=JIdPIJ||dt+dPIJ1dt+dPIJ2dt,
(31)

where dPIJ||, dPIJ1, and dPIJ2 represent the increments of momenta along eIJ, eIJ1, and eIJ2, respectively. By substituting Eq. (22) into Eq. (3b), we arrive at

dPIJ||dt=FIJ0tdtKIJ||ttVIJ||t+δFIJQ,||t,
(32a)
dPIJ1dt=0tdtKIJttVIJ1t+δFIJQ,1t,
(32b)
dPIJ2dt=0tdtKIJttVIJ2t+δFIJQ,2t.
(32c)

Here, VIJ1=VIJeIJ1eIJ1 and VIJ2=VIJeIJ2eIJ2, while δFIJQ,1=δFIJQeIJ1eIJ1 and δFIJQ,2=δFIJQeIJ2eIJ2.

FIG. 6.

Local coordinates assigned to each pair of CG particles. The unit vectors eIJ1 and eIJ2 are orthogonal to each other and eIJ.

FIG. 6.

Local coordinates assigned to each pair of CG particles. The unit vectors eIJ1 and eIJ2 are orthogonal to each other and eIJ.

Close modal

For each Eqs. (32a)–(32c), we consider the corresponding Markovian equations with auxiliary variables as shown in Eq. (4), necessitating the construction of the drift and diffusion matrices. The memory kernels KIJ||RIJ;t and KIJRIJ;t are divided into distance-dependent and time-dependent parts as

KIJ||RIJ;tΦIJ||RIJθIJ||t,
(33a)
KIJRIJ;tΦIJRIJθIJt,
(33b)

where θIJ||0=θIJ0=1. The time-dependent parts are approximated by a linear combination of exponentially damped oscillators given as30 

θIJ||,t=i=1ntκiexptτicosωitφi.
(34)

Equation (34) with nt = 4 and the fitting parameters listed in Table II faithfully represent the time-dependent parts of the memory kernels, as shown in Fig. 7(a). The distance-dependent parts ΦIJ||,RIJ are obtained through Eqs. (29) and (33) as

ΦIJ||,RIJ=γ||,RIJ0tcθIJ||,tdt,
(35)

which ensures accurate representations of the integral values of the memory kernels. Then, by fitting Eq. (6b) to the memory kernels using the Nelder–Mead method,52 the corresponding elements in the drift matrices in Eq. (4) are obtained as

Aps||=0,0.5894λ,0,0.6015λ,0,0.3917λ,0,0.3750λ(kg1/2s1),Asp||=Aps||T,Ass||=10.139.3260000009.32600000000025.7821.3600000021.3600000000011.2614.1500000014.150000000004.8453.3710000003.3710(ps1),Aps=0,0.8290μ,0,0.3432μ,0,0.02270μ,0,0.4442,μ(kg1/2s1),Asp=ApsT,Ass=15.6913.8200000013.8200000000019.8625.1500000025.150000000008.08287.8700000087.870000000009.3048.1910000008.1910(ps1),
(36)

where λ=ΦIJ||RIJ and μ=ΦIJRIJ. The method for constructing drift and diffusion matrices from the fitting functions is detailed in a previously published paper.30 Practically, care must be taken for K̃IJ=KIJ/M with M being the mass of a CG particle as indicated by Eqs. (3b) and (6a), meaning that Aps||, must be replaced by Aps||,/M in accordance with Eq. (6a). Once the matrices Ass||, are obtained, then the elements in the diffusion matrices Bs||, are automatically determined through the second FDT (7).

TABLE II.

Parameters for fitting Eq. (34) to the temporal parts of the memory kernels KIJ||, shown in Fig. 7(a).

Index iκiτi (fs)ωi (ps−1)φi
θIJ||t 0.4136 197.4 7.831 0.5741 
0.4536 77.58 17.03 0.6478 
0.1672 177.7 12.98 0.4092 
0.2021 412.8 2.345 0.8017 
θIJt 0.8312 127.4 11.37 0.6040 
0.1276 100.7 23.11 0.4059 
5.134 × 10−4 247.5 87.78 4.600 × 10−2 
0.2386 215.0 6.741 0.6040 
Index iκiτi (fs)ωi (ps−1)φi
θIJ||t 0.4136 197.4 7.831 0.5741 
0.4536 77.58 17.03 0.6478 
0.1672 177.7 12.98 0.4092 
0.2021 412.8 2.345 0.8017 
θIJt 0.8312 127.4 11.37 0.6040 
0.1276 100.7 23.11 0.4059 
5.134 × 10−4 247.5 87.78 4.600 × 10−2 
0.2386 215.0 6.741 0.6040 
FIG. 7.

The temporal parts of the memory kernels fitted by Eq. (34) for KIJ||, with nt = 4 and the parameters listed in Table II (a) and for CIJ||, with nt = 3 and the parameters listed in Table III (b).

FIG. 7.

The temporal parts of the memory kernels fitted by Eq. (34) for KIJ||, with nt = 4 and the parameters listed in Table II (a) and for CIJ||, with nt = 3 and the parameters listed in Table III (b).

Close modal

Likewise, when the memory kernel KIJ(t) is approximated by CIJ(t), their time-dependent parts are fitted by Eq. (34) with nt = 3 and the fitting parameters listed in Table III, as displayed in Fig. 7(b). Here, we neglect physically unnatural offsets in the long-time region of the temporal parts of the original data, which may be due to the time variations in inter-cluster relative positions.25 Subsequently, the elements in the drift matrices in Eq. (4) are given by

Aps||=0,0.2515λ,0,0.5584λ,0,0.8242λ(kg1/2s1),Asp||=Aps||T,Ass||=4.3972.54500002.545000000019.9320.20000020.20000000011.9711.74000011.740(ps1),Aps=0,0.2951μ,0,0.8350μ,0,0.5757μ(kg1/2s1),Asp=ApsT,Ass=5.3622.85800002.858000000025.6318.28000018.2800000008.11611.34000011.340(ps1).
(37)
TABLE III.

Parameters for fitting Eq. (34) to the temporal parts of the memory kernels CIJ||, shown in Fig. 7(b).

Index iκiτi (fs)ωi (ps−1)φi
θIJ||t 0.1200 454.8 1.281 1.043 
0.3424 100.4 17.57 0.5158 
0.7544 167.1 10.10 0.5352 
θIJt 0.2293 373.0 0.9893 1.217 
0.8910 78.05 13.04 0.7768 
0.3235 246.4 10.59 0.3660 
Index iκiτi (fs)ωi (ps−1)φi
θIJ||t 0.1200 454.8 1.281 1.043 
0.3424 100.4 17.57 0.5158 
0.7544 167.1 10.10 0.5352 
θIJt 0.2293 373.0 0.9893 1.217 
0.8910 78.05 13.04 0.7768 
0.3235 246.4 10.59 0.3660 

We perform the NMDPD simulations augmented with the auxiliary variables corresponding to KIJ and CIJ. The systems consist of 8000 CG particles in the periodic boxes of 158 × 158 × 158 Å3 in accordance with the MD system. The generalized Verlet scheme53 is employed for the numerical integration with a time step of Δt = 10 fs. The local coordinates assigned to each pair of CG particles are rotated at each time step.25 More specifically, eIJ1t+Δt and eIJ2t+Δt are updated by rotating eIJ1t and eIJ2t around eIJt×eIJt+Δt, which is uniquely determined from RIJ.

Figure 8 shows the dynamic properties in the MD, KIJ-based NMDPD, and CIJ-based NMDPD systems. The KIJ-based NMDPD model can successfully reproduce the VACF of the MD system, indicating an accurate representation of the short-time dynamics. In contrast, the CIJ-based NMDPD model fails to reproduce the VACF of the MD system, involving a noticeable undershoot after the exponential decay. This undershoot may be partly attributed to the undershoot of the memory kernel in the radial direction CIJ||t, as displayed in Fig. 4(a). Moreover, the mean square displacement (MSD) characterizing the long-time dynamics is accurately reproduced by the KIJ-based NMDPD model as shown in Fig. 8(b), while the CIJ-based NMDPD model underestimates the long-time diffusivity. Specifically, the self-diffusion coefficient in the KIJ-based NMDPD system is 4.20 × 10−10 m2/s, in excellent agreement with that in the MD system 4.23 × 10−10 m2/s, although that in the CIJ-based NMDPD system is 3.23 × 10−10 m2/s. These results indicate that the NMDPD model based on the exact memory kernels KIJ can accurately reproduce both the short- and long-time dynamics, while the approximate memory kernels CIJ are inappropriate for the present liquid system.

FIG. 8.

Dynamic properties in the MD, KIJ-based NMDPD, and CIJ-based NMDPD systems. The MD results are associated with the center-of-mass positions and velocities of the clusters. (a) Velocity autocorrelation functions (VACFs). (b) Mean square displacements (MSDs). The self-diffusion coefficients in the MD, KIJ-based NMDPD, and CIJ-based NMDPD systems are 4.23 × 10−10 m2/s, 4.20 × 10−10 m2/s, and 3.23 × 10−10 m2/s, respectively.

FIG. 8.

Dynamic properties in the MD, KIJ-based NMDPD, and CIJ-based NMDPD systems. The MD results are associated with the center-of-mass positions and velocities of the clusters. (a) Velocity autocorrelation functions (VACFs). (b) Mean square displacements (MSDs). The self-diffusion coefficients in the MD, KIJ-based NMDPD, and CIJ-based NMDPD systems are 4.23 × 10−10 m2/s, 4.20 × 10−10 m2/s, and 3.23 × 10−10 m2/s, respectively.

Close modal

Figure 9 shows the RDFs in the MD, KIJ-based NMDPD, and CIJ-based NMDPD systems. Both the CG models provide identical RDFs because static properties of the CG system depend primarily on the conservative force field.7 The RDFs in the CG systems differ significantly from that in the MD system. In addition, the pressures in the CG systems are 15.3 MPa, which is significantly smaller than that in the MD system, 46.4 MPa. These discrepancies imply that the many-body effect among CG particles should be effectively incorporated into the pair potential to accurately reproduce the static properties in the dense liquid system. As a matter of fact, the RDF in the MD system indicates significant overlapping among the clusters, highlighting the non-negligible many-body effect.

FIG. 9.

Radial distribution functions (RDFs) in the MD, KIJ-based NMDPD, and CIJ-based NMDPD systems. The MD results are associated with the center-of-mass positions of the clusters.

FIG. 9.

Radial distribution functions (RDFs) in the MD, KIJ-based NMDPD, and CIJ-based NMDPD systems. The MD results are associated with the center-of-mass positions of the clusters.

Close modal

To improve CG representations of the static properties, the pair CG potential is modified through the IBI described in Sec. II C. The MZ-guided pair potential shown in Fig. 2 is adopted as an initial guess, and the cutoff distance is set to be 22.32 Å according to that of the MZ-guided potential. The IBI procedure is implemented in the VOTCA package,54 where a series of CGMD simulations with a 2-fs time step is performed using GROMACS.55 In each iteration, the system is initially equilibrated for 60 ps followed by a production run of 440 ps.

Figure 10(a) depicts the original and modified pair potentials, where the modified one has a softer and longer-range repulsion core than the original one. The resultant RDF obtained from a KIJ-based NMDPD simulation employing the modified potential is displayed in Fig. 10(b), being an accurate representation of the MD counterpart. Moreover, the pressure in the KIJ-IBI-based NMDPD system is 46.3 MPa, showing excellent agreement with the MD pressure, 46.4 MPa. Most importantly, the dynamics properties characterized by the VACF and MSD are little influenced by the IBI correction of the pair potential, as shown in Figs. 10(c) and 10(d). Actually, the self-diffusion coefficient in the KIJ-IBI-based NMDPD system is 4.82 × 10−10 m2/s, still in good agreement with the MD value 4.23 × 10−10 m2/s, although it is slightly larger than that of the KIJ-based NMDPD model 4.20 × 10−10 m2/s because of less caging effects of the modified potential. These results indicate that the dynamic properties are primarily determined by the memory kernels, with negligible dependence on the form of the pair CG potential. Therefore, the MZ-guided construction of memory kernels with the subsequent IBI correction of a pair CG potential is a powerful technique to accurately represent both the static and dynamic properties of microscopic systems.

FIG. 10.

Modifying the pair CG potential through the IBI. The initial guess for the potential based on the MZ-guided one (solid line) and the modified potential (dashed line) (a). RDFs (b), VACFs (c), and MSDs (d) in the MD (solid line), KIJ-based NMDPD (dashed-dotted line), and KIJ-IBI-based NMDPD (dashed line) systems. The self-diffusion coefficients in the MD, KIJ-based NMDPD, and KIJ-IBI-based NMDPD systems are 4.23 × 10−10 m2/s, 4.20 × 10−10 m2/s, and 4.82 × 10−10 m2/s, respectively.

FIG. 10.

Modifying the pair CG potential through the IBI. The initial guess for the potential based on the MZ-guided one (solid line) and the modified potential (dashed line) (a). RDFs (b), VACFs (c), and MSDs (d) in the MD (solid line), KIJ-based NMDPD (dashed-dotted line), and KIJ-IBI-based NMDPD (dashed line) systems. The self-diffusion coefficients in the MD, KIJ-based NMDPD, and KIJ-IBI-based NMDPD systems are 4.23 × 10−10 m2/s, 4.20 × 10−10 m2/s, and 4.82 × 10−10 m2/s, respectively.

Close modal

Although the present study has investigated moderate coarse-graining of 10 LJ particles into a single cluster, the same methodology can be applied, in principle, to more aggressive coarse-graining, i.e., 100 LJ particles in a cluster. However, in that case, the many-body effect would be more pronounced27,47 and it would be difficult to separate the memory effect from the many-body effect. Hence, we plan a future study where we can systematically investigate the interplay of these two important effects in aggressively coarse-grained systems.

We have proposed a non-Markovian coarse-graining method, which takes advantage of the Mori–Zwanzig (MZ) formalism for extracting memory kernels together with the iterative Boltzmann inversion (IBI) for improving the representation of static properties. We have demonstrated that the MZ-guided non-Markovian dissipative particle dynamics model exploiting auxiliary variables (allowing for the replacement of a non-Markovian equation with a Markovian one in a higher dimensional space) can successfully reproduce the dynamic properties such as the velocity autocorrelation function and mean square displacement of the reference Lennard–Jones system, as long as the pairwise memory kernels are appropriately evaluated. More specifically, we have found that the memory kernels should be constructed based on the Volterra integral equation using force-velocity and velocity-velocity correlations, rather than based on the time-correlations of instantaneous force deviations from the mean force, which correspond to a zeroth-order approximation of the exact memory kernels. Furthermore, the IBI correction of the pair potential ensures the accurate representation of the radial distribution function and pressure, while it has little influence on the dynamic properties that depend primarily on the memory kernels. The present study indicates that constructing memory kernels based on the MZ formalism followed by the IBI correction of a pair potential is a promising method for obtaining coarse-grained models that can successfully represent both the static and dynamic properties of underlying microscopic systems.

This work was partly supported by JSPS KAKENHI Grant No. 16K18009. This work was also supported by the U.S. Army Research Laboratory and was accomplished under Cooperative Agreement No. W911NF-12-2-0023.

1.
C.
Peter
and
K.
Kremer
,
Soft Matter
5
,
4357
(
2009
).
2.
M. G.
Saunders
and
G. A.
Voth
,
Annu. Rev. Biophys.
42
,
73
(
2013
).
3.
M. P.
Allen
and
D. J.
Tildesley
,
Computer Simulation of Liquids
(
Clarendon Press
,
New York
,
1989
).
4.
D.
Frenkel
and
B.
Smit
,
Understanding Molecular Simulation: From Algorithms to Applications
, 2nd ed. (
Academic Press
,
San Diego
,
2001
).
5.
W. G.
Noid
,
J. Chem. Phys.
139
,
90901
(
2013
).
6.
P.
Español
and
P. B.
Warren
,
J. Chem. Phys.
146
,
150901
(
2017
).
7.
S.
Izvekov
and
G. A.
Voth
,
J. Phys. Chem. B
109
,
2469
(
2005
).
8.
S.
Izvekov
and
G. A.
Voth
,
J. Chem. Phys.
123
,
134105
(
2005
).
9.
D.
Reith
,
M.
Pütz
, and
F.
Müller-Plathe
,
J. Comput. Chem.
24
,
1624
(
2003
).
10.
V.
Agrawal
,
G.
Arya
, and
J.
Oswald
,
Macromolecules
47
,
3378
(
2014
).
11.
A. P.
Lyubartsev
and
A.
Laaksonen
,
Phys. Rev. E
52
,
3730
(
1995
).
12.
T.
Murtola
,
E.
Falck
,
M.
Karttunen
, and
I.
Vattulainen
,
J. Chem. Phys.
126
,
075101
(
2007
).
13.
P.
Español
and
I.
Zúñiga
,
Phys. Chem. Chem. Phys.
13
,
10538
(
2011
).
14.
H.
Lei
,
X.
Yang
,
Z.
Li
, and
G. E.
Karniadakis
,
J. Comput. Phys.
330
,
571
(
2017
).
15.
M. S.
Shell
,
J. Chem. Phys.
129
,
144108
(
2008
).
16.
A.
Chaimovich
and
M. S.
Shell
,
J. Chem. Phys.
134
,
94112
(
2011
).
17.
S. J.
Marrink
,
A. H.
de Vries
, and
A. E.
Mark
,
J. Phys. Chem. B
108
,
750
(
2004
).
18.
P. K.
Depa
and
J. K.
Maranas
,
J. Chem. Phys.
123
,
94901
(
2005
).
19.
R.
Zwanzig
,
Nonequilibrium Statistical Mechanics
(
Oxford University Press
,
New York
,
2001
).
20.
H.
Mori
,
Prog. Theor. Phys.
33
,
423
(
1965
).
21.
T.
Kinjo
and
S.
Hyodo
,
Phys. Rev. E
75
,
051109
(
2007
).
22.
E.
Darve
,
J.
Solomon
, and
A.
Kia
,
Proc. Natl. Acad. Sci. U. S. A.
106
,
10884
(
2009
).
23.
C.
Hijón
,
P.
Español
,
E.
Vanden-Eijnden
, and
R.
Delgado-Buscalioni
,
Faraday Discuss.
144
,
301
(
2010
).
24.
H.
Lei
,
B.
Caswell
, and
G. E.
Karniadakis
,
Phys. Rev. E
81
,
026704
(
2010
).
25.
Y.
Yoshimoto
,
I.
Kinefuchi
,
T.
Mima
,
A.
Fukushima
,
T.
Tokumasu
, and
S.
Takagi
,
Phys. Rev. E
88
,
043305
(
2013
).
26.
S.
Izvekov
,
J. Chem. Phys.
138
,
134106
(
2013
).
27.
S.
Izvekov
and
B. M.
Rice
,
J. Chem. Phys.
140
,
104104
(
2014
).
28.
Z.
Li
,
X.
Bian
,
B.
Caswell
, and
G. E.
Karniadakis
,
Soft Matter
10
,
8659
(
2014
).
29.
Z.
Li
,
X.
Bian
,
X.
Li
, and
G. E.
Karniadakis
,
J. Chem. Phys.
143
,
243128
(
2015
).
30.
Z.
Li
,
H. S.
Lee
,
E.
Darve
, and
G. E.
Karniadakis
,
J. Chem. Phys.
146
,
14104
(
2017
).
31.
P. J.
Hoogerbrugge
and
J. M. V. A.
Koelman
,
Europhys. Lett.
19
,
155
(
1992
).
32.
P.
Español
and
P. B.
Warren
,
Europhys. Lett.
30
,
191
(
1995
).
33.
R. D.
Groot
and
P. B.
Warren
,
J. Chem. Phys.
107
,
4423
(
1997
).
34.
M.
Ceriotti
,
G.
Bussi
, and
M.
Parrinello
,
J. Chem. Theory Comput.
6
,
1170
(
2010
).
35.
A. D.
Baczewski
and
S. D.
Bond
,
J. Chem. Phys.
139
,
044107
(
2013
).
36.
H.
Lei
,
N. A.
Baker
, and
X.
Li
,
Proc. Natl. Acad. Sci. U. S. A.
113
,
14183
(
2016
).
37.
T.
Koide
and
T.
Kodama
,
Phys. Rev. E
78
,
051107
(
2008
).
38.
R.
Zwanzig
,
J. Chem. Phys.
33
,
1338
(
1960
).
39.
R.
Kubo
,
Rep. Prog. Phys.
29
,
255
(
1966
).
40.
P.
Español
,
Phys. Rev. E
53
,
1572
(
1996
).
41.
X.
Li
,
Int. J. Numer. Methods Eng.
83
,
986
(
2010
).
42.
M.
Berkowitz
,
J. D.
Morgan
, and
J. A.
McCammon
,
J. Chem. Phys.
78
,
3256
(
1983
).
43.
G. H.
Golub
and
C. F.
Van Loan
,
Matrix Computations
, 4th ed. (
Johns Hopkins University Press
,
Baltimore
,
2012
).
44.
G. W.
Ford
,
J. T.
Lewis
, and
R. F.
O’Connell
,
Phys. Rev. A
37
,
4419
(
1988
).
45.
H. I.
Ingólfsson
,
C. A.
Lopez
,
J. J.
Uusitalo
,
D. H.
de Jong
,
S. M.
Gopal
,
X.
Periole
, and
S. J.
Marrink
,
Wiley Interdiscip. Rev.: Comput. Mol. Sci.
4
,
225
(
2014
).
46.
R. L.
Henderson
,
Phys. Lett. A
49
,
197
(
1974
).
47.
Z.
Li
,
X.
Bian
,
X.
Yang
, and
G. E.
Karniadakis
,
J. Chem. Phys.
145
,
044102
(
2016
).
48.
T.
Kinjo
and
S.
Hyodo
,
Mol. Simul.
33
,
417
(
2007
).
49.
M. B.
Liu
,
G. R.
Liu
, and
K. Y.
Lam
,
J. Comput. Appl. Math.
155
,
263
(
2003
).
50.
M.
Liu
,
P.
Meakin
, and
H.
Huang
,
Phys. Fluids
18
,
17101
(
2006
).
51.
L.
Gao
and
W.
Fang
,
J. Chem. Phys.
135
,
184101
(
2011
).
52.
W. H.
Press
,
S. A.
Teukolsky
,
W. T.
Vetterling
, and
B. P.
Flannery
,
Numerical Recipes: The Art of Scientific Computing
, 3rd ed. (
Cambridge University Press
,
New York
,
2007
).
53.
T.
Schlick
,
Molecular Modeling and Simulation: An Interdisciplinary Guide
, 2nd ed. (
Springer-Verlag
,
New York
,
2010
).
54.
V.
Rühle
,
C.
Junghans
,
A.
Lukyanov
,
K.
Kremer
, and
D.
Andrienko
,
J. Chem. Theory Comput.
5
,
3211
(
2009
).
55.
D.
Van Der Spoel
,
E.
Lindahl
,
B.
Hess
,
G.
Groenhof
,
A. E.
Mark
, and
H. J. C.
Berendsen
,
J. Comput. Chem.
26
,
1701
(
2005
).