Formal derivation of dissipative particle dynamics from first principles

Formal derivation of dissipative particle dynamics from first principles David Cubero*and Sophia N. Yaliraki† Department of Chemistry, South Kensington Campus, Imperial College London, London SW7 2AZ, United Kingdom 共Received 13 January 2005; revised manuscript received 20 May 2005; published 13 September 2005兲 We show that the Markovian approximation assumed in current particle-based coarse-grained techniques, like dissipative particle dynamics, is unreliable in situations in which sound plays an important role. As an example we solve analytically and numerically the dynamics of coarse-grained harmonic systems by using first principle methods, showing the presence of long-lived memory kernels. This effect raises questions about the connection of these approaches at their current form to molecular dynamics. DOI: 10.1103/PhysRevE.72.032101 PACS number共s兲: 05.40.⫺a, 02.70.Ns The processes underlying the behavior of large physical and biological macromolecular assemblies span a range of length and time scales that often require coupling of microscopic interactions with continuum-level properties that cannot be presently bridged by techniques like molecular dynamics. This has recently led to revisiting existing or developing new coarse-graining schemes. Particle-based methods 关1–5兴based on the dissipative particle dynamics 共DPD兲methodology have been used extensively in recent years as models for addressing such systems at the mesoscopic level 关6,7兴. One of the appeals of DPD is the capability of reproducing hydrodynamic behavior through the conservative, dissipative, and stochastic pairwise forces that are mass and momentum conserving. Manifestations of the continuum properties at the mesoscale have been extensively investigated. Nevertheless, this methodology is still lacking a clear connection with the underlying microscopic level. In Ref. 关4兴a derivation of DPD from the microscopic level was proposed. However, this derivation was based on a series of ad hoc assumptions about the form of the conservative, dissipative, and noise forces. In particular, it was assumed that the fluctuating forces are Markovian. A more fundamental approach was presented in Ref. 关5兴, in which the equation of motion for the dissipative particles was obtained by means of projection operators. As a result, the dissipative forces present memory kernels and the fluctuating forces are coupled to the dissipative terms via the fluctuationdissipation theorem. However, the DPD-like equations of motion are only obtained when it is assumed that the fluctuating forces are Markovian and the memory kernels collapse into the familiar dissipative forms, functions only of the instantaneous values of the coarse-grained variables. In Ref. 关5兴, this was justified by arguing that the coarse-grained variables are mesoscopic, and therefore slow when compared with the microscopic variables. Since this was a central requirement for obtaining the equations of motion, it merits verification. In this work we show that the Markovian approximation is not applicable in situations in which sound plays an important role. We first study the general case and the conditions for the failure of the Markovian approximation due to sound propagation. Then, as an example, we obtain the memory kernels of the coarse-grained dynamics of cubic harmonic lattices in several dimensions by using first principle methods. Harmonic lattices are paradigmatic for the kind of memory effect discussed here because elastic wave propagation is the only dissipation mechanism in these systems. We have recently used projector operators to formally derive the stochastic equations of motion for inhomogeneous harmonic lattices 关8兴. We use these results as a starting point in order to evaluate here for the first time the required memory kernels, analytically in one dimension and numerically in two and three dimensions. Consider a coarse-grained harmonic system by partitioning the volume into mesoscopic regions k. We define the coarse-grained particles by choosing the set of all center of mass coordinates Xkand velocities Vkof all particles inside each region kas the relevant variables. Then, by using Mori’s projection operator theory 关9兴, the following equation of motion for the coarse-grained variables is obtained 关8兴: Mk dVk共t兲 dt =−兺 l 冉 ⌳klXl共t兲+ 冕 0 t d ␶ ⌽kl共 ␶ 兲Vl共t− ␶ 兲 冊 +Rk共t兲, 共1兲 where Mkis the mass of each coarse-grained particle; ⌳kl is a conservative force matrix, ⌽kl共t兲is a memory kernel defining dissipation; and Rk共t兲is a random force satisfying 具Rk共t兲典=0 and the fluctuation-dissipation theorem 具Rk共t兲Rl共0兲典=kBT⌽kl共t兲, where kBis the Boltzmann constant and Tis the temperature. In these last equations the averages are taken with an initial distribution close to thermal equilibrium. Total momentum conservation implies that 兺l⌽kl共t兲 =兺lRl共t兲=0, so that if we decompose each random force as a sum of independent variables Rk共t兲=兺l⫽k ␰ kl, we can rewrite 共1兲in a form that resembles the equation of motion of a DPD system, *Present address: Física Teórica, Universidad de Sevilla, Apartado de Correos 1065, Sevilla 41080, Spain. †Electronic address: [email protected] PHYSICAL REVIEW E 72, 032101 共2005兲 1539-3755/2005/72共3兲/032101共4兲/$23.00 ©2005 The American Physical Society032101-1 Mk dVk共t兲 dt =−兺 l⫽k 再 ⌳kl关Xl共t兲−Xk共t兲兴 − ␰ kl共t兲 + 冕 0 t d ␶ ⌽kl共t− ␶ 兲关Vl共 ␶ 兲−Vk共 ␶ 兲兴 冎 .共2兲 However, only if the kernels are Dirac delta functions we obtain the DPD equation of motion. This would imply that the set of coarse-grained variables 兵Xk,Vk其is Markovian. To show that this is not the case when sound propagation is considered, let us focus on a simple one-dimensional solid. At the macroscopic level, the displacement field u共x,t兲关10兴 is defined as the deviation from equilibrium of the particles whose equilibrium site is x. If we neglect the microscopic fluctuations, we can express the mesoscopic variables as Xk共t兲=具Xk典eq+兰lk l共k+1兲dx u共x,t兲and Vk共t兲 =兰lk l共k+1兲dx ⳵ u共x,t兲/ ⳵ t, where lis the size of the coarsegrained regions and 具¯典eq denotes equilibrium averages. In the simplest case the displacement field verifies the wave equation ⳵ 2u/ ⳵ t2=v2 ⳵ 2u/ ⳵ x2, where vis the speed of sound in the material. The general solution to this equation is given by u共x,t兲=f共x−vt兲+g共x+vt兲, with a front f共x兲and g共x兲traveling to right and left, respectively. Note that the ratio ␶ s =l/vdefines a mesoscopic time scale because lis mesoscopic. Now consider that initially all coarse-grained particles are in their equilibrium positions with zero velocity except one, which we label as k=0. In Fig. 1共a兲we have plotted the corresponding fields u共x,0兲共solid line兲and 共 ⳵ u/ ⳵ t兲共x,0兲 共shaded area兲. Assume that the fluctuations are small enough so that the fronts are not destroyed. Then, Figs. 1共b兲and 1共c兲 show the profiles after a time ␶ s/2 and ␶ s, respectively. Now imagine that at the intermediate time t= ␶ s/2 we change the velocity field ⳵ u/ ⳵ tso that the set of coarse-grained variables 兵Xk,Vk其remains the same, as shown in Fig. 1共b兲共dashed line兲. Figure 1共d兲shows the profiles after a time ␶ s/2. We can clearly see that the values of the coarse-grained variables are very different from the former case shown in Fig. 1共c兲, even though they had the same values at t= ␶ s/2. Therefore, the knowledge of the values at a previous instant is not sufficient. Instead, we need to know the complete profiles of u and ⳵ u/ ⳵ tat a given time or, equivalently, the complete previous history of 兵Xk,Vk其to specify the present state. The arising memory kernels are the consequence of the perturbations traveling inside the coarse-grained particles by sound propagation. This phenomenon is due to the space discretization introduced by the coarse graining and we expect it to be present even if we consider fluctuations, higher dimensions, or take into account anharmonic effects, as long as the mesoscopic time scale associated with sound propagation is comparable with the relevant time scale of the mesoscopic particles. If ␶ sis much shorter than other characteristic times associated with distinct dissipation mechanisms present in the system, then the memory kernels due to sound propagation would collapse in the longer time scales and the previous discussion no longer applies. However, in situations in which sound propagation plays an important role, the coarse-grained variables must change in the same scale given by ␶ sin order to account properly for it, which implies the effect described above. As an example, let us consider a d-dimensional cubic harmonic lattice with only first neighbor interactions. Units are defined so that the mass of each oscillator, the force constants, and the equilibrium spacings are set to unity. We coarse-grain the system by grouping the oscillators into cubic clumps, each containing ndoscillators. It was shown in Ref. 关8兴that the conservative force ⌳kl is the inverse of the matrix 具XkXl典eq/kBT. In harmonic systems, the mesoscopic time scale that characterizes the dynamics of the coarsegrained particles is given by wave propagation only, which in our units is ␶ s=n. Therefore, by changing to t*=t/n, with Vk *=dXk/dt*=nVk, Eq. 共1兲becomes dVk * dt*=−兺 l 冉 ⌳kl *Xl+ 冕 0 t* d ␶ ⌽kl *共 ␶ 兲Vl *共t*− ␶ 兲 冊 +Rk *,共3兲 where ⌳kl *=n2−d⌳kl,⌽kl *=n2−d⌽kl, and Rk *=n2−dRk. Since we have set all internal parameters to unity, we expect that all variables in Eq. 共3兲, and, in particular, ⌳kl *and ⌽kl *, are on the order of O共1兲, even in the mesoscopic limit n→⬁. The analytic and simulation results we present later confirm this scaling in all three relevant physical dimensions, i.e., d=1, 2, and 3. To continue analytically, we restrict ourselves temporarily to the one-dimensional 共1-D兲case. Let us consider a finite chain of Noscillators with fixed ends and then take the thermodynamic limit 共N→⬁兲at the end of the calculations. The force matrix for the oscillators in this system is given by Aij=2 ␦ ij−1共 ␦ i,j+1+ ␦ i,j−1兲. Using the analytical expression for the inverse matrix Aij −1=min共i,j兲关1−max共i,j兲/共N+1兲兴 关11兴, we obtain 具XkXl典eq kBT=−n2−1 6n ␦ kl + 冉 关1+2min共k,l兲n−n兴 2 冊 ⫻ 冉 1−1+2max共k,l兲n−n 2共N+1兲 冊 .共4兲 Next, we invert this matrix in two steps. First we multiply it FIG. 1. Time evolution of the profiles of u共x,t兲共solid line兲and ⳵ u/ ⳵ t共x,t兲共shaded region兲in a one-dimensional system: 共a兲initial condition, 共b兲after a time ␶ s/2, 共c兲after a time ␶ s,共d兲after a time evolution ␶ s/2 from 共b兲, but with a different profile, ⳵ u/ ⳵ t, given by the dashed line in 共b兲. BRIEF REPORTS PHYSICAL REVIEW E 72, 032101 共2005兲 032101-2 by A, which leads to a tri-diagonal matrix Skl ⬅兺l⬘具XkXl⬘典eqAl⬘l/kBT=a ␦ ij+b共 ␦ i,j+1+ ␦ i,j−1兲, with a=1/3n +2n/3 and b=共n2−1兲/6n. Strictly speaking, the first and last columns of Sdo not follow that tri-diagonal form, but this is irrelevant in the thermodynamic limit. Now we invert the matrix Sby using the form Skl −1=C ␯ 兩k−l兩, which leads to C =1/冑a2−4b2and ␯ =−(a/b−冑共a/b兲2−4)/2. The fact that a/b⬎2 for any nguarantees that 兩 ␯ 兩⬍1 and consequently S−1 vanishes at the boundaries. The force matrix is then given by ⌳=A⫻S−1, resulting in ⌳kl =冑3 2+n2 冉 6n2 n2−1−2冑3n冑2+n2 n2−1 ␦ kl 冊 ⫻共−1兲兩k−l兩 冉 1+2n2−冑3n冑2+n2 n2−1 冊 兩k−l兩 .共5兲 In the mesoscopic limit 共n→⬁兲, it reduces to 关14兴 ⌳kl ⬃2冑3 n共3−冑3 ␦ kl兲共−1兲兩k−l兩共2−冑3兲兩k−l兩.共6兲 Note that since 兩2−冑3兩⬍1, Eq. 共6兲predicts an exponential decay of the conservative force, in contrast with the standard linear decay used in DPD simulations 关2兴. Let us now turn our attention to the memory kernels. The autocorrelation matrix Ckl *共t*兲=具Vk *共t*兲Vl *共0兲典eq satisfies the same equation as Vk *, Eq. 共3兲, but without the random force 关9兴. Taking the Laplace transforms ␾ 共s兲 =兰0 ⬁dt*exp共−st*兲⌽*共t*兲,c共s兲=兰0 ⬁dt*exp共−st*兲C*共t*兲,wearrive at 兺 l⬘ ␾ kl⬘共s兲cl⬘l共s兲=Ckl *共0兲−sckl共s兲−兺 l⬘ ⌳kl⬘ * scl⬘l共s兲,共7兲 where Ckl *共0兲= ␦ kln2kBT/Ml. The autocorrelation matrix can also be readily calculated in the mesoscopic limit using the macroscopic wave equation 关8兴, leading to Ckl *共t*兲=Cll *共0兲共1+ ␦ kl兲g共t*−兩k−l兩兲/2, 共8兲 where g共x兲=共1−兩x兩兲H共1−兩x兩兲 and H共x兲is the Heaviside unit step function. Therefore, we can compute the kernels by inverting the autocorrelation matrix ckl共s兲=Cll *共0兲 s2 冉 ␦ kl共e−s+s−1兲+e−s共1+兩k−l兩兲共es−1兲2 2共1− ␦ kl兲 冊 . 共9兲 We do that in two steps. First we multiply cby the intermediate matrix Rkl=cosh共s兲 ␦ kl−共1/2兲共 ␦ k,l+1+ ␦ k,l−1兲. The resulting matrix is tridiagonal and can be inverted using the same method we applied for the conservative force. The final result for the memory kernel is ␾ kl共s兲=s2关− ␰ 共s兲兴兩k−l兩 关s− sinh共s兲兴 ␮ 共s兲关sinh共s兲关1 − cosh共s兲兴 + ␦ kl ␮ 共s兲兴 −s ␦ kl −⌳kl * s,共10兲 where ␰ 共s兲=sinh共s兲−scosh共s兲+ ␮ 共s兲 s− sinh共s兲,共11兲 and ␮ 共s兲=冑s关cosh共s兲−1兴关s+scosh共s兲− 2 sinh共s兲兴.共12兲 The first piece of information we can extract from 共10兲is the friction coefficient that would be used in a DPD theory, 冕 0 ⬁ dt*⌽kl *共t*兲= lim s→0 ␾ kl共s兲=0, 共13兲 which vanishes here for all kand l. This result could be explained by the absence of real “dissipation” effects in harmonic systems like thermal conduction 关12兴. Additionally, due to the discontinuities in the time derivative of the autocorrelation function 共8兲, we expect that the memory kernels contain several Dirac delta functions, as we have also seen numerically 关8兴: ⌽kl *共t*兲=− ␦ 共t*−兩k−l兩+1兲− ␦ 共t*−兩k−l兩−1兲 +2 ␦ 共t*−兩k−l兩兲 + ␸ kl *共t*兲,共14兲 where ␸ kl *共t*兲is an otherwise smooth function of time. The discontinuities arise at the mesoscopic times t*=t/nthat the wave fronts of the perturbed coarse-grained particle 共containing noscillators, with n→⬁兲reach the boundaries of the coarse-grained particles as a consequence of the sound propagation. In principle, we could obtain ␸ kl *共t*兲by calculating the inverse Laplace transform of 共10兲. However, this is cumbersome and we actually prefer to calculate them numerically from molecular dynamics simulations, which also serves as a way to check the calculations. We present in Fig. 2 the first five memory kernels computed from a simulation with groups of n=10 000 oscillators. The highly oscillatory raw data has been coarse grained in time by applying the time convolution ⌽共t兲→兰d ␶ ⌽共 ␶ 兲 ␪ 共t− ␶ 兲, where ␪ 共t兲 =2g共2t/⌬兲/⌬, with ⌬*=⌬/n=0.01 关8兴. The numerical Laplace transform of these kernels is in complete agreement with the analytical prediction given by 共10兲. We can observe that the kernels are very slow decaying, having a significant magnitude at the end of the plot. In fact, we can use 共10兲to show that the absolute value of the kernel 兩 ␸ kl *共t*兲兩 is nonin- FIG. 2. Memory kernels ␸ kl *共t*兲in one dimension. BRIEF REPORTS PHYSICAL REVIEW E 72, 032101 共2005兲 032101-3 tegrable. This is a consequence of the fact that the analytic continuation of ␾ kl共s兲presents an infinite number of singularities on the imaginary axis 共s=i ␻ with real ␻ 兲. We have calculated numerically the memory kernels in higher dimensions. Figure 3 shows the first memory kernels of a two-dimensional system with square clumps containing 共400兲2oscillators. The data have been smoothed out with a mesoscopic time step of ⌬*=0.05. The behavior of the kernels is very similar to that of the one-dimensional case. The first two kernels 共solid and dashed lines兲display the same Dirac deltas 共conditioned in the figure by the convolution time step ⌬*兲, as predicted by Eq. 共14兲in the onedimensional case. Additionally, we present in Fig. 4 the results of a three-dimensional system with cubic clumps containing 共95兲3oscillators and ⌬*=0.1. Because the length of the clumps is too small to allow for a smooth representation of the kernels, we have also plotted in Fig. 4 the analogous results 共dotted and dash-dotted lines兲of a one-dimensional system with a similar clump size of 100 oscillators per group and ⌬*=0.1. The behavior of the kernels remains similar. The presence of the memory kernels, and the corresponding colored noise forces 关13兴, would render the simulation of the dissipative particles dynamics impractical due to computer memory limitations. A more promising approach would be to add new internal variables to the description so the process becomes Markovian. However, it is not clear if this can be done using a finite number of new internal variables. The thermodynamic variables 共local density and internal energy兲used in Refs. 关4,5兴would not account for the internal memory effect we have described in this work. This can be seen by considering a velocity perturbation traveling inside the coarse-grained particles. The center of mass, the local volume, and the internal energy of the coarse-grained particle would only be affected by the perturbation when it reaches the boundaries of the particle. On the other hand, a particle-based method could still be used if we regard them as virtual particles carrying energy, instead of representing portions of the system. In any case, the formalism used or the assumptions made would need to be closely related or checked at the microscopic level. D. C. thanks J. Casado-Pascual for helpful discussions. This work has been partially supported by a grant from the European Union and the EPSRC. 关1兴P. Hoogerbrugge and J. Koelman, Europhys. Lett. 19, 155 共1992兲. 关2兴P. Español and P. Warren, Europhys. Lett. 30, 191 共1995兲. 关3兴P. Español and M. Serrano, Phys. Rev. E 59, 6340 共1999兲. 关4兴E. G. Flekkøy and P. V. Coveney, Phys. Rev. Lett. 83, 1775 共1999兲. 关5兴P. Español, M. Serrano, and H. C. Öttinger, Phys. Rev. Lett. 83, 4542 共1999兲. 关6兴P. Prinsen, P. Warren, and M. Michels, Phys. Rev. Lett. 89, 148302 共2002兲. 关7兴M. Burke, R. Woscholski, and S. N. Yaliraki, Proc. Natl. Acad. Sci. U.S.A. 100, 13928 共2003兲. 关8兴D. Cubero and S. N. Yaliraki, J. Chem. Phys. 122, 034108 共2005兲. 关9兴R. Zwanzig, Nonequilibrium Statistical Mechanics 共Oxford University Press, New York, 2001兲. 关10兴N. W. Ashcroft and N. D. Mermin, Solid State Physics 共Saunders, Fort Worth, 1976兲. 关11兴R. Bird, C. F. Curtiss, R. Armstrong, and O. Hassager, Dynamics of Polymeric Liquids, Vol. 2: Kinetic Theory 共Wiley, New York, 1987兲. 关12兴G. Klein and I. Prigogine, Physica 共Amsterdam兲19,89 共1953兲. 关13兴P. Hänggi and P. Jung, Adv. Chem. Phys. 89, 239 共1995兲. 关14兴This force differs from that computed in P. Español, Phys. Rev. E53, 1572 共1996兲due to errors in the calculations. FIG. 3. Full memory kernels ⌽kl *共t*兲oscillating around the conservative force values ⌳kl *in two dimensions. For a fixed oscillator kchosen at the origin of the grid 共0, 0兲, the various curves correspond to different oscillators lalong the grid. FIG. 4. The same system as in Fig. 3, but in three dimensions. The dotted and the dash-dotted lines are the corresponding functions of a 1-D system, plotted here as a reference. BRIEF REPORTS PHYSICAL REVIEW E 72, 032101 共2005兲 032101-4