Mori-Zwanzig theory for dissipative forces in coarse-grained dynamics in the Markov limit

We derive alternative Markov approximations for the projected (stochastic) force and memory function in the coarse-grained (CG) generalized Langevin equation, which describes the time evolution of the center-of-mass coordinates of clusters of particles in the microscopic ensemble. This is done with the aid of the Mori-Zwanzig projection operator method based on the recently introduced projection operator [S. Izvekov, J. Chem. Phys. 138, 134106 (2013)]. The derivation exploits the “generalized additive fluctuating force” representation to which the projected force reduces in the adopted projection operator formalism. For the projected force, we present a first-order time expansion which correctly extends the static fluctuating force ansatz with the terms necessary to maintain the required orthogonality of the projected dynamics in the Markov limit to the space of CG phase variables. The approximant of the memory function correctly accounts for the momentum dependence in the lowest (second) order and indicates that such a dependence may be important in the CG dynamics approaching the Markov limit. In the case of CG dynamics with a weak dependence of the memory effects on the particle momenta, the expression for the memory function presented in this work is applicable to non-Markov systems. The approximations are formulated in a propagator-free form allowing their efficient evaluation from the microscopic data sampled by standard molecular dynamics simulations. A numerical application is presented for a molecular liquid (nitromethane). With our formalism we do not observe the “plateau-value problem” if the friction tensors for dissipative particle dynamics (DPD) are computed using the Green-Kubo relation. Our formalism provides a consistent bottom-up route for hierarchical parametrization of DPD models from atomistic simulations.

Figure 1

Center-of-mass RDFs, $g_{\mathrm{CoM}}\left(R_{IJ}\right)$; corresponding running coordination number, $n_{\mathrm{CoM}}\left(R_{IJ}\right)=4\pi \rho {\int }_0^{R_{IJ}}d{R}^{\prime }{R^{\prime }}^2{g}_{\mathrm{CoM}}\left(R^{\prime }\right)$, where $\rho$ is the particle number density; and pair contributions $u^{\mathrm{PMF}}\left(R_{IJ}\right)={\int }_{R_{IJ}}^{R_{\mathrm{cut}}^C}d{R}^{\prime }{f}^C\left(R^{\prime },\alpha \right)$ into the potential of mean force, $W^{\mathrm{PMF}}\left(R\right)\approx {\sum }_{I<J}{u}^{\mathrm{PMF}}\left(R_{IJ}\right)$, for the $N_{\mathrm{CG}}=1$ (a) and $N_{\mathrm{CG}}=2$ (b) models of liquid nitromethane from microscopic simulation (“Atm.” label). The microscopic $g_{\mathrm{CoM}}\left(R_{IJ}\right)$ and $n_{\mathrm{CoM}}\left(R_{IJ}\right)$ are compared to the respective functions obtained from frictionless MS-CG and MS-CG DPD simulations (“MS-CG, DPD” label). The insets depict the mapping for the respective CG model.

Figure 2

Correlations in the microscopic ensemble: normalized velocity autocorrelation function $c_{VV}\left(t\right)=\langle V_I\left(t\right)·V_I\left(0\right)\rangle /\langle V_I^2\left(0\right)\rangle$; two-particle velocity-velocity correlation function $c_{VV}\left(t,R_{IJ}\right)=\langle V_I\left(t\right)·V_J\left(0\right)\rangle /\langle V_J^2\left(0\right)\rangle$ at $R_{IJ}=0.5$ nm; normalized force-force correlation functions $c_{F\mathrm{\Delta }F}\left(t\right)=\langle F_I\left(t\right)·\mathrm{\Delta }{F}_I\left(0\right)\rangle /\langle F_I\left(0\right)·\mathrm{\Delta }{F}_I\left(0\right)\rangle$, $c_{\mathrm{\Delta }F\mathrm{\Delta }F}\left(t\right)=\langle \mathrm{\Delta }{F}_I\left(t\right)·\mathrm{\Delta }{F}_I\left(0\right)\rangle /\langle \mathrm{\Delta }{F}_I^2\left(0\right)\rangle$; and their respective running time integrals ${\overline{C}}_{F\mathrm{\Delta }F}\left(t\right)={\int }_0^{t}d\tau c_{F\mathrm{\Delta }F}\left(\tau \right)$ and ${\overline{C}}_{\mathrm{\Delta }F\mathrm{\Delta }F}\left(t\right)={\int }_0^{t}d\tau c_{\mathrm{\Delta }F\mathrm{\Delta }F}\left(\tau \right)$ (both quantities are plotted in units of 0.1 ps) for coarse-graining resolutions $N_{\mathrm{CG}}=1$ (a) and $N_{\mathrm{CG}}=2$ (b). The running time integrals describe (respectively) convergence of the cumulative frictions Eqs. (50), and (51) as functions of $\overline{t}$ because ${\overline{C}}_{F\mathrm{\Delta }F}\left(\overline{t}\right)\propto {\overline{\gamma }}_{II}^{\left|\right|}+2{\overline{\gamma }}_{II}^{\perp }$, ${\overline{C}}_{\mathrm{\Delta }F\mathrm{\Delta }F}\left(\overline{t}\right)\propto {\overline{\gamma }}_{II}^{\mathrm{\Delta }\left|\right|}+2{\overline{\gamma }}_{II}^{\mathrm{\Delta }\perp }$, respectively, and because of the rule given in Eq. (52).

