UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Bottom-up coarse-graining : theory, implementation, application Tran, Hoang 2018

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Item Metadata

Download

Media
24-ubc_2018_may_tran_hoang.pdf [ 2.45MB ]
Metadata
JSON: 24-1.0365974.json
JSON-LD: 24-1.0365974-ld.json
RDF/XML (Pretty): 24-1.0365974-rdf.xml
RDF/JSON: 24-1.0365974-rdf.json
Turtle: 24-1.0365974-turtle.txt
N-Triples: 24-1.0365974-rdf-ntriples.txt
Original Record: 24-1.0365974-source.json
Full Text
24-1.0365974-fulltext.txt
Citation
24-1.0365974.ris

Full Text

  Bottom-up coarse graining: Theory, implementation, application  by  Hoang Tran  B.A., Earlham College, 2014  A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF  MASTER OF SCIENCE in THE FACULTY OF GRADUATE AND POSTDOCTORAL STUDIES (Physics)  THE UNIVERSITY OF BRITISH COLUMBIA (Vancouver)  April 2018  ยฉ Hoang Tran, 2018ii  Abstract  Molecular dynamics simulations at the atomic scale are a powerful tool to study the structure and dynamics of biological systems. However, important processes, such as protein folding, are still too computationally expensive to be investigated with atomistic simulations. In these situations, different coarse-grained force fields have been developed to appropriately reproduce the static structure of the systems. These strategies can range from specially designed, ad hoc approaches to transferrable, bottom-up plans. Nevertheless, because of the reduction in number of degrees of freedom, the obtained dynamics typically suffer from inaccuracies, which may be not adjustable using trivial scaling. To correctly reintroduce these properties, one can apply the Mori-Zwanzig formalism, essentially injecting the memory effects into these coarse-graining schemes. In practice, because of the difficulty in algorithm designing and the limit on computing power, the memory is often assumed to be short-term, and therefore can be approximated as a delta function. Recently, some investigations have been done successfully on simple models with significantly long memory. This work aims to extend these successes on systems with more complex topology. In particular, we applied the Mori-Zwanzig formalism to the system of polyethylene chains. Using GROMACS united atom model, molecular dynamics simulation was performed to obtain the coarse-grained force field and the memory kernel of the nonbonded interaction. These elements allowed for a successful coarse-grained simulation, which captured both static and dynamical structures of the reference system.  iii  Lay summary  Molecular dynamics is a discipline of interest to biologists because it creates the bridge between theoretical formulation and experiment data. Full atomistic models have been developed and shown to give accurate predictions, but they take too much time to simulate many important processes using currently available computers. In these situations, coarse-grained simulations, which consider the overall motions of clusters of atoms and ignore less relevant details, are more efficient and useful. Nevertheless, the removal of less significant atoms produces observable consequences on the physics of these models, namely the pseudo-atoms have incorrect motion due to the lack of โ€œobstaclesโ€. In theory, it is possible to rectify these problems with the Mori-Zwanzig formalism. However, it is difficult to put this method in practice, so it has only been investigated recently, on a very simple system. In this work, we discuss the application of the Mori-Zwanzig formalism on a more complex model.  iv  Preface  The author did most of the tasks discussed in this work. The simulation of the coarse-grained model and the data analysis were done with his own implementation in C++. The atomistic simulation used for reference was done by Professor Mark Thachuk using GROMACS.  v  Table of Contents  Abstract......................................................................................................................................................................... ii Lay summary ............................................................................................................................................................... iii Preface ..........................................................................................................................................................................iv Table of Contents........................................................................................................................................................... v List of Figures ............................................................................................................................................................. vii List of Abbreviations ....................................................................................................................................................ix Acknowledgement ......................................................................................................................................................... x Chapter 1 Introduction............................................................................................................................................ 1 1.1. Integrator problem ......................................................................................................................................... 3 1.2. Memory problem ........................................................................................................................................... 4 1.3. Mean force problem ....................................................................................................................................... 5 1.4. Outline ........................................................................................................................................................... 5 Chapter 2 Theory .................................................................................................................................................... 6 2.1. Consistent coarse-graining ............................................................................................................................. 6 2.2. Mori-Zwanzig formalism ............................................................................................................................... 8 2.3. Memory as damped oscillators .................................................................................................................... 10 2.4. Integration of dissipative particle dynamics ................................................................................................ 12 Chapter 3 Implementation .................................................................................................................................... 14 3.1. Integrators .................................................................................................................................................... 14 3.1.1. Dissipative velocity-Verlet .................................................................................................................. 14 3.1.2. Lowe-Andersen .................................................................................................................................... 15 3.1.3. Impulsive leap-frog .............................................................................................................................. 15 3.1.4. Test result............................................................................................................................................. 16 3.1.4.1. Model systems ............................................................................................................................. 16 3.1.4.2. Physical quantities ....................................................................................................................... 17 3.1.4.3. Simulation conditions .................................................................................................................. 18 3.1.4.4. Result for Langevin dynamics ..................................................................................................... 18 vi  3.1.4.5. Result for dissipative particle dynamics ...................................................................................... 20 3.2. Memory ....................................................................................................................................................... 24 3.2.1. Standard memory ................................................................................................................................. 24 3.2.2. Test result for standard memory .......................................................................................................... 25 3.2.2.1. Model system ............................................................................................................................... 25 3.2.2.2. Result ........................................................................................................................................... 26 3.2.3. Auxiliary memory ................................................................................................................................ 30 3.2.4. Test result for auxiliary memory .......................................................................................................... 31 3.2.4.1. Model system ............................................................................................................................... 31 3.2.4.2. Result ........................................................................................................................................... 32 3.3. Mean force potential .................................................................................................................................... 35 3.3.1. Inverse Boltzmann method .................................................................................................................. 35 3.3.2. Test result............................................................................................................................................. 35 3.3.2.1. Model system ............................................................................................................................... 35 3.3.2.2. Result ........................................................................................................................................... 37 Chapter 4 Application .......................................................................................................................................... 39 4.1. Coarse-graining of bonded interactions ....................................................................................................... 39 4.1.1. Problem with coarse-graining of bonded interactions.......................................................................... 39 4.1.2. Solution ................................................................................................................................................ 40 4.2. Memory extraction ....................................................................................................................................... 46 4.3. Simulation model ......................................................................................................................................... 51 4.4. Result ........................................................................................................................................................... 53 Chapter 5 Conclusion ........................................................................................................................................... 55 Bibliography ................................................................................................................................................................ 57 Appendix A ................................................................................................................................................................. 61 Appendix B .................................................................................................................................................................. 68   vii  List of Figures  Fig. 1-1: The multitude of CG schemes. ........................................................................................................................ 2 Fig. 3-1: Time evolution of the average kinetic temperature ๐‘˜๐ต๐‘‡ for different integrators in LD model. .................. 19 Fig. 3-2: Time evolution of the mean square displacement for different integrators in LD model. ............................ 19 Fig. 3-3: The decay of velocity autocorrelation function over time for different integrators in LD model. ................ 20 Fig. 3-4: Time evolution of the average kinetic temperature ๐‘˜๐ต๐‘‡ for different integrators in DPD model. ................ 21 Fig. 3-5: Evolution of the mean square displacement for different integrators in DPD model. .................................. 21 Fig. 3-6: The decay of velocity autocorrelation function over time for different integrators in DPD model. ............. 22 Fig. 3-7: The radial distribution function ๐บ(๐‘…). .......................................................................................................... 23 Fig. 3-8: Coefficients for dissipative force (left) and random force (right), with some choices of ๐‘ and ๐‘›. ................ 26 Fig. 3-9: Average kinetic temperature ๐‘˜๐ต๐‘‡ as a function of cutoff length ๐‘š for the standard memory scheme. ........ 26 Fig. 3-10: Ratio between the kinetic temperature ๐‘˜๐ต๐‘‡ and the reference temperature ๐‘˜๐ต๐‘‡ as a function of ๐‘˜๐ต๐‘‡ in the standard memory scheme. ........................................................................................................................................... 27 Fig. 3-11: The relationship between ๐‘‡/โŒฉ๐‘‡โŸฉ and ๐‘˜๐ต๐‘‡ for different system sizes. ......................................................... 28 Fig. 3-12: Effect of the timestep value on โŸจ๐‘˜๐ต๐‘‡โŸฉ. ........................................................................................................ 29 Fig. 3-13: Time evolution of the average kinetic temperature ๐‘˜๐ต๐‘‡ for two different interpretations of ๐‘’๐ผ๐ฝ in Eq. (56). ..................................................................................................................................................................................... 29 Fig. 3-14: Time evolution of the average kinetic temperature ๐‘˜๐ต๐‘‡ in the auxiliary memory scheme. ........................ 33 Fig. 3-15: Time evolution of the mean square displacement for the auxiliary memory shame. .................................. 33 Fig. 3-16: Log-log plot of the above mean square displacement versus time for the auxiliary memory method. ....... 34 Fig. 3-17: The decay of velocity autocorrelation function over time. ......................................................................... 34 Fig. 3-18: CG of a polyethylene chain of 40 united atom monomers with the rate of 20 to 1. ................................... 36 Fig. 3-19: Inverted nonbonded potential and fitting function. ..................................................................................... 37 Fig. 3-20: Comparison between targeted distribution ๐‘”๐‘›๐‘(๐‘…) and the CG distribution ๐บ๐‘›๐‘(๐‘…) obtained from simulation using nonbonded Boltzmann-inverted potential with the CG ratio of 20 to 1 and fixed bond length. ....... 38 Fig. 4-1: Inverted bond potential and fitting function. ................................................................................................. 39 viii  Fig. 4-2: Comparison between targeted radial distribution ๐‘”(๐‘…) and the CG distribution ๐บ(๐‘…) for interactions of nonbonded pairs (square) and of pairs sharing a CG bond (circle), with the CG ratio of 20 to 1. ............................... 40 Fig. 4-3: Bond radial distribution from the thermalized simulation using Andersen thermostat with ๐‘‡ = 5.0๐พ. ....... 41 Fig. 4-4: Bond radial distribution from the thermalized simulation using Andersen thermostat with ๐‘‡ = 50.0๐พ (yellow triangle) and ๐‘‡ = 200.0๐พ (green circle). .................................................................................................................... 42 Fig. 4-5: Inverted potential and fitting function for the bonded interactions, with the CG ratio of 8 to 1................... 43 Fig. 4-6: Inverted nonbonded potential and fitting function for the CG ratio of 8 to 1. .............................................. 43 Fig. 4-7: Atomistic and CG distributions for the bonded interactions, with the CG ratio of 8 to 1. ............................ 44 Fig. 4-8: Atomistic and CG distributions for the bonded interactions for the CG ratio of 8 to 1. ............................... 45 Fig. 4-9: Two bonded CG particles. ............................................................................................................................ 47 Fig. 4-10: Decay rate of the radial dependence in the parallel (left) and orthogonal (right) memory kernel. ............. 48 Fig. 4-11: Temporal dependence at ๐‘… > 0.8๐‘›๐‘š of the parallel (left) and orthogonal (right) memory kernel. ........... 48 Fig. 4-12: Temporal correlation ๐พ๐ผ๐ฝ๐‘ฅ(๐‘…, ๐‘ก)/๐พ๐ผ๐ฝ๐‘ฅ(๐‘…, 0) and its fitted function for the parallel (left) and orthogonal (right) memory kernel. ............................................................................................................................................................ 49 Fig. 4-13: Valence angle. ............................................................................................................................................. 51 Fig. 4-14: CG of a polyethylene chain of 40 united atom monomers with the rate of 8 to 1. ..................................... 52 Fig. 4-15: Nonbonded distribution for the CG simulation using the auxiliary memory method. ................................ 53 Fig. 4-16: Nonbonded distribution for the CG simulation using the auxiliary memory method. ................................ 53 Fig. 4-17: Normalized VACF for the CG simulation using the auxiliary memory method. ........................................ 54   ix  List of Abbreviations  MD: Molecular dynamics CG: Coarse-graining MZ: Mori-Zwanzig EOM: Equation of motion LD: Langevin dynamics DPD: Dissipative particles dynamics COM: Center of mass FDT: Fluctuation-dissipation theorem VV: Velocity-Verlet DPD-VV: Dissipative velocity-Verlet Lowe: Lowe-Andersen Imp-LF: Impulsive leap-frog MSD: Mean square displacement VACF: Velocity autocorrelation function RDF: Radial distribution function NPT: Isothermalโ€“isobaric ensemble NVT: Canonical ensemble NVE: Microcanonical ensemble   x  Acknowledgement  I would like to thank Professor Mark Thachuk for his unreserved support and his never-ending patience, and Professor Joerg Rottler for his keen advices. 1  Chapter 1 Introduction    Computational techniques for the simulation of physical, chemical, or biological problems are based on a large variety of theoretical models. All these computational techniques can reach a typical time and length scale. Among them, classical molecular dynamics (MD) simulation has become standard for studying molecular systems, since it allows for the precise representation of molecular structure while still avoiding the complication of quantum physics [1]. The trajectory of the particles is computed by numerically solving equations of classical mechanics (most commonly Newtonโ€™s equations of motion) for an N-body interacting system, where potential energy is defined by known molecular mechanics force fields. With the current CPU speed, atomistic MD simulations that can be done in tolerable real time can support tens of thousands of individual objects, and up to a few nanoseconds for the simulation times [2]. However, many important chemical and biological mechanisms involve much larger numbers of atoms and much longer durations. For example, a typical protein folding process occurs in microseconds and some polymer chains contain millions of atoms [3]. In addition, the interactions between the particles are long-range, so the force computation scales as O(N2) and is usually the most CPU-demanding evaluation of the whole simulation. This inefficiency can be mitigated partly by reducing the algorithmic complexity (Verlet neighbor list for Lennard-Jones force, particle-mesh method for electrostatic interaction) or improving the computational resources (using machines specifically adapted to MD simulations or combining a lot of general-purpose CPUs with concurrent computing). Notwithstanding, one can still only simulate a small, limited fraction of these important processes. This is a major reason why development of coarse-graining (CG) methods is needed [4]. Coarse-graining is the process of merging multiple particles and describing their overall motion using fewer variables. The first CG models were proposed practically at the birth of MD, but they only became popular in the last decade [5]. There are good reasons for this increasing role. First, as mentioned above, despite rapid advances in computing power, CG approaches are still the only appropriate choices when the system is relatively large or the timescale is relatively long. Second, atomistic studies produce vast amount of data, some of which are necessary for the understanding of the mechanism in question, but the majority is irrelevant [6]. By focusing on key details, while averaging over more trivial elements, CG provides crucial conceptual advantages with respect to atomistic models. Notwithstanding, only recently did the computer become fast enough for MD simulations to support sufficiently large number of particles in a reasonably long duration, so that the result is interesting and verifiable [2]. Thus, for many processes, CG schemes are now in the โ€œsweet spotโ€ where there is adequate power to enable acceptable resolutions, but not adequate to do so with all-atomistic techniques. Consequently, CG methods have quickly grown in popularity in recent years and are now recognized as an equally valuable tool for the study of molecules. 2  In the last twenty years, many CG models have been introduced, targeting different, specific types of systems. They can assume various levels of reduction, ranging from united-atom, with the hydrogen atoms being absorbed into their connected, larger partners, to very crude representations, where tens of these groups are treated as a single interaction site [7].  Preferably, the elements in a CG cluster are related through some predetermined, fixed links, such as bonds or angles, but some schemes allow more extreme combinations. As an example, Fig. 1-1 exhibits some of the possible designs: a typical CG group represents a whole molecule (A), or parts of one (B, C, D). However, there are models that combine multiple independent molecules (E) [8], or even modify the grouping on-the-fly [9, 10]. While CG strategies are many, in essence, we can classify them into two major philosophies: bottom-up (structure-based) and top-down (thermodynamic-based). In the first approach, the CG force fields are extracted from reference atomistic simulations using approximation schemes such as Boltzmann inversion, inverse Monte Carlo, or force matching. On the other hand, top-down CG starts with simple, ad hoc potential form, then optimizes the parameters in an iterative procedure to reproduce key experimental data, especially static properties [8]. Although the bottom-up approaches can capture more of the fine details of the interaction, the top-down method is usually easier to employ and transfer [6]. On the other hand, bottom-up schemes enable a multiscale-modelling pipeline, where one uses more-detailed systems as the basis for the force field of the coarser models [2]. Although the distinction is quite intuitive, it should be emphasized that this is just a simple guideline, not a hard boundary, as an increasing number of models integrate both paradigms [6, 11, 12].   Fig. 1-1: The multitude of CG schemes. Either way, the dynamic properties of the CG simulation, such as diffusion coefficients or time correlations, all suffer inaccuracies compared to that of the reference system, which can be an actual, physical system or a higher-resolution model [13]. This problem arises because of the elimination of degrees of freedom that should appear in the form of thermal noise: intuitively, with less collisions, the particles will move more freely and display different dynamic behaviors. These disparities can be as simple as a scaling of the system, or much more complicated and 3  irreducible. To emulate the correct dynamical properties, one can employ the Mori-Zwanzig (MZ) formalism [14] on the CG model, essentially injecting a memory effect, in which the systemโ€™s past trajectory influences its current motion. Mathematically, the CG procedure is basically a projection operator from a larger vector space (the reference system) into a smaller space (the CG system). In the process, one must sacrifice some information, of which the reduction is inevitably displayed in the dynamics [15]. The MZ formalism deals with this loss of information through the orthogonal operator and shows that it amounts to a new equation of motion (EOM): with  โŒฉ๐น๐‘–๐‘—โŒชโƒ—โƒ—โƒ—โƒ—โƒ—โƒ— โƒ—โƒ—  being the average of the instantaneous force ๐น๐‘–๐‘—โƒ—โƒ—โƒ—โƒ— , ๐‘ฒ๐‘–๐‘—(๐‘ก โˆ’ ๐‘ ) being the pairwise memory kernel, which relates the past velocities with the dissipation at the current time, and ๐›ฟ๐น๐‘–๐‘—โƒ—โƒ— โƒ—โƒ— โƒ—โƒ—  โƒ— the random component, representing the unresolved degrees of freedom during the CG procedure. The first term results from the application of the CG procedure, while the last two are the consequence of the MZ formalism. In principle, their explicit form can all be determined from the reference system [16], but one will encounter multiple practical obstacles when trying to design the numerical implementation. These hindrances contributed to why this scheme has not been investigated until recently, even though the MZ formalism has been known for more than forty years.  1.1. Integrator problem For the discerning eyes, the three terms on the right-hand side of Eq. (1) brings to mind the familiar Langevin dynamics (LD) [17], with the first addend corresponding to the conservative force, the second to the friction force and the third being the random force. In fact, as long as the timescale of the memory and of the change in momentum are well-separated, we can apply the Markovian approximation [18]. This condition reduces Eq. (1) into the EOM of dissipative particle dynamics (DPD), a pairwise generalization of LD equation. Originally proposed during the 90s in the soft matter community [19], DPD is a solution to correct the aforementioned discrepancy in the CG dynamics by introducing fluid viscosity effects and replicating the missing internal degrees of freedom. In the most generic form, the EOM in DPD can be written as The conservative force ๐น๐‘–๐‘—๐ถโƒ—โƒ— โƒ—โƒ—   acts as the source of the thermodynamic structure of the system, while the dissipative term ๐น๐‘–๐‘—๐ทโƒ—โƒ—โƒ—โƒ—  โƒ— and the random term ๐น๐‘–๐‘—๐‘…โƒ—โƒ— โƒ—โƒ—   together form a thermostat to reproduce the dynamic properties [20]. Although the viscosity effects can also be reproduced using LD, through the employment of pairwise dissipative and random forces, DPD describes how the hydrodynamic interactions between the particles affect the systemโ€™s approach to equilibrium, and thus maintains the correct canonical distribution. However, the finite integration of Eq. (2) proves to be troublesome, as the friction force between the particles at a particular time ๐‘ก is an explicit function of their current relative velocity ๐‘ฃ๐‘–๐‘—โƒ—โƒ— โƒ—โƒ—  (๐‘ก). To calculate that velocity, the  ๐‘š ๐‘Ž๐‘–โƒ—โƒ—  โƒ— = โˆ‘ (โŒฉ๐น๐‘–๐‘—โƒ—โƒ—โƒ—โƒ— โŒช +1๐‘˜๐ต๐‘‡โˆซ ๐‘ฒ๐‘–๐‘—(๐‘ก โˆ’ ๐‘ )๐‘ฃ๐‘–๐‘—โƒ—โƒ— โƒ—โƒ—  (๐‘ )๐‘‘๐‘ ๐‘ก0+ ๐›ฟ๐น๐‘–๐‘—โƒ—โƒ—โƒ—โƒ— )๐‘–โ‰ ๐‘— , (1)  ๐‘š ๐‘Ž๐‘–โƒ—โƒ—  โƒ— = โˆ‘ ๐น๐‘–๐‘—๐ถโƒ—โƒ— โƒ—โƒ—  + ๐น๐‘–๐‘—๐ทโƒ—โƒ—โƒ—โƒ—  โƒ— + ๐น๐‘–๐‘—๐‘…โƒ—โƒ— โƒ—โƒ—  ๐‘–โ‰ ๐‘—  . (2) 4  standard MD integration strategy requires the forces ๐น๐ผ๐ฝโƒ—โƒ—โƒ—โƒ—  โƒ—(๐‘ก) to be known beforehand, including the dissipative force ๐น๐ผ๐ฝ๐ทโƒ—โƒ—โƒ—โƒ—  โƒ—(๐‘ก). Therefore, ๐‘ฃ๐‘–๐‘—โƒ—โƒ— โƒ—โƒ—  (๐‘ก) and ๐น๐ผ๐ฝ๐ทโƒ—โƒ—โƒ—โƒ—  โƒ—(๐‘ก) are circularly dependent on each other [1], and a self-consistent approach is required to solve the DPD EOM exactly. Thus, we cannot simply commandeer the usual MD packages and utilize their algorithms in DPD simulations. It should be noted here that the conclusion may or may not be applicable to Eq. (1), since most of the velocities involved are determined from the trajectorial history. On the other hand, the historical velocities are all propagated with the same integrator, thus it is arguable that any possible inaccuracy simply compounds instead of mitigating. Recent works using the MZ formalism have not addressed the integrator issue [18, 21], so it is unclear whether it will be a problem. Nevertheless, by the point this study was conducted, most of the popular simulation programs had not incorporated, or incorporated poorly, a correct DPD integrator [22, 23]. Thus, we took it upon ourselves to implement different DPD algorithms available in the literature and assess their accuracy/performance.  1.2. Memory problem One of the main reasons for the lack of interest in the MZ formalism is the requirement for a memory kernel. Studies using the MZ formalism have been sporadic over the years, and the calculation is often simplified by applying the Markovian assumption, where the memories at different time are uncorrelated. In short, the features of interest evolve on a time scale much larger than the correlation time of the memory, so the memory integration can be approximated as a time-independent friction tensor and the EOM turns into a classic stochastic differential equation [14]. This solution works well in some applications, but obviously it is not universally true and fails in many others. For instance, they lead to the โ€œplateau-value problemโ€, where the memory tensor does not go to zero at infinite time, so some ad hoc estimations should be introduced [15]. In these situations, the employment of a memory function is imperative. However, it is very costly, computation-wise, if one tries the direct, naรฏve implementation, where the simulation remembers its past trajectory: since a history must exist for every pair of particles, it grows as O(N2). While this implementation is algorithmically similar to that of the Verlet neighbor list, it demands a recalculation at every timestep, unlike the latter. This, combined with the fact that the memory can be significant over tens or hundreds of timesteps, increases the CPU time by a few orders of magnitude and partially neutralizes the advantages of CG schemes. Recently, Li et al. have shown that if the memory vanishes at infinity, it can be approximated as a series of dampened harmonic oscillators [24]. These oscillators absorb the memory effect into their time-evolution, so the program does not need to keep the history anymore and the time complexity reduces to linear. Li et al. have successfully employed this technique in the star-polymer melts system, where each polymer is one pseudo-atom. In our study, we investigate whether the methodology is applicable to more a complex system of polyethylene chains.  5  1.3. Mean force problem The problem of deriving the mean force potential from a given atomistic system and the corresponding CG scheme is well-known and can be solved with many different methods [2]. In fact, as the aim of this study is not identifying a new CG force field, we can simply borrow from the literature a CG model that has been investigated thoroughly. Nevertheless, to reconstruct the memory kernel, we still need to perform the atomistic simulation ourselves, so we decided that having our own CG machinery is preferable. Since the model of interest contains polyethylene molecules, which cannot be converted into singular, non-bonded CG particles, we implement a solution that can derive both non-bond and bond mean force potentials, including that of the bond angles.  1.4. Outline In this study, we try to replicate and enhance the results done by Li et al., who applied successfully the MZ formalism on a system of star polymer melts [13, 21]. Chapter 2 discusses in more detail the theoretical foundations of the various problems we have mentioned earlier. In Chapter 3, we describe the implementation of different approaches to these questions, and the degree of success for each. Finally, in Chapter 4, we choose the best-working solution and apply it to a polyethylene system. Our probe provides some improvements over the works of Liโ€™s group, as they did not examine the validity of the MZ formalism on bonded CG particles.  6  Chapter 2 Theory    2.1. Consistent coarse-graining While this study mainly deals with the CG procedure of an all-atom model, the following framework can be applied to any two models for a given physical system. After all, the idea of bonded connections already contains a notion of CG, since the true atoms have no such structure and only interact with each other through their quantum fields. Henceforth, we shall refer to the higher-resolution model as atomistic, its interaction sites as atoms, whereas the lower-resolution model is called CG, with its sites being particles. Finally, we reserve lower-case symbols for atomistic variables, while CG variables are capitalized.  This section largely follows the derivation done by Noid et al. [2]. Now, let us consider an atomistic system consisting of ๐‘› structureless atoms with mass ๐‘š๐‘–. We define ๐œ‘(?ฬ‚?) = ๐œ‘(๐‘Ÿ๐‘›โƒ—โƒ—โƒ—โƒ— , ๐‘๐‘›โƒ—โƒ— โƒ—โƒ—  ) to be the atomistic phase space, where ๐‘Ÿ๐‘›โƒ—โƒ—โƒ—โƒ— = (๐‘Ÿ1โƒ—โƒ—โƒ—  , โ€ฆ , ๐‘Ÿ๐‘›โƒ—โƒ—  โƒ—) and ๐‘๐‘›โƒ—โƒ— โƒ—โƒ—  = (๐‘1โƒ—โƒ—  โƒ—, โ€ฆ , ๐‘๐‘›โƒ—โƒ—โƒ—โƒ— ) are the canonical positions and momenta of its atoms, respectively. The Hamiltonian is with ๐‘ข(๐‘Ÿ๐‘›โƒ—โƒ—โƒ—โƒ— ) being the atomistic multi-body potential, and the corresponding equilibrium probability density in the canonical ensemble is with Next, we divide the atomistic system into ๐‘ particles, with ๐‘ < ๐‘›, using some mapping operator ๐‘ท. Obviously, in exchange for computational performance, we need to reduce the number of variables, necessarily discarding some information. It is possible to recover this loss partially by embedding the particles with extra data, such as rotation, but usually we want to employ the same formalism. Therefore, the CG particles should also be structure-less:  โ„Ž(?ฬ‚?) = โ„Ž( ๐‘Ÿ๐‘›โƒ—โƒ—โƒ—โƒ— , ๐‘๐‘›โƒ—โƒ— โƒ—โƒ—  ) =  โˆ‘๐‘๐‘–22๐‘š๐‘–๐‘–+ ๐‘ข( ๐‘Ÿ๐‘›โƒ—โƒ—โƒ—โƒ— ) , (3)  ๐‘๐œ‘(?ฬ‚?) = ๐‘๐‘Ÿ(๐‘Ÿ๐‘›โƒ—โƒ—โƒ—โƒ— )๐‘๐‘(๐‘๐‘›โƒ—โƒ— โƒ—โƒ—  ) , (4)  ๐‘๐‘Ÿ(๐‘Ÿ๐‘›โƒ—โƒ—โƒ—โƒ— ) โˆ exp (โˆ’๐›ฝ๐‘ข( ๐‘Ÿ๐‘›โƒ—โƒ—โƒ—โƒ— )) , (5)  ๐‘๐‘(๐‘๐‘›โƒ—โƒ— โƒ—โƒ—  ) โˆ exp (โˆ’๐›ฝ โˆ‘๐‘๐‘–22๐‘š๐‘–๐‘–) . (6) 7  where ๐‘…๐‘โƒ—โƒ— โƒ—โƒ—  โƒ— = (๐‘…1โƒ—โƒ—โƒ—โƒ— , โ€ฆ , ๐‘…๐‘โƒ—โƒ—โƒ—โƒ—  โƒ—) and ๐‘ƒ๐‘โƒ—โƒ— โƒ—โƒ—  โƒ— = (๐‘ƒ1โƒ—โƒ—  โƒ—, โ€ฆ , ๐‘ƒ๐‘โƒ—โƒ— โƒ—โƒ—  ) are the CG coordinates and momenta. Although ๐‘ท can still be quite general, in practice, we typically choose the CG particle to be the center of mass (COM) of its constituent atoms, with the corresponding phase space ฮฆ(?ฬ‚?) = ฮฆ(๐‘…๐‘โƒ—โƒ— โƒ—โƒ—  โƒ—, ๐‘ƒ๐‘โƒ—โƒ— โƒ—โƒ—  โƒ—) being determined by [13, 25] where ๐ผ denotes the particle. Because Eq. (7) becomes right-unique in this case, i.e. each atomistic state maps to only one CG state, the equilibrium probability density ๐‘ƒ๐œ‘ of a given COM state (๐‘…๐‘โƒ—โƒ— โƒ—โƒ—  โƒ—, ๐‘ƒ๐‘โƒ—โƒ— โƒ—โƒ—  โƒ—) is simply the probabilistic sum for all microscopic states (๐‘Ÿ๐‘›โƒ—โƒ—โƒ—โƒ— , ๐‘๐‘›โƒ—โƒ— โƒ—โƒ—  ) that are projected to that CG state, i.e. [14] However, we acquire this CG equilibrium distribution from the reference atomistic simulation, not the CG simulation itself. In line with the traditional CG schemes, we want the equilibrium averages obtained from the time evolution of the new CG system to be consistent with the atomistic ones, that is with being the Hamiltonian of the CG system. Using Eqs. (4), (9) and (10), we separate the variables and find the consistency conditions for the CG Hamiltonian into [2] Equation (12) shows how we can resolve the CG potential of mean force ๐‘ˆ(๐‘…๐‘โƒ—โƒ— โƒ—โƒ—  โƒ—)   from the atomistic potential ๐‘ข(๐‘Ÿ๐‘›โƒ—โƒ—โƒ—โƒ— ) and the mapping operator ๐‘ท. We shall discuss the actual calculation in Section 3.3. For now, it should only be noted  ๐‘ท: (๐‘Ÿ๐‘›โƒ—โƒ—โƒ—โƒ— , ๐‘๐‘›โƒ—โƒ— โƒ—โƒ—  ) โ†ฆ (๐‘…๐‘โƒ—โƒ— โƒ—โƒ—  โƒ—, ๐‘ƒ๐‘โƒ—โƒ— โƒ—โƒ—  โƒ—) , (7)  ๐‘€๐ผ๐‘…๐ผโƒ—โƒ—โƒ—โƒ— = โˆ‘๐‘š๐‘–๐‘Ÿ๐‘– โƒ—๐‘–โˆˆ๐ผ , ๐‘ƒ๐ผโƒ—โƒ—  โƒ— = โˆ‘ ๐‘๐‘–โƒ—โƒ—โƒ—  ๐‘–โˆˆ๐ผ , ๐‘€๐ผ = โˆ‘๐‘š๐‘–๐‘–โˆˆ๐ผ , (8)  ๐‘ƒ๐œ‘(๐‘…๐‘โƒ—โƒ— โƒ—โƒ—  โƒ—, ๐‘ƒ๐‘โƒ—โƒ— โƒ—โƒ—  โƒ—) = โˆซ๐‘‘๐‘Ÿ๐‘›โƒ—โƒ—โƒ—โƒ— โˆซ ๐‘‘๐‘๐‘›โƒ—โƒ— โƒ—โƒ—  ๐‘๐œ‘(๐‘Ÿ๐‘›โƒ—โƒ—โƒ—โƒ— , ๐‘๐‘›โƒ—โƒ— โƒ—โƒ—  )๐›ฟ (๐‘ท(๐‘Ÿ๐‘›โƒ—โƒ—โƒ—โƒ— ) โˆ’ ๐‘…๐‘โƒ—โƒ— โƒ—โƒ—  โƒ—) ๐›ฟ (๐‘ท(๐‘๐‘›โƒ—โƒ— โƒ—โƒ—  ) โˆ’ ๐‘ƒ๐‘โƒ—โƒ— โƒ—โƒ—  โƒ—) . (9)  ๐‘ƒ๐œ‘(?ฬ‚?) = ๐‘ƒฮฆ(?ฬ‚?) โˆ exp (โˆ’๐›ฝ โˆ‘๐‘ƒ๐ผ22๐‘š๐ผ๐ผโˆ’ ๐›ฝ๐‘ˆ ( ๐‘…๐‘โƒ—โƒ— โƒ—โƒ—  โƒ—)) , (10)  ๐ป(?ฬ‚?) = โˆ‘๐‘ƒ๐ผ22๐‘š๐ผ๐ผ+ ๐‘ˆ ( ๐‘…๐‘โƒ—โƒ— โƒ—โƒ—  โƒ—) , (11)  ๐‘ƒ๐‘… (๐‘…๐‘โƒ—โƒ— โƒ—โƒ—  โƒ—) โˆ exp (โˆ’๐›ฝ๐‘ˆ ( ๐‘…๐‘โƒ—โƒ— โƒ—โƒ—  โƒ—)) โˆ โˆซexp (โˆ’๐›ฝ๐‘ข( ๐‘Ÿ๐‘›โƒ—โƒ—โƒ—โƒ— )) ๐›ฟ (๐‘ท(๐‘Ÿ๐‘›โƒ—โƒ—โƒ—โƒ— ) โˆ’ ๐‘…๐‘โƒ—โƒ— โƒ—โƒ—  โƒ—) ๐‘‘๐‘Ÿ๐‘›โƒ—โƒ—โƒ—โƒ—  , (12)  ๐‘ƒ๐‘ƒ (๐‘ƒ๐‘โƒ—โƒ— โƒ—โƒ—  โƒ—) โˆ exp (โˆ’๐›ฝ โˆ‘๐‘ƒ๐ผ22๐‘š๐ผ๐ผ) โˆ โˆซexp (โˆ’๐›ฝ โˆ‘๐‘๐‘–22๐‘š๐‘–๐‘–) ๐›ฟ (๐‘ท(๐‘๐‘›โƒ—โƒ— โƒ—โƒ—  ) โˆ’ ๐‘ƒ๐‘โƒ—โƒ— โƒ—โƒ—  โƒ—) ๐‘‘๐‘๐‘›โƒ—โƒ— โƒ—โƒ—   . (13) 8  that the computed function is unique up to an additive constant. On the other hand, our specific choice of  ๐‘ท satisfies Eq. (13) automatically. The reader is invited to explore the justification for the above claims in the work on multiscale CG by Noid et al. [2].  2.2. Mori-Zwanzig formalism Previously, we mentioned that the CG procedure results in unavoidable information loss. Our choice of CG variables, determined by Eq. (8), preserves the static structure of the atomistic system, hence the disparity can only manifest in the dynamics. In this section, we shall clarify our assertion and its implications. However, the actual derivation is long and complicated, so the bulk of the calculations is presented in Appendix A. First, a distinction must be made between the CG variables obtained from the CG simulation and their corresponding value computed directly using the atomistic trajectory. As the potential of mean force ๐‘ˆ(๐‘…๐‘โƒ—โƒ— โƒ—โƒ—  โƒ—) is an ensemble average, the force derived from this potential differs from the instantaneous total force acting on the actual COM, so the two time evolutions diverge from each other. Mathematically, we can express the time evolution of the (๐‘…๐‘โƒ—โƒ— โƒ—โƒ—  โƒ—, ๐‘ƒ๐‘โƒ—โƒ— โƒ—โƒ—  โƒ—) along the atomistic path as [14] where ฮฉ is the CG phase space density, ๐’ซ is the ensemble-averaged projection, ๐’ฌ โ‰ก 1 โˆ’ ๐’ซ and is the atomistic Liouville operator. The first term on the right-hand side of Eq. (14) can be resolved completely from the projection operator ๐’ซ and the atomistic Hamiltonian. Intuitively, we can see that it is the time evolution of ฮฉ acquired from the CG simulation. On the other hand, ๐’ฌ is the complementary operator to ๐’ซ, thus the second term represents the data lost in the CG projection. From Eq. (14), we can derive the EOM of the particles as [14] with ๐‘ฒ๐ผ๐ฝ being the pairwise memory kernel and  ๐›ฟ๐น๐ผ๐“ โƒ—โƒ—โƒ—โƒ—  โƒ— the deviation of the total atomistic force from the mean force acting on ๐ผ, namely  ๐‘‘๐‘‘๐‘กฮฉ = ๐‘–๐‘™ฮฉ = ๐’ซ๐‘–๐‘™ฮฉ + ๐’ฌ๐‘–๐‘™ฮฉ , (14)  ๐‘–๐‘™ = โˆ’โˆ‘ [๐œ•โ„Ž๐œ•๐‘๐‘–โƒ—โƒ—โƒ—  ๐œ•๐œ•๐‘Ÿ๐‘– โƒ—โˆ’๐œ•โ„Ž๐œ•๐‘Ÿ๐‘– โƒ—๐œ•๐œ•๐‘๐‘–โƒ—โƒ—โƒ—  ]๐‘–  (15)  ๐‘…๐ผโƒ—โƒ—โƒ—โƒ— ฬ‡ =๐‘ƒ๐ผโƒ—โƒ—  โƒ—๐‘€๐ผ , (16)  ๐‘ƒ๐ผโƒ—โƒ—  ฬ‡โƒ— = โˆ’๐œ•๐œ•๐‘…๐ผโƒ—โƒ—โƒ—โƒ— ๐‘ˆ (๐‘…๐‘โƒ—โƒ— โƒ—โƒ—  โƒ—) โˆ’ โˆ‘โˆซ ๐‘ฒ๐ผ๐ฝ(๐‘ก โˆ’ ๐‘ )๐‘ก0๐‘ƒ๐ฝโƒ—โƒ—  โƒ—(๐‘ )๐‘€๐ฝ๐‘‘๐‘ ๐ฝ + ๐›ฟ๐น๐ผ๐“ โƒ—โƒ—โƒ—โƒ—  โƒ—(๐‘ก) , (17) 9  Equation (17) reminds us of the EOM for DPD, where the three terms on the right-hand side corresponds to the mean, friction, and fluctuating forces. Theoretically, all three terms can be determined completely with the CG variables. The mean force is calculated from the N-body potential of mean force. The friction force depends on the history of the system because its integrand involves past momenta. Finally, the fluctuating force is not necessarily random, since it is the atomistic fluctuation, which has a known equilibrium distribution in ๐œ‹(๐‘Ÿ๐‘›โƒ—โƒ—โƒ—โƒ— , ๐‘๐‘›โƒ—โƒ— โƒ—โƒ—  ). Essentially, while we cannot accurately predict the microscopic configuration (๐‘Ÿ๐‘›โƒ—โƒ—โƒ—โƒ— , ๐‘๐‘›โƒ—โƒ— โƒ—โƒ—  ) from a given CG state, the causal relation in the atomistic space restricts the possible CG events, manifesting through the temporal correlation in the friction force, which requires where โŸจโ€ฆ โŸฉ is defined as the thermal average over the atomistic canonical ensemble: Evidently, this is the generalized formulation of the second fluctuation-dissipation theorem (FDT) for non-Markovian fluctuations [13]. In principle, assessing Eq. (17) directly allows us to calculate the particlesโ€™ motion. However, such evaluation is nontrivial because the force field depends intrinsically on the ๐‘-body configuration of the CG system [14]. To make Eq. (17) practical, we shall employ three assumptions. First, we approximate the non-bonded interactions as sums of pairwise-additive forces, i.e. Second, we take which implies the fluctuating forces are random for different pair, to get: Third, we assume that the forces between neighboring particles depend only on their COM positions, essentially neglecting the ๐‘-body correlation:  ๐›ฟ๐น๐ผ๐“ โƒ—โƒ—โƒ—โƒ—  โƒ— = (โˆ’โˆ‘๐œ•๐œ•๐‘Ÿ๐‘– โƒ—๐‘ข(๐‘Ÿ๐‘›โƒ—โƒ—โƒ—โƒ— )๐‘–โˆˆ๐ผ) โˆ’ (โˆ’๐œ•๐œ•๐‘…๐ผโƒ—โƒ—โƒ—โƒ— ๐‘ˆ (๐‘…๐‘โƒ—โƒ— โƒ—โƒ—  โƒ—)) โ‰ก ๐น๐ผ๐‘ƒโƒ—โƒ—โƒ—โƒ—  โƒ— โˆ’ ๐น๐ผโƒ—โƒ—  โƒ— . (18)  โŒฉ๐›ฟ๐น๐ผ๐“ โƒ—โƒ—โƒ—โƒ—  โƒ—(๐‘ก) โŠ— ๐›ฟ๐น๐ฝ๐“ โƒ—โƒ—โƒ—โƒ—  โƒ—(0)โŒช = ๐‘˜๐ต๐‘‡๐‘ฒ๐ผ๐ฝ(๐‘ก) , (19)  โŸจ๐‘Ž(?ฬ‚?)โŸฉ = โˆซ ๐‘Ž(?ฬ‚?)๐‘๐œ‘(?ฬ‚?)๐‘‘?ฬ‚? (20)  ๐น๐ผโƒ—โƒ—  โƒ— โ‰ˆ โˆ‘ ๐น๐ผ๐ฝโƒ—โƒ—โƒ—โƒ—  โƒ—๐ฝโ‰ ๐ผ , ๐›ฟ๐น๐ผ๐’ฌโƒ—โƒ—โƒ—โƒ—  โƒ— โ‰ˆ โˆ‘ ๐›ฟ๐น๐ผ๐ฝ๐’ฌโƒ—โƒ—โƒ—โƒ—  โƒ—๐ฝโ‰ ๐ผ . (21)  โŒฉ๐›ฟ๐น๐ผ๐ฝ๐“ โƒ—โƒ—โƒ—โƒ—  โƒ—(๐‘ก) โŠ— ๐›ฟ๐น๐พ๐ฟ๐“ โƒ—โƒ—โƒ—โƒ—โƒ—โƒ— (0)โŒช = 0 for ๐ผ๐ฝ โ‰  ๐พ๐ฟ , (22)  โŒฉ๐›ฟ๐น๐ผ๐“ โƒ—โƒ—โƒ—โƒ—  โƒ—(๐‘ก) โŠ— ๐›ฟ๐น๐ฝ๐“ โƒ—โƒ—โƒ—โƒ—  โƒ—(0)โŒช = โŒฉโˆ‘ ๐น๐ผ๐พโƒ—โƒ— โƒ—โƒ—  โƒ—(๐‘ก)๐พโ‰ ๐ผโŠ— โˆ‘๐น๐ผ๐พโƒ—โƒ— โƒ—โƒ—  โƒ—(0)๐พโ‰ ๐ผโŒช = โŒฉ๐›ฟ๐น๐ผ๐ฝ๐“ โƒ—โƒ—โƒ—โƒ—  โƒ—(๐‘ก) โŠ— ๐›ฟ๐น๐ผ๐ฝ๐“ โƒ—โƒ—โƒ—โƒ—  โƒ—(0)โŒช (23) 10  Thus, we arrive at the pairwise form of Eq. (17) where ๐‘…๐ผ๐ฝโƒ—โƒ— โƒ—โƒ—  โƒ— = ๐‘…๐ผโƒ—โƒ—โƒ—โƒ— โˆ’ ๐‘…๐ฝโƒ—โƒ—โƒ—โƒ— , ๐น๐ผ๐ฝ๐‘ƒโƒ—โƒ—โƒ—โƒ—  โƒ— = โˆ‘ ๐น๐‘–๐‘—โƒ—โƒ—โƒ—โƒ— ๐‘–โˆˆ๐ผ,๐‘—โˆˆ๐ฝ  is the pairwise instantaneous force between particles ๐ผ and ๐ฝ, ๐น๐ผ๐ฝโƒ—โƒ—โƒ—โƒ—  โƒ— and ๐›ฟ๐น๐ผ๐ฝ๐’ฌโƒ—โƒ—โƒ—โƒ—  โƒ— are its mean and fluctuating force, and the matrix is the pairwise memory kernel of dissipation. Our next step is to compute each term on the right-hand side of Eq. (25). The pairwise approximation for the CG force field is a well-known problem in MD, so we only need to find a workable solution from the available literature. We shall explain our method of choice in Section 3.3. On the other hand, the non-conservative forces depend intrinsically on the systemโ€™s setup and their evaluation greatly affects the simulationโ€™s success. Thus, in Section 3.2.1, we shall take some tests using a simple, a priori form of ๐‘ฒ๐ผ๐ฝ to verify the accuracy of our implementation, while leaving the evaluation of the actual memory kernel to Chapter 4.  2.3. Memory as damped oscillators At a glance, it is quite straightforward to discretize Eq. (25) as with ๐‘€ being the cutoff length for the memory kernel, that is The term ๐น๐ผ๐ฝ๐’ฌโƒ—โƒ—โƒ—โƒ—  โƒ— has been replaced with ๐น๐ผ๐ฝ๐‘…โƒ—โƒ—โƒ—โƒ—  โƒ— to denote that the random force is now truly random, drawing from a known distribution. We can determine this distribution and the explicit form of the memory kernel from the reference atomistic simulation. Thus, it seems once all the coefficients are found, the CG simulation will run smoothly.  However, as we said in the introduction, ๐‘€ must be rather large for any system of interest; otherwise, it implies a significant separation between the time scale for the random force and momentum, so we can simply extract the ๐‘‰๐ผ๐ฝโƒ—โƒ—โƒ—โƒ—  โƒ— from the inner sum in Eq. (27) and apply the Markovian approximation, reducing the MZ EOM to the standard DPD:  ๐น๐ผ๐ฝโƒ—โƒ—โƒ—โƒ—  โƒ— =  ๐น๐ผ๐ฝโƒ—โƒ—โƒ—โƒ—  โƒ—(๐‘…๐ผ๐ฝโƒ—โƒ— โƒ—โƒ—  โƒ—), ๐›ฟ๐น๐ผ๐ฝ๐’ฌโƒ—โƒ—โƒ—โƒ—  โƒ— = ๐›ฟ๐น๐ผ๐ฝ๐’ฌโƒ—โƒ—โƒ—โƒ—  โƒ—(๐‘…๐ผ๐ฝโƒ—โƒ— โƒ—โƒ—  โƒ—) . (24)  ๐‘‘๐‘‘๐‘ก๐‘ƒ๐ผโƒ—โƒ—  โƒ—(๐‘ก) = โˆ‘ [๐น๐ผ๐ฝโƒ—โƒ—โƒ—โƒ—  โƒ—(๐‘…๐ผ๐ฝโƒ—โƒ— โƒ—โƒ—  โƒ—) โˆ’ โˆซ ๐‘ฒ๐ผ๐ฝ๐‘ก0(๐‘ก โˆ’ ๐‘ )๐‘‰๐ผ๐ฝโƒ—โƒ—โƒ—โƒ—  โƒ—(๐‘ )๐‘‘๐‘  + ๐›ฟ๐น๐ผ๐ฝ๐’ฌโƒ—โƒ—โƒ—โƒ—  โƒ—(๐‘…๐ผ๐ฝโƒ—โƒ— โƒ—โƒ—  โƒ—)] ๐ฝโ‰ ๐ผ, (25)  ๐‘ฒ๐ผ๐ฝ(๐‘ก) = ๐›ฝ โŒฉ๐›ฟ๐น๐ผ๐ฝ๐“ โƒ—โƒ—โƒ—โƒ—  โƒ—(๐‘ก) โŠ— ๐›ฟ๐น๐ผ๐ฝ๐“ โƒ—โƒ—โƒ—โƒ—  โƒ—(0)โŒช (26)  ฮ”๐‘ƒ๐ผโƒ—โƒ—  โƒ—ฮ”๐‘ก= โˆ‘ [โˆ’๐œ•๐œ•๐‘…๐ผ๐ฝโƒ—โƒ— โƒ—โƒ—  โƒ—๐‘ˆ (๐‘…๐‘โƒ—โƒ— โƒ—โƒ—  โƒ—) โˆ’ ๐›ฝ โˆ‘ ๐‘ฒ๐ผ๐ฝ(๐‘šฮ”๐‘ก)๐‘‰๐ผ๐ฝโƒ—โƒ—โƒ—โƒ—  โƒ—(๐‘ก โˆ’ ๐‘šฮ”๐‘ก)ฮ”๐‘ก๐‘€๐‘š=0+๐น๐ผ๐ฝ๐‘…โƒ—โƒ—โƒ—โƒ—  โƒ—(๐‘ก)โˆšฮ”๐‘ก] ๐ฝโ‰ ๐ผ , (27)  ๐‘ฒ๐ผ๐ฝ(๐‘šฮ”๐‘ก) = 0 for ๐‘š > ๐‘€ . (28) 11  in which ๐›พ๐ผ๐ฝ is the strength of the friction force. Moreover, the N-body function is very difficult to work with [2, 13], so Eq. (27) still needs some further decomposition, such as pairwise decomposition, to be practical. Since a pairwise memory kernel scales as ๐‘‚(๐‘€๐‘2), a direct application of Eq. (17) would consume a lot of the computerโ€™s memory and run inefficiently due to the large number of cache misses [26]. There exists an alternative approach that mitigates the complexity of the MZ formalism: instead of keeping track of all the relevant historical terms, we supplement the system with multiple series of ๐‘˜ auxiliary variables ๐‘ ๐ผ๐ฝโƒ—โƒ— โƒ—โƒ—  =(๐‘ 1, ๐‘ 2, โ€ฆ ๐‘ ๐‘˜), whose coupling with the relative momentum ๐‘ƒ๐ผ๐ฝโƒ—โƒ—โƒ—โƒ—  โƒ— mimics the memory effect. The coupled equation is given as [21] where ๐‘จ๐‘๐‘  = โˆ’๐‘จ๐‘ ๐‘๐‘‡  and ๐œ‰๐ผ๐ฝโƒ—โƒ— โƒ—โƒ—   is a vector of ๐‘˜ Gaussian random numbers with โŒฉ๐œ‰๐ผ๐ฝโƒ—โƒ— โƒ—โƒ—  (๐‘ก) โŠ— ๐œ‰๐พ๐ฟโƒ—โƒ— โƒ—โƒ—โƒ—โƒ—  (0)โŒช = ๐›ฟ(๐‘ก)๐›ฟ๐ผ๐ฝ,๐พ๐ฟ. This is the so-called Ornstein-Uhlenbeck process [24]. By assuming again that the memory length is finite, we can take ๐‘ ๐ผโƒ—โƒ—โƒ—  (โˆ’โˆž) = 0โƒ—  and an ansatz for ๐‘ ๐ผ๐ฝโƒ—โƒ— โƒ—โƒ—   as Using Eq. (31), we eliminate the auxiliary terms from Eq. (30) and arrive at in which We can see that the form of the first equation functions similarly to the MZ EOM of Eq. (17). In Appendix B, we shall present the proof of equivalence. For now, it should be noted that Eq. (32) reduces to the pairwise-decomposed Eq. (25) with some specific choices of ๐‘จ and ๐‘ฉ. In particular, the relation is sufficient to satisfy the second FDT. We can always solve Eq. (34) to find ๐‘ฉ๐‘ ๐‘  using Cholesky factorization if (๐‘จ๐‘ ๐‘  + ๐‘จ๐‘ ๐‘ ๐‘‡ ) is positive-semidefinite [24]. In fact, when ๐‘จ๐‘ ๐‘  is real and has complex eigenvalues with positive real part, ๐‘ฒ(๐‘ก) becomes a linear combination of damped harmonic oscillators [21]. Thus, for any decaying memory kernel  โˆ‘ โŒฉ๐›ฟ๐น๐ผ๐ฝ๐‘…โƒ—โƒ—โƒ—โƒ—  โƒ—(๐‘ก โˆ’ ๐‘šฮ”๐‘ก) โŠ— ๐›ฟ๐น๐ผ๐ฝ๐‘…โƒ—โƒ—โƒ—โƒ—  โƒ—(๐‘šฮ”๐‘ก)โŒช ๐‘‰๐ผ๐ฝโƒ—โƒ—โƒ—โƒ—  โƒ—(๐‘ก โˆ’ ๐‘šฮ”๐‘ก)ฮ”๐‘ก๐‘€๐‘š=0|๐‘€โ‰ช1โ‰ˆ ๐›พ๐ผ๐ฝ๐‘‰๐ผ๐ฝโƒ—โƒ—โƒ—โƒ—  โƒ—(๐‘ก) , (29)  (๐‘ƒ๐ผ๐ฝโƒ—โƒ—โƒ—โƒ—  โƒ—ฬ‡๐‘ ๐ผ๐ฝโƒ—โƒ— โƒ—โƒ—  ฬ‡) = (๐น๐ผ๐ฝโƒ—โƒ—โƒ—โƒ—  โƒ—0โƒ— ) โˆ’ (0 ๐‘จ๐‘๐‘ ๐‘จ๐‘ ๐‘ ๐‘จ๐‘ ๐‘ ) (๐‘ƒ๐ผ๐ฝโƒ—โƒ—โƒ—โƒ—  โƒ—๐‘ ๐ผ๐ฝโƒ—โƒ— โƒ—โƒ—  ) + (0 00 ๐‘ฉ๐‘ ๐‘ ) (0โƒ— ๐œ‰๐ผ๐ฝโƒ—โƒ— โƒ—โƒ—  ) , (30)  ๐‘ ๐ผ๐ฝโƒ—โƒ— โƒ—โƒ—  (๐‘ก) = โˆซ ๐‘’โˆ’(๐‘กโˆ’๐‘กโ€ฒ)๐‘จ๐‘ ๐‘  (โˆ’๐‘จ๐‘ ๐‘๐‘ƒ๐ผ๐ฝโƒ—โƒ—โƒ—โƒ—  โƒ—(๐‘กโ€ฒ) + ๐‘ฉ๐‘ ๐‘ ๐œ‰๐ผ๐ฝโƒ—โƒ— โƒ—โƒ—  (๐‘กโ€ฒ)) ๐‘‘๐‘กโ€ฒ๐‘กโˆž . (31)  ๐‘ƒ๐ผ๐ฝโƒ—โƒ—โƒ—โƒ—  โƒ—ฬ‡ (๐‘ก) = ๐น๐ผ๐ฝโƒ—โƒ—โƒ—โƒ—  โƒ—(๐‘ก) โˆ’ โˆซ ๐‘ฒ(๐‘ก โˆ’ ๐‘กโ€ฒ)๐‘ƒ๐ผ๐ฝโƒ—โƒ—โƒ—โƒ—  โƒ—(๐‘กโ€ฒ)๐‘‘๐‘กโ€ฒ๐‘กโˆ’โˆž+ ๐œ๐ผ๐ฝโƒ—โƒ— โƒ—โƒ—  (๐‘ก) , (32)  ๐‘ฒ(๐‘ก) = โˆ’๐‘จ๐‘๐‘ ๐‘’โˆ’๐‘ก๐‘จ๐‘ ๐‘ ๐‘จ๐‘ ๐‘ , ๐œ๐ผ๐ฝโƒ—โƒ— โƒ—โƒ—  (๐‘ก) = โˆซ ๐‘จ๐‘๐‘ ๐‘’โˆ’(๐‘กโˆ’๐‘กโ€ฒ)๐‘จ๐‘ ๐‘ ๐‘ฉ๐‘ ๐‘ ๐œ‰๐ผ๐ฝโƒ—โƒ— โƒ—โƒ—  (๐‘กโ€ฒ)๐‘‘๐‘กโ€ฒ๐‘กโˆ’โˆž . (33)  ๐‘ฉ๐‘ ๐‘ ๐‘ฉ๐‘ ๐‘ ๐‘‡ = ๐‘˜๐ต๐‘‡(๐‘จ๐‘ ๐‘  + ๐‘จ๐‘ ๐‘ ๐‘‡ ) (34) 12  in Eq. (25), we can approximate it with a series of damped oscillators and convert the EOM to Eq. (30). Section 3.2.3 will clarify the practical implementation of this method.  2.4. Integration of dissipative particle dynamics So, we have seen that the EOM resulting from applying the MZ formalism is closely related to the DPD equation, of which the general form was presented in Eq. (2). In practice, the forces are typically assumed to be central with the following functional form [20] where ๐‘Ÿ๐ผ๐ฝ = |๐‘Ÿ๐ผโƒ—โƒ— โˆ’ ๐‘Ÿ๐ฝโƒ—โƒ— |, ๐‘’๐ผ๐ฝโƒ—โƒ— โƒ—โƒ—  = ๐‘Ÿ๐ผ๐ฝโƒ—โƒ—โƒ—โƒ— /๐‘Ÿ๐ผ๐ฝ, ๐‘ฃ๐ผ๐ฝโƒ—โƒ—โƒ—โƒ—  โƒ— = ๐‘ฃ๐ผโƒ—โƒ—  โƒ— โˆ’ ๐‘ฃ๐ฝโƒ—โƒ—  โƒ—, ๐›พ and ๐œŽ are the strength of the dissipative and random forces, respectively. ๐œ‰๐ผ๐ฝ are standard normal random variables, uncorrelated for different ๐ผ๐ฝ pairs and different times, with ๐œ‰๐ผ๐ฝ = ๐œ‰๐ฝ๐ผ. The justification for the last property is the conservation of momentum, while the consequence of the former shall be explored later in this section. By convention, ๐น๐ผ๐ฝ๐ถโƒ—โƒ—โƒ—โƒ—  โƒ— can include any force term appropriate to the model, as long as it is conservative and independent of the other two. By contrast, the dissipative and random forces are coupled through a fluctuation-dissipation relation, which gives [20] in our particular case. This formula ensures the simulation samples the canonical distribution at the correct temperature. Once we have specified the explicit form for each force term, we can determine their value at the current point in phase space and numerically propagate our position to the next timestep. At the lowest order of approximation, the velocity-Verlet (VV) integration is the most popular scheme, due to simplicity, stability, and efficiency [27]. The shortened implementation of the VV algorithm is [28] However, the calculation of the second equation is nontrivial: ๐‘ฃ (๐‘ก + ฮ”๐‘ก) appears on both sides, explicitly on the left-hand side and implicitly inside ๐น (๐‘ก + ฮ”๐‘ก). Since the friction force in DPD always contain some dependency on velocity, it is thus impossible to accurately integrate Eq. (2) with the VV algorithm [1]. In order to obtain the correct result, a self-consistent solution is required; i.e. the value for ๐‘ฃ (๐‘ก + ฮ”๐‘ก) calculated on the left-hand side must be the  ๐น๐ผ๐ฝ๐ถโƒ—โƒ—โƒ—โƒ—  โƒ— = ๐น๐ผ๐ฝ๐ถ(๐‘Ÿ๐ผ๐ฝ)๐‘’๐ผ๐ฝโƒ—โƒ— โƒ—โƒ—   , ๐น๐ผ๐ฝ๐ทโƒ—โƒ—โƒ—โƒ—  โƒ— = โˆ’๐›พ(๐‘Ÿ๐ผ๐ฝ)(๐‘ฃ๐ผ๐ฝโƒ—โƒ—โƒ—โƒ—  โƒ— โ‹… ๐‘’๐ผ๐ฝโƒ—โƒ— โƒ—โƒ—  )๐‘’๐ผ๐ฝโƒ—โƒ— โƒ—โƒ—   , ๐น๐ผ๐ฝ๐‘…โƒ—โƒ—โƒ—โƒ—  โƒ— = ๐œŽ(๐‘Ÿ๐ผ๐ฝ)๐œ‰๐ผ๐‘—๐‘’๐ผ๐ฝโƒ—โƒ— โƒ—โƒ—   , (35)  ๐›พ(๐‘Ÿ๐ผ๐ฝ) =12๐›ฝ[๐œŽ(๐‘Ÿ๐ผ๐ฝ)]2 (36)  ๐‘ฅ (๐‘ก + ฮ”๐‘ก) = ๐‘ฅ (๐‘ก) + ๐‘ฃ (๐‘ก)ฮ”๐‘ก +12๐‘š๐น (๐‘ก)ฮ”๐‘ก2 , ๐‘ฃ (๐‘ก + ฮ”๐‘ก) = ๐‘ฃ (๐‘ก) +12๐‘š(๐น (๐‘ก) + ๐น (๐‘ก + ฮ”๐‘ก)) ฮ”๐‘ก . (37) 13  same as the one used to find ๐น (๐‘ก + ฮ”๐‘ก). This self-consistency guarantees time-reversibility, which in turn is necessary for detailed balance and a correct canonical sampling. Moreover, although the discussion above applies to the VV algorithm, most of the known MD integration schemes experience the same problem [20], while the self-consistent ones, such as the predictor-corrector method [1], are too clunky for our purpose. Therefore, alternative integration strategies are in order. We shall explore some of them in Section 3.1. There is another subtlety in the DPD force calculation. As a force describing the rapidly fluctuating interactions between the atoms, ๐น ๐ผ๐ฝ๐‘… is assumed to have a zero mean and to be uncorrelated with itself for any time scale relevant to the study [29], that is with ๐ถ being a scaling constant. Thus, we can divide the integral of ๐น ๐ผ๐ฝ๐‘… into smaller terms, each spanning over some time ๐›ฟ๐‘ก and drawing from the same distribution: Applying the central limit theorem on the right-hand side, we discover that the integral on the left-hand side also has a mean of zero and follows a distribution with a variance proportional to ๐‘ก๐›ฟ๐‘ก. This means when propagating the velocity, the impulse due to the random force must be scaled with โˆš๐›ฟ๐‘ก:    โŒฉ๐น ๐ผ๐ฝ๐‘…(๐‘ก)โŒช = 0 , โŒฉ๐น ๐ผ๐ฝ๐‘…(0) โ‹… ๐น ๐ผ๐ฝ๐‘…(๐‘ก)โŒช = ๐ถ๐›ฟ(๐‘ก) , (38)  โˆซ ๐น ๐ผ๐ฝ๐‘…(๐‘ )๐‘‘๐‘ ๐‘ก0= โˆ‘ โˆซ ๐น ๐ผ๐ฝ๐‘…(๐‘š๐›ฟ๐‘ก)๐‘‘๐‘ (๐‘š+1)๐›ฟ๐‘ก๐‘š๐›ฟ๐‘ก๐‘ก/๐›ฟ๐‘กโˆ’1๐‘š=0 . (39)  ฮ”๐‘ฃ๐ผโƒ—โƒ—  โƒ— =1๐‘€๐ผ(โˆ‘๐น๐ผ๐ฝ๐ถโƒ—โƒ—โƒ—โƒ—  โƒ—ฮ”๐‘ก + ๐น๐ผ๐ฝ๐ทโƒ—โƒ—โƒ—โƒ—  โƒ—ฮ”๐‘ก + ๐น๐ผ๐ฝ๐‘…โƒ—โƒ—โƒ—โƒ—  โƒ—๐ผโ‰ ๐ฝโˆšฮ”๐‘ก) . (40) 14  Chapter 3 Implementation    3.1. Integrators As we have seen previously, there are two complications with the practical application of a numerical DPD integrator. The second one is a superficial problem that can be solved by just scaling the random force with the appropriate factor. However, the first problem is much more fundamental; we cannot solve it with a simple insertion of the new forces in the force calculation step. In this section, we shall investigate the behavior of different DPD algorithms. The test we perform is very similar to what has been done by Nikunen et al. [20]. The reinvestigation is necessary as we have to implement our code, due to a lack of correct DPD integrators in popular MD programs.  3.1.1. Dissipative velocity-Verlet The dissipative velocity-Verlet (DPD-VV) integrator uses a very simple modification of the standard VV version, namely [30] in which we use ๐น๐ถโƒ—โƒ— โƒ—โƒ—  (๐‘ก) = ๐น๐ถโƒ—โƒ— โƒ—โƒ—  (?โƒ—? (๐‘ก)) and ๐น๐‘…โƒ—โƒ—โƒ—โƒ— (๐‘ก) = ๐น๐‘…โƒ—โƒ—โƒ—โƒ— (?โƒ—? (๐‘ก)). Clearly, without the last calculation, Eq. (41) just describes the traditional VV algorithm in DPD. The difference happens in the last line, where the dissipative force is recomputed using the velocity at the new time with ?โƒ—? (๐‘ก + ฮ”๐‘ก)  adjusted accordingly. In short, this is an iterative approximation to the true solution of the DPD EOM. Obviously, we should not expect that the updated velocity in DPD-VV is always a good approximation, since it only goes over a single iteration. To achieve self-consistency, we can repeat the last line of Eq. (41), in which ๐‘‰โ€ฒโƒ—โƒ—โƒ—โƒ— (๐‘ก + ฮ”๐‘ก) is replaced by the velocity obtained from the previous iteration, until we reach the desired level of convergence. We denote this variant as self-consistent velocity-Verlet (SC-VV).  ๐‘‰๐ผโƒ—โƒ—  โƒ— (๐‘ก +ฮ”๐‘ก2) = ๐‘‰๐ผโƒ—โƒ—  โƒ—(๐‘ก) +12๐‘€๐ผ[๐น๐ผ๐ถโƒ—โƒ— โƒ—โƒ—  (๐‘ก)ฮ”๐‘ก + ๐น๐ผ๐ทโƒ—โƒ—โƒ—โƒ—  โƒ— (๐‘…๐ผโƒ—โƒ—โƒ—โƒ— (๐‘ก), ๐‘‰๐ผโƒ—โƒ—  โƒ—(๐‘ก)) ฮ”๐‘ก + ๐น๐ผ๐‘…โƒ—โƒ—โƒ—โƒ—  โƒ—(๐‘ก)โˆšฮ”๐‘ก] , ๐‘…๐ผโƒ—โƒ—โƒ—โƒ— (๐‘ก + ฮ”t) = ๐‘…๐ผโƒ—โƒ—โƒ—โƒ— (๐‘ก) + ๐‘‰๐ผโƒ—โƒ—  โƒ— (๐‘ก +ฮ”๐‘ก2) ฮ”๐‘ก , ๐‘‰๐ผโ€ฒโƒ—โƒ—โƒ—โƒ— (๐‘ก + ฮ”๐‘ก) = ๐‘‰๐ผโƒ—โƒ—  โƒ— (๐‘ก +ฮ”๐‘ก2) +12๐‘€๐ผ[๐น๐ผ๐ถโƒ—โƒ— โƒ—โƒ—  (๐‘ก + ฮ”๐‘ก)ฮ”๐‘ก + ๐น๐ผ๐ทโƒ—โƒ—โƒ—โƒ—  โƒ— (๐‘…๐ผโƒ—โƒ—โƒ—โƒ— (๐‘ก + ฮ”๐‘ก), ๐‘‰๐ผโƒ—โƒ—  โƒ— (๐‘ก +ฮ”๐‘ก2))ฮ”๐‘ก+ ๐น๐ผ๐‘…โƒ—โƒ—โƒ—โƒ—  โƒ—(๐‘ก + ฮ”๐‘ก)โˆšฮ”๐‘ก] , ๐‘‰๐ผโƒ—โƒ—  โƒ—(๐‘ก + ฮ”๐‘ก) = ๐‘‰๐ผโ€ฒโƒ—โƒ—โƒ—โƒ— (๐‘ก + ฮ”๐‘ก) +12๐‘€๐ผ๐น๐ผ๐ทโƒ—โƒ—โƒ—โƒ—  โƒ— (๐‘…๐ผโƒ—โƒ—โƒ—โƒ— (๐‘ก + ฮ”๐‘ก), ๐‘‰๐ผโ€ฒโƒ—โƒ—โƒ—โƒ— (๐‘ก + ฮ”๐‘ก)) ฮ”๐‘ก , (41) 15  Nevertheless, in practice, DPD-VV has been shown to give acceptable results while still being computationally efficient [20].  3.1.2. Lowe-Andersen Like DPD-VV, this method (Lowe) is also a modified version of the standard VV integrator, but utilizes a very different idea [31]. where ๐‘’๐ผ๐ฝโƒ—โƒ— โƒ—โƒ—   is the radial unit vector between ๐ผ and ๐ฝ, ๐œ‡๐ผ๐ฝ =๐‘€๐ผ๐‘€๐ฝ๐‘€๐ผ+๐‘€๐ฝ the reduced mass of ๐ผ and ๐ฝ, ๐œ‰๐ผ๐ฝโƒ—โƒ— โƒ—โƒ—   a vector of normalized Gaussian random numbers, independent for each ๐ผ๐ฝ pair, and ๐‘ƒ the selection probability. We need to highlight the sequential, i.e. non-parallelizable, characteristic of the indented block in Eq. (42): the velocity for the current pair of particles must be updated before moving to the next pair, which implies the values of ฮ”๐‘‰๐ผ๐ฝโƒ—โƒ—โƒ—โƒ—  โƒ— are not independent of each other. However, for each pair of particles, we execute the chain of actions with independent probability ๐‘ƒ, meaning it is possible that one pair always goes through the treatment while another never. Furthermore, ๐‘ƒ โ‰ช 1 weakly couples the system to the heat bath, while thermalization happens at every timestep in the limiting case of ๐‘ƒ = 1. Therefore, its dynamical properties can be adjusted by the appropriate choice of ๐‘ƒ. Finally, we note that Lowe only has the DPD version; for an LD system, it is simplified to the Andersen thermostat [32], hence the name.  3.1.3. Impulsive leap-frog This section presents the impulsive leap-frog algorithm (Imp-LF) as described by Goga et al. [33]. Conceptually, it is the leap-frogged version of the Andersen thermostat, that is    ๐‘‰๐ผโƒ—โƒ—  โƒ— (๐‘ก +ฮ”๐‘ก2) = ๐‘‰๐ผโƒ—โƒ—  โƒ— (๐‘ก โˆ’ฮ”๐‘ก2) +1๐‘€๐ผ๐น๐ผ๐ถโƒ—โƒ— โƒ—โƒ—  (๐‘ก)ฮ”๐‘ก , ๐‘…๐ผโƒ—โƒ—โƒ—โƒ— (๐‘ก + ฮ”๐‘ก) = ๐‘…๐ผโƒ—โƒ—โƒ—โƒ— (๐‘ก) + ๐‘‰๐ผโƒ—โƒ—  โƒ— (๐‘ก +ฮ”๐‘ก2)ฮ”๐‘ก , ๐‘‰๐ผโƒ—โƒ—  โƒ— (๐‘ก +ฮ”๐‘ก2) = ๐‘‰๐ผโƒ—โƒ—  โƒ— (๐‘ก +ฮ”๐‘ก2) +1๐‘€๐ผ๐น๐ผ๐ถโƒ—โƒ— โƒ—โƒ—  (๐‘ก + ฮ”๐‘ก)ฮ”๐‘ก . For every pair ๐ผ, ๐ฝ, do the following sequentially with probability ๐‘ƒ: ฮ”๐‘‰๐ผ๐ฝโƒ—โƒ—โƒ—โƒ—  โƒ— (๐‘ก +ฮ”๐‘ก2) = ๐œ‰๐ผ๐ฝโƒ—โƒ— โƒ—โƒ—  (๐‘ก)โˆš12๐›ฝ๐œ‡๐ผ๐ฝโˆ’12๐‘‰๐ผ๐ฝโƒ—โƒ—โƒ—โƒ—  โƒ— (๐‘ก +ฮ”๐‘ก2) , ๐‘‰๐ผโƒ—โƒ—  โƒ— (๐‘ก +ฮ”๐‘ก2) = ๐‘‰๐ผโƒ—โƒ—  โƒ— (๐‘ก +ฮ”๐‘ก2) +๐œ‡๐ผ๐ฝ๐‘š๐ผฮ”๐‘‰๐ผ๐ฝโƒ—โƒ—โƒ—โƒ—  โƒ— (๐‘ก +ฮ”๐‘ก2) , ๐‘‰๐ฝโƒ—โƒ—  โƒ— (๐‘ก +ฮ”๐‘ก2) = ๐‘‰๐ฝโƒ—โƒ—  โƒ— (๐‘ก +ฮ”๐‘ก2) โˆ’๐œ‡๐ผ๐ฝ๐‘š๐ฝฮ”๐‘‰๐ผ๐ฝโƒ—โƒ—โƒ—โƒ—  โƒ— (๐‘ก +ฮ”๐‘ก2) , (42) 16  For the uncorrelated system, the velocity reduction factor ๐‘“ can be chosen to be anything between 0 and 1, and may differ for different particles.  In this case, the limiting condition ๐‘“ = 1 reduces the algorithm to the Andersen thermostat. In addition, the Imp-LF algorithm for the LD system is very similar to Lowe. However, to provide the correct dynamics, ๐‘“ must be equal to ๐‘’๐‘ฅ๐‘(โˆ’12๐›ฝฮ”๐‘ก) and the value of ๐‘ƒ is chosen according to the actual form of the dissipative force.  Compared to the two integrators we have seen, Imp-LF employs a much more convoluted modification, so it is not immediately clear how this algorithm yields dynamically equivalent characteristics. Nevertheless, by applying Taylorโ€™s expansion on ๐‘“ = ๐‘’๐‘ฅ๐‘(โˆ’12๐›ฝฮ”๐‘ก) , we can at least confirm that in the small ฮ”๐‘ก limit, the resulting impulse conforms to the requirement mentioned at the end of Section 2.4, i.e. ฮ”๐‘‰๐ผ๐ฝโƒ—โƒ—โƒ—โƒ—  โƒ— โ‰ˆ โˆ’12๐›ฝ๐‘‰๐ผ๐ฝโƒ—โƒ—โƒ—โƒ—  โƒ—ฮ”๐‘ก +1โˆš๐œ‡๐ผ๐ฝ๐œ‰๐ผ๐ฝโˆš๐›ฅ๐‘ก satisfies the scaling condition of Eq. (40).  3.1.4. Test result 3.1.4.1. Model systems We tested the integrators on two models that capture some basic properties of a fluid of identical particles. The first one is an LD system without the conservative force ๐น๐ผ๐ถ = 0, so that the EOM becomes  ๐‘‰๐ผโƒ—โƒ—  โƒ— (๐‘ก +๐›ฅ๐‘ก2) = ๐‘‰๐ผโƒ—โƒ—  โƒ— (๐‘ก โˆ’๐›ฅ๐‘ก2) +1๐‘€๐ผ๐น๐ผ๐ถโƒ—โƒ— โƒ—โƒ—  (๐‘ก)๐›ฅ๐‘ก . For the LD system: ๐‘‰๐ผโƒ—โƒ—  โƒ— (๐‘ก +๐›ฅ๐‘ก2) = (1 โˆ’ ๐‘“)๐‘‰๐ผโƒ—โƒ—  โƒ— (๐‘ก โˆ’๐›ฅ๐‘ก2) + โˆš๐‘“(2 โˆ’ ๐‘“)๐‘˜๐ต๐‘‡๐‘€๐ผ๐œ‰ . For the DPD system: For every pair ๐ผ, ๐ฝ, do the following sequentially with probability ๐‘ƒ: ๐›ฅ๐‘‰๐ผ๐ฝโƒ—โƒ—โƒ—โƒ—  โƒ—(๐‘ก) = โˆ’(๐‘“ โˆ’ 1)๐‘‰๐ผ๐ฝโƒ—โƒ—โƒ—โƒ—  โƒ— (๐‘ก +๐›ฅ๐‘ก2) + โˆš1 โˆ’ ๐‘“2๐›ฝ๐œ‡๐ผ๐ฝ๐œ‰๐ผ๐ฝโƒ—โƒ— โƒ—โƒ—  (๐‘ก) , ๐‘‰๐ผโƒ—โƒ—  โƒ— (๐‘ก +๐›ฅ๐‘ก2) = ๐‘‰๐ผโƒ—โƒ—  โƒ— (๐‘ก +๐›ฅ๐‘ก2) +๐œ‡๐ผ๐ฝ๐‘š๐ผ๐›ฅ๐‘‰๐ผ๐ฝโƒ—โƒ—โƒ—โƒ—  โƒ—(๐‘ก) , ๐‘‰๐ฝโƒ—โƒ—  โƒ— (๐‘ก +๐›ฅ๐‘ก2) = ๐‘‰๐ฝโƒ—โƒ—  โƒ— (๐‘ก +๐›ฅ๐‘ก2) โˆ’๐œ‡๐ผ๐ฝ๐‘š๐ฝ๐›ฅ๐‘‰๐ผ๐ฝโƒ—โƒ—โƒ—โƒ—  โƒ—(๐‘ก) . ๐‘…๐ผโƒ—โƒ—โƒ—โƒ— (๐‘ก + ฮ”๐‘ก) = ๐‘…๐ผโƒ—โƒ—โƒ—โƒ— (๐‘ก) + ๐‘‰๐ผโƒ—โƒ—  โƒ— (๐‘ก +ฮ”๐‘ก2)ฮ”๐‘ก . (43)  ๐‘‘๐‘‘๐‘ก๐‘ƒ๐ผโƒ—โƒ—  โƒ—(๐‘ก) = โˆ’๐›พ๐‘‰๐ผโƒ—โƒ—  โƒ—(๐‘ก) + ๐œŽ๐œ‰๐ผ(๐‘ก) , (44) 17  with ๐›พ and ๐œŽ being the strength of the damping and random forces, having a constant value as the particles do not interact with each other. Thus, this simulates ๐‘ independent particles, essentially describing the motion of large, non-interacting particles in a solvent. This simple system provides a quick preliminary test for the correctness of our implementation. Once the basic expected results had been verified, we tried the integrators on the DPD model, with the conservative force being zero again and the dissipative force applied to each pairwise velocity ๐‘‰๐ผ๐ฝโƒ—โƒ—โƒ—โƒ—  โƒ— for in which ๐‘…๐ถ is the cutoff distance for the correlation function and ๐›พ the scale factor. Thus, the EOM for this model is The soft-repulsive dissipative force is a very common choice in DPD simulations and has been employed in various investigations for integration schemes [1, 34, 20], thus allowing comparable results.  3.1.4.2. Physical quantities Our survey was chiefly concerned with the average kinetic temperature whose convergence is one of the main conditions for reliable simulations in the canonical ensemble. Here, ๐‘‘ = 3 is the dimension of the system. However, having the correct temperature only validates the macroscopic properties of the system. In order to verify the accuracy of its microscopic behaviors, we need to consider the mean square displacement (MSD) and the velocity autocorrelation function ๐ถ(๐‘ก) For the LD system, theoretically [29], we expect that the MSD grows linearly at large ๐‘ก, with while the VACF decays exponentially over time:  ๐น๐ทโƒ—โƒ—โƒ—โƒ—  โƒ— = โˆ‘ ๐น๐ผ๐ฝ๐ทโƒ—โƒ—โƒ—โƒ—  โƒ—๐ฝโ‰ ๐ผ= โˆ‘โˆ’๐›พ(๐‘…๐ผ๐ฝ)๐‘‰๐ผ๐ฝโƒ—โƒ—โƒ—โƒ—  โƒ—(๐‘ก)๐ฝโ‰ ๐ผ , (45)  ๐›พ(๐‘…๐ผ๐ฝ) = {๐›พ (1 โˆ’๐‘…๐ผ๐ฝ๐‘…๐ถ)2   for ๐‘…๐ผ๐ฝ โ‰ค ๐‘…๐‘0   otherwise , (46)  ๐‘‘๐‘‘๐‘ก๐‘ƒ๐ผโƒ—โƒ—  โƒ— = โˆ‘ โˆ’๐›พ๐ฝโ‰ ๐ผ (๐‘…๐ผ๐ฝ)๐‘‰๐ผ๐ฝโƒ—โƒ—โƒ—โƒ—  โƒ—(๐‘ก) + ๐œŽ(๐‘…๐ผ๐ฝ)๐œ‰๐ผ๐ฝโƒ—โƒ— โƒ—โƒ—  (๐‘ก) . (47)  โŒฉ๐‘˜๐ต๐‘‡โŒช =1๐‘ โˆ’ 1๐‘€๐‘‘โˆ‘ ๐‘‰๐ผ2๐ผ , (48)  ๐‘€๐‘†๐ท(๐‘ก) =1๐‘โˆ‘(๐‘‹๐ผ(๐‘ก) โˆ’ ๐‘‹๐ผ(0))2๐ผ , ๐ถ(๐‘ก) =1๐‘โˆ‘โŸจ๐‘‰๐ผโƒ—โƒ—  โƒ—(๐‘ก) โ‹… ๐‘‰๐ผโƒ—โƒ—  โƒ—(0)โŸฉ๐ผ . (49)  ๐‘€๐‘†๐ท โ†’2๐‘‘2๐‘˜๐ต๐‘‡๐›พ๐‘ก , (50) 18  On the other hand, we do not have the algebraic expression for these two quantities in the DPD model, but we do know they exhibit the same behavior mentioned above, namely linearity for the MSD and exponential decay for the VCAF, under the condition [35] Finally, for the DPD system, we also compute the radial distribution function (RDF) of which the result has been investigated by Nikunen et al. [20]. Naturally, the RDF is unity for a system with no conservative force, because it is structureless.  3.1.4.3. Simulation conditions This experiment employed the DPD reduced units, i.e. ๐‘… = ๐‘…โˆ—/๐‘…๐ถ, ๐‘‡ = ๐‘‡โˆ—/๐›พ and ๐‘ก = ๐‘กโˆ—๐‘…๐ถ(๐‘€/๐›พ)1/2, with the asterisk denoting values in real units. The simulations were all performed with ๐‘ = 4096 identical, non-interacting particles of ๐‘€ = 1.0, in periodic cubic boxes of size 10 ร—  10 ร—  10, giving a number density of ๐œŒ โ‰ˆ 4.1. The energy scale is defined by setting the temperature to unity ๐‘˜๐ต๐‘‡ = 1.0. The damping factor is chosen to be 0.18 in both cases, so the strength of the random force can be evaluated using Eq. (36), specifically ๐œŽ = 0.6 for the LD system and ๐œŽ(๐‘…๐ผ๐ฝ) = 0.6(1 โˆ’ ๐‘…๐ผ๐ฝ/๐‘…๐ถ)๐ป(๐‘…๐ถ โˆ’ ๐‘…๐ผ๐ฝ) for the DPD system, with ๐ป(๐‘ฅ) being the Heaviside step function.  3.1.4.4. Result for Langevin dynamics In Fig. 3-1, we plot the average kinetic temperature for the DPD-VV, Imp-LF and standard VV integrators for the LD system. Obviously, as the Lowe integrator reduces to the Andersen thermostat in this case, we shall not consider its result. We can see that all three integrators give temperature fluctuating around the correct value of ๐‘˜๐ต๐‘‡ =1.0. However, the VV simulation produced thermal fluctuations visibly larger than those by the DPD integrators, hinting at a divergence of the VV result from the other two.  ๐ถ(๐‘ก) =๐‘˜๐ต๐‘‡๐‘€๐‘’โˆ’๐›พ๐‘ก/๐‘€ . (51)  ๐›พ๐œŒโˆ’1๐‘‘๐‘‘โˆš๐‘˜๐ต๐‘‡/๐‘€โ‰ช 1 . (52)  ๐บ(๐‘…) =1๐‘(๐‘ โˆ’ 1)โˆ‘โŒฉ๐›ฟ(๐‘…๐ผ๐ฝ โˆ’ ๐‘…)โŒช4๐œ‹๐‘…2๐ผ,๐ฝโ‰ ๐ผ , (53) 19   Fig. 3-1: Time evolution of the average kinetic temperature โŒฉ๐’Œ๐‘ฉ๐‘ปโŒช for different integrators in LD model. As mentioned in Section 3.1.2, the Lowe integrator only has the pairwise form, so its result is not presented here. In fact, the discrepancy of the standard VV integrator becomes obvious in Fig. 3-2, which shows the computed MSD. For each result, we fit the values from ๐‘ก = 20.0 with linear regression. We find that both DPD integrators yield very identical MSD lines, the asymptotic behavior of which conforms with the theoretical formula, namely ๐‘€๐‘†๐ท โ†’ 100.0๐‘ก for our specified choice of variables, while the corresponding result of the standard VV algorithm clearly deviates from the other two.  Fig. 3-2: Time evolution of the mean square displacement for different integrators in LD model. Using linear regression, we fit the data from ๐‘ก = 20.0 of the Imp-LF and VV integrators to two lines, whose equations are displayed on the graph. 20   Fig. 3-3: The decay of velocity autocorrelation function over time for different integrators in LD model. Unlike the previous plot, the trendline is simply here as visual help. Finally, Fig. 3-3 plots ln(๐ถ(๐‘ก)/๐ถ(0)) against time. Here, the two trendlines are not fitted to the data and are simply there for visual support, unlike Fig. 3-2. Akin to the previous result, the two DPD integrators generate data compatible with the theory, which dictates the logarithm of the VACF to be linear in time with a slope of โˆ’0.18. Likewise, the behavior of the VV algorithm is displayed again to highlight the discrepancy of the standard version, thus justifying the necessity of DPD integrators. Moreover, we observe that the fluctuations at the tail of the two plots are expected too, as the VACF approaches 0 and starts to be dominated by thermal noises.  3.1.4.5. Result for dissipative particle dynamics Next, we consider the DPD model, where we took ๐น๐ผ๐ฝ๐ถ = 0. Akin to the above discussion, our interest was focused on the three quantities โŒฉ๐‘˜๐ต๐‘‡โŒช, MSD and VACF as a function of time. Here, we must remember that the Lowe and Imp-LF integrators are fundamentally different than DPD-VV, because the dynamics of the former cases can be controlled by choosing the appropriate collision probability ๐‘ƒ, while the latterโ€™s effective friction rate depends solely on ๐›พ(๐‘…๐ผ๐ฝ). Thus, the value of ๐‘ƒ was chosen such that the simulations for Lowe and Imp-LF gave an effective friction rate similar to that of DPD-VV, namely ๐‘ƒ = 0.037 for Lowe and ๐‘ƒ = 2๐›พ(๐‘…๐ผ๐ฝ)/๐›พ(0) for Imp-LF [33]. Fig. 3-4 exhibits the average kinetic temperature, which should fluctuate around ๐‘˜๐ต๐‘‡ = 1.0, for all three integrators. We see that by having a similar friction rate, the three simulations yielded thermal fluctuations of very similar size. 21   Fig. 3-4: Time evolution of the average kinetic temperature โŒฉ๐’Œ๐‘ฉ๐‘ปโŒช for different integrators in DPD model. This fact is further substantiated by comparing the MSD and the logarithm of the VACF, which are plotted in Fig. 3-5 and Fig. 3-6, with the shown trendlines fitted to the DPD-VV data. We observe that the data are very similar to each other and display the anticipated behavior, namely the MSD are linear in time for large ๐‘ก, while the VACF decay exponentially at the expected rate of 0.1 [35]. We also want to note that this rate is significantly different than the LD result, even though the strength of the friction force in the DPD model is about the same, and less on average. This occurs due to the fact that the force is now pairwise, meaning there are now more collisions per timestep, which alters the mean free path and the effective friction rate.  Fig. 3-5: Evolution of the mean square displacement for different integrators in DPD model. Using linear regression, we find the fit for the data from ๐‘ก = 20.0 of the DPD-VV integrators. 22   Fig. 3-6: The decay of velocity autocorrelation function over time for different integrators in DPD model. The trendline is fitted to the DPD-VV result, with the y-intercept fixed at 0.0. Finally, we compared the radial distribution of DPD-VV, Lowe and Imp-LF, of which the result for the first two algorithms has been obtained by Nikunenโ€™s group [20]. We can confirm by eye the expected unity of ๐บ(๐‘…) at large distance.  Furthermore, we notice the anomalies at small ๐‘…๐ผ๐ฝ, namely ๐บ(๐‘…) for DPD-VV and Lowe both become larger than one near ๐‘… = 0, and Loweโ€™s distribution has a small dip around ๐‘… = 0.2, while DPD-VVโ€™s does not. These deviations from the ideal gas limit have both been observed by Nikunen et al., but the group did not offer any explanation other than they are simulation artifacts. We can only partially explain the first behavior by the fact that the particles in our models do not collide and can pass through each other, which disturbs the distribution at small distance. Nevertheless, we verify our implementation produced results similar to Nikunen et al., so we conclude that our propagators are working correctly.   23   Fig. 3-7: The radial distribution function ๐‘ฎ(๐‘น).   24  3.2. Memory With the integrators tested and proven to work correctly, we can now work on the implementation of the memory kernel, the core of our work. In this section, we follow the simulation model used by Liโ€™s group for the standard and auxiliary memory [13, 21].  3.2.1. Standard memory As we have discussed earlier, the objective of this section is not the retrieval of the memory kernel from an atomistic simulation, but the verification of our implementation for the memory kernel. For now, we assume that the relevant values have been found and focus on the computation of the EOM in Eq. (25).To facilitate the computation of the pairwise force between two particles, we want to express it in terms of the radial and orthogonal directions. For a system without any preferred direction, we expect that the fluctuating force ๐›ฟ๐น๐ผ๐ฝ๐’ฌโƒ—โƒ—โƒ—โƒ—  โƒ— has no bias in the orthogonal plane. Therefore, we can decompose ๐›ฟ๐น๐ผ๐ฝ๐’ฌโƒ—โƒ—โƒ—โƒ—  โƒ— into in which ๐›ฟ๐น๐ผ๐ฝ๐’ฌโˆฅโƒ—โƒ— โƒ—โƒ— โƒ—โƒ—   and ๐›ฟ๐น๐ผ๐ฝ๐’ฌโŠฅโƒ—โƒ— โƒ—โƒ— โƒ—โƒ—  โƒ— are the radial and orthogonal components of ๐›ฟ๐น๐ผ๐ฝ๐‘„โƒ—โƒ—โƒ—โƒ—  โƒ—. The memory kernel can also be approximated as where ๐พ๐ผ๐ฝโˆฅ  and ๐พ๐ผ๐ฝโŠฅ are presumed to be known. Discretizing and decomposing Eq. (25), we finally find with ๐‘š being the cutoff length for the memory kernel, ๐‘‰๐ผ๐ฝโˆฅโƒ—โƒ—โƒ—โƒ—  = (๐‘‰๐ผ๐ฝโƒ—โƒ—โƒ—โƒ—  โƒ— โ‹… ๐‘’๐ผ๐ฝโƒ—โƒ— โƒ—โƒ—  )๐‘’๐ผ๐ฝโƒ—โƒ— โƒ—โƒ—  ,  ๐‘‰๐ผ๐ฝโŠฅโƒ—โƒ—โƒ—โƒ—  โƒ— = ๐‘‰๐ผ๐ฝโƒ—โƒ—โƒ—โƒ—  โƒ— โˆ’ ๐‘‰๐ผ๐ฝโˆฅโƒ—โƒ—โƒ—โƒ—   the parallel and orthogonal velocity between ๐ผ and ๐ฝ, ๐œ‰ and ๐ธ Gaussian random number and random antisymmetric matrix, both of which are independent for different particles pairs and at different times [13]. The coefficients ๐›ผ๐ผ๐ฝโˆฅ  and ๐›ผ๐ผ๐ฝโŠฅ  are used to produce random forces that satisfy the generalized FDT:  ๐›ฟ๐น๐ผ๐ฝ๐’ฌโƒ—โƒ—โƒ—โƒ—  โƒ— = (๐‘’๐ผ๐ฝโƒ—โƒ— โƒ—โƒ—  โŠ— ๐‘’๐ผ๐ฝโƒ—โƒ— โƒ—โƒ—  ) โ‹… ๐›ฟ๐น๐ผ๐ฝ๐’ฌโƒ—โƒ—โƒ—โƒ—  โƒ— + (1 โˆ’ ๐‘’๐ผ๐ฝโƒ—โƒ— โƒ—โƒ—  โŠ— ๐‘’๐ผ๐ฝโƒ—โƒ— โƒ—โƒ—  ) โ‹… ๐›ฟ๐น๐ผ๐ฝ๐‘„โƒ—โƒ—โƒ—โƒ—  โƒ— = ๐›ฟ๐น๐ผ๐ฝ๐’ฌโˆฅโƒ—โƒ— โƒ—โƒ— โƒ—โƒ—  + ๐›ฟ๐น๐ผ๐ฝ๐’ฌโŠฅโƒ—โƒ— โƒ—โƒ— โƒ—โƒ—  โƒ— , (54)  ๐‘ฒ๐ผ๐ฝ = (๐‘’๐ผ๐ฝโƒ—โƒ— โƒ—โƒ—  โŠ— ๐‘’๐ผ๐ฝโƒ—โƒ— โƒ—โƒ—  )๐พ๐ผ๐ฝโˆฅ + (1 โˆ’ ๐‘’๐ผ๐ฝโƒ—โƒ— โƒ—โƒ—  โŠ— ๐‘’๐ผ๐ฝโƒ—โƒ— โƒ—โƒ—  )๐พ๐ผ๐ฝโŠฅ , (55)  ฮ”ฮ”๐‘ก๐‘ƒ๐ผโƒ—โƒ—  โƒ—(๐‘ก) = โˆ‘[๐น๐ผ๐ฝ๐ถ(๐‘…๐ผ๐ฝ)๐‘’๐ผ๐ฝโƒ—โƒ— โƒ—โƒ—  ๐ฝโ‰ ๐ผโˆ’ ๐›ฅ๐‘ก โˆ‘ (๐พ๐ผ๐ฝโˆฅ (๐‘…๐ผ๐ฝ(๐‘ก โˆ’ ๐‘›๐›ฅ๐‘ก), ๐‘›๐›ฅ๐‘ก)๐‘‰๐ผ๐ฝโˆฅโƒ—โƒ—โƒ—โƒ—  (๐‘ก โˆ’ ๐‘›๐›ฅ๐‘ก) + ๐พ๐ผ๐ฝโŠฅ(๐‘…๐ผ๐ฝ(๐‘ก โˆ’ ๐‘›๐›ฅ๐‘ก), ๐‘›๐›ฅ๐‘ก)๐‘‰๐ผ๐ฝโŠฅโƒ—โƒ—โƒ—โƒ—  โƒ—(๐‘ก โˆ’ ๐‘›๐›ฅ๐‘ก))๐‘š๐‘›=0+ โˆš๐›ฅ๐‘ก โˆ‘(๐›ผ๐ผ๐ฝโˆฅ (๐‘…๐ผ๐ฝ(๐‘ก โˆ’ ๐‘›๐›ฅ๐‘ก), ๐‘›๐›ฅ๐‘ก)๐œ‰๐ผ๐ฝ(๐‘›๐›ฅ๐‘ก)๐‘’๐ผ๐ฝโƒ—โƒ— โƒ—โƒ—  + ๐›ผ๐ผ๐ฝโŠฅ(๐‘…๐ผ๐ฝ(๐‘ก โˆ’ ๐‘›๐›ฅ๐‘ก), ๐‘›๐›ฅ๐‘ก)๐ธ(๐‘›๐›ฅ๐‘ก) โ‹… ๐‘’๐ผ๐ฝโƒ—โƒ— โƒ—โƒ—  )๐‘š๐‘›=0] , (56)  ๐›ฝ โˆ‘ ๐›ผ๐ผ๐ฝ๐‘ฅ (๐‘…๐ผ๐ฝ, ๐‘›๐›ฅ๐‘ก)๐›ผ๐ผ๐ฝ๐‘ฅ (๐‘…๐ผ๐ฝ, (๐‘› + ๐‘ )๐›ฅ๐‘ก)๐‘šโˆ’๐‘ ๐‘›=0= ๐พ๐ผ๐ฝ๐‘ฅ(๐‘…๐ผ๐ฝ, ๐‘ ๐›ฅ๐‘ก) , (57) 25  for ๐‘ฅ โˆˆ {โˆฅ, โŠฅ}. We can now run the simulation by applying on Eq. (56) any of the integration schemes discussed previously. However, there is one last crucial problem we must address. Up until now, we have ignored the fact that ๐‘’๐ผ๐ฝโƒ—โƒ— โƒ—โƒ—   is a function of time and may vary greatly during our period of interest, as the particlesโ€™ trajectory evolves. This conveniently allows us to write Eq. (55) and (56) in terms of pure radial and orthogonal components, avoiding the complication due to the cross-terms. However, we should expect that these terms are not negligible in general. Unfortunately, all of the works we have consulted do not share this concern; in fact, they seem to imply that the radial unit vector appearing in Eq. (54), (55) and (56) is the current unit vector, evaluated at time ๐‘ก [13, 21, 15]. This simplification is somewhat puzzling, since conceptually, the colored noise generated in the last two terms of Eq. (56) is the residual effects of the fluctuating microscopic forces in the past, which may not be parallel or orthogonal to the unit vector at the present. In Section 3.2.2, we shall also investigate whether this distinction creates any significant discrepancy.  3.2.2. Test result for standard memory 3.2.2.1. Model system As mentioned before, the determination of the memory kernel is rather elaborate, as it requires a reference MD simulation, but the actual values are of little importance right now. Consequentially, for the propagatorsโ€™ tests, we employed memory kernels with simple functional forms. Specifically, for the standard memory, we assumed the temporal and spatial dependence of the kernel are separable, and the random coefficients were in which ๐‘ฅ โˆˆ {โˆฅ, โŠฅ}, 0 < ๐‘ โ‰ค 1 and ๐œŽ(๐‘…๐ผ๐ฝ) has been defined previously in Eq. (46). Equation (57) then required This choice of the memory kernel is beneficial because it vanishes at large time ๐‘ก with the appropriate value for ๐‘ and the cutoff length ๐‘š, while still being easily calculable. Clearly, the larger ๐‘š is, the smaller ๐›ผ๐ผ๐ฝ๐‘ฅ (๐‘šฮ”๐‘ก) and ๐พ๐ผ๐ฝ๐‘ฅ(๐‘šฮ”๐‘ก) become. Fig. 3-8 shows ๐›ผ๐ผ๐ฝ๐‘ฅ /๐œŽ(๐‘…๐ผ๐ฝ) and ๐พ๐ผ๐ฝ๐‘ฅ/๐›ฝ๐œŽ(๐‘…๐ผ๐ฝ)  for two combinations of ๐‘ and ๐‘š, i.e. ๐‘ = 0.7, ๐‘š = 20 and ๐‘ =0.836, ๐‘š = 40. We can see that, for each ๐‘, the chosen ๐‘š is about the lower limit where the vanishing condition still holds, meaning ๐พ๐ผ๐ฝ๐‘ฅ(๐‘›ฮ”๐‘ก) โ‰ช ๐พ๐ผ๐ฝ๐‘ฅ(0) and ๐›ผ๐ผ๐ฝ๐‘ฅ (๐‘›ฮ”๐‘ก) โ‰ช ๐›ผ๐ผ๐ฝ๐‘ฅ (0).  ๐›ผ๐ผ๐ฝ๐‘ฅ (๐‘…๐ผ๐ฝ, ๐‘›๐›ฅ๐‘ก) = ๐‘๐‘›๐œŽ(๐‘…๐ผ๐ฝ) , (58)  ๐พ๐ผ๐ฝ๐‘ฅ(๐‘…๐ผ๐ฝ, ๐‘›๐›ฅ๐‘ก) =๐›ฝ๐‘๐‘›๐‘2 โˆ’ 1(๐‘2(๐‘š+1โˆ’๐‘›) โˆ’ 1)๐œŽ2(๐‘…๐ผ๐ฝ) . (59) 26    Fig. 3-8: Coefficients for dissipative force (left) and random force (right), with some choices of ๐’„ and ๐’. Lastly, for the simulation conditions, all tests used DPD reduced units and were performed on a system of ๐‘ identical particles with ๐‘€ = 1.0 and ๐น๐ผ๐ฝ๐ถ = 0, in periodic cubic boxes of size 10 ร—  10 ร—  10.  3.2.2.2. Result  Fig. 3-9: Average kinetic temperature โŒฉ๐’Œ๐‘ฉ๐‘ปโŒช as a function of cutoff length ๐’Ž for the standard memory scheme. Here, we used ๐‘ = 1000, ๐‘‡ = 3.0, ๐‘ = 0.7 and employed DPD-VV. Unfortunately, even the most basic benchmark, the average temperature, could not be achieved. Our first test was run with ๐‘ = 0.7 and ๐‘š = 20, whose data point is second from the left on Fig. 3-9. We discovered that the kinetic temperature converged nicely, but to a value slightly less than the expected temperature. Naturally, we thought this happened because the memory is insufficient, so we extended the memory cutoff. This means that the memory kernel varies with each case, but the random force is roughly identical, since ๐›ผ๐ผ๐ฝ๐‘ฅ (๐‘› > 20) are relatively very small. Bizarrely, 27  we found that the longer ๐‘š is, the further the average kinetic temperature โŒฉ๐‘‡โŒช drifts from the expected temperature, as we can see from Fig. 3-9. Apparently, our thermostat produced a converging temperature, but, for some reason, violated the FDT. Next, we realized the FDT was also violated by simply changing the absolute value of the reference temperature ๐‘‡. Fig. 3-10 plots โŒฉ๐‘‡โŒช/๐‘‡ for different ๐‘‡, with ๐‘ = 0.836, ๐‘š = 40 and two different system sizes (green dotted line for ๐‘ = 125 and purple dash line for ๐‘ = 256). As in Fig. 3-9, we find โŒฉ๐‘‡โŒช does not converge to the expected temperature ๐‘‡, indicating that the thermostat somehow satisfied the FDT in a wrong way. In the same figure, we also show the corresponding results from simulations with ๐‘ = 0.7, which reinforce our observation, that there is a very strange, but in some sense well-behaved, relationship between the average temperature โŒฉ๐‘˜๐ต๐‘‡โŒช and the reference ๐‘˜๐ต๐‘‡. Once more, the relation hints that we were having an elementary trouble, as it both overestimated and underestimated ๐‘˜๐ต๐‘‡, eliminating the possibility of a systematic error.  Fig. 3-10: Ratio between the kinetic temperature โŒฉ๐’Œ๐‘ฉ๐‘ปโŒช and the reference temperature ๐’Œ๐‘ฉ๐‘ป as a function of ๐’Œ๐‘ฉ๐‘ป in the standard memory scheme. For both ๐‘, the system size was chosen so that the density was constant ๐œŒ = 0.4. The simulations were run with ฮ”๐‘ก = 0.0005 and ๐‘š = 40. Finally, we remark that the proportion is independent of system size. This is confirmed by Fig. 3-11, which shows the ratio โŒฉ๐‘‡โŒช/๐‘‡ for another series of tests with changing ๐‘, but constant ๐œŒ = 0.4, ๐‘ = 0.7 and ๐‘š = 20. The result supports that the identified behavior is not a finite boundary effect: the detected trend actually becomes smoother with larger system size.  28   Fig. 3-11: The relationship between โŒฉ๐‘ปโŒช/๐‘ป and ๐’Œ๐‘ฉ๐‘ป for different system sizes. Here, the simulations were all executed with ฮ”๐‘ก = 0.01, ๐‘ = 0.7 and ๐‘š = 20. We also notice two other trends. First, we want to inspect the effect of the timestep ฮ”๐‘ก on the earlier observations. In fact, the previous two figures actually used two different values for ฮ”๐‘ก: the simulations in Fig. 3-11 were executed with ฮ”๐‘ก = 0.01, while Fig. 3-10 had ฮ”๐‘ก = 0.0005. We isolate the data for ๐‘ = 0.7 in Fig. 3-12 and conclude that the timestep ฮ”๐‘ก affects the results quantitatively, but does not alter the qualitative relationship between โŒฉ๐‘‡โŒช/๐‘‡ and ๐‘˜๐ต๐‘‡. Second, Fig. 3-13 examines the difference of reading ๐‘’๐ผ๐ฝโƒ—โƒ— โƒ—โƒ—   as the past vector or the present vector. Namely, the EOM for the present interpretation is with ๐‘‰๐ผ๐ฝโˆฅโƒ—โƒ—โƒ—โƒ—  (๐‘ก โˆ’ ๐‘›๐›ฅ๐‘ก) = (๐‘‰๐ผ๐ฝโƒ—โƒ—โƒ—โƒ—  โƒ—(๐‘ก โˆ’ ๐‘›๐›ฅ๐‘ก) โ‹… ๐‘’๐ผ๐ฝโƒ—โƒ— โƒ—โƒ—  (๐‘ก)) ๐‘’๐ผ๐ฝโƒ—โƒ— โƒ—โƒ—  (๐‘ก), while the EOM for the past case is obtained by replacing all instances of ๐‘’๐ผ๐ฝโƒ—โƒ— โƒ—โƒ—  (๐‘ก) with ๐‘’๐ผ๐ฝโƒ—โƒ— โƒ—โƒ—  (๐‘ก โˆ’ ๐‘›ฮ”๐‘ก) in Eq. (60). Evidently, for the simulation conditions of ๐‘ = 0.7, ๐‘š = 20, ๐‘ =1000 and ๐‘˜๐ต๐‘‡ = 1.0, there is no appreciable discrepancy between the two choices: both average kinetic temperatures converge to the same, incorrect value of 0.91. We also see that the thermal fluctuations in the present interpretation are larger than, but do not significantly differ from the historical ones.  ฮ”ฮ”๐‘ก๐‘ƒ๐ผโƒ—โƒ—  โƒ—(๐‘ก) = โˆ‘[๐น๐ผ๐ฝ๐ถ(๐‘…๐ผ๐ฝ)๐‘’๐ผ๐ฝโƒ—โƒ—โƒ—โƒ—  โƒ—(๐‘ก)๐ฝโ‰ ๐ผโˆ’ ๐›ฅ๐‘ก โˆ‘ (๐พ๐ผ๐ฝโˆฅ (๐‘…๐ผ๐ฝ(๐‘ก โˆ’ ๐‘›๐›ฅ๐‘ก), ๐‘›๐›ฅ๐‘ก)๐‘‰๐ผ๐ฝโˆฅโƒ—โƒ— โƒ—โƒ— (๐‘ก โˆ’ ๐‘›๐›ฅ๐‘ก) + ๐พ๐ผ๐ฝโŠฅ(๐‘…๐ผ๐ฝ(๐‘ก โˆ’ ๐‘›๐›ฅ๐‘ก), ๐‘›๐›ฅ๐‘ก)๐‘‰๐ผ๐ฝโŠฅโƒ—โƒ—โƒ—โƒ—  โƒ—(๐‘ก โˆ’ ๐‘›๐›ฅ๐‘ก))๐‘š๐‘›=0+ โˆš๐›ฅ๐‘ก โˆ‘ (๐›ผ๐ผ๐ฝโˆฅ (๐‘…๐ผ๐ฝ(๐‘ก โˆ’ ๐‘›๐›ฅ๐‘ก), ๐‘›๐›ฅ๐‘ก)๐œ‰๐ผ๐ฝ(๐‘›๐›ฅ๐‘ก)๐‘’๐ผ๐ฝโƒ—โƒ—โƒ—โƒ—  โƒ—(๐‘ก) + ๐›ผ๐ผ๐ฝโŠฅ(๐‘…๐ผ๐ฝ(๐‘ก โˆ’ ๐‘›๐›ฅ๐‘ก), ๐‘›๐›ฅ๐‘ก)๐ธ(๐‘›๐›ฅ๐‘ก) โ‹… ๐‘’๐ผ๐ฝโƒ—โƒ—โƒ—โƒ—  โƒ—(๐‘ก))๐‘š๐‘›=0] , (60) 29   Fig. 3-12: Effect of the timestep value on โŸจ๐’Œ๐‘ฉ๐‘ปโŸฉ. The simulations all used ๐‘ = 0.7.  Fig. 3-13: Time evolution of the average kinetic temperature โŒฉ๐’Œ๐‘ฉ๐‘ปโŒช for two different interpretations of ๐’†๐‘ฐ๐‘ฑโƒ—โƒ—โƒ—โƒ—  โƒ— in Eq. (56).  In short, our attempt at implementing the standard memory is a confounding blunder. We spent a lot of resources to investigate the problem, testing many different suspicions. For example, we have discussed earlier our test on the different interpretation of Eq. (56). We also inspected whether the various integration methods mentioned in Section 3.1 produced any variation. Furthermore, as the memory was uninitialized (i.e. equal to zero) at the beginning of the simulation, we suspected this initial condition might somehow compound and affect the steady-state behavior, so we attempted a run using the DPD EOM initially, recording the memory along the way, then switched to the MZ EOM. Nevertheless, no matter what test we tried, the resulting characteristics stayed the same. Finally, we 30  tested if we could reduce the MZ EOM to a standard DPD EOM, namely taking ๐‘š = 0. However, this was also not possible as the functional form of our memory kernel cannot satisfy the FDT in this case1.   The obvious pathological bias of โŒฉ๐‘˜๐ต๐‘‡โŒช on the reference temperature ๐‘˜๐ต๐‘‡ implies the thermostat is a complete failure, which means the simulation code may have some algorithmic problems. Curiously, the available literature seems to treat the actual implementation of the memory kernel as a trivial task, since we could not find any detailed discussion. On the other hand, the stability of the thermostat, the independence of the observed behavior from the periodic boundary conditions (and arguably the timestep) hint at more fundamental issues, for example, wrong interpretation of the mathematical foundation. We strongly suspect that our interpretation of Eq. (56) or (57), which we copied directly from the paper by Li et al. [13], is incorrect.  3.2.3. Auxiliary memory This lack of success motivated us to try the auxiliary memory method, with which we obtained much better results. Similar to the previous section, we want to decompose Eq. (30) into pairwise terms, with components in the radial and orthogonal direction [21]: with the parallel and perpendicular components of ๐‘จ and ๐‘ฉ satisfying the second FDT independently. Again, we assume that their functional form has been found. We shall demonstrate in Section 4.2 how their value can be obtained from the memory kernel. In the meantime, we want to discuss the various issues in applying Eq. (61) in the simulation.  First, Eq. (61) is an EOM for the pairwise momenta, not for the momentum of each particle. Intuitively, the decomposition is reasonable, since each ๐‘ƒ๐ผ๐ฝโƒ—โƒ—โƒ—โƒ—  โƒ— can be considered as a degree of freedom. However, it is difficult to design an algorithm that can determine the momenta of ๐‘ particles from the ๐‘(๐‘ โˆ’ 1) constraints found by solving Eq. (61). Thus, we want to use a form where each individual momentum appears explicitly: where ๐‘’๐ผ๐ฝโŠฅโƒ—โƒ— โƒ—โƒ—  =๐‘ƒ๐ผ๐ฝโŠฅโƒ—โƒ— โƒ—โƒ— โƒ—โƒ— ๐‘ƒ๐ผ๐ฝโŠฅ is the unit vector parallel to the orthogonal momentum.                                                           1 It is a bit nuanced than simply taking ๐‘š = 0 in Eq. (58) and (59) then applying the standard FDT relation in Eq. (36). Instead, we must use the generalized formula (Eq. (19)) under the limit ๐‘š โ†’ 0.  (๐‘ƒ๐ผ๐ฝโˆฅฬ‡๐‘ ๐ผ๐ฝโˆฅโƒ—โƒ— โƒ—โƒ—  ฬ‡) = (๐น๐ผ๐ฝ๐ถโƒ—โƒ—โƒ—โƒ—  โƒ—0โƒ— ) โˆ’ (0 (๐‘จ๐ผ๐ฝโˆฅ )๐‘๐‘ (๐‘จ๐ผ๐ฝโˆฅ )๐‘ ๐‘(๐‘จ๐ผ๐ฝโˆฅ )๐‘ ๐‘ )(๐‘ƒ๐ผ๐ฝโˆฅ๐‘ ๐ผ๐ฝโˆฅโƒ—โƒ— โƒ—โƒ—  ) + (0 00 (๐‘ฉ๐ผ๐ฝโˆฅ )๐‘ ๐‘ ) (0โƒ— ๐œ‰๐ผ๐ฝโˆฅโƒ—โƒ— โƒ—โƒ—  ) , (๐‘ƒ๐ผ๐ฝโŠฅฬ‡๐‘ ๐ผ๐ฝโŠฅโƒ—โƒ— โƒ—โƒ—  ฬ‡) = โˆ’ (0 (๐‘จ๐ผ๐ฝโŠฅ )๐‘๐‘ (๐‘จ๐ผ๐ฝโŠฅ )๐‘ ๐‘(๐‘จ๐ผ๐ฝโŠฅ )๐‘ ๐‘ )(๐‘ƒ๐ผ๐ฝโŠฅ๐‘ ๐ผ๐ฝโŠฅโƒ—โƒ— โƒ—โƒ—  ) + (0 00 (๐‘ฉ๐ผ๐ฝโŠฅ )๐‘ ๐‘ ) (0โƒ— ๐œ‰๐ผ๐ฝโŠฅโƒ—โƒ— โƒ—โƒ—  ) , (61)  ๐‘‘๐‘‘๐‘ก๐‘ƒ๐ผโƒ—โƒ—  โƒ— = โˆ‘(๐‘ƒ๐ผ๐ฝโˆฅฬ‡๐‘’๐ผ๐ฝโƒ—โƒ— โƒ—โƒ—  + ๐‘ƒ๐ผ๐ฝโŠฅฬ‡๐‘’๐ผ๐ฝโŠฅโƒ—โƒ— โƒ—โƒ—  )๐ฝโ‰ ๐ผ , (62) 31   Subsequently, we need to find the finite-time propagator for Eq. (62). We must apply the Trotter splitting here, as the direct, naรฏve discretization of Eq. (61) produces incorrect behavior [24]. Here, we shall present the VV form of the finite integrator [24]: where ๐‘ฅ โˆˆ {โˆฅ, โŠฅ} and Again, we note that the generalized FDT is satisfied by the second formula in Eq. (64).  3.2.4. Test result for auxiliary memory 3.2.4.1. Model system  Except for the functional form of the memory, we employed the same conditions as the tests for the standard memory, namely, same system size and no conservative force. On the other hand, the auxiliary memory matrices were appropriated directly from the work by Li et al. [21]2, so we could easily compare the results. Thus, for each parallel degree of freedom, i.e. pair momentum ๐‘ƒ๐ผ๐ฝโˆฅ , we have eight auxiliary variables, corresponding to four damped oscillators, while the orthogonal ones only have four auxiliary variables. The entries for the matrices ๐‘จ are                                                             2 Their standard memory is much more complicated, so we did not use the available function and opted for a simpler one.  ๐‘ƒ๐ผโƒ—โƒ—  โƒ— (๐‘ก +ฮ”๐‘ก2) = ๐‘ƒ๐ผโƒ—โƒ—  โƒ—(๐‘ก) +ฮ”๐‘ก2โˆ‘๐น๐ผ๐ฝ๐ถโƒ—โƒ—โƒ—โƒ—  โƒ—๐ฝโ‰ ๐ผ ,  ๐‘ƒ๐ผ๐ฝโˆฅ = (๐‘ƒ๐ผโƒ—โƒ—  โƒ— (๐‘ก +ฮ”๐‘ก2) โˆ’ ๐‘ƒ๐ฝโƒ—โƒ—  โƒ— (๐‘ก +ฮ”๐‘ก2)) โ‹… ๐‘’๐ผ๐ฝโƒ—โƒ— โƒ—โƒ—   and ๐‘ƒ๐ผ๐ฝโŠฅ = ๐‘ƒ๐ผ๐ฝโƒ—โƒ—โƒ—โƒ—  โƒ— โˆ’ ๐‘ƒ๐ผ๐ฝโˆฅ๐‘’๐ผ๐ฝโƒ—โƒ— โƒ—โƒ—   , (๐‘ƒ๐ผ๐ฝ๐‘ฅ๐‘ ๐ผ๐ฝ๐‘ฅโƒ—โƒ— โƒ—โƒ—  ) = [๐‘ป๐‘ฅ โˆ’ (1 ๐ŸŽ๐ŸŽ ๐ŸŽ)] (๐‘ƒ๐ผ๐ฝ๐‘ฅ๐‘ ๐ผ๐ฝ๐‘ฅโƒ—โƒ— โƒ—โƒ—  ) + ๐‘บ๐‘ฅ (0โƒ— ๐œ‰๐ผ๐ฝ๐‘ฅโƒ—โƒ— โƒ—โƒ—  ) , ๐‘ƒ๐ผโƒ—โƒ—  โƒ—(๐‘ก + ฮ”๐‘ก) = ๐‘ƒ๐ผโƒ—โƒ—  โƒ— (๐‘ก +ฮ”๐‘ก2) +ฮ”๐‘ก2โˆ‘๐น๐ผ๐ฝ๐ถโƒ—โƒ—โƒ—โƒ—  โƒ—๐ฝโ‰ ๐ผ+ โˆ‘(๐‘ƒ๐ผ๐ฝโˆฅ๐‘’๐ผ๐ฝโƒ—โƒ— โƒ—โƒ—  + ๐‘ƒ๐ผ๐ฝโŠฅ๐‘’๐ผ๐ฝโŠฅโƒ—โƒ— โƒ—โƒ—  )๐ฝโ‰ ๐ผ , ๐‘…๐ผโƒ—โƒ—โƒ—โƒ— (๐‘ก + ฮ”๐‘ก) = ๐‘…๐ผโƒ—โƒ—โƒ—โƒ— (๐‘ก) +ฮ”๐‘ก๐‘€๐ผ๐‘ƒ๐ผโƒ—โƒ—  โƒ— (๐‘ก +ฮ”๐‘ก2) , (63)  ๐‘ป๐‘ฅ = ๐‘’โˆ’ฮ”๐‘ก๐‘จ๐‘ฅ , ๐›ฝ๐‘บ๐‘ฅ(๐‘บ๐‘ฅ)๐‘‡ = ๐Ÿ โˆ’ ๐‘ป๐‘ฅ(๐‘ป๐‘ฅ)๐‘‡ . (64) 32  and (๐‘จ๐ผ๐ฝ๐‘ฅ )๐‘ ๐‘= โˆ’(๐‘จ๐ผ๐ฝ๐‘ฅ )๐‘๐‘ , with being the spatial dependence of the memory kernel, in which ๐‘…๐ถโˆฅ  and ๐‘…๐ถโŠฅ are the cutoff distances. The corresponding matrices ๐‘ฉ can then be computed from the FDT relation of Eq. (34). This particular form of the matrices was chosen because, applying Eq. (33), it yields a memory kernel ๐‘ฒ(๐‘ก) equal to a sum of damped harmonic oscillators, which is very convenient for the fitting process [21]. Lastly, for the simulation conditions, we ignored the conservative force again. The simulation was run in a box of 10 ร— 10 ร— 10, with density ๐œŒ = 0.4 and ๐‘˜๐ต๐‘‡ = 3.0.  3.2.4.2. Result Fig. 3-14 displays the average kinetic temperature of the simulation using the auxiliary memory scheme. We confirm the correct temperature is finally obtained, but we cannot comment on the thermal fluctuations, since Li et al. did not provide their corresponding result.  (๐‘จ๐ผ๐ฝโŠฅ )๐‘๐‘ = [0 0.907๐œ‡ 0 0.421๐œ‡] , (๐‘จ๐ผ๐ฝโŠฅ )๐‘ ๐‘ = [51.69 38.78 0 0โˆ’38.78 0 0 00 0 15.00 12.720 0 โˆ’12.72 0] , (๐‘จ๐ผ๐ฝโˆฅ )๐‘๐‘ = [0 0.280๐œ† 0 0.611๐œ† 0 0.658๐œ† 0 0.339๐œ†] , (๐‘จ๐ผ๐ฝโˆฅ )๐‘ ๐‘ =[       19.66 60.55 0 0 0 0 0 0โˆ’60.55 0 0 0 0 0 0 00 0 26.30 23.87 0 0 0 00 0 โˆ’23.87 0 0 0 0 00 0 0 0 23.55 39.44 0 00 0 0 0 โˆ’39.44 0 0 00 0 0 0 0 0 14.29 7.150 0 0 0 0 0 โˆ’7.15 0 ]        , (65)  ๐œ‡ = 30.57(1 + 2.62๐‘…๐ผ๐ฝ ๐‘…๐ถโŠฅโ„ )0.5(1 โˆ’ ๐‘…๐ผ๐ฝ ๐‘…๐ถโŠฅโ„ )1.31 , ๐œ† = 38.06(1 + 1.78๐‘…๐ผ๐ฝ ๐‘…๐ถโˆฅโ„ )0.5(1 โˆ’ ๐‘…๐ผ๐ฝ ๐‘…๐ถโˆฅโ„ )1.31 (66) 33   Fig. 3-14: Time evolution of the average kinetic temperature โŸจ๐’Œ๐‘ฉ๐‘ปโŸฉ in the auxiliary memory scheme.  Fig. 3-15: Time evolution of the mean square displacement for the auxiliary memory shame. Here, we used ๐‘ = 1000 and ๐œŒ = 0.4. 34   Fig. 3-16: Log-log plot of the above mean square displacement versus time for the auxiliary memory method. On the other hand, the obtained MSD and VACF are shown in Fig. 3-15 and Fig. 3-17, respectively, which exhibit similar features when compared with the results of Li et al. [21]. Specifically, the log-log plot of the MSD has a slope of 2.0 in the ballistic regime, then 1.0 at large ๐‘ก (highlighted in red in Fig. 3-16), Furthermore, the VACF creates a negative โ€œtroughโ€ at small ๐‘ก then fluctuates around zero at larger time. We also note that compared to the figures by Liโ€™s group [21], the analogous patterns of our MSD and VACF occur at slightly earlier time, which can be attributed to a smaller system size, since we did not want to spend too much resources on basic testing.   Fig. 3-17: The decay of velocity autocorrelation function over time.  These developments confirm our implementation of the auxiliary memory and its respective integrator work as intended.  35  3.3. Mean force potential In this section, we present the result for our implementation of Boltzmann inversion, a popular method to derive coarse-grained (CG) potentials due to its straightforward nature and general applicability.  3.3.1. Inverse Boltzmann method As its name implies, the technique inverts the identity where ๐‘”(๐‘…) is the radial distribution function, to calculate the potential of mean force ๐‘ˆ(๐‘…). The practical procedure is very direct: Using the atomistic trajectory, we count the number of neighboring particle that fall within the distance ๐‘… โˆ’ ฮ”/2 and ๐‘… + ฮ”/2 for each particle ๐ผ. Then the discrete ๐‘”(๐‘…) is computed by with ๐ป(๐‘ฅ) being the Heaviside step function, the average taken over the ensemble, and so The result can be further refined by running a CG simulation with the newly found values for ๐‘ˆ(๐‘…), computing the CG radial distribution ๐บ(๐‘…) and updating the potential with where ๐›ผ(๐‘…) is a scaling factor to prevent large fluctuations. This approach can be applied iteratively until ๐บ(๐‘…) matches the distribution from the atomistic simulation [36]. However, in the preliminary experiment, we only employed the single Boltzmann inversion, since we deemed the result was already acceptable.  3.3.2. Test result 3.3.2.1. Model system For the test of the inverse Boltzmann method, we wanted to compare with the result by Padding and Briels [25]. This procedure was chosen because our final objective was to apply the MZ formalism on a system of polyethylene chains and find the functional forms for the CG EOM, whose conservative part had been found in the paper. Again, we need to remark that, obviously, we could simply borrow the force field from there, but we still had to run the atomistic simulation to calculate the memory kernel anyway, so we just wanted to have our own complete CG program.  ๐‘”(๐‘…) = ๐‘’โˆ’๐›ฝ๐‘ˆ(๐‘…) , (67)  ๐‘”(๐‘…) =1๐‘(๐‘ โˆ’ 1)โˆ‘โŒฉ๐ป(ฮ”/2 โˆ’ |๐‘…๐ผ๐ฝ โˆ’ ๐‘…|)โŒช4๐œ‹๐‘…2๐ผ,๐ฝ , (68)  ๐‘ˆ(๐‘…) = โˆ’๐‘˜๐ต๐‘‡ ln๐‘”(๐‘…) . (69)  ๐‘ˆ(๐‘…) = ๐‘ˆ(๐‘…) โˆ’ ๐›ผ(๐‘…)๐‘˜๐ต๐‘‡ ln [๐บ(๐‘…)๐‘”(๐‘…)] , (70) 36  Thus, our initial conditions were similar to those used by Padding and Briels and we used the GROMACS MD units, namely nanometer, picosecond, Kelvin, and atomic mass unit. To obtain the atomistic data, we simulated polyethylene chains of 40 monomers per chain at ๐‘‡ = 450.0๐พ using GROMACS united atom force field. As mentioned in Section 1.3, we simply wish to obtain some results for reference, not perfectly realistic data, so we employed the united atom model, where each monomer is one indivisible pseudo-atom [23]. However, the chains are long and easy to entangle, so the system can be stuck in a configuration. To avoid this problem, we prepared ten identical boxes of size much larger than the targeted value, namely ๐ฟ = 5.383๐‘›๐‘š, slowly squeezed them with an NPT thermostat until they reached the expected size and pressure. Because of the randomness in the heat bath, the boxes arrive at very different final configurations. Then, we ran them in NVT simulation and averaged the result. For the CG process, we combined the monomers with the rate of 20 to 1 (see Fig. 3-18, where the monomers are displayed in gray spheres, while the red circles represent the groupings). and computed the various ๐‘”(๐‘…) corresponding to different components of the conservative force. Here, we must distinguish the difference between the bonded and nonbonded pair. The particles are considered nonbonded if they do not share a bonded structure, even if they are on the same molecule. For example, in standard MD [23], the first and fifth particles in a chain are nonbonded because they have no common bond or angle. In this pilot test, we shall consider the nonbonded interaction only.  Fig. 3-18: CG of a polyethylene chain of 40 united atom monomers with the rate of ๐Ÿ๐ŸŽ to ๐Ÿ.   37  3.3.2.2. Result  Fig. 3-19: Inverted nonbonded potential and fitting function. Once the RDF had been obtained, we calculated the potential using Eq. (69), then fitted the discrete data with an analytical function. Fig. 3-19 displays both the RDF resulting from the atomistic simulation and the fitting function, which we chose to be Gaussian, i.e. in which ๐‘… is measured in ๐‘›๐‘š. The fitting parameters were found to be  ๐›ผ = 7.062๐‘˜๐ฝ๐‘š๐‘œ๐‘™โˆ’1, ๐›ฝ = 4.015๐‘›๐‘šโˆ’2 and ๐›พ = โˆ’0.0926. Finally, we ran an NVE, i.e. canonical ensemble, simulation on the CG system, using fixed bond length and initial conditions identical to those of the reference simulation, then the nonbonded CG distribution ๐บ๐‘›๐‘(๐‘…) was extracted and compared with ๐‘”๐‘›๐‘(๐‘…) in Eq. (68). As we can see, the two lines agree reasonably. Moreover, the discrepancy can be attributed to the failure of our analytical estimation to capture the finer shape of ๐‘ˆ๐‘›๐‘(๐‘…) around 0.4 < ๐‘… < 0.7. In fact, ๐บ๐‘›๐‘(๐‘…) follows closely the approximated function, constructed by applying the definition of the RDF to Eq. (71) and drawn with solid cyan line in the figure. Here, it should be emphasized that our test only employed the single Boltzmann inversion, so the achieved result is reasonable. Lastly, we note that the atomistic distribution is very similar to that of Padding and Briels [25], so our machinery for the Boltzmann inversion seems to be working as intended.     ๐‘ˆ๐‘›๐‘(๐‘…) = ๐›ผ๐‘’โˆ’๐›ฝ๐‘…2+๐›พ ๐‘˜๐ฝ๐‘š๐‘œ๐‘™โˆ’1 , (71) 38   Fig. 3-20: Comparison between targeted distribution ๐’ˆ๐’๐’ƒ(๐‘น) and the CG distribution ๐‘ฎ๐’๐’ƒ(๐‘น) obtained from simulation using nonbonded Boltzmann-inverted potential with the CG ratio of ๐Ÿ๐ŸŽ to ๐Ÿ and fixed bond length. The solid line is the approximated function of the atomistic distribution, constructed from the fitting potential in Eq. (70).   39  Chapter 4 Application    4.1. Coarse-graining of bonded interactions To stay consistent with the terminology in the previous section, we shall use the adjective โ€œbondedโ€ to describe anything related to the systemโ€™s topological structures, which can be bonds, valence angles or dihedral angles in standard MD [23]. So, for example, the bonded interactions include the interaction between two particles sharing a bond, and three particles sharing a bond angle, and four particles sharing a torsion angle.  4.1.1. Problem with coarse-graining of bonded interactions  Fig. 4-1: Inverted bond potential and fitting function. In the previous chapter, we ran a trial NVE simulation on the polyethylene system, in which the dimers had fixed bond length, verified the consistency of the nonbonded distributions, and claimed that our implementation of Boltzmann inversion was working properly. However, the test failed spectacularly with regard to the bonded interactions, which we had intentionally avoided earlier. Following the same procedure described in Section 3.3.2, we inverted the mean force potential for the bond between two immediate particles and found the fitting polynomial  ๐‘ˆ๐‘(๐‘…) = 11.87๐‘…6 โˆ’ 88.04๐‘…5 + 270.8๐‘…4 โˆ’ 438.1๐‘…3 + 391.2๐‘…2 โˆ’ 184.0๐‘… + 36.87 ๐‘˜๐ฝ๐‘š๐‘œ๐‘™โˆ’1 . (72) 40  Again, the distance ๐‘… is measured in ๐‘›๐‘š. Both the discrete potential of mean force and its analytical fit are displayed in Fig. 4-1. We then performed an NVE simulation of the CG system, where the particles could interact through either ๐‘ˆ๐‘›๐‘(๐‘…) (Eq. (71)) or ๐‘ˆ๐‘(๐‘…), and their initial conditions were found by applying Eq. (8). The bond RDF ๐บ๐‘(๐‘…) was then extracted and compared with the targeted one ๐‘”๐‘(๐‘…) in Fig. 4-2. We see that the result is simply incorrect: the CG distribution has a very different shape, extending very far into forbidden regions, where the potential is large. This is clearly disconcerting, as the same figure shows the nonbonded distribution ๐บ๐‘›๐‘(๐‘…) still exhibits the same behavior that we observed in Section 3.3.2.2, namely the small deviation from the reference ๐‘”๐‘›๐‘(๐‘…) and the high similarity to the distribution resulting from the fitting nonbonded potential ๐‘ˆ๐‘›๐‘(๐‘…). Additionally, the corresponding approximated distribution for the CG bond correlates strongly with ๐‘”๐‘(๐‘…), but not at all with ๐บ๐‘(๐‘…). The two contrasting observations hint that this complication does not arise from our Boltzmann inversion process, but from something else more fundamental about the bonded interactions.  Fig. 4-2: Comparison between targeted radial distribution ๐’ˆ(๐‘น) and the CG distribution ๐‘ฎ(๐‘น) for interactions of nonbonded pairs (left) and of pairs sharing a CG bond (right), with the CG ratio of ๐Ÿ๐ŸŽ to ๐Ÿ. The solid lines are the approximated function of the atomistic distribution.  4.1.2. Solution After many more failed tests, including ones with all possible interactions in standard MD (nonbonded, bond, angle, torsion), we finally realized that maybe something other than our implementation of the bonded interactions went wrong. First, we considered a system consisting of 500 dimers, which can interact through the CG potentials in Eq. (71) and (72), in a simulation box of size 10 ร— 10 ร— 10๐‘›๐‘š. This gives an average spacing between the monomers of (๐‘‰/๐‘)โˆ’1/3 = 1.5๐‘›๐‘š, about the point of minimum ๐‘ˆ๐‘(๐‘…). We then ran a simulation using the Andersen thermostat with ๐‘‡ = 5.0๐พ and examined the bond radial distribution in Fig. 4-3. Obviously, at this low density, the nonbonded distribution is not correct, so we would not analyze it. Evidently, the resulting data for ๐บ๐‘(๐‘…) agree quite well with 41  the targeted values, except for the region of ๐‘… < 0.2๐‘›๐‘š, where our sixth-order polynomial deviates significantly from the inverted potential. This shows that, at least independently, our code for the bond force works as intended.  Fig. 4-3: Bond radial distribution from the thermalized simulation using Andersen thermostat with ๐‘ป = ๐Ÿ“. ๐ŸŽ๐‘ฒ. The data for targeted distribution ๐‘”๐‘(๐‘…) are represented by the dashed blue line, while the solid line is the approximated function of the atomistic distribution. However, when we ran similar simulations with higher temperatures, i.e. ๐‘‡ = 50.0๐พ and ๐‘‡ = 200.0๐พ, the generated data are incorrect: from Fig. 4-4, we see that they are still different from that of the simulation in Section 4.1.1, but clearly wrong nonetheless. Moreover, the deviations become progressively apparent as the temperature increases. This explains why ๐บ๐‘(๐‘…) diverges remarkably in Fig. 4-2: at ๐‘‡ = 450.0๐พ, the bond structure is completely obscured by the thermal noises. In other words, the bond potential of Eq. (72) is too soft, so we can only obtain the correct pair correlation by exploring the phase-space at low ๐‘‡. In fact, the paper by Padding and Briels [25] specifically addressed this problem and proposed a solution in terms of an extra potential preventing the crossing of the bonds, essentially creating a multi-body potential. In short, with a CG rate of 20 to 1, it is not possible to completely capture the static structure of the polyethylene system by using the standard MD force fields. 42   Fig. 4-4: Bond radial distribution from the thermalized simulation using Andersen thermostat with ๐‘ป = ๐Ÿ“๐ŸŽ. ๐ŸŽ๐‘ฒ (yellow triangle) and ๐‘ป = ๐Ÿ๐ŸŽ๐ŸŽ. ๐ŸŽ๐‘ฒ (green circle). The targeted distribution ๐‘”๐‘(๐‘…) is represented in dashed blue line. Subsequently, we reduced the CG ratio to 8 to 1, meaning the CG chains now have 5 monomers. So far, as a dimer only contains two particles, we had only considered the bond force, but the new chain is too long to ignore the angle force. The formalism to evaluate the functional form of the new interaction is similar to the one discussed in Section 3.3.1, except the angular distribution must be normalized with a factor of 1/ sin ๐œƒ, i.e. with ๐‘๐‘Ž being the total number of angles in the system. Fig. 4-5 shows the inverted data for both bonded potentials and their fitting function, which we chose to be polynomial: where ๐‘ˆ๐‘Ž(๐œƒ) is the angle potential, with ๐œƒ in radian.  ๐‘”๐‘Ž(๐œƒ) =1๐‘๐‘Žโˆ‘โŒฉ๐ป(ฮ”/2 โˆ’ |๐œƒ๐‘– โˆ’ ๐œƒ|)โŒชsin ๐œƒ๐‘–=1 , (73)  ๐‘ˆ๐‘(๐‘…) = 1918๐‘…4 โˆ’ 4971๐‘…3 + 4763๐‘…2 โˆ’ 2022๐‘… + 328.6 ๐‘˜๐ฝ๐‘š๐‘œ๐‘™โˆ’1 , ๐‘ˆ๐‘Ž(๐œƒ) = 8.835๐œƒ6 โˆ’ 103.0๐œƒ5 + 482.8๐œƒ4 โˆ’ 1158๐œƒ3 + 1489๐œƒ2 โˆ’ 974.5๐œƒ + 265.4 ๐‘˜๐ฝ๐‘š๐‘œ๐‘™โˆ’1, (74) 43     Fig. 4-5: Inverted potential and fitting function for the bonded interactions, with the CG ratio of ๐Ÿ– to ๐Ÿ. On the other hand, the nonbonded potential, shown in Fig. 4-6, possesses a behavior characteristic of the Lennard-Jones potential [37], which is where ๐œ– and ๐‘…๐‘š are the depth of the potential well and the distance where ๐‘ˆ๐‘›๐‘(๐‘…) reaches the minimum. In our case, they have the value of  ๐œ– = 0.690๐‘˜๐ฝ๐‘š๐‘œ๐‘™โˆ’1  and ๐‘…๐‘š = 0.570๐‘›๐‘š. In fact, we also expect to find the wavelike structure at the tail of the inverted function: as our setup is far from the dilute limit, the pair correlation function for the nonbonded interaction is influenced by the indirect effects from the neighboring particles [1].  Fig. 4-6: Inverted nonbonded potential and fitting function for the CG ratio of ๐Ÿ– to ๐Ÿ.  ๐‘ˆ๐‘›๐‘(๐‘…) = ๐œ– [(๐‘…๐‘š๐‘…)12โˆ’ 2(๐‘…๐‘š๐‘…)6] , (75) 44  Next, we prepared a CG system following the same procedure described in Section 3.3.2, except each particle is now mapped from eight atoms and can act on the others through three different potentials, then performed an NVE simulation. The reader can find more detailed explanation on the implementation of the new interactions in Section 0. For now, we want to discuss the resulting distribution for the bonded (Fig. 4-7) and nonbonded (Fig. 4-8) interactions. We find the bonded distributions match quite well to the targeted values, save for the regions where the fitting polynomials fail to capture the inverted potentialsโ€™ subtler features, such as the curve around 0.5 < ๐‘… < 0.75 for ๐‘ˆ๐‘(๐‘…) or the potential well near ๐œƒ = 0.7 for ๐‘ˆ๐‘Ž(๐œƒ). On the other hand, ๐บ๐‘›๐‘(๐‘…) only agrees qualitatively, but not quantitatively with ๐‘”๐‘›๐‘(๐‘…), since the repulsive portion of the Lennard-Jones fit is much harder than the raw ๐‘ˆ๐‘›๐‘(๐‘…). Moreover, even if the wanted potential is Lennard-Jones, the minimum is not really at ๐‘… = 0.57๐‘›๐‘š: it has been shifted from its true value due to, again, the influence from the systemโ€™s density. However, for the purpose of this chapter, these finer details can be disregarded, so we conclude that the ratio of 8 to 1 provides an acceptable CG model for the polyethylene system.   Fig. 4-7: Atomistic and CG distributions for the bonded interactions, with the CG ratio of ๐Ÿ– to ๐Ÿ.  45   Fig. 4-8: Atomistic and CG distributions for the CG ratio of ๐Ÿ– to ๐Ÿ. We do not provide the approximated function here since the correct form must be modified using the density, in an expansion akin to the virial equation.        46  4.2. Memory extraction This section describes the memory extraction process, mostly following the work done by Liโ€™s group [13, 21]. At a glance, the definition of the memory kernel ๐‘ฒ๐ผ๐ฝ(๐‘ก) = ๐›ฝ โŒฉ๐›ฟ๐น๐ผ๐ฝ๐“ โƒ—โƒ—โƒ—โƒ—  โƒ—(๐‘ก) โŠ— ๐›ฟ๐น๐ผ๐ฝ๐“ โƒ—โƒ—โƒ—โƒ—  โƒ—(0)โŒช seems to offer a straightforward method to compute the kernel: once the mean force field has been determined, one only has to calculate, directly from the atomistic trajectory, ๐›ฟ๐น๐ผ๐ฝโƒ—โƒ—โƒ—โƒ—  โƒ—(๐‘ก) = ๐น๐ผ๐ฝ๐‘ƒโƒ—โƒ—โƒ—โƒ—  โƒ—(๐‘ก) โˆ’ ๐น๐ผ๐ฝโƒ—โƒ—โƒ—โƒ—  โƒ—(๐‘ก) and find its temporal autocorrelation function. Regrettably, the true fluctuating microscopic force ๐›ฟ๐น๐ผ๐ฝ๐‘„โƒ—โƒ—โƒ—โƒ—  โƒ—(๐‘ก) is not the same as ๐›ฟ๐น๐ผ๐ฝโƒ—โƒ—โƒ—โƒ—  โƒ—(๐‘ก) for long memory. The explanation for the difference is quite nuanced, so we shall elaborate it at the end of Appendix A, but here we need to know that the correlation ๐‘ช๐ผ๐ฝ(๐‘ก) = ๐›ฝโŒฉ๐›ฟ๐น๐ผ๐ฝโƒ—โƒ—โƒ—โƒ—  โƒ—(๐‘ก) โŠ— ๐›ฟ๐น๐ผ๐ฝโƒ—โƒ—โƒ—โƒ—  โƒ—(0)โŒช is just the zeroth-order approximation of the memory kernel. Instead, ๐‘ฒ๐ผ๐ฝ can be calculated exactly through This equation is found by multiplying ๐‘‰๐ผ๐ฝโƒ—โƒ—โƒ—โƒ—  โƒ—(0) on both sides of Eq. (25) and taking the ensemble average. Either way, we must assume ๐‘ฒ๐ผ๐ฝ is a function of ๐‘…๐ผ๐ฝโƒ—โƒ— โƒ—โƒ—  โƒ— only, as we cannot know the atomistic configuration of ๐ผ and ๐ฝ from the CG simulation. Unfortunately, this means that not all autocorrelated values can be used: the calculation can only be done with the ones where the two particles have the same relative position. Thus, in practice, we divide the intermolecular distance ๐‘…๐ผ๐ฝ into many bins of size ฮ”, preferably with the same value used previously for the Boltzmann inversion, then compute ๐‘ฒ๐ผ๐ฝ(๐‘…, ๐‘ก) by averaging over all ๐ผ๐ฝ pairs staying within ๐‘… โˆ’ ฮ”/2 < ๐‘…๐ผ๐ฝ < ๐‘… + ฮ”/2 during the period ๐‘ก, i.e. where we have discretized Eq. (76) and applied the decomposition similar to Eq. (55). The parallel and orthogonal factors of the velocity are ๐‘‰๐ผ๐ฝโˆฅ (๐‘ ) = ๐‘‰๐ผ๐ฝโƒ—โƒ—โƒ—โƒ—  โƒ—(๐‘ ) โ‹… ๐‘’๐ผ๐ฝโƒ—โƒ— โƒ—โƒ—  (๐‘ ) and ๐‘‰๐ผ๐ฝโŠฅ(๐‘ ) = ๐‘‰๐ผ๐ฝโƒ—โƒ—โƒ—โƒ—  โƒ—(๐‘ ) โˆ’ ๐‘‰๐ผ๐ฝโˆฅ (๐‘ )๐‘’๐ผ๐ฝโƒ—โƒ— โƒ—โƒ—  (๐‘ ), for โˆ€๐‘  < ๐‘ก, and the components of ๐›ฟ๐น๐ผ๐ฝโƒ—โƒ—โƒ—โƒ—  โƒ— are defined identically. The two formulae in Eq. (77) allow the memory kernel to be computed recursively, with the base cases being ๐พ๐ผ๐ฝ๐‘ฅ(๐‘…, 0): for ๐‘ฅ โˆˆ {โˆฅ, โŠฅ} and the average taken over all pairs ๐ผ๐ฝ with |๐‘…๐ผ๐ฝ(0) โˆ’ ๐‘…| < ฮ”/2. Next, we must make the distinction between different interaction types: it is reasonable to expect that the memory kernel should behave differently for a  โŒฉ๐›ฟ๐น๐ผ๐ฝโƒ—โƒ—โƒ—โƒ—  โƒ—(๐‘ก) โŠ— ๐‘‰๐ผ๐ฝโƒ—โƒ—โƒ—โƒ—  โƒ—(0)โŒช = โˆซ ๐‘ฒ๐ผ๐ฝ(๐‘ก โˆ’ ๐‘ )โŒฉ๐‘‰๐ผ๐ฝโƒ—โƒ—โƒ—โƒ—  โƒ—(๐‘ ) โŠ— ๐‘‰๐ผ๐ฝโƒ—โƒ—โƒ—โƒ—  โƒ—(0)โŒช๐‘‘๐‘ ๐‘ก0 . (76)  [โŒฉ๐›ฟ๐น๐ผ๐ฝโˆฅ (๐‘ก)๐‘‰๐ผ๐ฝโˆฅ (0)โŒช = โˆ‘ ๐พ๐ผ๐ฝโˆฅ (๐‘…, ๐‘ก โˆ’ ๐‘šฮ”๐‘ก)โŒฉ๐‘‰๐ผ๐ฝโˆฅ (๐‘šฮ”๐‘ก)๐‘‰๐ผ๐ฝโˆฅ (0)โŒชฮ”๐‘ก๐‘ก/ฮ”๐‘ก๐‘š]|๐‘…๐ผ๐ฝ(๐‘šฮ”๐‘ก<๐‘ก)โˆ’๐‘…|<ฮ”/2 , [โŒฉ๐›ฟ๐น๐ผ๐ฝโŠฅ(๐‘ก)๐‘‰๐ผ๐ฝโŠฅ(0)โŒช = โˆ‘ ๐พ๐ผ๐ฝโŠฅ(๐‘…, ๐‘ก โˆ’ ๐‘šฮ”๐‘ก)โŒฉ๐‘‰๐ผ๐ฝโŠฅ(๐‘šฮ”๐‘ก)๐‘‰๐ผ๐ฝโŠฅ(0)โŒชฮ”๐‘ก๐‘ก/ฮ”๐‘ก๐‘š]|๐‘…๐ผ๐ฝ(๐‘šฮ”๐‘ก<๐‘ก)โˆ’๐‘…|<ฮ”/2 , (77)  ๐พ๐ผ๐ฝ๐‘ฅ(๐‘…, 0) =1ฮ”๐‘กโŒฉ๐›ฟ๐น๐ผ๐ฝ๐‘ฅ(0)๐‘‰๐ผ๐ฝ๐‘ฅ(0)โŒชโŒฉ๐‘‰๐ผ๐ฝ๐‘ฅ(0)๐‘‰๐ผ๐ฝ๐‘ฅ(0)โŒช  , (78) 47  different potential. However, at the point of writing of this work, we have not been able to clarify a meaningful method to determine the force between bonded particles, even when they are only connected through some CG bonds. As an example, let us examine the atomistic angle ๐œƒ lying in between two particles ๐ผ and ๐ฝ in Fig. 4-9: it is unclear whether ๐น๐‘Ÿ๐‘Žโƒ—โƒ— โƒ—โƒ—   and ๐น๐‘ ๐‘Žโƒ—โƒ— โƒ—โƒ—   should be included in the internal force of ๐ผ or considered to be the external forces acting on ๐ฝ by ๐ผ. Therefore, in this chapter, we shall only evaluate the nonbonded memory kernel.  Fig. 4-9: Two bonded CG particles. To evaluate the functional form of the memory kernel, we made one last approximation that the temporal and spatial dependence of the kernel are separable, i.e. Indeed, we found that ๐พ๐ผ๐ฝ๐‘ฅ(๐‘…, ๐‘ก) at different times all decay exponentially. Moreover, from Fig. 4-10, which plots the logarithm of ๐พ๐ผ๐ฝ๐‘ฅ(๐‘…, ๐‘ก) for various ๐‘ก, we note that the decay rate does not change much over time, supporting our stated assumption. Disregarding the data for ๐‘… < 0.4๐‘›๐‘š, where the nonbonded force is large, making poor sampling, and averaging over 0 < ๐‘ก < 0.3, we received a slope of 7.38๐‘›๐‘šโˆ’1 for the radial and 7.36๐‘›๐‘šโˆ’1 for the perpendicular component of ๐พ๐ผ๐ฝ.  ๐พ๐ผ๐ฝ๐‘ฅ(๐‘…, ๐‘ก) = ฮ˜(๐‘…)ฮค(๐‘ก) . (79) 48    Fig. 4-10: Decay rate of the radial dependence in the parallel (left) and orthogonal (right) memory kernel.  Fig. 4-11: Temporal dependence at ๐‘น > ๐ŸŽ. ๐Ÿ–๐’๐’Ž of the parallel (left) and orthogonal (right) memory kernel.  On the other hand, the temporal function ฮค(๐‘ก) has more fluctuations and its form depends highly on the value of ๐‘…. In particular, Fig. 4-11 shows that, for ๐‘… > 0.8๐‘›๐‘š, the ratio ๐พ๐ผ๐ฝ๐‘ฅ(๐‘…, ๐‘ก)/๐พ๐ผ๐ฝ๐‘ฅ(๐‘…, 0) may not only stop decaying but can also increase again. However, we must recognize that, at this large intermolecular distance, the absolute value of the memory kernel is relatively small, meaning the fluctuating force ๐›ฟ๐น๐ผ๐ฝ๐’ฌโƒ—โƒ—โƒ—โƒ—  โƒ— is minuscule and its effect on the particlesโ€™ motion negligible. Consequently, we decided to average over 0.4 < ๐‘… < 0.7, where ๐พ๐ผ๐ฝ๐‘ฅ(๐‘…, ๐‘ก)/๐พ๐ผ๐ฝ๐‘ฅ(๐‘…, 0) behaves quite homogenously, and approximated the ratio as a series of damped harmonic oscillators. Using three oscillators for the fitting process, we found the parallel component to be   while the orthogonal part is  ๐พ๐ผ๐ฝโˆฅ (๐‘…, ๐‘ก)๐พ๐ผ๐ฝโˆฅ (๐‘…, 0)= 14.4๐‘’โˆ’6.74๐‘ก cos(0.0866๐‘ก โˆ’ 1.56) + 0.919๐‘’โˆ’42.6๐‘ก cos(21.2๐‘ก โˆ’ 1.11)+ 0.474๐‘’โˆ’11.2๐‘ก cos(18.4๐‘ก โˆ’ 0.547) , (80) 49  The entries for the auxiliary matrices corresponding to each oscillator ๐ด(๐‘ก) = ๐‘ exp(โˆ’๐‘ž๐‘ก) cos(๐‘Ÿ๐‘ก โˆ’ ๐‘ ) could then be computed with  Fig. 4-12: Temporal correlation ๐‘ฒ๐‘ฐ๐‘ฑ๐’™ (๐‘น, ๐’•)/๐‘ฒ๐‘ฐ๐‘ฑ๐’™ (๐‘น, ๐ŸŽ) and its fitted function for the parallel (left) and orthogonal (right) memory kernel. Fig. 4-12 illustrates the radial and perpendicular temporal correlations and their fitting function. We note that the raw data do not extend very far in ๐‘ก, an unfortunate consequence from the compactness and the high diffusion rate of our setup: from Fig. 4-6, we concluded that the spatial resolution could be at most 0.02๐‘›๐‘š, giving reliable sampling until around ๐‘ก = 0.4๐‘๐‘ . Thus, we demanded the oscillators to decay quickly to obtain a reasonable kernel. Also, we observe the derivative of the memory kernel in both directions at ๐‘ก = 0 is zero, which is expected according to Li et al. [21]. To force the same condition on the damped oscillators, we employed the constraint tan ๐‘  = ๐‘ž/๐‘Ÿ. In our test, the coupling matrices were found to be    ๐พ๐ผ๐ฝโˆฅ (๐‘…, ๐‘ก)๐พ๐ผ๐ฝโˆฅ (๐‘…, 0)= 5.94๐‘’โˆ’4.00๐‘ก cos(0.0791๐‘ก โˆ’ 1.55) + 1.30๐‘’โˆ’33.0๐‘ก cos(21.8๐‘ก โˆ’ 0.195)+ 0.182๐‘’โˆ’7.73๐‘ก cos(17.4๐‘ก โˆ’ 0.418) . (81)  (๐‘จ๐ผ๐ฝ๐‘ฅ )๐‘๐‘ = [โˆš๐‘2(cos ๐‘  โˆ’๐‘Ÿ๐‘žsin ๐‘ ) โˆš๐‘2(cos ๐‘  +๐‘Ÿ๐‘žsin ๐‘ )] , (๐‘จ๐ผ๐ฝ๐‘ฅ )๐‘ ๐‘ = [2๐‘ž โˆš๐‘Ÿ2 + ๐‘ž2โˆ’โˆš๐‘Ÿ2 + ๐‘ž2 0] , (๐‘จ๐ผ๐ฝ๐‘ฅ )๐‘ ๐‘= (๐‘จ๐ผ๐ฝ๐‘ฅ )๐‘๐‘  . (82) 50  with     (๐‘จ๐ผ๐ฝโˆฅ )๐‘๐‘ = [ 0 1.838๐œ† 0 4.616๐œ† 0 2.368๐œ†] , (๐‘จ๐ผ๐ฝโˆฅ )๐‘ ๐‘ =[     13.48 6.745 0 0 0 0โˆ’6.745 0 0 0 0 00 0 85.24 47.59 0 00 0 โˆ’47.59 0 0 00 0 0 0 22.43 21.590 0 0 0 โˆ’21.59 0 ]      , (๐‘จ๐ผ๐ฝโŠฅ )๐‘๐‘ = [ 0 1.413๐œ‡ 0 4.060๐œ‡ 0 1.966๐œ‡] , (๐‘จ๐ผ๐ฝโŠฅ )๐‘ ๐‘ =[     7.998 4.000 0 0 0 0โˆ’4.000 0 0 0 0 00 0 65.93 39.52 0 00 0 โˆ’39.52 0 0 00 0 0 0 15.46 19.050 0 0 0 โˆ’19.05 0 ]      , (83)  ๐œ†2 = 8.224 ร— 105๐‘’โˆ’7.38๐‘… , ๐œ‡2 = 2.591 ร— 105๐‘’โˆ’7.36๐‘… . (84) 51  4.3. Simulation model For the simulation of the CG polyethylene system, we shall ignore the effect from the torsion angle; thus, particles that have 3 neighbors in between on a chain are considered nonbonded. Consequentially, our conservative force field possesses three terms [25]: in which ๐น๐ผ๐ฝ๐‘›๐‘ is the nonbonded force, ๐น๐ผ๐ฝ๐‘  the bond force and ๐น๐ผ๐‘Ž the angle force acting on the particle ๐ผ. The first sum runs over all particles not bonded to ๐ผ, the second sum goes over all immediate neighbors of ๐ผ, while the third over angles having ๐ผ as one of their vertices. Assuming the form of potentials is known, the nonbonded and bond forces are where ๐‘ˆ(๐‘…๐ผ๐ฝ) is the corresponding potential. On the other hand, when ๐ผ is not the central particle of the angle (see Fig. 4-13), the angle force is calculated by [27] with ๐‘ˆ๐‘Ž(๐œƒ) being the angle potential, ๐œƒ = ๐ผ๐ฝ?ฬ‚? โ‰ก arccos(๐‘’๐ฝ๐ผโƒ—โƒ— โƒ—โƒ—  โ‹… ๐‘’๐ฝ๐พโƒ—โƒ— โƒ—โƒ—  โƒ—) the angle, ๐‘’๐ผโƒ—โƒ—โƒ—  = ๐‘›๐‘œ๐‘Ÿ๐‘š (๐‘’๐ฝ๐ผโƒ—โƒ— โƒ—โƒ—  ร— (๐‘’๐ฝ๐ผโƒ—โƒ— โƒ—โƒ—  ร— ๐‘’๐ฝ๐พโƒ—โƒ— โƒ—โƒ—  โƒ—)) the normalized vector in the plane of ๐ผ๐ฝ?ฬ‚?, orthogonal to ๐‘…๐ฝ๐ผโƒ—โƒ— โƒ—โƒ—  โƒ—. The angle force acting on the center can be found by applying Newtonโ€™s third law. The form of all three potentials can be found in Section 4.1.2.  Fig. 4-13: Valence angle. ๐‘’๐ฝ๐ผโƒ—โƒ— โƒ—โƒ—   and ๐‘’๐ฝ๐พโƒ—โƒ—โƒ—โƒ—โƒ—โƒ—  are normalized vectors pointing from the center particle ๐ฝ.  ๐น๐ผ๐ถโƒ—โƒ— โƒ—โƒ—  = โˆ‘ ๐น๐ผ๐ฝ๐‘›๐‘๐ฝโ‰ ๐ผ(๐‘…๐ผ๐ฝ)๐‘’๐ผ๐ฝโƒ—โƒ— โƒ—โƒ—  + โˆ‘ ๐น๐ผ๐ฝ๐‘(๐‘…๐ผ๐ฝ)๐‘’๐ผ๐ฝโƒ—โƒ— โƒ—โƒ—  ๐ฝ+ โˆ‘ ๐น๐ผ๐‘Žโƒ—โƒ— โƒ—โƒ—  (๐œƒ)๐œƒโˆ‹๐ผ ,  (85)  ๐น๐ผ๐ฝโƒ—โƒ—โƒ—โƒ—  โƒ—(๐‘…๐ผ๐ฝ) = โˆ’๐œ•๐‘ˆ(๐‘…๐ผ๐ฝ)๐œ•๐‘…๐ผ๐ฝ๐‘’๐ผ๐ฝโƒ—โƒ— โƒ—โƒ—   , (86)  ๐น๐ผ๐‘Žโƒ—โƒ— โƒ—โƒ—  (๐œƒ) = โˆ’๐œ•๐‘ˆ๐‘Ž(๐œƒ)๐œ•๐œƒ๐‘’๐ผโƒ—โƒ—โƒ—  ๐‘…๐ฝ๐ผ , (87) 52  Combining with the memory kernel ๐พ๐ผ๐ฝ๐‘ฅ(๐‘…, ๐‘ก) in Section 4.2, we obtain the EOM for the CG simulation: with ๐‘ฅ โˆˆ {โˆฅ, โŠฅ} and ๐‘ ๐ผ๐ฝ๐‘ฅโƒ—โƒ— โƒ—โƒ—   being the corresponding component of the auxiliary variables for the ๐ผ๐ฝ pair. The integration scheme for the auxiliary memory is presented in Section 3.2.3. For the simulation conditions, we prepared an atomistic system of 125 polyethylene chains with length of 40 monomers, followed the thermalization procedure in Section 3.3.2.1, then applied a CG ratio of 8 to 1 (see Fig. 4-14, with the empty circles being the CG particles). While the two end particles of the resulting molecule should behave somewhat differently from the rest of the body, we considered all particles to be of the same type, with mass of 120.28๐‘ข. The simulation was performed in a periodic box of size 5.38 ร— 5.38 ร— 5.38๐‘›๐‘š, at ๐‘‡ = 450.0๐พ and ฮ”๐‘ก =0.001๐‘๐‘ . The cutoff distance for both the nonbonded interaction and the memory kernel was chosen at 1.0๐‘›๐‘š.  Fig. 4-14: CG of a polyethylene chain of 40 united atom monomers with the rate of ๐Ÿ– to ๐Ÿ.       ๐‘‘๐‘‘๐‘ก๐‘ƒ๐ผโƒ—โƒ—  โƒ— = โˆ‘ (๐น๐ผ๐ฝ๐‘›๐‘(๐‘…๐ผ๐ฝ)๐‘’๐ผ๐ฝโƒ—โƒ— โƒ—โƒ—  โˆ’ [(๐‘จ๐ผ๐ฝโˆฅ )๐‘๐‘ ๐‘ ๐ผ๐ฝโˆฅโƒ—โƒ— โƒ—โƒ—  ] ๐‘’๐ผ๐ฝโˆฅโƒ—โƒ— โƒ—โƒ—  โˆ’ [(๐‘จ๐ผ๐ฝโŠฅ )๐‘๐‘ ๐‘ ๐ผ๐ฝโŠฅโƒ—โƒ— โƒ—โƒ—  ] ๐‘’๐ผ๐ฝโŠฅโƒ—โƒ— โƒ—โƒ—  )๐ผ๐ฝ ๐‘›๐‘œ๐‘›๐‘๐‘œ๐‘›๐‘‘๐‘’๐‘‘ +โˆ‘๐น๐ผ๐ฝ๐‘(๐‘…๐ผ๐ฝ)๐‘’๐ผ๐ฝโƒ—โƒ— โƒ—โƒ—  ๐ฝ+ โˆ‘๐น๐ผ๐‘Žโƒ—โƒ— โƒ—โƒ—  (๐œƒ)๐œƒโˆ‹๐ผ , ๐‘‘๐‘‘๐‘ก๐‘ ๐ผ๐ฝ๐‘ฅโƒ—โƒ— โƒ—โƒ—  = โˆ’(๐‘จ๐ผ๐ฝ๐‘ฅ )๐‘ ๐‘๐‘ƒ๐ผ๐ฝ๐‘ฅโƒ—โƒ— โƒ—โƒ—  โˆ’ (๐‘จ๐ผ๐ฝ๐‘ฅ )๐‘ ๐‘ ๐‘ ๐ผ๐ฝ๐‘ฅโƒ—โƒ— โƒ—โƒ—  + (๐‘ฉ๐ผ๐ฝ๐‘ฅ )๐‘ ๐‘ ๐œ‰๐ผ๐ฝ๐‘ฅโƒ—โƒ— โƒ—โƒ—   , (88)  53  4.4. Result  We obtained a very reasonable result for the MZ simulation. First, let us inspect the distribution for the nonbonded (Fig. 4-15) and bonded (Fig. 4-16) interactions. In the same figures, we also show the corresponding data extracted from the atomistic (blue dot) and from the NVE (red dash) simulation in Section 4.1.2. Again, we observe that the MZ distributions exhibit the rough features of the atomistic distribution, except the area where the fitting potentials do not match the inverted ones. Moreover, the resulting ๐บ๐‘›๐‘(๐‘…), ๐บ๐‘(๐‘…), ๐บ๐‘Ž(๐œƒ) closely follow those of the NVE simulation, and, by proxy, the approximated function constructed from the fitting potential. This means as far as we are concerned, the MZ simulation replicated the correct static structure of the CG polyethylene system.   Fig. 4-15: Nonbonded distribution for the CG simulation using the auxiliary memory method. The blue dotted line is the atomistic distribution, while the red dashed one represents the NVE data.  Fig. 4-16: Nonbonded distribution for the CG simulation using the auxiliary memory method. 54   Next, we want to examine the dynamical property. Fig. 4-17 displays ๐ถ(๐‘ก)/๐ถ(0) from all three simulations: atomistic, NVE and MZ. The VACFs had to be scaled since under the mapping in Eq. (8), many degrees of freedom were eliminated, reducing the effective temperature in the mapped system. Looking at the figure, we find that the MZ and atomistic normalized VACFs match each other quite well, including an inflection point near ๐‘ก = 0.3๐‘๐‘ . However, the former has a zero derivative at ๐‘ก = 0, a feature not shared with the latter. Apparently, this is a side-effect of the memory method: It persisted over various trial memory kernels ๐พ๐ผ๐ฝ๐‘ฅ(๐‘ก), even when ๐‘‘๐พ๐ผ๐ฝ๐‘ฅ(๐‘ก)/๐‘‘๐‘ก|๐‘ก=0โ‰  0. Nevertheless, this is a successful implementation of the memory method, as we can see that the NVE VACF behaves very differently from the rest, decaying at a much slower rate. This is another expected consequence of the CG projection, which decreased the systemโ€™s thermal fluctuations and the effective friction rate. Moreover, we had not included the kernel for the bonded interactions, so this level of agreement is very reasonable.  Fig. 4-17: Normalized VACF for the CG simulation using the auxiliary memory method. Similar to the graphs above, this also displays the data from the atomistic and the NVE simulation.   55  Chapter 5 Conclusion    Previously, Li et al. had successfully employed the MZ projection to recover the correct dynamics for a system of CG star-polymer melts, where each CG particle represented one polymer [13]. In this work, we investigated the applicability of the MZ formalism on the more complex polyethylene system, which included bonded structures between the particles. While the formalism itself is quite straightforward, there were many subtle problems needed to be resolved before the practical application could be realized. First, the MZ EOM is formally very similar to the DPD EOM. Thus, we wanted to confirm whether a correct DPD integration scheme was necessary. However, at the point of writing, many popular simulation packages did not contain a working DPD integrator, so we decided to build our own version. While the implementation for the integrator was successful, we were not so lucky with our MZ endeavor: our simulations using the standard memory could not generate the correct temperature, no matter what we tried. Nevertheless, our systemโ€™s temperature still converged, albeit to a wrong value, and it clearly depended on the trial memory kernel we used. We believed this violation of the FDT was a consequence of our erroneous interpretation of the MZ EOM, but we were not able to resolve the issue. Fortunately, there is an alternative approach, which substitutes the direct memory with a series of damped oscillators. The behavior of these extra variables imitates the effect of the memory, in a fashion akin to the Nosรฉ-Hoover heat bath [38]. We managed to implement the auxiliary memory correctly. Finally, we discovered that the CG process of the polyethylene system was not as simple as it seemed: it turned out that a naรฏve decomposition of the multi-body CG potential into standard MD components, namely nonbonded, bond and angle potentials, could not work with an aggressive CG strategy, such as a rate of 20 monomers to 1. Thus, in the end, we decided to use a rate of 8 to 1, whose potential replicated the systemโ€™s distributions reasonably. Next, we constructed a CG model for the polyethylene system by applying the MZ formalism on the MD trajectory. The atomistic data were obtained by simulating a system of 125 polyethylene chains of 40 monomers, using a GROMACS united atom force field. Then, we coarsened the chains with a rate of 8 to 1 and retrieved the CG potentials by employing the Boltzmann inversion. Once the potentials were obtained, we evaluated the fluctuating forces, defined as the difference between the sum of the atomistic forces and the mean force on the particle. These forces allowed us to compute the memory kernel and its corresponding coupling oscillators. Since we did not have a correct method to extract the kernel for the bonded interactions, we only calculated the nonbonded one, but our final simulation demonstrated that it was sufficient to correct the general dynamics of the system. Besides, while we could not provide an exact number on the gain in speed, as our implementation is certainly not as optimized as GROMACS, we estimated the speedup to be at least one order of magnitude by employing our own code for both the atomistic and 56  the CG runs. Thus, the MZ formalism achieved both targets outlined in Chapter 1, namely producing the accurate dynamics and improving the simulationโ€™s performance. In the future, we plan to derive a formalism to evaluate the bonded memory kernel. Previously in Chapter 4, we discussed one aspect of the problem, which examines the ownership of the atomistic angle forces in a CG bond (Fig. 4-9). The issue can be solved if we simply consider the shared angle ๐œƒ as an intermediary for the interactions between the two particles ๐ผ and ๐ฝ: thus, ๐น๐‘Ÿ๐‘Žโƒ—โƒ— โƒ—โƒ—   is the atomistic force on ๐ฝ by ๐ผ, while (๐น๐‘™๐‘Žโƒ—โƒ— โƒ—โƒ—   + ๐น๐‘๐‘Žโƒ—โƒ— โƒ—โƒ—  ) is the reaction by ๐ฝ on ๐ผ. Due to Newtonโ€™s third law, this clarification guarantees the conservation of momentum in the system. However, this question alludes to a more fundamental complication, that the valence (or torsion) angle interaction involves more than two bodies, so it is unclear how we can obtain a formula to determine its respective memory kernel. One possible approach is to apply to the standard MD decomposition of the force field on Eq. (17), the most general MZ EOM.      57  Bibliography    [1]  D. Frenkel and B. Smit, Understanding molecular simulation: from algorithms to applications, vol. 1, Academic press, 2001.  [2]  W. G. Noid, J. W. Chu, G. S. Ayton, V. Krishna, S. Izvekov, G. A. Voth and H. C. Andersen, "The multiscale coarse-graining method. I. A rigorous bridge between atomistic and coarse-grained models," The Journal of chemical physics, vol. 128, no. 24, p. 244114, 2008.  [3]  M. Neri, C. Anselmi, M. Cascella, A. Maritan and P. Carloni, "Coarse-grained model of proteins incorporating atomistic detail of the active site," Physical review letters, vol. 95, no. 21, p. 218102, 2005.  [4]  S. Kmiecik, D. Gront, M. Kolinski, L. Wieteska, A. E. Dawid and A. Kolinski, "Coarse-grained protein models and their applications," Chemical Reviews, vol. 116, no. 14, pp. 7898-7936, 2016.  [5]  B. J. Alder and T. E. Wainwright, "Studies in molecular dynamics. I. General method," The Journal of Chemical Physics, vol. 31, no. 2, pp. 459-466, 1959.  [6]  W. G. Noid, "Perspective: Coarse-grained models for biomolecular systems," The Journal of chemical physics, vol. 139, no. 9, p. 09B201_1, 2013.  [7]  S. J. Marrink, H. J. Risselada, S. Yefimov, D. P. Tieleman and A. H. De Vries, "The MARTINI force field: coarse grained model for biomolecular simulations," The journal of physical chemistry B, vol. 111, no. 27, pp. 7812-7824, 2007.  [8]  L. Monticelli, S. K. Kandasamy, X. Periole, R. G. Larson, D. P. Tieleman and S. J. Marrink, "The MARTINI coarse-grained force field: extension to proteins," Journal of chemical theory and computation, vol. 4, no. 5, pp. 819-834, 2008.  [9]  J. S. Andreu, C. Calero, J. Camacho and J. Faraudo, "On-the-fly coarse-graining methodology for the simulation of chain formation of superparamagnetic colloids in strong magnetic fields," Physical Review E, vol. 85, no. 3, p. 036709, 2012.  [10]  W. H. Wang, C. Schuฬˆtte and L. Delle Site, "Adaptive resolution simulation (AdResS): a smooth thermodynamic and structural transition from atomistic to coarse grained resolution and vice versa in a grand canonical fashion," Journal of chemical theory and computation, vol. 8, no. 8, pp. 2878-2887, 2012.  58  [11]  H. I. Ingรณlfsson, C. A. Lopez, J. J. Uusitalo, D. H. de Jong, S. M. Gopal, X. Periole and S. J. Marrink, "The power of coarse graining in biomolecular simulations," Wiley Interdisciplinary Reviews: Computational Molecular Science, vol. 4, no. 3, pp. 225-248, 2014.  [12]  R. Chudoba, J. Heyda and J. Dzubiella, "Temperature-Dependent Implicit-Solvent Model of Polyethylene Glycol in Aqueous Solution," Journal of chemical theory and computation, vol. 13, no. 12, pp. 6317-6327, 2017.  [13]  . Z. Li, X. Bian, X. Li and G. E. Karniadakis, "Incorporation of memory effects in coarse-grained modeling via the Mori-Zwanzig formalism," The Journal of chemical physics, vol. 143, no. 24, p. 243128, 2015.  [14]  C. Hijรณn, P. Espaรฑol, E. Vanden-Eijnden and R. Delgado-Buscalioni, "Moriโ€“Zwanzig formalism as a practical computational tool," Faraday discussions, vol. 144, pp. 301-322, 2010.  [15]  Y. Yoshimoto, I. Kinefuchi, T. Mima, A. Fukushima, T. Tokumasu and S. Takagi, "Bottom-up construction of interaction models of non-Markovian dissipative particle dynamics," Physical Review E, vol. 88, no. 4, p. 043305, 2013.  [16]  T. Kinjo and S. A. Hyodo, "Equation of motion for coarse-grained simulation based on microscopic description," Physical review E, vol. 75, no. 5, p. 051109, 2007.  [17]  Garcรญa-Palacios, J. Luis and F. J. Lรกzaro, "Langevin-dynamics study of the dynamical properties of small magnetic particles," Physical Review B, vol. 58, no. 22, p. 14937, 1998.  [18]  S. Izvekov, "Mori-Zwanzig theory for dissipative forces in coarse-grained dynamics in the Markov limit," Physical Review E, vol. 95, no. 1, p. 013303, 2017.  [19]  P. J. Hoogerbrugge and J. M. V. A. Koelman, "Simulating microscopic hydrodynamic phenomena with dissipative particle dynamics," EPL (Europhysics Letters), vol. 19, no. 3, p. 155, 1992.  [20]  P. Nikunen, M. Karttunen and I. Vattulainen, "How would you integrate the equations of motion in dissipative particle dynamics simulations?," Computer Physics Communications, vol. 153, no. 3, pp. 407-423, 2003.  [21]  Z. Li, H. S. Lee, E. Darve and G. E. Karniadakis, "Computing the non-Markovian coarse-grained interactions derived from the Moriโ€“Zwanzig formalism in molecular systems: Application to polymer melts," The Journal of chemical physics, vol. 146, no. 1, p. 014104, 2017.  [22]  "Re: [lammps-users] The temperature keeps rising with NVE, pair_style dpd, to simulate an amphiphilic polymeric system," 07 Jul 2016. [Online]. Available: http://lammps.sandia.gov/threads/msg63653.html. [Accessed 02 Feb 2018]. 59  [23]  M. J. Abraham, D. van der Spoel, E. Lindahl, D. Hess and and the GROMACS development team, "GROMACS User Manual version 5.0.4," 2014. [Online]. Available: ftp://ftp.gromacs.org/pub/manual/manual-5.0.4.pdf. [Accessed 02 Feb 2018]. [24]  M. Ceriotti, G. Bussi and M. Parrinello, "Colored-noise thermostats ร  la carte," Journal of Chemical Theory and Computation, pp. 1170-1180, 2010.  [25]  J. T. Padding and W. J. Briels, "Uncrossability constraints in mesoscopic polymer melt simulations: non-Rouse behavior of C 120 H 242," The Journal of Chemical Physics, vol. 115, no. 6, pp. 2846-2859, 2001.  [26]  J. A. Anderson, C. D. Lorenz and A. Travesset, "General purpose molecular dynamics simulations fully implemented on graphics processing units," Journal of Computational Physics, vol. 227, no. 10, pp. 5342-5359, 2008.  [27]  M. P. Allen and D. J. Tildesley, Computer simulation of liquids, Oxford university press, 2017.  [28]  W. C. Swope, H. C. Andersen, P. H. Berens and K. R. Wilson, "A computer simulation method for the calculation of equilibrium constants for the formation of physical clusters of molecules: Application to small water clusters," The Journal of Chemical Physics, vol. 76, no. 1, pp. 637-649, 1982.  [29]  R. Zwanzig, Nonequilibrium statistical mechanics, Oxford University Press, 2001.  [30]  G. Besold, I. Vattulainen, M. Karttunen and J. M. Polson, "Towards better integrators for dissipative particle dynamics simulations," Physical Review E, vol. 62, no. 6, p. R7611, 2000.  [31]  C. P. Lowe, "An alternative approach to dissipative particle dynamics," EPL (Europhysics Letters), vol. 47, no. 2, p. 145, 1999.  [32]  H. C. Andersen, "Molecular dynamics simulations at constant pressure and/or temperature," The Journal of chemical physics, vol. 72, no. 4, pp. 2384-2393, 1980.  [33]  N. Goga, A. J. Rzepiela, A. H. De Vries, S. J. Marrink and H. J. C. Berendsen, "Efficient algorithms for Langevin and DPD dynamics," Journal of chemical theory and computation, vol. 8, no. 10, pp. 3637-3649, 2012.  [34]  R. D. Groot and P. B. Warren, "Dissipative particle dynamics: Bridging the gap between atomistic and mesoscopic simulation," The Journal of chemical physics, vol. 107, no. 11, pp. 4423-4435, 1997.  [35]  P. Espanol and M. Serrano, "Dynamical regimes in the dissipative particle dynamics model," Physical Review E, vol. 59, no. 6, p. 6340, 1999.  [36]  T. C. Moore, C. R. Iacovella and C. McCabe, "Derivation of coarse-grained potentials via multistate iterative Boltzmann inversion," The Journal of chemical physics, vol. 140, no. 22, p. 06B606_1, 2014.  60  [37]  J. E. Jones, "On the determination of molecular fields.โ€”II. From the equation of state of a gas," Proceedings of the Royal Society of London A, vol. 106, no. 738, pp. 463-477, 1924.  [38]  W. G. Hoover, "Canonical dynamics: equilibrium phase-space distributions," Physical review A, vol. 31, no. 3, p. 1695, 1985.      61  Appendix A First, let us introduce some new functions of interest. Since we want to investigate the dynamical aspect of a CG mapping, we now care about the trajectory of the system. Let ๐›พ = ?ฬ‚?(๐‘ก) and ฮ“๐‘ƒ โ‰ก ๐‘ท(๐›พ) denote the trajectory of the atoms and of their COM in the atomistic simulation, respectively, with ๐‘ท being the mapping operator in Eq. (8). Obviously, ๐›พ is determined completely by the atomistic Hamiltonian โ„Ž( ๐‘Ÿ๐‘›โƒ—โƒ—โƒ—โƒ— , ๐‘๐‘›โƒ—โƒ— โƒ—โƒ—  ), that is while ฮ“๐‘ƒ can be deduced from ๐›พ using Eq. (8).  However, we are not interested in the time-evolution of the ฮ“๐‘ƒ, but in a CG EOM that replicates similar statistical properties. Namely, for any trajectory ฮ“ created using the CG EOM, ๐ด(ฮ“) has the same distribution as ๐ด(ฮ“๐‘ƒ), where ๐ด being an observable in CG phase space. There are different ways of deriving the MZ EOM [14], but in this appendix, we shall follow the work of Kinjo and Hyodo [16], which produces the desired formula through the time evolution of the phase space density function. Let us define the COM phase space density with ?ฬ‚? being the instantaneous coordinates of the atomistic system in the COM phase space ฮฆ. This definition allows us to express an observable of the COM as meaning ๐ด(๐‘ท(?ฬ‚?)) only depends on the atomistic coordinates {?ฬ‚?} through the mapping ๐‘ท in ฮฉ๐‘ƒ(?ฬ‚?). Then, we assume we can divide an atomistic observable ๐‘Ž(?ฬ‚?) into where ๐‘Ž๐’ซ(๐‘ท(?ฬ‚?)) can be expanded by the basis set {ฮฉ๐‘ƒ(?ฬ‚?)} as We can think of the operator ๐’ซ as the projection from ๐œ‘ into ฮฆ, while ๐’ฌ = 1 โˆ’ ๐’ซ is the operator orthogonal to ๐’ซ. Next, we define the scalar product โŒฉโ€ฆ ,โ€ฆ โŒช as the correlation function of the canonical ensemble, i.e. with ๐‘๐œ‘(?ฬ‚?) โˆ ๐‘’โˆ’๐›ฝโ„Ž  being the equilibrium probability density. Then, we apply another assumption  ๐‘‘๐‘‘๐‘ก๐‘Ÿ๐‘›โƒ—โƒ—โƒ—โƒ— (๐‘ก) =๐œ•โ„Ž๐œ•๐‘๐‘›โƒ—โƒ— โƒ—โƒ—   , ๐‘‘๐‘‘๐‘ก๐‘๐‘›โƒ—โƒ— โƒ—โƒ—  (๐‘ก) = โˆ’๐œ•โ„Ž๐œ•๐‘Ÿ๐‘›โƒ—โƒ—โƒ—โƒ—  , (89)  ฮฉ๐‘ƒ(?ฬ‚?) = ๐›ฟ(๐‘ท(?ฬ‚?) โˆ’ ?ฬ‚?) = ๐›ฟ (๐‘ท(๐‘Ÿ๐‘›โƒ—โƒ—โƒ—โƒ— ) โˆ’ ๐‘…๐‘โƒ—โƒ— โƒ—โƒ—  โƒ—) ๐›ฟ (๐‘ท(๐‘๐‘›โƒ—โƒ— โƒ—โƒ—  ) โˆ’ ๐‘ƒ๐‘โƒ—โƒ— โƒ—โƒ—  โƒ—) , (90)  ๐ด(๐‘ท(?ฬ‚?)) = โˆซ๐ด(?ฬ‚?)๐›ฟ(๐‘ท(?ฬ‚?) โˆ’ ?ฬ‚?)๐‘‘?ฬ‚? = โˆซ๐ด(?ฬ‚?)ฮฉ๐‘ƒ(?ฬ‚?)๐‘‘?ฬ‚? , (91)  ๐‘Ž(?ฬ‚?) = ๐‘Ž๐’ซ(๐‘ท(?ฬ‚?)) + ๐‘Ž๐’ฌ(?ฬ‚?) โ‰ก ๐’ซ๐‘Ž(?ฬ‚?) + ๐’ฌ๐‘Ž(?ฬ‚?) , (92)  ๐‘Ž๐’ซ(๐‘ท(?ฬ‚?)) = โˆซ๐ถ(?ฬ‚?)ฮฉ๐‘ƒ(?ฬ‚?)๐‘‘?ฬ‚? . (93)  โŒฉ๐‘Ž(?ฬ‚?), ๐‘(?ฬ‚?)โŒช = โˆซ๐‘Ž(?ฬ‚?)๐‘(?ฬ‚?)๐‘๐œ‘(?ฬ‚?)๐‘‘?ฬ‚? , (94) 62  and take the scalar product of ๐‘Ž(?ฬ‚?) with ฮฉ๐‘ƒ(?ฬ‚?) to get which gives an expression for the projected part of ๐‘Ž(?ฬ‚?) where โŒฉฮฉ๐‘ƒ(?ฬ‚?)โŒช = โˆซ ๐›ฟ(๐‘ท(?ฬ‚?) โˆ’ ?ฬ‚?)๐‘๐œ‘(?ฬ‚?)๐‘‘?ฬ‚? = ๐‘ƒ๐œ‘(?ฬ‚?) is the equilibrium distribution function of the mapped COM. Equation (97) shows the projection into the CG space is the thermal average of all atomistic states that correspond to the same CG state. We note that the two employed assumptions are reasonable: Eq. (93) simply means ๐‘Ž๐’ซ is the part of ๐‘Ž that can be described completely using the COM variables, while Eq. (95) implies the time scales of motion of the atoms and their corresponding CG particles are separable, or equivalently, Now, let us apply this formalism to the EOM of the phase space density of the COM. The time-evolution of ฮฉ๐‘ƒ(?ฬ‚?) along the atomistic trajectory is where the second equality occurs by grouping the atoms by their CG particle, ๐‘–๐‘™ is the atomistic Liouville operator and ๐น๐‘–โƒ—โƒ—  is the conservative force acting on the atom ๐‘–. Using the linearity of ๐‘ท and the property of the Dirac delta function, the first term in the last line becomes with ๐น๐ผ๐‘ƒโƒ—โƒ—โƒ—โƒ—  โƒ— = โˆ‘ ๐น๐‘–โƒ—โƒ— ๐‘–โˆˆ๐ผ  being total atomistic force acting on the particle ๐ผ. We can also transform the second term similarly to obtain the definition for the CG Liouville operator:  โŒฉ๐‘Ž๐’ฌ(?ฬ‚?), ฮฉ๐‘ƒ(?ฬ‚?)โŒช = 0 , (95)  โŒฉฮฉ๐‘ƒ(?ฬ‚?), ๐‘Ž(?ฬ‚?)โŒช = โŒฉฮฉ๐‘ƒ(?ฬ‚?), ๐‘Ž๐’ซ(๐‘ท(?ฬ‚?))โŒช = โˆซ๐›ฟ(๐‘ท(?ฬ‚?) โˆ’ ?ฬ‚?)๐‘Ž๐’ซ(๐‘ท(?ฬ‚?))๐‘๐œ‘(?ฬ‚?)๐‘‘?ฬ‚?= ๐‘Ž๐’ซ(?ฬ‚?)โˆซ๐›ฟ(๐‘ท(?ฬ‚?) โˆ’ ?ฬ‚?)๐‘๐œ‘(?ฬ‚?)๐‘‘?ฬ‚? , (96)  ๐‘Ž๐’ซ(?ฬ‚?) =โŒฉฮฉ๐‘ƒ(?ฬ‚?), ๐‘Ž(?ฬ‚?)โŒชโŒฉฮฉ๐‘ƒ(?ฬ‚?)โŒช=โˆซ ๐›ฟ(๐‘ท(?ฬ‚?) โˆ’ ?ฬ‚?)๐‘Ž(?ฬ‚?)๐‘๐œ‘(?ฬ‚?)๐‘‘?ฬ‚?โˆซ ๐›ฟ(๐‘ท(?ฬ‚?) โˆ’ ?ฬ‚?)๐‘๐œ‘(?ฬ‚?)๐‘‘?ฬ‚? , (97)  โŒฉ๐‘Ž๐’ฌ(?ฬ‚?), ๐‘๐’ซ(๐‘ท(?ฬ‚?))โŒช = 0 (98)      (๐‘‘๐‘‘๐‘ก)๐›พฮฉ๐‘ƒ(?ฬ‚?) = ๐‘–๐‘™ฮฉ๐‘ƒ(?ฬ‚?) โ‰ก โˆ’โˆ‘ [๐œ•โ„Ž๐œ•๐‘Ÿ๐‘– โƒ—โ‹…๐œ•๐œ•๐‘๐‘–โƒ—โƒ—โƒ—  โˆ’๐œ•โ„Ž๐œ•๐‘๐‘–โƒ—โƒ—โƒ—  โ‹…๐œ•๐œ•๐‘Ÿ๐‘– โƒ—] ฮฉ๐‘ƒ(?ฬ‚?)๐‘–= โˆ‘ โˆ‘[๐น๐‘–โƒ—โƒ— โ‹…๐œ•๐œ•๐‘๐‘–โƒ—โƒ—โƒ—  +๐‘๐‘–โƒ—โƒ—โƒ—  ๐‘š๐‘–โ‹…๐œ•๐œ•๐‘Ÿ๐‘– โƒ—] ๐›ฟ(๐‘ท(๐‘๐‘›โƒ—โƒ— โƒ—โƒ—  ) โˆ’ ๐‘ƒ๐ผโƒ—โƒ—  โƒ—)๐›ฟ(๐‘ท(๐‘Ÿ๐‘›โƒ—โƒ—โƒ—โƒ— ) โˆ’ ๐‘…๐ผโƒ—โƒ—โƒ—โƒ— )๐‘–โˆˆ๐ผ๐ผ= โˆ‘ โˆ‘[๐›ฟ(๐‘ท(๐‘Ÿ๐‘›โƒ—โƒ—โƒ—โƒ— ) โˆ’ ๐‘…๐ผโƒ—โƒ—โƒ—โƒ— )๐น๐‘–โƒ—โƒ— โ‹…๐œ•๐œ•๐‘๐‘–โƒ—โƒ—โƒ—  ๐›ฟ(๐‘ท(๐‘๐‘›โƒ—โƒ— โƒ—โƒ—  ) โˆ’ ๐‘ƒ๐ผโƒ—โƒ—  โƒ—) + ๐›ฟ(๐‘ท(๐‘๐‘›โƒ—โƒ— โƒ—โƒ—  ) โˆ’ ๐‘ƒ๐ผโƒ—โƒ—  โƒ—)๐‘๐‘–โƒ—โƒ—โƒ—  ๐‘š๐‘–โ‹…๐œ•๐œ•๐‘Ÿ๐‘– โƒ—๐›ฟ(๐‘ท(๐‘Ÿ๐‘›โƒ—โƒ—โƒ—โƒ— ) โˆ’ ๐‘…๐ผโƒ—โƒ—โƒ—โƒ— )]๐‘–โˆˆ๐ผ๐ผ , (99)  โˆ‘ [โˆ’๐›ฟ(๐‘ท(๐‘Ÿ๐‘›โƒ—โƒ—โƒ—โƒ— ) โˆ’ ๐‘…๐ผโƒ—โƒ—โƒ—โƒ— ) (โˆ‘๐น๐‘–โƒ—โƒ— ๐‘–โˆˆ๐ผ) โ‹…๐œ•๐œ•๐‘ƒ๐ผโƒ—โƒ—  โƒ—๐›ฟ(๐‘ท(๐‘๐‘›โƒ—โƒ— โƒ—โƒ—  ) โˆ’ ๐‘ƒ๐ผโƒ—โƒ—  โƒ—)]๐ผ= โˆ’โˆ‘๐น๐ผ๐‘ƒโƒ—โƒ—โƒ—โƒ—  โƒ— โ‹…๐œ•๐œ•๐‘ƒ๐ผโƒ—โƒ—  โƒ—ฮฉ๐‘ƒ(?ฬ‚?)๐ผ , (100)  (๐‘‘๐‘‘๐‘ก)๐›พฮฉ๐‘ƒ(?ฬ‚?) = โˆ’โˆ‘ [๐น๐ผ๐‘ƒโƒ—โƒ—โƒ—โƒ—  โƒ— โ‹…๐œ•๐œ•๐‘ƒ๐ผโƒ—โƒ—  โƒ—+๐‘ƒ๐ผ๐‘ƒโƒ—โƒ—โƒ—โƒ—  โƒ—๐‘€๐ผโ‹…๐œ•๐œ•๐‘…๐ผโƒ—โƒ—โƒ—โƒ— ] ฮฉ๐‘ƒ(?ฬ‚?)๐ผโ‰ก ๐‘–๐ฟฮฉ๐‘ƒ(?ฬ‚?) , (101) 63  in which ๐‘ƒ๐ผ๐‘ƒโƒ—โƒ—โƒ—โƒ—  โƒ— = โˆ‘ ๐‘๐‘–โƒ—โƒ—โƒ—  ๐‘–โˆˆ๐ผ = โˆ‘ ๐œ•โ„Ž/๐œ•๐‘Ÿ๐‘– โƒ—๐‘–โˆˆ๐ผ  is the total momentum of particle ๐ผ. From Eq. (97), we find the projected part of the time-evolution to be We want to derive the explicit form for  The first scalar product is Using ๐œ“(?ฬ‚?) = ๐‘’โˆ’๐›ฝโ„Ž/๐‘ and integration by parts, the inner integral can be rewritten as which gives since ๐œ•๐œ•๐‘…๐ผโƒ—โƒ— โƒ—โƒ— โ€ฒ๐›ฟ(?ฬ‚?โ€ฒ โˆ’ ?ฬ‚?) = โˆ’๐œ•๐œ•๐‘…๐ผโƒ—โƒ— โƒ—โƒ— ๐›ฟ(?ฬ‚?โ€ฒ โˆ’ ?ฬ‚?). Similarly, we have  ๐‘“๐’ซ(?ฬ‚?) โ‰ก ๐’ซ๐‘–๐ฟฮฉ๐‘ƒ(?ฬ‚?) =โŒฉฮฉ๐‘ƒ(?ฬ‚?โ€ฒ), ๐‘–๐ฟฮฉ๐‘ƒ(?ฬ‚?)โŒชโŒฉฮฉ๐‘ƒ(?ฬ‚?โ€ฒ)โŒช=โˆซ ๐›ฟ(๐‘ท(?ฬ‚?) โˆ’ ?ฬ‚?โ€ฒ)๐‘–๐ฟ๐›ฟ(๐‘ท(?ฬ‚?) โˆ’ ?ฬ‚?)๐‘๐œ‘(?ฬ‚?)๐‘‘?ฬ‚?โˆซ ๐›ฟ(๐‘ท(?ฬ‚?) โˆ’ ?ฬ‚?โ€ฒ)๐‘๐œ‘(?ฬ‚?)๐‘‘?ฬ‚? . (102)  โŒฉฮฉ๐‘ƒ(?ฬ‚?โ€ฒ), ๐‘–๐ฟ๐‘ƒฮฉ๐‘ƒ(?ฬ‚?)โŒช = โŸจฮฉ๐‘ƒ(?ฬ‚?โ€ฒ), โˆ’โˆ‘[๐น๐ผ๐‘ƒโƒ—โƒ—โƒ—โƒ—  โƒ— โ‹…๐œ•๐œ•๐‘ƒ๐ผโƒ—โƒ—  โƒ—+๐‘ƒ๐ผ๐‘ƒโƒ—โƒ—โƒ—โƒ—  โƒ—๐‘€๐ผโ‹…๐œ•๐œ•๐‘…๐ผโƒ—โƒ—โƒ—โƒ— ] ฮฉ๐‘ƒ(?ฬ‚?)IโŸฉ= โˆ’โˆ‘๐œ•๐œ•๐‘ƒ๐ผโƒ—โƒ—  โƒ—โ‹… โŒฉ๐น๐ผ๐‘ƒโƒ—โƒ—โƒ—โƒ—  โƒ—ฮฉ๐‘ƒ(?ฬ‚?โ€ฒ), ฮฉ๐‘ƒ(?ฬ‚?)โŒช๐ผโˆ’ โˆ‘๐œ•๐œ•๐‘…๐ผโƒ—โƒ—โƒ—โƒ—  โ‹… โŒฉ๐‘ƒ๐ผ๐‘ƒโƒ—โƒ—โƒ—โƒ—  โƒ—๐‘€๐ผฮฉ๐‘ƒ(?ฬ‚?โ€ฒ), ฮฉ๐‘ƒ(?ฬ‚?)โŒช๐ผ . (103)  โŒฉ๐น๐ผ๐‘ƒโƒ—โƒ—โƒ—โƒ—  โƒ—ฮฉ๐‘ƒ(?ฬ‚?โ€ฒ), ฮฉ๐‘ƒ(?ฬ‚?)โŒช = โˆซ(โˆ’โˆ‘๐œ•โ„Ž๐œ•๐‘Ÿ๐‘– โƒ—๐‘–โˆˆ๐ผ) ๐›ฟ(๐‘ท(?ฬ‚?) โˆ’ ?ฬ‚?โ€ฒ)๐›ฟ(๐‘ท(?ฬ‚?) โˆ’ ?ฬ‚?)๐‘๐œ‘(?ฬ‚?)๐‘‘?ฬ‚?= โˆซ๐›ฟ (๐‘ท(๐‘๐‘›โƒ—โƒ— โƒ—โƒ—  ) โˆ’ ๐‘ƒ๐‘โƒ—โƒ— โƒ—โƒ—  โƒ—โ€ฒ) ๐›ฟ (๐‘ท(๐‘๐‘›โƒ—โƒ— โƒ—โƒ—  ) โˆ’ ๐‘ƒ๐‘โƒ—โƒ— โƒ—โƒ—  โƒ—) ๐‘‘๐‘๐‘›โƒ—โƒ— โƒ—โƒ—   ร— โˆซ๐›ฟ (๐‘ท(๐‘Ÿ๐‘›โƒ—โƒ—โƒ—โƒ— ) โˆ’ ๐‘…๐‘โƒ—โƒ— โƒ—โƒ—  โƒ—โ€ฒ) ๐›ฟ (๐‘ท(๐‘Ÿ๐‘›โƒ—โƒ—โƒ—โƒ— ) โˆ’ ๐‘…๐‘โƒ—โƒ— โƒ—โƒ—  โƒ—) (โˆ’โˆ‘๐œ•โ„Ž๐œ•๐‘Ÿ๐‘– โƒ—๐‘–โˆˆ๐ผ)๐‘๐œ‘(?ฬ‚?)๐‘‘๐‘Ÿ๐‘›โƒ—โƒ—โƒ—โƒ—  . (104)      โˆซ ๐›ฟ (๐‘ท(๐‘Ÿ๐‘›โƒ—โƒ—โƒ—โƒ— ) โˆ’ ๐‘…๐‘โƒ—โƒ— โƒ—โƒ—  โƒ—โ€ฒ) ๐›ฟ (๐‘ท(๐‘Ÿ๐‘›โƒ—โƒ—โƒ—โƒ— ) โˆ’ ๐‘…๐‘โƒ—โƒ— โƒ—โƒ—  โƒ—)1๐›ฝโˆ‘๐œ•๐‘๐œ‘(?ฬ‚?)๐œ•๐‘Ÿ๐‘– โƒ—๐‘–โˆˆ๐ผ๐‘‘๐‘Ÿ๐‘›โƒ—โƒ—โƒ—โƒ— = โˆ’1๐›ฝโˆซ๐‘๐œ‘(?ฬ‚?) (โˆ‘๐œ•๐œ•๐‘Ÿ๐‘– โƒ—๐‘–โˆˆ๐ผ) [๐›ฟ (๐‘ท(๐‘Ÿ๐‘›โƒ—โƒ—โƒ—โƒ— ) โˆ’ ๐‘…๐‘โƒ—โƒ— โƒ—โƒ—  โƒ—โ€ฒ) ๐›ฟ (๐‘ท(๐‘Ÿ๐‘›โƒ—โƒ—โƒ—โƒ— ) โˆ’ ๐‘…๐‘โƒ—โƒ— โƒ—โƒ—  โƒ—)] ๐‘‘๐‘Ÿ๐‘›โƒ—โƒ—โƒ—โƒ— =1๐›ฝโˆซ๐‘๐œ‘(?ฬ‚?) (๐œ•๐œ•๐‘…๐ผโƒ—โƒ—โƒ—โƒ— โ€ฒ+๐œ•๐œ•๐‘…๐ผโƒ—โƒ—โƒ—โƒ— ) [๐›ฟ (๐‘ท(๐‘Ÿ๐‘›โƒ—โƒ—โƒ—โƒ— ) โˆ’ ๐‘…๐‘โƒ—โƒ— โƒ—โƒ—  โƒ—โ€ฒ) ๐›ฟ (๐‘ท(๐‘Ÿ๐‘›โƒ—โƒ—โƒ—โƒ— ) โˆ’ ๐‘…๐‘โƒ—โƒ— โƒ—โƒ—  โƒ—)] ๐‘‘๐‘Ÿ๐‘›โƒ—โƒ—โƒ—โƒ—  , (105)  โŒฉ๐น๐ผ๐‘ƒโƒ—โƒ—โƒ—โƒ—  โƒ—ฮฉ๐‘ƒ(?ฬ‚?โ€ฒ), ฮฉ๐‘ƒ(?ฬ‚?)โŒช =1๐›ฝ(๐œ•๐œ•๐‘…๐ผโƒ—โƒ—โƒ—โƒ— โ€ฒ+๐œ•๐œ•๐‘…๐ผโƒ—โƒ—โƒ—โƒ— )โˆซ๐›ฟ(๐‘ท(?ฬ‚?) โˆ’ ?ฬ‚?โ€ฒ)๐›ฟ(๐‘ท(?ฬ‚?) โˆ’ ?ฬ‚?)๐‘๐œ‘(?ฬ‚?)๐‘‘?ฬ‚?=1๐›ฝ(๐œ•๐œ•๐‘…๐ผโƒ—โƒ—โƒ—โƒ— โ€ฒ+๐œ•๐œ•๐‘…๐ผโƒ—โƒ—โƒ—โƒ— ) ๐›ฟ(?ฬ‚?โ€ฒ โˆ’ ?ฬ‚?)๐‘ƒ๐œ‘(?ฬ‚?)=1๐›ฝ๐›ฟ(?ฬ‚?โ€ฒ โˆ’ ?ฬ‚?)๐œ•๐œ•๐‘…๐ผโƒ—โƒ—โƒ—โƒ— ๐‘ƒ๐œ‘(?ฬ‚?) , (106) 64  Thus, the projected portion of the CG Liouville equation is Next, we consider the time evolution for the rest of the CG Liouville equation, namely ๐’ฌ๐‘–๐ฟฮฉ๐‘ƒ, along ๐›พ In contrast to the projected term, the orthogonal term can only be described completely using the full atomistic phase space coordinates, so its time evolution depends explicitly on ๐›พ and therefore it is a function of ๐‘ก. Hence, we shall relabel it as ๐‘„๐‘–๐ฟฮฉ๐‘ƒ(?ฬ‚?) โ‰ก ๐‘“๐’ฌ(?ฬ‚?, ๐‘ก), where we define ๐‘ก = 0 as the point in time when the projection takes place3. Then, the formal solution of Eq. (126) is where we have applied the following identity for linear operators: Now, we consider in which the second equality occurs due to the Hermicity of ๐‘–๐‘™, the fourth equality due to Eq. (98). The last line only involves thermal averages, so it does not depend on time and thus is constant under the effect of ๐‘’โˆ’๐‘–๐‘™(๐‘กโˆ’๐‘ ). Then, Eq. (110) is equivalent to                                                           3 We did not make this clear earlier, but the projection ๐’ซ is necessarily a projection at a fixed time. This, however, does not affect the result about ๐’ซ in the previous paragraphs, since by the ergodic hypothesis, thermal averages in a simulation do not depend on time. Nevertheless, ๐‘„ is not an average, so ๐‘ก is relevant here.  โŒฉ๐‘ƒ๐ผ๐‘ƒโƒ—โƒ—โƒ—โƒ—  โƒ—๐‘€๐ผฮฉ๐‘ƒ(?ฬ‚?โ€ฒ), ฮฉ๐‘ƒ(?ฬ‚?)โŒช = โˆ’1๐›ฝ๐›ฟ(?ฬ‚?โ€ฒ โˆ’ ?ฬ‚?)๐œ•๐œ•๐‘ƒ๐ผโƒ—โƒ—  โƒ—๐‘ƒ๐œ‘(?ฬ‚?) . (107)  ๐‘“๐’ซ(?ฬ‚?) = โˆ’1๐›ฝ1๐‘ƒ๐œ‘(?ฬ‚?โ€ฒ)โˆ‘ [๐œ•๐œ•๐‘ƒ๐ผโƒ—โƒ—  โƒ—โ‹… (๐›ฟ(?ฬ‚?โ€ฒ โˆ’ ?ฬ‚?)๐œ•๐œ•๐‘…๐ผโƒ—โƒ—โƒ—โƒ— ๐‘ƒ๐œ‘(?ฬ‚?)) โˆ’๐œ•๐œ•๐‘…๐ผโƒ—โƒ—โƒ—โƒ—  โ‹… (๐›ฟ(?ฬ‚?โ€ฒ โˆ’ ?ฬ‚?)๐œ•๐œ•๐‘ƒ๐ผโƒ—โƒ—  โƒ—๐‘ƒ๐œ‘(?ฬ‚?))]๐ผ . (108)  (๐‘‘๐‘‘๐‘ก)๐›พ๐’ฌ๐‘–๐ฟฮฉ๐‘ƒ(?ฬ‚?) = โˆ’๐‘–๐‘™๐’ฌ๐‘–๐ฟฮฉ๐‘ƒ(?ฬ‚?) . (109)  ๐‘“๐’ฌ(?ฬ‚?, ๐‘ก) = ๐‘’โˆ’๐‘–๐‘™๐‘ก๐‘“๐’ฌ(?ฬ‚?, 0) = ๐‘’โˆ’(๐’ซ+๐’ฌ)๐‘–๐‘™๐‘ก๐‘“๐’ฌ(?ฬ‚?, 0)= ๐‘’โˆ’๐’ฌ๐‘–๐‘™๐‘ก๐‘“๐’ฌ(?ฬ‚?, 0) โˆ’ โˆซ ๐‘’โˆ’๐‘–๐‘™(๐‘กโˆ’๐‘ )๐’ซ๐‘–๐‘™๐‘’โˆ’๐’ฌ๐‘–๐‘™๐‘ ๐‘“๐’ฌ(?ฬ‚?, 0)๐‘‘๐‘ ๐‘ก0 , (110)  ๐‘’(๐’œ+โ„ฌ)๐‘ก = ๐‘’โ„ฌ๐‘ก + โˆซ ๐‘’(๐’œ+โ„ฌ)(๐‘กโˆ’๐‘ )๐’œ๐‘’โ„ฌ๐‘ ๐‘‘๐‘ ๐‘ก0 . (111)  โˆ’๐’ซ๐‘–๐‘™๐‘’โˆ’๐’ฌ๐‘–๐‘™๐‘ ๐‘“๐’ฌ(?ฬ‚?, 0) =โˆ’โŸจฮฉ๐‘ƒ(?ฬ‚?โ€ฒ), ๐‘–๐‘™๐‘’โˆ’๐’ฌ๐‘–๐‘™๐‘ ๐‘“๐’ฌ(?ฬ‚?, 0)โŸฉโŸจฮฉ๐‘ƒ(?ฬ‚?โ€ฒ)โŸฉ=โŸจ๐‘–๐‘™ฮฉ๐‘ƒ(?ฬ‚?โ€ฒ), ๐‘’โˆ’๐’ฌ๐‘–๐‘™๐‘ ๐‘“๐’ฌ(?ฬ‚?, 0)โŸฉโŸจฮฉ๐‘ƒ(?ฬ‚?โ€ฒ)โŸฉ=  โŸจ๐‘–๐ฟฮฉ๐‘ƒ(?ฬ‚?โ€ฒ), ๐‘’โˆ’๐’ฌ๐‘–๐‘™๐‘ ๐’ฌ๐‘–๐ฟฮฉ๐‘ƒ(?ฬ‚?)โŸฉโŸจฮฉ๐‘ƒ(?ฬ‚?โ€ฒ)โŸฉ=โŸจ๐’ฌ๐‘–๐ฟฮฉ๐‘ƒ(?ฬ‚?โ€ฒ), ๐‘’โˆ’๐’ฌ๐‘–๐‘™๐‘ ๐’ฌ๐‘–๐ฟฮฉ๐‘ƒ(?ฬ‚?)โŸฉโŸจฮฉ๐‘ƒ(?ฬ‚?โ€ฒ)โŸฉ=1๐‘ƒ๐œ‘(?ฬ‚?โ€ฒ)โŸจ๐‘“๐’ฌ(?ฬ‚?โ€ฒ, 0), ๐‘’โˆ’๐’ฌ๐‘–๐‘™๐‘ ๐‘“๐’ฌ(?ฬ‚?, 0)โŸฉ , (112) 65   We have now derived the explicit form for the EOM of the CG phase space density function: To obtain the time evolution of the particle ๐ฝ, we need to multiply the former EOM with ๐‘ƒ๐ฝโƒ—โƒ—  โƒ— then integrate over ฮฆ. For example, using integration by parts, the left-hand side of Eq. (101) is while the projected term becomes Here, we remark this expression yields the CG conservative force by using the consistency conditions in Eq. (12) and (13), i.e. Next, for the orthogonal term, we note that  ๐‘“๐’ฌ(?ฬ‚?, ๐‘ก) = ๐‘’โˆ’๐’ฌ๐‘–๐‘™๐‘ก๐‘“๐’ฌ(?ฬ‚?, 0) +1๐‘ƒ๐œ‘(?ฬ‚?โ€ฒ)โˆซ โŸจ๐‘“๐’ฌ(?ฬ‚?โ€ฒ, 0), ๐‘’โˆ’๐’ฌ๐‘–๐‘™๐‘ ๐‘“๐’ฌ(?ฬ‚?, 0)โŸฉ๐‘‘๐‘ ๐‘ก0 . (113)  (๐‘‘๐‘‘๐‘ก)๐›พฮฉ๐‘ƒ(?ฬ‚?) = ๐‘“๐’ซ(?ฬ‚?) + ๐‘“๐’ฌ(?ฬ‚?, ๐‘ก)= โˆ’1๐›ฝ1๐‘ƒ๐œ‘(?ฬ‚?โ€ฒ)โˆ‘[๐œ•๐œ•๐‘ƒ๐ผโƒ—โƒ—  โƒ—โ‹… (๐›ฟ(?ฬ‚?โ€ฒ โˆ’ ?ฬ‚?)๐œ•๐œ•๐‘…๐ผโƒ—โƒ—โƒ—โƒ— ๐‘ƒ๐œ‘(?ฬ‚?)) +๐œ•๐œ•๐‘…๐ผโƒ—โƒ—โƒ—โƒ—  โ‹… ((?ฬ‚?โ€ฒ โˆ’ ?ฬ‚?)๐œ•๐œ•๐‘ƒ๐ผโƒ—โƒ—  โƒ—๐‘ƒ๐œ‘(?ฬ‚?))]๐ผ+1๐‘ƒ๐œ‘(?ฬ‚?โ€ฒ)โˆซ โŸจ๐‘“๐’ฌ(?ฬ‚?โ€ฒ, 0), ๐‘’โˆ’๐’ฌ๐‘–๐‘™๐‘ ๐‘“๐’ฌ(?ฬ‚?, 0)โŸฉ๐‘‘๐‘ ๐‘ก0+ ๐‘’โˆ’๐’ฌ๐‘–๐‘™๐‘ก๐‘“๐’ฌ(?ฬ‚?, 0) . (114)  โˆซ๐‘ƒ๐ฝโƒ—โƒ—  โƒ— (๐‘‘๐‘‘๐‘ก)๐›พฮฉ๐‘ƒ(?ฬ‚?)๐‘‘?ฬ‚? = โˆซ(๐‘‘๐‘‘๐‘ก๐‘ƒ๐ฝโƒ—โƒ—  โƒ—) ฮฉ๐‘ƒ(?ฬ‚?)๐‘‘?ฬ‚? =๐‘‘๐‘‘๐‘ก๐‘ƒ๐ฝโƒ—โƒ—  โƒ— . (115)      โˆซ๐‘ƒ๐ฝโƒ—โƒ—  โƒ—๐‘“๐’ซ(?ฬ‚?)๐‘‘?ฬ‚? = โˆซ๐‘ƒ๐ฝโƒ—โƒ—  โƒ—๐’ซ๐‘–๐ฟฮฉ๐‘ƒ(?ฬ‚?)๐‘‘?ฬ‚?= โˆ’1๐›ฝโˆซ๐‘ƒ๐ฝโƒ—โƒ—  โƒ—1๐‘ƒ๐œ‘(?ฬ‚?โ€ฒ)โˆ‘ [๐œ•๐œ•๐‘ƒ๐ผโƒ—โƒ—  โƒ—โ‹… (๐›ฟ(?ฬ‚?โ€ฒ โˆ’ ?ฬ‚?)๐œ•๐œ•๐‘…๐ผโƒ—โƒ—โƒ—โƒ— ๐‘ƒ๐œ‘(?ฬ‚?)) +๐œ•๐œ•๐‘…๐ผโƒ—โƒ—โƒ—โƒ—  โ‹… ((?ฬ‚?โ€ฒ โˆ’ ?ฬ‚?)๐œ•๐œ•๐‘ƒ๐ผโƒ—โƒ—  โƒ—๐‘ƒ๐œ‘(?ฬ‚?))]๐ผ๐‘‘?ฬ‚?=1๐›ฝโˆซ1๐‘ƒ๐œ‘(?ฬ‚?โ€ฒ)(๐œ•๐œ•๐‘ƒ๐ฝโƒ—โƒ—  โƒ—โ‹… ๐‘ƒ๐ฝโƒ—โƒ—  โƒ—) ๐›ฟ(?ฬ‚?โ€ฒ โˆ’ ?ฬ‚?)๐œ•๐œ•๐‘…๐ฝโƒ—โƒ—โƒ—โƒ— ๐‘ƒ๐œ‘(?ฬ‚?) ๐‘‘?ฬ‚?=1๐›ฝ1๐‘ƒ๐œ‘(?ฬ‚?โ€ฒ)๐œ•๐œ•๐‘…๐ฝโ€ฒโƒ—โƒ—โƒ—โƒ— ๐‘ƒ๐œ‘(?ฬ‚?โ€ฒ) =1๐›ฝ๐œ•๐œ•๐‘…๐ฝโƒ—โƒ—โƒ—โƒ— ln ๐‘ƒ๐œ‘(?ฬ‚?) . (116)  1๐›ฝ๐œ•๐œ•๐‘…๐ฝโƒ—โƒ—โƒ—โƒ— ln ๐‘ƒ๐œ‘(?ฬ‚?) =1๐›ฝ๐œ•๐œ•๐‘…๐ฝโƒ—โƒ—โƒ—โƒ— ln ๐‘ƒฮฆ(?ฬ‚?) =1๐›ฝ๐œ•๐œ•๐‘…๐ฝโƒ—โƒ—โƒ—โƒ— ln [๐‘ƒ๐‘… (๐‘…๐‘โƒ—โƒ— โƒ—โƒ—  โƒ—) ๐‘ƒ๐‘ƒ (๐‘ƒ๐‘โƒ—โƒ— โƒ—โƒ—  โƒ—)]=1๐›ฝ๐œ•๐œ•๐‘…๐ฝโƒ—โƒ—โƒ—โƒ— ln ๐‘ƒ๐‘… (๐‘…๐‘โƒ—โƒ— โƒ—โƒ—  โƒ—) = โˆ’๐œ•๐œ•๐‘…๐ฝโƒ—โƒ—โƒ—โƒ— ๐ป(?ฬ‚?) โ‰ก ๐น๐ฝโƒ—โƒ—  โƒ— . (117)  โˆซ๐‘ƒ๐ฝโƒ—โƒ—  โƒ—๐‘–๐ฟฮฉ๐‘ƒ(?ฬ‚?)๐‘‘?ฬ‚? = โˆ’โˆซโˆ‘ ๐‘ƒ๐ฝโƒ—โƒ—  โƒ— [๐น๐ผ๐‘ƒโƒ—โƒ—โƒ—โƒ—  โƒ— โ‹…๐œ•๐œ•๐‘ƒ๐ผโƒ—โƒ—  โƒ—+๐‘ƒ๐ผ๐‘ƒโƒ—โƒ—โƒ—โƒ—  โƒ—๐‘€๐ผโ‹…๐œ•๐œ•๐‘…๐ผโƒ—โƒ—โƒ—โƒ— ] ๐›ฟ(?ฬ‚?โ€ฒ โˆ’ ?ฬ‚?)๐‘‘?ฬ‚?๐ผ= โˆซ๐น๐ฝ๐‘ƒโƒ—โƒ—โƒ—โƒ—  โƒ— โ‹… (๐œ•๐œ•๐‘ƒ๐ฝโƒ—โƒ—  โƒ—๐‘ƒ๐ฝโƒ—โƒ—  โƒ—) ๐›ฟ(?ฬ‚?โ€ฒ โˆ’ ?ฬ‚?)๐‘‘?ฬ‚? = ๐น๐ฝ๐‘ƒโƒ—โƒ—โƒ—โƒ—  โƒ— , (118) 66  so solving the same computation as in Eq. (115) and (116) for ๐‘“๐’ฌ(?ฬ‚?โ€ฒ, 0) gives which we have defined earlier as the atomistic fluctuating force. Moreover, by applying the consistency conditions to Eq. (108), we have therefore in which we define ๐›ฟ๐น๐ฝ๐’ฌโƒ—โƒ—โƒ—โƒ—  โƒ—(๐‘ก) โ‰ก ๐‘’โˆ’๐’ฌ๐‘–๐‘™๐‘ก๐›ฟ๐น๐ฝ๐’ฌโƒ—โƒ—โƒ—โƒ—  โƒ—(0). Thus, we obtain the MZ EOM:  Finally, we want to clarify the difference between ๐›ฟ๐น๐ผ๐“ โƒ—โƒ—โƒ—โƒ—  โƒ—(๐‘ ) and ๐›ฟ๐น๐ผโƒ—โƒ—  โƒ—(๐‘ ) and fulfill our promise in the beginning of Section 4.2. Since ๐›ฟ๐น๐ผโƒ—โƒ—  โƒ—(๐‘ )  โ‰ก [๐น๐ผ๐‘ƒโƒ—โƒ—โƒ—โƒ—  โƒ— โˆ’ ๐น๐ผโƒ—โƒ—  โƒ—]๐‘ก=๐‘  is the deviation of the total atomistic force from the mean force, evaluated from the atomistic trajectory at time ๐‘ , it evolves using the atomistic Liouville propagator, that is On the other hand, we have  โˆซ๐‘ƒ๐ฝโƒ—โƒ—  โƒ—๐‘“๐’ฌ(?ฬ‚?, 0)๐‘‘?ฬ‚? = โˆซ๐‘ƒ๐ฝโƒ—โƒ—  โƒ—๐’ฌ๐‘–๐ฟฮฉ๐‘ƒ(?ฬ‚?)๐‘‘?ฬ‚?|๐‘ก=0= [โˆซ๐‘ƒ๐ฝโƒ—โƒ—  โƒ—๐‘–๐ฟฮฉ๐‘ƒ(?ฬ‚?)๐‘‘?ฬ‚? โˆ’ โˆซ๐‘ƒ๐ฝโƒ—โƒ—  โƒ—๐’ซ๐‘–๐ฟฮฉ๐‘ƒ(?ฬ‚?)๐‘‘?ฬ‚?]๐‘ก=0= [๐น๐ฝ๐‘ƒโƒ—โƒ—โƒ—โƒ—  โƒ— โˆ’ ๐น๐ฝโƒ—โƒ—  โƒ—]๐‘ก=0โ‰ก ๐›ฟ๐น๐ฝ๐’ฌโƒ—โƒ—โƒ—โƒ—  โƒ—(0) , (119)  ๐’ซ๐‘–๐ฟฮฉ๐‘ƒ(?ฬ‚?) = โˆ’1๐›ฝ1๐‘ƒฮฆ(?ฬ‚?โ€ฒ)โˆ‘[๐œ•๐œ•๐‘ƒ๐ผโƒ—โƒ—  โƒ—โ‹… (๐›ฟ(?ฬ‚?โ€ฒ โˆ’ ?ฬ‚?)๐œ•๐œ•๐‘…๐ผโƒ—โƒ—โƒ—โƒ— ๐‘ƒฮฆ(?ฬ‚?)) โˆ’๐œ•๐œ•๐‘…๐ผโƒ—โƒ—โƒ—โƒ—  โ‹… (๐›ฟ(?ฬ‚?โ€ฒ โˆ’ ?ฬ‚?)๐œ•๐œ•๐‘ƒ๐ผโƒ—โƒ—  โƒ—๐‘ƒฮฆ(?ฬ‚?))]๐ผ= โˆ’1๐›ฝ1๐‘ƒฮฆ(?ฬ‚?โ€ฒ)โˆ‘[๐œ•๐œ•๐‘…๐ผโƒ—โƒ—โƒ—โƒ— ๐‘ƒฮฆ(?ฬ‚?) โ‹…๐œ•๐œ•๐‘ƒ๐ผโƒ—โƒ—  โƒ—๐›ฟ(?ฬ‚?โ€ฒ โˆ’ ?ฬ‚?) โˆ’๐œ•๐œ•๐‘ƒ๐ผโƒ—โƒ—  โƒ—๐‘ƒฮฆ(?ฬ‚?) โ‹…๐œ•๐œ•๐‘…๐ผโƒ—โƒ—โƒ—โƒ—  ๐›ฟ(?ฬ‚?โ€ฒ โˆ’ ?ฬ‚?)]๐ผ=๐‘ƒฮฆ(?ฬ‚?)๐‘ƒฮฆ(?ฬ‚?โ€ฒ)โˆ‘[๐น๐ผโƒ—โƒ—  โƒ— โ‹…๐œ•๐œ•๐‘ƒ๐ผโƒ—โƒ—  โƒ—๐›ฟ(?ฬ‚?โ€ฒ โˆ’ ?ฬ‚?) +๐‘ƒ๐ผโƒ—โƒ—  โƒ—๐‘€๐ผโ‹…๐œ•๐œ•๐‘…๐ผโƒ—โƒ—โƒ—โƒ—  ๐›ฟ(?ฬ‚?โ€ฒ โˆ’ ?ฬ‚?)]๐ผ , (120)      โˆซ ๐‘ƒ๐ฝโƒ—โƒ—  โƒ—1๐‘ƒ๐œ‘(?ฬ‚?โ€ฒ)โŸจ๐‘“๐’ฌ(?ฬ‚?โ€ฒ, 0), ๐‘’โˆ’๐’ฌ๐‘–๐‘™๐‘ ๐‘“๐’ฌ(?ฬ‚?, 0)โŸฉ๐‘‘?ฬ‚?=1๐‘ƒฮฆ(?ฬ‚?โ€ฒ)โŸจ๐’ฌ๐‘–๐ฟฮฉ๐‘ƒ(?ฬ‚?โ€ฒ), ๐‘’โˆ’๐’ฌ๐‘–๐‘™๐‘  โˆซ๐‘ƒ๐ฝโƒ—โƒ—  โƒ— ๐‘“๐’ฌ(?ฬ‚?, 0)๐‘‘?ฬ‚?โŸฉ =1๐‘ƒฮฆ(?ฬ‚?โ€ฒ)โŸจ(1 โˆ’ ๐’ซ)๐‘–๐ฟฮฉ๐‘ƒ(?ฬ‚?โ€ฒ), ๐‘’โˆ’๐’ฌ๐‘–๐‘™๐‘ ๐›ฟ๐น๐ฝ๐’ฌโƒ—โƒ—โƒ—โƒ—  โƒ—(0)โŸฉ=1๐‘ƒฮฆ(?ฬ‚?โ€ฒ)โˆ‘ โŸจ[(๐‘ƒฮฆ(?ฬ‚?โ€ฒ)๐‘ƒฮฆ(?ฬ‚?โ€ฒโ€ฒ)๐น๐ผโ€ฒโƒ—โƒ—  โƒ— โˆ’ ๐น๐ผ๐‘ƒโ€ฒโƒ—โƒ— โƒ—โƒ—โƒ—โƒ— ) โ‹…๐œ•๐œ•๐‘ƒ๐ผโ€ฒโƒ—โƒ—  โƒ—+ (๐‘ƒฮฆ(?ฬ‚?โ€ฒ)๐‘ƒฮฆ(?ฬ‚?โ€ฒโ€ฒ)๐‘ƒ๐ผโ€ฒโƒ—โƒ—  โƒ—๐‘€๐ผโˆ’๐‘ƒ๐ผ๐‘ƒโ€ฒโƒ—โƒ—โƒ—โƒ— โƒ—โƒ— ๐‘€๐ผ) โ‹…๐œ•๐œ•๐‘…๐ผโ€ฒโƒ—โƒ—โƒ—โƒ—  ] ๐›ฟ(?ฬ‚?โ€ฒโ€ฒ โˆ’ ?ฬ‚?โ€ฒ), ๐›ฟ๐น๐ฝ๐’ฌโƒ—โƒ—โƒ—โƒ—  โƒ—(๐‘ )โŸฉ๐ผ= โˆ’๐›ฝ โˆ‘โŒฉ๐›ฟ๐น๐ผ๐“ โƒ—โƒ—โƒ—โƒ—  โƒ—(๐‘ ) โŠ— ๐›ฟ๐น๐ฝ๐“ โƒ—โƒ—โƒ—โƒ—  โƒ—(0)โŒช๐ผ๐‘ƒ๐ผโ€ฒโƒ—โƒ—  โƒ—๐‘€๐ผ , (121)  ๐‘‘๐‘‘๐‘ก๐‘ƒ๐ฝโƒ—โƒ—  โƒ— = ๐น๐ฝโƒ—โƒ—  โƒ— โˆ’ ๐›ฝ โˆ‘ โˆซ โŒฉ๐›ฟ๐น๐ผ๐“ โƒ—โƒ—โƒ—โƒ—  โƒ—(๐‘ก โˆ’ ๐‘ ) โŠ— ๐›ฟ๐น๐ฝ๐“ โƒ—โƒ—โƒ—โƒ—  โƒ—(0)โŒช๐‘ƒ๐ผโƒ—โƒ—  โƒ—(๐‘ )๐‘€๐ผ๐‘‘๐‘ ๐‘ก0๐ผ+ ๐›ฟ๐น๐ผ๐“ โƒ—โƒ—โƒ—โƒ—  โƒ—(๐‘ก) . (122)  ๐›ฟ๐น๐ผโƒ—โƒ—  โƒ—(๐‘ ) = ๐‘’โˆ’๐‘–๐‘™๐‘ ๐›ฟ๐น๐ผโƒ—โƒ—  โƒ—(0) = ๐‘’โˆ’๐‘–๐‘™๐‘  [๐น๐ฝ๐‘ƒโƒ—โƒ—โƒ—โƒ—  โƒ— โˆ’ ๐น๐ฝโƒ—โƒ—  โƒ—]๐‘ก=0= ๐‘’โˆ’๐‘–๐‘™๐‘ ๐›ฟ๐น๐ฝ๐’ฌโƒ—โƒ—โƒ—โƒ—  โƒ—(0) . (123)  ๐›ฟ๐น๐ผ๐“ โƒ—โƒ—โƒ—โƒ—  โƒ—(๐‘ ) = ๐‘’โˆ’๐’ฌ๐‘–๐‘™๐‘ ๐›ฟ๐น๐ผ๐’ฌโƒ—โƒ—โƒ—โƒ—  โƒ—(0) = ๐›ฟ๐น๐ผโƒ—โƒ—  โƒ—(๐‘ ) + ๐‘’๐’ซ๐‘–๐‘™๐‘ ๐›ฟ๐น๐ผ๐’ฌโƒ—โƒ—โƒ—โƒ—  โƒ—(0) . (124) 67  Basically, while ๐›ฟ๐น๐ผโƒ—โƒ—  โƒ— and ๐›ฟ๐น๐ฝ๐’ฌโƒ—โƒ—โƒ—โƒ—  โƒ— are identical at ๐‘ก = 0, their future values propagate according to different operator, so ๐‘ฒ๐ผ๐ฝ(๐‘ก) and ๐‘ช๐ผ๐ฝ(๐‘ก) are not equal to each other. For short memory, however, ๐’ฌ โ‰ˆ 1, so the memory kernel can be approximated as ๐›ฝโŒฉ๐›ฟ๐น๐ผโƒ—โƒ—  โƒ—(๐‘ก โˆ’ ๐‘ ) โŠ— ๐›ฟ๐น๐ฝโƒ—โƒ—  โƒ—(0)โŒช.    68  Appendix B  In this appendix, we want to prove the equivalence between the explicit and the auxiliary memory formulation. Specifically, under the condition ๐‘ฉ๐‘ ๐‘ ๐‘ฉ๐‘ ๐‘ ๐‘‡ = ๐‘˜๐ต๐‘‡(๐‘จ๐‘ ๐‘  + ๐‘จ๐‘ ๐‘ ๐‘‡ ), the time-correlation โŒฉ๐œ๐ผ๐ฝโƒ—โƒ— โƒ—โƒ—  (๐‘ก) โŠ— ๐œ๐พ๐ฟโƒ—โƒ—โƒ—โƒ—โƒ—โƒ— (0)โŒช also satisfies an FDT relation similar to Eq. (26), meaning with ๐œ๐ผ๐ฝโƒ—โƒ— โƒ—โƒ—  (๐‘ก) and ๐‘ฒ(๐‘ก) defined in Eq. (33). By definition, we have By remembering that the average is just another integral and the matrices ๐‘จ, ๐‘ฉ are all constant, we can reorder the above equation into which is using the property of the random vector. Now let us consider the matrix ๐’(๐‘ก) = โˆซ ๐‘‘๐‘กโ€ฒ๐‘’๐‘กโ€ฒ๐‘จ๐‘ ๐‘ ๐‘’๐‘กโ€ฒ๐‘จ๐‘ ๐‘ ๐‘‡๐‘กโˆ’โˆž: the derivative of ๐’ is in which ๐‘ฐ is the identity matrix. Since ๐’(0) is constant, we have Thus, under the condition of Eq. (34), we find    โŒฉ๐œ๐ผ๐ฝโƒ—โƒ— โƒ—โƒ—  (๐‘ก) โŠ— ๐œ๐พ๐ฟโƒ—โƒ—โƒ—โƒ—โƒ—โƒ— (0)โŒช = ๐‘˜๐ต๐‘‡๐›ฟ๐ผ๐ฝ,๐พ๐ฟ๐‘ฒ(๐‘ก) , (125)  โŒฉ๐œ๐ผ๐ฝโƒ—โƒ— โƒ—โƒ—  (๐‘ก) โŠ— ๐œ๐พ๐ฟโƒ—โƒ—โƒ—โƒ—โƒ—โƒ— (0)โŒช = โŒฉโˆซ ๐‘จ๐‘๐‘ ๐‘’โˆ’(๐‘กโˆ’๐‘กโ€ฒ)๐‘จ๐‘ ๐‘ ๐‘ฉ๐‘ ๐‘ ๐œ‰๐ผ๐ฝโƒ—โƒ— โƒ—โƒ—  (๐‘กโ€ฒ)๐‘‘๐‘กโ€ฒ๐‘กโˆ’โˆžโŠ— โˆซ ๐‘จ๐‘๐‘ ๐‘’tโ€ฒโ€ฒ๐‘จ๐‘ ๐‘ ๐‘ฉ๐‘ ๐‘ ๐œ‰๐พ๐ฟโƒ—โƒ— โƒ—โƒ—โƒ—โƒ—  (๐‘กโ€ฒโ€ฒ)๐‘‘๐‘กโ€ฒโ€ฒ0โˆ’โˆžโŒช . (126)  โˆซ ๐‘‘๐‘กโ€ฒโ€ฒ โˆซ ๐‘‘๐‘กโ€ฒ๐‘จ๐‘๐‘ ๐‘’โˆ’(๐‘กโˆ’๐‘กโ€ฒ)๐‘จ๐‘ ๐‘ ๐‘ฉ๐‘ ๐‘ โŒฉ๐œ‰๐ผ๐ฝโƒ—โƒ— โƒ—โƒ—  (๐‘กโ€ฒ) โŠ— ๐œ‰๐พ๐ฟโƒ—โƒ— โƒ—โƒ—โƒ—โƒ—  (๐‘กโ€ฒโ€ฒ)โŒช๐‘ฉ๐‘ ๐‘ ๐‘‡ ๐‘’๐‘กโ€ฒโ€ฒ๐‘จ๐‘ ๐‘ ๐‘‡๐‘จ๐‘๐‘ ๐‘‡๐‘กโˆ’โˆž0โˆ’โˆž , (127)      โˆซ ๐‘‘๐‘กโ€ฒโ€ฒ๐‘จ๐‘๐‘ ๐‘’โˆ’๐‘ก๐‘จ๐‘ ๐‘  [โˆซ ๐‘‘๐‘กโ€ฒ๐‘’๐‘กโ€ฒ๐‘จ๐‘ ๐‘ ๐›ฟ๐ผ๐ฝ,๐พ๐ฟ๐›ฟ(๐‘กโ€ฒ โˆ’ ๐‘กโ€ฒโ€ฒ)๐‘กโˆ’โˆž] ๐‘ฉ๐‘ ๐‘ ๐‘ฉ๐‘ ๐‘ ๐‘‡ ๐‘’๐‘กโ€ฒโ€ฒ๐‘จ๐‘ ๐‘ ๐‘‡๐‘จ๐‘๐‘ ๐‘‡0โˆ’โˆž= ๐›ฟ๐ผ๐ฝ,๐พ๐ฟ๐‘จ๐‘๐‘ ๐‘’โˆ’๐‘ก๐‘จ๐‘ ๐‘  [โˆซ ๐‘‘๐‘กโ€ฒโ€ฒ๐‘’๐‘กโ€ฒโ€ฒ๐‘จ๐‘ ๐‘ ๐‘ฉ๐‘ ๐‘ ๐‘ฉ๐‘ ๐‘ ๐‘‡ ๐‘’๐‘กโ€ฒโ€ฒ๐‘จ๐‘ ๐‘ ๐‘‡0โˆ’โˆž] ๐‘จ๐‘๐‘ ๐‘‡  , (128)  ๐‘‘๐‘‘๐‘ก๐’(๐‘ก) = ๐‘จ๐‘ ๐‘ ๐’(๐‘ก) + ๐’(๐‘ก)๐‘จ๐‘ ๐‘ ๐‘‡ + ๐‘ฐ , (129)  โˆ’๐‘ฐ = ๐‘จ๐‘ ๐‘ ๐’(0) + ๐’(0)๐‘จ๐‘ ๐‘ ๐‘‡ = โˆซ ๐‘‘๐‘กโ€ฒ๐‘’๐‘กโ€ฒ๐‘จ๐‘ ๐‘ ๐‘จ๐‘ ๐‘ ๐‘’๐‘กโ€ฒ๐‘จ๐‘ ๐‘ ๐‘‡0โˆ’โˆž+ โˆซ ๐‘‘๐‘กโ€ฒ๐‘’๐‘กโ€ฒ๐‘จ๐‘ ๐‘ ๐‘จ๐‘ ๐‘ ๐‘‡ ๐‘’๐‘กโ€ฒ๐‘จ๐‘ ๐‘ ๐‘‡0โˆ’โˆž= โˆซ ๐‘‘๐‘กโ€ฒ๐‘’๐‘กโ€ฒ๐‘จ๐‘ ๐‘ (๐‘จ๐‘ ๐‘  + ๐‘จ๐‘ ๐‘ ๐‘‡ )๐‘’๐‘กโ€ฒ๐‘จ๐‘ ๐‘ ๐‘‡0โˆ’โˆž . (130)  โŒฉ๐œ๐ผ๐ฝโƒ—โƒ— โƒ—โƒ—  (๐‘ก) โŠ— ๐œ๐พ๐ฟโƒ—โƒ—โƒ—โƒ—โƒ—โƒ— (0)โŒช = โˆ’๐‘˜๐ต๐‘‡๐›ฟ๐ผ๐ฝ,๐พ๐ฟ๐‘จ๐‘๐‘ ๐‘’โˆ’๐‘ก๐‘จ๐‘ ๐‘ ๐‘จ๐‘๐‘ ๐‘‡ = ๐‘˜๐ต๐‘‡๐›ฟ๐ผ๐ฝ,๐พ๐ฟ๐‘ฒ(๐‘ก) . (131) 

Cite

Citation Scheme:

        

Citations by CSL (citeproc-js)

Usage Statistics

Share

Embed

Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                        
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            src="{[{embed.src}]}"
                            data-item="{[{embed.item}]}"
                            data-collection="{[{embed.collection}]}"
                            data-metadata="{[{embed.showMetadata}]}"
                            data-width="{[{embed.width}]}"
                            async >
                            </script>
                            </div>
                        
                    
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:
https://iiif.library.ubc.ca/presentation/dsp.24.1-0365974/manifest

Comment

Related Items