Inhomogeneous multiscale dynamics in harmonic lattices

Inhomogeneous multiscale dynamics in harmonic lattices David Cubero and Sophia N. Yaliraki Citation: The Journal of Chemical Physics 122, 034108 (2005); doi: 10.1063/1.1829253 View online: http://dx.doi.org/10.1063/1.1829253 View Table of Contents: http://scitation.aip.org/content/aip/journal/jcp/122/3?ver=pdfcov Published by the AIP Publishing Articles you may be interested in Effect of hydrodynamic correlations on the dynamics of polymers in dilute solution J. Chem. Phys. 138, 144902 (2013); 10.1063/1.4799877 Fluctuating hydrodynamics for multiscale simulation of inhomogeneous fluids: Mapping all-atom molecular dynamics to capillary waves J. Chem. Phys. 135, 044111 (2011); 10.1063/1.3615719 Structural dynamics of supercooled water from quasielastic neutron scattering and molecular simulations J. Chem. Phys. 134, 144508 (2011); 10.1063/1.3578472 Lattice dynamics in the half-space. Energy transport equation J. Math. Phys. 51, 083301 (2010); 10.1063/1.3451109 Comment on “Comment on ‘Constant temperature molecular dynamics simulations by means of a stochastic collision model. II. The harmonic oscillator’ [J. Chem. Phys. 104, 3732 (1996)]” [J. Chem. Phys. 106, 1646 (1997)] J. Chem. Phys. 120, 4991 (2004); 10.1063/1.1644801 This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 150.214.182.194 On: Mon, 01 Jun 2015 13:21:49 Inhomogeneous multiscale dynamics in harmonic lattices David Cubero and Sophia N. Yalirakia) Department of Chemistry, South Kensington Campus, Imperial College, London SW7 2AZ, United Kingdom 共Received 10 June 2004; accepted 14 October 2004; published online 3 January 2005兲 We use projection operators to address the coarse-grained multiscale problem in harmonic systems. Stochastic equations of motion for the coarse-grained variables, with an inhomogeneous level of coarse graining in both time and space, are presented. In contrast to previous approaches that typically start with thermodynamic averages, the key element of our approach is the use of a projection matrix chosen both for its physical appeal in analogy to mechanical stability theory and for its algebraic properties. We show that thermodynamic equilibrium can be recovered and obtain the fluctuation dissipation theorem a posteriori. All system-specific information can be computed from a series of feasible molecular dynamics simulations. We recover previous results in the literature and show how this approach can be used to extend the quasicontinuum approach and comment on implications for dissipative particle dynamics type of methods. Contrary to what is assumed in the latter models, the stochastic process of all coarse-grained variables is not necessarily Markovian, even though the variables are slow. Our approach is applicable to any system in which the coarse-grained regions are linear. As an example, we apply it to the dynamics of a single mesoscopic particle in the infinite one-dimensional harmonic chain. © 2005 American Institute of Physics. 关DOI: 10.1063/1.1829253兴 I. INTRODUCTION The processes underlying the properties of materials and biological assemblies often span a range of length and time scales. In understanding and predicting their behavior, it would be desirable to start from an atomistic description which could be capable of exhibiting the continuum and macroscopic behavior of the system. Molecular dynamics 共MD兲alone is presently incapable of bridging scales of such orders of magnitude. Coarse-graining schemes that tackle different aspects of this broad problem are currently an active area of research. Of particular interest is the possibility to follow the system dynamically and not just obtain its thermodynamic properties. A number of ‘‘mesoscale’’ dissipative particle dynamics 共DPD兲type methods, for example, introduce particlelike variables in continuous equations and associate with them conservative, dissipative, and random forces.1,2 The origin and appropriate functional form of these forces is still not fully resolved.3,4 A different approach that aims to avoid these equations altogether relies on exploiting through efficient numerical integration techniques a large number of short MD simulations.5Another level of complexity is introduced when the system needs to be inhomogeneously coarse grained. One of the most successful methods under this category, the quasicontinuum,6involves no explicit time nor finite temperature. A generalization has been attempted but the dynamics have been introduced ad hoc.7One of the main remaining challenges in such methods remains the reflections at boundaries between regions of different level of description.8Coarse-grained dynamics that can systematically describe the system at different scales in both time and space remain a challenge. A general scheme that prescribes the form of the equations to be used based on an equivalent footing or framework is currently lacking. We have used the projection operator approach of Mori and Zwanzig9to obtain the dynamical equations of an inhomogeneous multiscale system. To illustrate this, we solve exactly the multiscale inhomogeneous problem for harmonic systems. In this regard, harmonic systems, apart from their obvious adequacy as a first approximation for solids, provide an excellent framework to probe these ideas. The projector operator approach has been used extensively in the past. One of the main drawbacks has been the evaluation of the formal expressions obtained by the theory. By using an algebraic projection and in analogy to mechanical stability, we obtain functional forms of the forces that can be computed by a series of short molecular dynamics simulations. The key point of our approach is that we avoid altogether the first step of averaging and obtain the thermodynamic equilibrium properties a posteriori. The mechanical analogy has been first introduced in the literature by Deutch and Silbey for a single particle in a lattice.10 In that respect, our results can be considered a generalization of their work. We also recover as a special case the results of Adelman and Doll for atom/solid surface scattering in harmonic lattices.11 Additionally, the formalism provides a connection to the quasicontinuum method with explicit time and finite temperature. Finally, we comment on the form of the random forces and memory kernels that are usually assumed in other mesoscale dissipative dynamics type of methods. The paper is organized as follows: In Sec. II we introduce the system and the coarse-grained scheme. In Sec. III we describe the equations for the coarse-grained variables using the Mori formalism. In Sec. IV we introduce a different projection matrix and rederive the theory using that inner a兲Electronic mail: [email protected] THE JOURNAL OF CHEMICAL PHYSICS 122, 034108 共2005兲 122, 034108-10021-9606/2005/122(3)/034108/9/$22.50 © 2005 American Institute of Physics This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 150.214.182.194 On: Mon, 01 Jun 2015 13:21:49 product, which allows us to obtain closed expressions for the results in Sec. IV and a numerical procedure to compute them in Sec. V. In Sec. VI we propose a simulation procedure to study the dynamics. We apply the theory developed in the previous sections to study the coarse-grained dynamics of a single mesoscopic particle in the infinite onedimensional 共1D兲harmonic chain in Sec. VII. Finally, Sec. VIII provides a short summary and conclusions. II. SYSTEM DETAILS AND COARSE-GRAINING SCHEME We consider a harmonic system of dimension dconsisting of Nparticles of mass mi. We will assume that the interaction between any two particles only depends on the distance between them. The Hamiltonian of this system is H⫽兺 ␮ ⫽1 d 冉 兺 i 1 2mivi, ␮ 2⫹1 2兺 i,jxi, ␮ Aijxj, ␮ 冊 ,共1兲 where xi, ␮ and vi, ␮ are the ␮ components of the deviation from the equilibrium position and velocity of the particle i, respectively, and Aij is a symmetric matrix that satisfies the stability condition 兺 iAij⫽0. 共2兲 The variables of interest are the coordinates and velocities of the center of mass of each of the coarse-grained regions of the system Xk, ␮ ⫽1 Mk兺 i k mixi, ␮ ,Vk, ␮ ⫽1 Mk兺 i k mivi, ␮ ,共3兲 where Mk⫽兺i kmiis the total mass of the coarse-grained region or particle k. The number of oscillators in each region is allowed to vary in accordance with the desired coarsegraining level. For example, if we are interested in keeping the atomistic level in one part of the system, the corresponding regions contain only one oscillator and the coarsegrained variables are just the coordinates and velocities of the original oscillators. At the same time, not every oscillator need to be related to a coarse-grained particle. We may be interested in coarse-graining completely a part of the system comprising many oscillators. In that case no index kis associated with such region. We present in Fig. 1 an example of an inhomogeneous coarse-graining scheme. III. MORI THEORY The equations of motion for these variables can be obtained in a straightforward manner by using the Mori projection operator formalism.9The equations are obtained by projecting onto the subspace spanned by the relevant variables. The projection operator is determined by the choice of the inner product (B,C) in the Hilbert space of all functions B(x,v) and C(x,v) of phase space coordinates.9It is customary to choose canonical equilibrium averages (B,C) ⫽ 具 BC 典 eq , where 具 B共x,v兲 典 eq⫽ 兰 dNxdNvB共x,v兲e⫺H/kBT 兰 dNxdNve⫺H/kBT,共4兲 kBbeing the Boltzmann constant and Tthe temperature. We then obtain the following linear generalized Langevin equations: dXk共t兲 dt ⫽Vk共t兲,共5兲 Mk dVk共t兲 dt ⫽⫺兺 l 冋 ⌳klXl共t兲⫹ 冕 0 td ␶ ⌽kl共 ␶ 兲Vl共t⫺ ␶ 兲 册 ⫹Rk共t兲,共6兲 where ⌳kl is a renormalized matrix force, ⌽kl(t) a memory kernel, and Rk(t) a force which contains the information of the nonrelevant variables we have left out. We have dropped the spatial labels ␮ in these expressions for clarity, because the Hamiltonian we are using here does not provide any coupling between the different components. The explicit expressions for the unknown quantities in Eq. 共6兲involve complicated functions containing formal projection operators. However, the following information can be obtained relatively easily 共see Chap. 8 of Ref. 9兲: the force Rk(t) satisfies 具 Rk共t兲 典 ⫽0共7兲 and 具 Rk共t兲Rl共0兲 典 ⫽kBT⌽kl共t兲,共8兲 where 具¯典denotes averages over a statistical ensemble of the initial conditions that is close to equilibrium. Equation 共8兲is called the non-Markovian fluctuation-dissipation theorem. It guarantees that the system admits the thermodynamic equilibrium solution. The advantage of this approach is that it provides us with equations of motion where the forces are decomposed in three parts. The first two are functions of the relevant variables X(t) and V(t) alone, and can be identified with conservative and dissipative type force terms, respectively. The force Rk(t) can be considered as a random force. Therefore, the equations can be cast inside the theory of stochastic processes. While the physical interpretation of these equations is clear, the main problem resides in finding both closed forms of these expressions for ⌳kl and ⌽kl(t) as well as ways of actually computing them in practice. This is not an easy task FIG. 1. Example of a coarse-graining scheme in a 2D harmonic system. 034108-2 D. Cubero and S. N. Yaliraki J. Chem. Phys. 122, 034108 (2005) This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 150.214.182.194 On: Mon, 01 Jun 2015 13:21:49 even in the case of linear systems. For example, with this projection method, we find the matrix ⌳kl as the inverse of the matrix 具 Xl⬘Xl 典 eq /kBT, i.e., 兺 l⬘ ⌳kl⬘ 具 Xl⬘Xl 典 eq kBT⫽兺 l⬘ 具 XkXl⬘ 典 eq kBT⌳l⬘l⫽ ␦ kl .共9兲 This inverse matrix is generally difficult to compute. Special care is required as can be seen by noting that the individual elements 具 Xl 2 典 eq do not exist in the thermodynamic limit (N→⬁) because of translational invariance. IV. MORI THEORY WITH A SPECIAL INNER PRODUCT A. The mechanical analog: Random forces The expressions in the preceding section have been obtained by using equilibrium averages as the inner product. Although with this approach we obtained straightaway Eqs. 共5兲–共9兲, it is not clear how one can proceed further. A better insight is gained if we use instead an inner product with the following properties: 共xi, ␮ ,xj, ␯ 兲⫽共vi, ␮ ,vj, ␯ 兲⫽ ␦ ij ␦ ␮ ␯ mi ⫺1,共xi, ␮ ,vj, ␯ 兲⫽0. 共10兲 Using this inner product, the projection operator formalism produces the equivalent equations of the preceding section. Specifically, we recover Eq. 共5兲and a modified equation for the time derivative of V(t), Eq. 共6兲that now involves a memory kernel related to the time convolution of the position X(t). Integrating by parts this integral we arrive at Mk dVk共t兲 dt ⫽⫺兺 l 冋 AklXl共t兲⫹ 冕 0 td ␶␤ kl共 ␶ 兲Vl共t⫺ ␶ 兲 ⫹ ␤ kl共t兲Xl共0兲 册 ⫹R ˆk共t兲,共11兲 where we have written Akl⫽兺i k兺j lAij. We now need to compute the random force R ˆk(t) and the memory kernel function ␤ kl(t). The latter is connected to the former by the projection operators formalism as ␤ kl共t兲⫽⫺ 冕 0 td ␶ 关LR ˆk共 ␶ 兲,Xl共0兲兴Ml,共12兲 where Lis the Liouville9operator of the system. We later relate each quantity to those of the preceding section 关see Eq. 共6兲兴. The random force is in general a complicated function of the initial coordinates and velocities of the system. In harmonic systems, however, it can be expressed as a linear function of the initial conditions 共see Sec. 8.5 of Ref. 9兲. This will further allow us to establish a mechanical force analogy with a reference mechanical system in which the coarse-grained variables are held fixed. We start by assuming R ˆk共t兲⫽兺 iqi k共t兲mixi共0兲⫹ri k共t兲mivi共0兲.共13兲 Consistency requires that the time dependent coefficients are given by ri k(t)⫽ 兰 0 td ␶ qi k( ␶ ) and the set of equations miq ¨i k共t兲⫽⫺兺 jA ˆijqj k共t兲,共14兲 miqi k共0兲⫽⫺兺 j k 关共I⫺P ˆ兲A兴ij,q ˙i k共0兲⫽0, 共15兲 where the dot denotes time derivative, A ˆ⫽共1⫺P ˆ兲A共1⫺P ˆt兲共16兲 and P ˆis a projection matrix defined by P ˆij⫽mi Mk共i兲 ␦ k共i兲,k共j兲,共17兲 where k(i) denotes the index of the group in which the oscillator ibelongs and P ˆtthe transposed matrix. Note that if x⫽(x1,x2,...) is the vector formed by all oscillators’ coordinates, P ˆtx⫽Xgives a vector whose components are the centers of mass of the coarse-grained particles. In addition, if ior jbelongs to a region that has been totally coarse grained, then P ˆij⫽0. Thus, if both iand j belong to such regions, then A ˆij⫽Aij. It is easy to show that the new force matrix A ˆis symmetric and satisfies the stability condition 共2兲. Therefore, qk(t), and thus the random forces R ˆk(t), are expressed in terms of the mechanical problem 共14兲–共15兲. In fact, using the algebra we present below, we can show that the random forces can be further written as R ˆk共t兲⫽⫺兺 i k 冉 兺 jAij关x ˆj共t兲⫹v ˆj共t兲兴 冊 ,共18兲 where x ˆ(t) and v ˆ(t) are the solution of the mechanical problem 共14兲with the initial conditions x ˆ(0)⫽(I⫺P ˆt)x(0) and v ˆ(0)⫽(I⫺P ˆt)v(0). This is a generalization of the results presented in Ref. 10 for a single particle in a lattice. The mechanical problem generated by A ˆis related to the original mechanical problem of the underlying lattice but instead the center of mass of the coarse-grained particles are fixed, as illustrated by the fact that the following property is obeyed at all times: P ˆtqk共t兲⫽P ˆMqk共t兲⫽0, 共19兲 where Mij⫽ ␦ ijmiis the mass matrix. Note that if initially there is no disorder in the coarsegrained lattice, i.e., there are only collective initial deviations from equilibrium, then R ˆk(t)⫽0 at all times. This is a crucial property of linear systems that will allow us to compute the relevant magnitudes of the coarse-grained theory from simple simulations. The solution of Eqs. 共14兲–共15兲can be formally written as qk共t兲⫽Re共ei⍀t兲qk共0兲,共20兲 where ⍀is the frequency matrix given by ⍀2⫽M⫺1A ˆand i is the imaginary unit. Equation 共19兲is a consequence of A ˆP ˆt⫽0, which implies ⍀P ˆt⫽0. Thus, the subspace generated by P ˆtcorresponds to eigenvectors of A ˆwith zero eigenvalues, each one corresponding to every coarse-grained par- 034108-3 Inhomogeneous multiscale dynamics J. Chem. Phys. 122, 034108 (2005) This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 150.214.182.194 On: Mon, 01 Jun 2015 13:21:49 ticle. As a consequence, neither ⍀nor A ˆpossess an inverse. However, since both matrices are diagonalizable (A ˆis symmetric and ⍀is Hermitian with the scalar product xy ⫽xtMy), we can define the pseudoinverse ⍀⫺1so that ⍀⫺1⍀⫽⍀⍀⫺1⫽I⫺P ˆt,⍀⫺1Pt⫽0. 共21兲 Then, the pseudoinverse A ˆ⫺1⫽(⍀⫺1)2M⫺1turns out to be symmetric. B. Memory kernel We have obtained the random forces through the transformation to an equivalent mechanical system. We can now tackle the kernel ␤ kl(t), which using Eqs. 共13兲and 共12兲, can be written as ␤ kl共t兲⫽兺 j 冕 0 td ␶ 冕 0 ␶ d ␶ ⬘qj k共 ␶ ⬘兲Ai l共22兲 where Ai l⫽兺j lAij. In the limit t→⬁, each qi k(t) decays to zero as the initial perturbation 共15兲is propagated to the boundaries of the system by elastic waves. However, since Eq. 共22兲involves a double time integral, in general ␤ kl(t) does not decay to zero with time. A better representation is obtained by calculating the integrals in Eq. 共22兲using Eq. 共20兲, which leads to ␤ kl共t兲⫽⌽kl共t兲⫺⌽kl共0兲,共23兲 where ⌽kl共t兲⫽⫺兺 jAj l关A ˆ⫺1Mqk共t兲兴j共24兲 with the following properties 兺 l⌽kl共t兲⫽0 and ⌽kl共t兲⫽⌽lk共t兲.共25兲 This kernel is expected to vanish as t→⬁in the thermodynamic limit (N→⬁), in which case the boundaries are removed. Note that we have used the same letter to name this kernel and the one we obtained in Sec. III using equilibrium averages as the inner product. We will show in Sec. IVC that they are indeed the same. We can now return to the generalized Langevin equation 共11兲, which becomes Mk dVk共t兲 dt ⫽⫺兺 l关Akl⫺⌽kl共0兲兴Xl共t兲 ⫺兺 l 冕 0 td ␶ ⌽kl共 ␶ 兲Vl共t⫺ ␶ 兲 ⫺兺 l⌽kl共t兲Xl共0兲⫹R ˆk共t兲.共26兲 Equation 共26兲is the generalization to the present arbitrary inhomogeneous coarse graining of the results of Adelman and Doll for atom/solid surface scattering in harmonic solids.11 Note that with the identifications ⌳kl⫽Akl⫺⌽kl共0兲共27兲 and Rk共t兲⫽⫺兺 l⌽kl共t兲Xl共0兲⫹R ˆk共t兲,共28兲 Eq. 共26兲is just Eq. 共6兲. Finally, note that Eqs. 共25兲and 共2兲imply the following expected property of the conservative force 兺 l⌳kl⫽0. 共29兲 C. Fluctuation-dissipation theorems Using the explicit expressions 共13兲and 共24兲, and the integration formula for Gaussian distributions, we readily obtain 具 R ˆk共t兲R ˆl共0兲 典 A ˆ⫽kBT⌽kl共t兲,共30兲 where 具 ¯ 典 A ˆdenotes averages over the initial conditions with the canonical distribution ␳ ˆ(x,v)⬀exp(⫺H ˆ/kBT), where H ˆhas the same form as the Hamiltonian of Eq. 共1兲, but with A ˆinstead of A. This is the version of the fluctuation dissipation theorem that appears in Ref. 11. We show in Appendix A that the averages taken with 具 ¯ 典 A ˆin Eq. 共30兲are the same as the constrained averages 具 ¯ 典 X0⫽0defined by the canonical distribution with the center of mass of the ‘‘coarse-grained’’ particles X0fixed at their equilibrium values (X0⫽0), providing a physical meaning to the fluctuation-dissipation equation 共30兲. Nevertheless, note that Eq. 共30兲does not hold for general constrained averages 具 B共x,v兲 典 X0⫽ 兰 dNxdNvB共x,v兲e⫺H/kBT ␦ 共P ˆtx⫺X0兲 兰 dNxdNve⫺H/kBT ␦ 共P ˆtx⫺X0兲,共31兲 with X0different from their equilibrium values. For our purposes it is more convenient to invoke the standard form of the fluctuation-dissipation theorem, involving the random force Rk(t) 具 Rk共t兲Rl共0兲 典 eq⫽kBT⌽kl共t兲,共32兲 where the averages are taken with the full equilibrium distribution using A. We can prove Eq. 共32兲straightaway by using the following algebraic identity: A ˆ⫺1⫽A ˆ⫺1AP ˆtA⫺1P ˆAA ˆ⫺1⫹A ˆ⫺1AP ˆtA⫺1共I⫺P ˆ兲 ⫹共I⫺P ˆt兲A ˆ⫺1AP ˆtA⫺1⫹共I⫺P ˆt兲A⫺1共I⫺P ˆ兲,共33兲 where A⫺1is the pseudoinverse of Adefined in similar fashion to Eq. 共21兲by A⫺1A⫽AA⫺1⫽I⫺PT.(PT)ij⫽1/Nis the projector on the translational mode 关given by Eq. 共2兲兴, the only zero frequency mode in regular lattices. Using a similar algebra it can be readily shown that the conservative force ⌳kl defined by Eq. 共27兲is the same as the force we obtained in Sec. III, defined by Eq. 共9兲. Alternatively, we can prove this by taking the limit T→0 in Eqs. 共6兲 and 共26兲. Since Vk(t)→0 in that limit and both equations are exact for any initial condition, they must contain necessarily the same conservative force. Finally, note that if we write Eq. 共5兲as 034108-4 D. Cubero and S. N. Yaliraki J. Chem. Phys. 122, 034108 (2005) This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 150.214.182.194 On: Mon, 01 Jun 2015 13:21:49 冕 0 td ␶ Vl共t⫺ ␶ 兲⫽Xl共t兲⫺Xl共0兲,共34兲 multiply this equation by an arbitrary constant ␥ kl and add it to Eq. 共6兲we obtain Mk dVk共t兲 dt ⫽⫺兺 l 冋 ⌳ ¯ klXl共t兲⫹ 冕 0 td ␶ ⌽ ¯ kl共 ␶ 兲Vl共t⫺ ␶ 兲 册 ⫹R ¯ k共t兲,共35兲 where ⌳ ¯ kl⫽⌳kl⫹ ␥ kl ,⌽ ¯ kl⫽⌽kl⫺ ␥ kl , and R ¯ k⫽兺l ␥ klXl ⫹Rk. Therefore, by changing ␥ kl we can obtain an infinite number of different equations of motion, all with a random force R ¯ satisfying the fluctuation-dissipation formula 共30兲 共though with ⌽instead of ⌽ ¯ ). However, the corresponding kernel ⌽ ¯ kl at long times converges to ⫺ ␥ kl , resulting in an artificial nondecaying memory effect. Furthermore, R ¯ k共as well as R ˆk) does not satisfy the standard fluctuationdissipation formula 共32兲. In fact, in general the fluctuations of those random forces diverge in the unconstrained thermodynamic equilibrium due to translational invariance. Therefore, we will refer in the following to the Langevin equation given by Eq. 共6兲or equivalently Eq. 共26兲. V. A FEW IMPORTANT RESULTS The expression 共24兲we obtained for the memory kernel is still too formal for most practical situations. In this section we present a few exact results that will facilitate its calculation both analytically and numerically. A. Velocity autocorrelation function Let us consider the velocity autocorrelation matrix Ckl共t兲⫽ 具 Vk共t兲Vl共0兲 典 eq .共36兲 Multiplying Eq. 共6兲by Vl(0) and taking equilibrium averages we obtain a set of differential equations without noise Mk dCkl共t兲 dt ⫽⫺兺 l⬘ 冕 0 td ␶ 关⌳kl⬘⫹⌽kl⬘共t⫺ ␶ 兲兴Cl⬘l共 ␶ 兲, 共37兲 which together with the initial condition Ckl(0) ⫽ ␦ klMkkBTdetermine Ckl(t). Alternatively, if we know Ckl(t) we can derive ⌽kl(t) from these equations. In Appendix B we take advantage of this relationship to show that if Ckl(t) presents a discontinuity in its derivative at time t, then the memory kernel presents a Dirac-delta singularity at the same time t. To be more precise, ⌽kl共t兲⫽⫺ 冋 C ˙kl共0兲 Cll共0兲 ␦ 共t兲⫹兺 r ⌬C ˙kl共tr兲 Cll共0兲 ␦ 共t⫺tr兲 册 ⫹ ␸ kl共t兲,共38兲 where ⌬C ˙kl(tr)⫽C ˙kl(tr ⫹)⫺C ˙kl(tr ⫺), trdenotes the times where C ˙kl(t) is discontinuous, and ␸ kl(t) is an otherwise smooth function. This property will become useful when we consider coarse graining in the time scale in Sec. VII. The above results are general and not particular to linear systems. However, in harmonic systems there is a simpler way to compute the autocorrelation matrix. In fact, Ckl共t兲⫽Vk共t兲Vl共0兲,共39兲 where Vk(t) is the velocity of the coarse-grained particle kin the system starting from the following initial conditions: xi共0兲⫽0 and vi共0兲⫽ ␦ k共i兲,lVl共0兲,共40兲 with Vl(0)2⫽MlkBT. This method exploits the fact that the noise Rk(t) vanishes with this particular choice of the initial conditions, as shown in Sec. IVA. B. Numerical calculation of the kernel and conservative forces Note that Eqs. 共15兲and 共19兲imply that if the group k contains only one oscillator, and it is connected through the force matrix Awith other single-oscillator groups, then qi k(t)⫽0 for all iand hence ⌽kl(t)⫽Rk(t)⫽0 at all times. Thus, in a region where we have kept the atomistic description, the generalized Langevin equation 共6兲reduces to the original Newton’s equation. Therefore, the problem is reduced to computing the conservative forces ⌳kl and the kernels ⌽kl(t) for those coarsegrained particles that belong to actual coarse-grained regions or in the proximities of them. Because of Eq. 共27兲, we just need to calculate the memory kernel ⌽kl(t). This can be done by solving Eq. 共14兲numerically. Then, by using Eq. 共22兲we can obtain ␤ kl(t). This involves a double numerical integration. The long time limit of ␤ kl(t) gives ⌽kl(0), from which we can determine ⌳kl and ⌽kl(t). Alternatively, we can use the more direct method of Cai et al.8In this approach we run a molecular dynamics simulation of the harmonic system starting from the following initial condition: xi共0兲⫽ ⑀ ␦ k共i兲land vi共0兲⫽0. 共41兲 Next, the center of mass of the coarse-grained particles is kept fixed by means of an external force on each oscillator Fk(i). Therefore, Eq. 共26兲becomes 兺 i k Fk共i兲共t兲⫺关⌳kl⫹⌽kl共t兲兴 ⑀ ⫽0. 共42兲 By determining the external force required to keep the center of mass fixed we can obtain the conservative forces and the memory kernel. This method is also applicable to situations in which the harmonic character is a first-order approximation of a more complex system. The methods above provide us with a numerical representation of ⌽kl(t) at a high accuracy. This is numerically ‘‘exact’’ in the microscopic time scale. However, the kernel is likely to display a highly oscillatory behavior when we look at it on a coarser time scale. In order to obtain a reasonable smooth function in the mesoscopic time scale we need to coarse grain in time further. We can do that as follows: let us call ␦ and ⌬Ⰷ ␦ the basic time steps in the microscopic and mesoscopic time scales, respectively. One im- 034108-5 Inhomogeneous multiscale dynamics J. Chem. Phys. 122, 034108 (2005) This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 150.214.182.194 On: Mon, 01 Jun 2015 13:21:49 mediate approach would be to calculate the coarse-grained kernel ⌽ ¯ kl(t) at the center of the mesoscopic interval jas a simple time integral: ⌽ ¯ 共tj兲⫽1 ⌬ 冕 tj⫺⌬/2 tj⫹⌬/2d ␶ ⌽共 ␶ 兲,共43兲 with tj⫽⌬(j⫹1/2). However, the set of data ⌽ ¯ (tj) is going to depend strongly on the particular choice we have made of the exact location of the mesoscopic times tj. For example, if we look at ⌽ ¯ (tj⫹ ␦ ), defined as in Eq. 共43兲but integrating from tj⫺/⌬/2⫹ ␦ to tj⫹/⌬/2⫹ ␦ instead, we may obtain a very different value due to the high oscillatory behavior of the raw data. Nevertheless, we can define instead the coarsegrained kernel at tjas the average of all these possible values inside that time step ⌬. This is equivalent to calculating the following convolution: ⌽ ¯ 共t兲⫽ 冕 d ␶ ⌽共 ␶ 兲 ␪ 共t⫺ ␶ 兲,共44兲 where ␪ 共t兲⫽2 ⌬ 冉 1⫺ 兩 t 兩 2 ⌬ 冊 H 冉 1⫺ 兩 t 兩 2 ⌬ 冊 ,共45兲 and H(x) is the Heaviside unit-step function. By applying Eq. 共44兲two times, we are able to reproduce the analytical results for the example we present in Sec. VII. VI. AN ‘‘ AB INITIO ’’ METHOD FOR THE SIMULATION OF THE DYNAMICS Once we know ⌽kl(t) numerically we can simulate the dynamics by using a DPD-like algorithm. If we assume that the stochastic process Rk(t) is Gaussian 共which is justified as long as we consider small deviations from equilibrium兲, then it is determined by the first moments: 共7兲and 共8兲. In this case, we can exploit the properties 共25兲to generate the set of correlated variables Rk(t) by using a set of independent Gaussian variables ␰ kl⫽⫺ ␰ lk , so that Rk(t)⫽兺l⫽k ␰ kl and 具 ␰ kl共t兲 ␰ k⬘l⬘共0兲 典 eq⫽kBT⌽kl共t兲共 ␦ kl⬘ ␦ lk⬘⫺ ␦ kk⬘ ␦ ll⬘兲.共46兲 This is the essential ingredient used in DPD to preserve momentum conservation. In fact, since Eqs. 共29兲and 共25兲imply ⌳kk⫽⫺兺 l⫽k⌳kl and ⌽kk共t兲⫽⫺兺 l⫽k⌽kl共t兲,共47兲 respectively, we can write Eq. 共6兲in the form Mk dVk共t兲 dt ⫽兺 l⫽k 再 ⫺⌳kl关Xl共t兲⫺Xk共t兲兴 ⫺ 冕 0 td ␶ ⌽kl共t⫺ ␶ 兲关Vl共 ␶ 兲⫺Vk共 ␶ 兲兴⫹ ␰ kl共t兲 冎 , 共48兲 which very much resembles the DPD equations of motion. However, the first question that arises is whether we can consider ⌽kl(t) as Dirac-delta functions, as customary in DPD. It is frequently argued that they can be considered so because the coarse-grained variables are slow compared to the atomic time scales. This would imply that the full stochastic process 兵 X(t),V(t) 其 would be Markovian in the relevant time scale. We will show next that this should not be assumed. VII. AN EXAMPLE: COARSE-GRAINING IN THE 1D HARMONIC CHAIN Let us illustrate the above results with a simple example: the dynamics of a single coarse-grained particle or blob formed by nconsecutive oscillators in the one-dimensional infinite harmonic chain. For simplicity we set all intrinsic parameters of the chain to unity 共the mass of each oscillator, the elastic constant, and the equilibrium spacing between them兲. Let us denote with Xand Vthe relevant variables 共now scalar magnitudes兲of the single blob. Due to translational invariance, the conservative force ⌳on the blob vanishes. This can be shown by using 具 X2 典 →⬁in Eq. 共9兲. Therefore, the generalized Langevin equation 共6兲becomes ndV共t兲 dt ⫽⫺ 冕 0 td ␶ ⌽共 ␶ 兲V共t⫺ ␶ 兲⫹R共t兲.共49兲 We now use the method proposed in Sec. VA to compute the memory kernel. The autocorrelation function can be obtained by studying the dissipation of an initial velocity perturbation consisting of all oscillators at rest at the equilibrium positions except the ones in the blob, which start instead with velocity V0. The time dependent velocity V(t)of the coarse-grained particle can be computed directly from the exact solution of the infinite 1D harmonic chain, which is expressed in terms of Bessel functions, V共t兲⫽V0兺 i⫽1 n 兺 j⫽1 n J2 兩 i⫺j 兩 共2t兲,共50兲 or from the macroscopic displacement field u(x,t),12 which verifies the wave equation ⳵ 2u ⳵ t2⫽ ⳵ 2u ⳵ x2.共51兲 This equation can be solved easily by standard methods, and the velocity is obtained by integrating ⳵ u(x,t)/ ⳵ tover the blob’s region. Both methods provide the same result for large n, V共t兲⬃V0共1⫺t/n兲H共1⫺t/n兲.共52兲 Note that V(t) changes in a time scale of order ⬃n.We already see from Eq. 共52兲that the velocity autocorrelation decay is not exponential but linear, which is an indication that the stochastic process is non-Markovian. We can solve Eq. 共49兲by using Laplace transforms.13 This is done in Appendix C. By changing to the proper time scale t*⫽t/n, with ⌽*(t*)⫽n⌽(nt*), we obtain 冕 0 ⬁dt*⌽*共t*兲⫽2, 共53兲 and 034108-6 D. Cubero and S. N. Yaliraki J. Chem. Phys. 122, 034108 (2005) This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 150.214.182.194 On: Mon, 01 Jun 2015 13:21:49 ⌽*共t*兲⫽ ␦ 共t*兲⫺ ␦ 共t*⫺1兲⫹et*⫺H共t*⫺1兲et*⫺1共t*⫹1兲 ⫹兺 k⫽2 ⬁ H共t*⫺k兲et*⫺k共⫺1兲k共t*⫺k兲k⫺2 共k⫺2兲! ⫻ 冉 共t*⫺k兲2 k共k⫺1兲⫹2共t*⫺k兲 k⫺1⫹1 冊 .共54兲 Note that the Dirac deltas appear as a consequence of the time scale change. At the microcospic time scale the autocorrelation function 共50兲does not present discontinuities in its time derivative. We have plotted in Fig. 2 the memory kernel 共without the Dirac-delta functions兲given by Eq. 共54兲. We have also calculated numerically this kernel by using the method of Sec. VB. Both results are indistinguishable in the plot. It is clear that the kernel has not decayed completely in the time scale of V(t)(t*⬃1), and thus we cannot approximate it with a single Dirac delta and the random term as white noise, even in the limit n→⬁. This phenomenon is general and has implications beyond this particular example. The reason for this long-lived kernels is that the time scale of Vk(t) is given by the dissipation, which in harmonic systems is due to wave propagation, and the decay time scale of the kernels is on the same time scale. This is clear from Eqs. 共14兲–共15兲. In the mesoscopic limit A ˆtends to A, and qi k(t) tends to the solution of the mechanical problem consisting of an initial perturbation around the edge of the coarse-grained region k. Thus, ⌽kl(t) will not appreciably decay at least until the elastic waves have removed the perturbation, which will be about the time it takes for the sound to cross the mesoscopic region, i.e., a mesoscopic time. This is a memory effect that occurs inside the coarse-grained particles, causing the whole process to be non-Markovian. Therefore, we cannot consider the kernels or the autocorrelation function of the random forces Dirac-delta functions, such as in DPD. In Ref. 14, Espan ˜ ol justified the DPD method by studying the equation of motion of an infinite number of coarsegrained particles in the one-dimensional harmonic chain. In contrast to the example discussed above, now the infinite harmonic chain is partitioned into an infinite number of blobs, each one formed by nconsecutive oscillators, which is more closely related to the customary coarse-graining scheme in DPD. However, this work relies on the Markovian assumption for the dynamics, which we have shown to be unjustified. We present in Fig. 3 the scaled memory kernel ⌽3 *⫽n⌽3connecting two coarse-grained particles kand l separated by two coarse-grained particles, i.e., 兩 k⫺l 兩 ⫽3, as computed from a molecular dynamics simulation with blobs of size n⫽10000. As before, the raw data has been coarse grained in time with a mesoscopic time step ⌬*⫽⌬/n ⫽0.01. It is clearly seen that the kernel is on the mesoscopic time scale t*⫽t/n, as expected. In fact, by applying the theory of Sec. VA we obtain that this kernel contains three Dirac-delta functions at different mesoscopic times, 2 ␦ (t*⫺3)⫺ ␦ (t*⫺2)⫺ ␦ (t*⫺4) 共recognizable in Fig. 3 as lines going out of scale兲, around t*⫽3, which is the time it takes the sound to propagate from one blob to the other. We report the corresponding analytical expressions of this problem and the genelarization to higher dimensions elsewhere. VIII. CONCLUSIONS We have applied projection operators to the coarsegrained multiscale problem in harmonic systems with an inhomogeneous level of coarse graining. The customary approach of using equilibrium averages as inner product produces the right generalized Langevin equation, but the explicit expressions are difficult to calculate. Using an alternative inner product we have been able to provide explicit expressions for the conservative forces, the memory kernels and the random forces, in terms of a mechanical analog. Based on these expressions, we have shown that the memory kernels and the velocity autocorrelation functions can be computed in linear systems from a single molecular dynamics simulation. These results represent a generalization of previous analytical work on harmonic lattices involving the derivation of a series of algebraic properties. In addition, a method that resembles DPD has been proposed and is applicable to any system in which the coarsegrained region is linear. This can also be seen as a natural extension of the quasicontinuum method to account for finite temperature and dynamics. Moreover, we have shown that the memory kernels and also the autocorrelation function of the random forces used in the simulations should not be considered as Dirac-delta functions, as customary in most approaches, without justification. Instead, the proposed simulations would need to include the computation of the time FIG. 2. Scaled memory kernel ⌽*⫽n⌽as a function of t*⫽t/n共without the Dirac deltas at t*⫽0, 1兲for a coarse-grained particle in a 1D harmonic chain. FIG. 3. Memory kernel connecting two coarse-grained particles separated by two blobs in the 1D harmonic chain. 034108-7 Inhomogeneous multiscale dynamics J. Chem. Phys. 122, 034108 (2005) This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 150.214.182.194 On: Mon, 01 Jun 2015 13:21:49 convolutions involving memory kernels. In order to account for the fluctuation-dissipation theorem, the random forces need to be generated using a colored noise algorithm. Otherwise, one risks that the system simulated does not tend to the thermodynamic equilibrium. Methods addressing this issue exist in the literature,15,16 however, further work is required to test their efficiency in these particular problems. ACKNOWLEDGMENT This work was funded in part by the EU growth project SENTIMATS under Contract No. G5RD-CT-2001. APPENDIX A: CONSTRAINED AVERAGES Starting from Eq. 共13兲we obtain 具 R ˆk共t兲R ˆl共0兲 典 A ˆ⫽qk共t兲tM 具 xxt 典 A ˆMql共0兲 ⫽qk共t兲tM共1⫺P ˆt兲 具 xxt 典 A ˆ共1⫺P ˆ兲Mql共0兲 ⫽qk共t兲tM 具 共x⫺X兲共x⫺X兲t 典 A ˆMql共0兲. 共A1兲 Let us now consider the variable change from xto center mass coordinates X⫽P ˆtxplus a number of relative coordinates which we will represent using vectorial notation as ␰ . Each variable ␰ ican be defined as the difference between the component xiof the oscillator iand a fixed reference oscillator inside the same group. Then, the quantity B(x) ⫽(x⫺X)(x⫺X)tas a function of 兵 X, ␰ 其 does not depend on the center of mass coordinates X, i.e., it is a function of ␰ alone: B(x⫹␭X)⫽B(x) for all ␭and x. We show next that for any function with the same property, 具 B(x) 典 A ˆ ⫽ 具 B(x) 典 X0⫽0. Equation 共16兲implies xtAx⫽xtA ˆx⫺XtAX⫹2xtAX,共A2兲 being xtA ˆxa function of ␰ alone. Therefore, using this expression and changing variables we obtain 具 B共x兲 典 X0⫽0⫽ 兰 dNxB共x兲exp共⫺xtAx/2kBT兲 ␦ 共P ˆtx兲 兰 dNxexp共⫺xtAx/2kBT兲 ␦ 共P ˆtx兲 ⫽ 兰 dNxB共x兲exp共⫺xtA ˆx/2kBT兲 ␦ 共P ˆtx兲 兰 dNxexp共⫺xtA ˆx/2kBT兲 ␦ 共P ˆtx兲 ⫽ 兰 dXd ␰ B共x兲exp共⫺xtA ˆx/2kBT兲 ␦ 共X兲 兰 dXd ␰ exp共⫺xtA ˆx/2kBT兲 ␦ 共X兲 ⫽ 兰 d ␰ B共x兲exp共⫺xtA ˆx/2kBT兲 兰 d ␰ exp共⫺xtA ˆx/2kBT兲 ⫽ 兰 dXd ␰ B共x兲exp共⫺xtA ˆx/2kBT兲 兰 dXd ␰ exp共⫺xtA ˆx/2kBT兲 ⫽ 具 B共x兲 典 A ˆ.共A3兲 APPENDIX B: SINGULARITIES IN THE MEMORY KERNEL Let us consider the Fourier transform of the autocorrelation function ckl共 ␻ 兲⫽ 冕 0 ⬁dte⫺i ␻ tCkl共t兲共B1兲 and analogously for the memory kernel ␾ kl( ␻ ). In this representation, Eq. 共37兲becomes 兺 l⬘ ␾ kl⬘共 ␻ 兲cl⬘l共 ␻ 兲⫽Mk关Ckl共0兲⫺i ␻ ckl共 ␻ 兲兴 ⫺兺 l⬘ ⌳kl⬘ cl⬘l共 ␻ 兲 i ␻ .共B2兲 Assume C ˙kl(t) is discontinuous at the times t ⫽t1,t2,.... Then, the asymptotic behavior of ckl( ␻ ) in the limit ␻ →⬁can be obtained by integration by parts17 c共 ␻ 兲⬃C共0兲 i ␻ ⫺1 ␻ 2 冋 C ˙共0兲⫹兺 r⌬C ˙共tr兲e⫺i ␻ tr 册 ,共B3兲 where ⌬C ˙(tr)⫽C ˙(tr ⫹)⫺C ˙(tr ⫺). Inserting this expression in Eq. 共B2兲we obtain 兺 l⬘ ␾ kl⬘共 ␻ 兲cl⬘l共 ␻ 兲⬃⫺Mk i ␻ 冋 C ˙共0兲⫹兺 r⌬C ˙共tr兲e⫺i ␻ tr 册 . 共B4兲 On the other hand, since ␾ 共 ␻ 兲should be bounded in the limit ␻ →⬁, we also have 兺 l⬘ ␾ kl⬘共 ␻ 兲cl⬘l共 ␻ 兲⬃ ␾ kl共 ␻ 兲Mk Cll共0兲 i ␻ .共B5兲 Combining both expressions ␾ kl共 ␻ 兲⬃⫺1 Cll共0兲 冋 C ˙共0兲⫹兺 r⌬C ˙共tr兲e⫺i ␻ tr 册 ,共B6兲 in the limit ␻ →⬁. The inverse transform is given by Eq. 共38兲. APPENDIX C: CALCULATION OF THE MEMORY KERNEL The Laplace transform of the autocorrelation function V(s) is obtained from Eq. 共52兲as 034108-8 D. Cubero and S. N. Yaliraki J. Chem. Phys. 122, 034108 (2005) This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 150.214.182.194 On: Mon, 01 Jun 2015 13:21:49