Figure 3

Power spectral densities of particle velocity, ${\tilde{c}}_{VV}\left(\lambda \right)$; projected force, ${\tilde{c}}_{\delta F\delta F}\left(\lambda \right)={\tilde{c}}_{F\mathrm{\Delta }F}\left(\lambda \right)$; and projected force in the static fluctuation approximation [Eq. (26)], ${\tilde{c}}_{\mathrm{\Delta }F\mathrm{\Delta }F}\left(\lambda \right)$. Power spectral densities are obtained by Fourier transform $\tilde{c}\left(\lambda \right)={\int }_0^{\infty }dtcos\left(2\pi \lambda t\right)c\left(t\right)$ of corresponding microscopic correlation functions shown in Fig. 2.

Figure 4

Parallel ${\overline{\gamma }}^{\left|\right|}\left(R_{IJ}\right)$ (a) and normal ${\overline{\gamma }}^{\perp }\left(R_{IJ}\right)$ (b) frictions in the Markov approximation for the $N_{\mathrm{CG}}=1$ coarse graining using Eq. (50) (lines labeled by $F\mathrm{\Delta }F$) and static fluctuation formula Eq. (51) (lines labeled by $\mathrm{\Delta }F\mathrm{\Delta }F$). The frictions from Ref. [27] (lines labeled by $VV$), which are obtained using velocity-velocity and force-velocity correlations, are included for a comparison. The $F\mathrm{\Delta }F$ and $\mathrm{\Delta }F\mathrm{\Delta }F$ curves were obtained by applying the linear low-pass Gaussian filter [27, 28] $\overline{\gamma }\left(R\right)\approx \mathrm{\Delta }R{\sum }_{l=1}^{N_{\mathrm{bin}}}\overline{\gamma }\left(R_l\right)G\left(R_l-R\right)$ to the respective $R_{IJ}$ histograms of bin size $\mathrm{\Delta }R=0.01$ nm, shown in cyan. We used the filter $G\left(R\right)=1/\left(\sqrt{2\pi }\sigma \right)e^{-R^2/\phantom{R^{2}2{\sigma }^2}2{\sigma }^2}H\left(1-\left|R\right|/\phantom{\left|R\right|{\sigma }^{\prime }}{\sigma }^{\prime }\right)$, where the $H\left(t\right)$ is the Heaviside function and $\sigma$ is the bandwidth parameter which controls the smoothness of the convoluted frictions. For a chosen value of $\sigma$ (we used $\sigma =8\mathrm{\Delta }R$), the value of ${\sigma }^{\prime }$ is selected to have $G\left({\sigma }^{\prime }\right)=\varepsilon$, where the tolerance $\varepsilon ={10}^{-3}$.

Figure 5

Same data as in Fig. 4 for the $N_{\mathrm{CG}}=2$ coarse graining.

Figure 6

Comparison of normalized velocity autocorrelation functions $c_{VV}\left(t\right)$ and their power spectra ${\tilde{c}}_{VV}\left(\lambda \right)$ from microscopic simulation (“Atm.” label), MS-CG simulations in which dissipative forces are absent (“MS-CG” label), and MS-CG DPD simulations with frictions ${\overline{\gamma }}^{\alpha }\left(R_{IJ}\right)$ in approximation equation (50) (“DPD” label). For $N_{\mathrm{CG}}=2$, the MS-CG DPD results using frictions ${\overline{\gamma }}^{VV\alpha }\left(R_{IJ}\right)$ from Ref. [27] are shown (lines labeled by “DPD, $VV$”). Definitions of $c_{VV}\left(t\right)$ and ${\tilde{c}}_{VV}\left(\lambda \right)$ are given in captions (respectively) of Figs. 2 and 3.