{'fraction_non_alphanumeric': 0.02558281377066956, 'fraction_numerical': 0.007251287611818921, 'mean_word_length': 4.759758001561281, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Effective thermostat induced by coarse-graining of SPC water

We investigate how the transport properties of a united atoms fluid with a dissipative particle dynamics thermostat depend on the functional form and magnitude of both the conservative and the stochastic interactions. We demonstrate how the thermostat strongly affects the hydrodynamics, especially diffusion, viscosity, and local escape times. As model system we use SPC water, from which projected trajectories are used to determine the effective interactions in the united atoms model. The simulation results support our argument that the thermostat should be viewed as an integral part of the coarse-grained dynamics, rather than a tool for approaching thermal equilibrium. As our main result we show that the united atoms model with the adjusted effective interactions approximately reproduce the diffusion constant and the viscosity of the underlying detailed SPC water model.


I. INTRODUCTION
In molecular dynamics, the addition of a thermostat is usually motivated as a representation of interactions with the surroundings, important primarily for the thermal equilibration of the system. The choice of thermostat generally depends on whether the system is characterized by constant temperature and volume (e.g. the Nosé-Hoover 1,2,3 thermostat) or constant temperature and pressure (e.g. the Andersen 4 or Parrinello-Rahman 5 thermostats). Under the standard assumption that the system is sufficiently chaotic (and therefore mixing), averages over the thermodynamic equilibrium ensemble can be calculated as time averages, provided that the system is allowed to equilibrate before measuring the time average. This allows for measuring, e.g., the pressure and heat capacity of the system, and to investigate complex phenomena such as phase transitions. Due to the fluctuation dissipation theorem, thermodynamic properties defined through the partition function are unaffected by the choice of thermostat (for further discussion, see section II).
In contrast to the equilibrium properties, transport processes are intimately tied to the trajectories as they depend on the auto-correlation of velocities and forces through the Green-Kubo relations 6,7 . The thermostat changes how the trajectories approach equilibrium, and as a consequence the transport properties change as well. This is usually considered to be a problem; especially when running non-equilibrium (NEMD) simulations to measure transport properties, where the system is brought to a stationary, out of equilibrium, state by an external force field. The standard way of minimizing this effect is to make the system as large as possible, and the interaction with the thermostat as weak as possible.
In this article we take a different view on the role of the thermostat: Rather than minimizing its effects, we view the thermostat as an integrated part of the dynamics. The thermostat is a mesoscopic representation of the exchange of energy between the coarse-grained degrees of freedom and the microscopic, rapidly fluctuating, degrees of freedom. As a result we expect the effective coarse-grained dynamics to be of Langevin type. To demonstrate this point explicitly, we study a coarse-grained model of SPC water. We use MD simulations of SPC water as our microscopic system and define the projected dynamics as the center of mass motion (united atoms) of the SPC particles. We demonstrate that it is possible to choose the effective conservative and stochastic interactions for the united atoms model, such that the equilibrium and transport properties are close to those of the projected microscopic dynamics.
The motivation for viewing the thermostat as a natural part of the system's dynamics lies in the Mori-Zwanzig theory on projection operators 8,9,10,11,12 . In short, the theory states that given a microscale dynamics, a lower dimensional representation can be attained through a projection of the phase space (e.g. the map from the atomic coordinates to the center of mass of the molecules), where fast degrees of freedom can either give rise to Markovian (white) noise and dissipation, or be eliminated due to averaging 11 . Which of these two scenarios that best describe the system at hand depends on the exchange of energy (heat) between the coarse-grained degrees of freedom and the degrees lost in the coarse-graining procedure. If the energy exchange is substantial, then the motion of the coarse-grained particles will not be smooth or deterministic, but rather evolve according to a stochastic (Langevin type) differential equation. We will show that our test system is described by the latter case. Naturally, how faithfully a particular coarse-grained model represents the underlying dynamics depends strongly on the choice of projection, the mixing time of the fast degrees of freedom, and the details of the coupling between fast and slow dynamics.
In principle, given a specific projection it is possible to derive the effective coarse-grained dynamics using the Mori-Zwanzig theory. In practice, however, the direct approach is neither computationally feasible, nor is the resulting coarsegrained system typically represented as a particle based simulation. To make the method practically useful we use the united atoms coarse-graining as an ansatz. It is important that the forces in the coarse-grained model respects the fundamental mechanical symmetries of the system, especially Galilean invariance, which implies conservation of linear and angular momentum. Therefore, we restrict ourselves to central forces that obey Newton's third law, and where the force between two particles depends only on the distance between the particles. If we additionally assume that the stochastic component of the pair-wise forces is statistically independent, we arrive at the dissipative particle dynamics (DPD) model, introduced in 1992 by Hoogerbrugge and Koelman 13 as a simulation technique for complex hydrodynamic phenomena. In the DPD method, the thermostat is represented explicitly as dissipative and stochastic forces 13,14 . DPD can now be considered a standard method for mesoscopic simulation. It has been used to study, for example, complex fluids 15 , spontaneous self-assembly of amphiphilic molecules into bilayered membranes 16 , vesicles 17,18 , and hydrodynamics 19 .
We elaborate briefly on the direct consequences of choosing a Galilean invariant thermostat that faithfully characterize the local transport of linear and angular momentum. First, sufficiently close to equilibrium, the system obeys the classical result of the asymptotic t −d/2 decay of the velocity autocorrelation (d is the dimensionality of the system) 20 . In addition, the interactions give rise to hydrodynamic modes in the fluid 20,21 , which lead to the Navier-Stokes equations on the macroscopic level 14 . Galilean invariance is especially important for instance in NEMD measurements of the viscosity of complex fluids 15,22 . If the equations of motion are not Galilean invariant the thermostat causes a thermodynamic screening (see, e.g. Ref. 22 and references therein). A striking example of the importance of correct hydrodynamic interactions is the block co-polymer melts studied by Groot and Madden 23,24 , which pass through a disordered meta-stable phase before organizing themselves in a hexagonal pattern. Using DPD allowed the system to go between the two regions, whereas a similar model without hydrodynamic interactions failed to pass through the barrier. Hence, the local interactions between the coarse-grained molecules are important for the observed macroscopic state. Furthermore, the folding pathways of many proteins and peptides are sensitive to the interactions with the surrounding water, and exhibit meta-stable intermediate states separated by kinetic barriers 25,26 .
In the standard approach to DPD the forces are usually not derived from the microscopic dynamics (for exceptions, see e.g. Ref. 27 and references therein). Instead, generic and simple functional forms, such as constant or linearly decreasing up to a cut-off radius, are used. Moreover, the same functional form is usually used for both conservative and stochastic forces. It should be emphasized that this choice is guided by maximizing simplicity rather than strict physical arguments (the generic repulsive nature of the conservative force is, however, consistent with the effective interactions between the center of mass of clusters of particles 28 ). To account for transport properties it has become common practice to rescale time in order to match either diffusion or viscosity. A serious and well documented pitfall of the DPD method in this version is that it cannot be tuned to match both the diffusion rate and the viscosity simultaneously 14,29,30,31 . In an attempt to overcome this problem, Junghans et al. 32 recently presented a study where a transversal component, orthogonal to the central force, was added to the DPD interactions. It was shown that this gives freedom to tune both the viscosity and the diffusion rate independently. However, it is clear that with this ansatz the conservation of angular momentum is no longer manifest, and we should therefore expect that the hydrodynamic behavior of the system is not necessarily faith-fully represented. In the current study we limit the interactions in the DPD dynamics to central forces but change the magnitude and support of the interactions. We demonstrate that we can find parameter settings where both diffusion and viscosity are consistent between the coarse-grained and the microscopic model.
The remainder of the article is organized as follows. First, we introduce DPD as a thermostat and show how the radial dependence of the stochastic force influences the approach to equilibrium. Second, we review how to estimate the conservative force from the center of mass motion of SPC water molecules in atomistic MD simulations. Third, in the results section, we compare transport properties of the coarse-grained model to those of SPC water. Finally, we conclude with a discussion and summarize the results.

II. THE DPD THERMOSTAT
In this section we review the DPD thermostat, and illustrate how the structure of the stochastic force gives rise to Galilean invariance of the dissipative force. We also show how the structure of the noise affects the system's relaxation towards equilibrium.
In its simplest form, the equations of motion for a DPD model, with particles positioned at r i , with velocities v i and momenta p i , can be written as a system of Langevin equationṡ where F C i j , F D i j and F S i j are the conservative, dissipative, and stochastic forces between particles i and j. Both the conservative and non-conservative interactions in DPD are modeled by central forces obeying Newton's third law, ensuring that linear and angular momentum are conserved 13 . In DPD, the stochastic force between particles i and j take the form where r i j is the distance between particles i and j, e i j is the unit vector pointing from j to i, k B is Boltzmann's constant, and T is the temperature in Kelvin. The scalar function ω(r i j ) describe how the stochastic force depends on the distance between the particles, and ζ i j is interpreted as a symmetric Gaussian white noise term with mean zero and covariance where δ i j and δ (t) are the Kronecker and Dirac delta functions, respectively. At thermal equilibrium the system is distributed according to the canonical ensemble where Z is the normalization term for the distribution. The equilibrium ensemble must be invariant under the equations of motion. Since H is a constant of motion for Hamiltonian dynamics, the ensemble is invariant under this part of the dynamics (this is true for any ensemble where the probability of finding the system in a given micro-state depends only on the energy), but in order for the dissipative and random forces to conserve the equilibrium, we need to choose the dissipative force F D i such that the combined contributions from the dissipative and stochastic forces cancel when acting on the equilibrium distribution: where L is the Fokker-Planck operator of Eq. (1), and A i j is a 3 × 3 matrix given by the covariance of the total forces on particles i and j. The equilibrium Fokker-Planck equation (5) is commonly referred to as a fluctuation-dissipation relation. For the DPD model, this was first analyzed by Español and Warren 14 . Since it must hold for all points in phase-space, the only possible solution for the dissipative force is For the DPD model the force covariance is given by where ⊗ denotes an outer product. Inserting Eq. (7) into Eq. (6) we can write the dissipative force on a particle as a sum of pair-wise dissipative forces: Note that since F D i depends only on the velocity differences between interacting particles, it is manifestly Galilean invariant, and it is clear from the derivation above that this is a direct consequence of the covariance property of the stochastic forces [c.f. Eq. (3)], which in turn stems from the assumption that Newton's third law applies.
A subtle point is that there is no one-to-one relation between the stochastic forces and A, the total force covariance matrix: Though F S i j appears in the Langevin equation, the dynamics depends only on A (as is seen from the Fokker-Plank equation for the system), and for each choice of the force covariance matrix there are infinitely many choices of F S i j that obey Eq. (7). For instance, the standard formulation of the DPD equations of motion contains N(N − 1)/2 independent random variables, where N is the number of particles, but it is possible to find a representation using no more than 3N independent random variables. In order to use this in the simulations, however, we would need to calculate the square root of the matrix A in each time step. Although both representations are equally valid, the standard form of the DPD equations of motion is much more efficient.
The dissipation-fluctuation relation [c.f. Eq. (5)] asserts that the systems thermal equilibrium is a fixed-point for the Fokker-Planck equation, or equivalently an invariant measure of Eq. (1). In order to better understand the effect of the stochastic forces on the path to thermal equilibrium it is illuminating to study the time-evolution of the entropy of a nonequilibrium ensemble f (x, p) over the phase-space. Following Green 33 , we can express the difference in entropy of the ensemble f to the equilibrium distribution f eq as which is negative for all ensemble distributions with the same support as the equilibrium distribution, and is zero if and only if the two distributions are equal. With strictly Hamiltonian dynamics, the entropy is constant in time. Intuitively, this is because the internal energy of the system needs to change in order for the ensemble to approach the equilibrium distribution, but the Hamiltonian conserves the energy. With the addition of dissipative and random forces it can be shown that entropy difference is a negative Lyapunov function on the space of ensemble distributions f : which clearly is positive. In agreement with the second law of thermodynamics, the entropy continues to increase until f = f eq . The main point of Eq. (10) is that it shows explicitly how the relaxation towards the equilibrium distribution depends on the structure of the noise. As we have seen above, the conservative and dissipative forces follow from the Langevin equation and the shape of the equilibrium distribution. Hence, for the purpose of capturing the transport properties of the system, only the shape of the stochastic forces remains to be specified.
Eqs. (1-3, 8) together establishes the general form of the DPD dynamics. Both the conservative force F C i j , or equivalently the corresponding scalar potential, and the scalar function ω(r) depend on the particular system of interest and need to be determined to obtain the correct DPD model. In practice, this is the difficult part of DPD, and also the rationale behind the often used heuristic approach for deciding the interactions. In the following section we discuss a more systematic approach for deriving the conservative interaction from a microscopic dynamics.

III. CONSERVATIVE POTENTIAL RECONSTRUCTION
Several methods exist for deriving potentials from given equilibrium radial distribution functions (RDFs) 27,34,35,36,37,38 (alternatively, direct time averaging over the fast degrees of freedom can be used, see e.g. Refs. 28,39). The methods for reconstructing effective potentials from the RDF rely on the result by Henderson 40 , that two pairwise potentials resulting in the same RDF cannot differ by more than an additive constant. The importance of this theorem lies in the one-to-one correspondence between pairwise central force and the radial distribution function.
In order to determine the effective potential corresponding to a given RDF, we use the inverse Monte Carlo method developed by Lyubartsev and Laaksonen 35 . This method starts from a discretized Hamiltonian of the system, which corresponds to using a stepwise constant potential, Φ α . S α denotes the number of particle pairs separated by a distance in the range r α to r α+1 , where r 0 = 0, r 1 = dr and r α = α · dr.
The average of S α is directly connected to the radial distribution function, g(r α ), by the relation where N is the number of particles, L the side length of the simulation box and V α the volume of the spherical shell between radii r α and r α+1 , Using a start potential, Φ α , normally chosen as the potential of mean force, Φ (0) α = −k B T ln g (r α ), a Monte Carlo (MC) simulation is made, where at each time step a random particle is moved a small distance in a random direction. From the simulation, the quantity S α is measured, and the difference from the desired value, S * α (given directly form the wanted RDF), is calculated, Also measured is the correlation in the number of particle pairs at different distances in the system, S α S γ , and the covariance matrix, S α S γ − S α S γ . By taking these averages in the canonical ensemble, each entry in the covariance matrix can be related to the partial derivative of S α with respect to the potential Φ γ , through As Lyubartsev and Laaksonen 35 suggests, this can be used in an expansion of S α , from which an estimate of the error in the current potential, i.e. in Φ (0) α , is given by ∆Φ α . Updating the potential and restarting the MC simulation with the new potential will result in a better approximation of the RDF, and after some iterations the method converges on a potential giving rise to the desired RDF. Difficulties with convergence in the method can normally be overcome by not changing the potential as much as described by ∆Φ, but instead chose the new potential as 1]. For further discussion on the efficiency of the inverse Monte Carlo method, see e.g. Ref. 41.

IV. SIMULATIONS
The procedure and technical details of our simulations are outlined below, corresponding to the schematics in Fig. 1. From a detailed MD simulation of SPC water, coarse grained conservative forces are constructed using the inverse Monte Carlo method. By also adding a DPD thermostat we wish to demonstrate how the thermostat can be tuned to obtain good agreement between the transport properties of the coarsegrained and the microscopic levels of description. It is important to notice that the added thermostat does not affect the RDF, and therefore the conservative forces can be determined independently of the stochastic interactions. This is the rationale behind decomposing the coarse-graining procedure into two steps, as shown in Fig. 1.

A. SPC water
Simulations using the Simple Point Charge (SPC) model of water are performed using the molecular dynamics package Gromacs. The simulation is set up with 2180 water molecules, target temperature 298 K and target pressure 1 bar (NPT ensemble with Berendsen thermostat). Electrostatic interactions are modeled using the reaction field approach for separation distances larger than 1.4 nm.

B. United atoms
As a coarse-grained description of the same system, we use spherically symmetric point particles representing the center of mass of the SPC molecules, and assume pairwise interactions between the coarse-grained particles. This type of coarse-graining is usually referred to as United Atoms (UA). The interaction potential between the particles is obtained by applying the inverse MC method (section III) to RDF data for the center of mass of the SPC molecules. Figure 2 shows the coarse-grained potential as the piecewise constant approximation given by the inverse MC method. In the coarse-grained simulations, the particle density and temperature are the same as in the original SPC simulation, but we now use the NVT ensemble (Berendsen thermostat), as this makes comparison with simulations using a DPD thermostat easier. Simulations of SPC water with the NVT ensemble and box size 4.06 nm (which was the average box size using the NPT ensemble) do not show any differences compared to the NPT simulations, so we conclude that the choice of ensemble is not important for the analysis.

C. United atoms with DPD-thermostat (UA-DPD)
The UA-DPD simulations are performed using a standard velocity Verlet integration of the DPD equations of motion. This method gives second-order accuracy in the time-step for the conservative part, but is only first-order for the stochastic part 15 . In the case of unbounded conservative interactions for small distances, the accuracy of the integration is mainly determined by the conservative part 22 ; therefore, the integration error in the stochastic forces is of less concern. The simulation set-up is identical to the United Atoms simulation, with the addition of dissipative and stochastic forces. We examine the simple example where ω(r) is linearly decreasing: The parameter σ defines the strength of the dissipative force, and r c is the cutoff radius. Note that ω(r) is also set to zero for r ≤ r 0 . The reason for this is that for the rare occasions when the inter-particle distance goes below r 0 = 0.277 nm, we do not want a stochastic force capable of pushing the particles closer together. This would result in uncontrolled magnitude of conservative forces, unless the time step is significantly lowered.
In the simulations we vary the parameters σ and r c to examine the effect on the diffusion and viscosity values. Tuning the parameters, we try if it is possible to obtain diffusion and viscosity values consistent with those of the SPC simulation.

D. Measurements of diffusion and viscosity
The diffusion coefficients for SPC water and the UA model without DPD thermostat are obtained directly from Gromacs. In the UA-DPD simulations, the diffusion coefficients are calculated from the mean square displacement, using the Einstein relation The viscosities for SPC and UA are measured with the NEMD method implemented in Gromacs 42 , while the viscosity of UA-DPD is measured applying a Poiseuille flow method 43 .

A. Self-diffusion and viscosity
For the UA-DPD, we investigated the dependency of the diffusion and viscosity on the cutoff radius r c . By maintaining a constant σ for different values of r c , we obtained the results shown in Figure 3. Both transport properties are strongly dependent on the value of r c : While the diffusion decreases with increasing r c , the viscosity increases. The diffusion of SPC water molecules was measured to D = 4.11 × 10 −9 m 2 /s. Similarly, the viscosity was measured to η = 0.42 cP. As these values represents the correct dynamics of the SPC water molecules and hence the coarse-grained representation, a coarse-grained model should be able to reproduce these numbers. The simplest approach, the UA model with a simple Berendsen thermostat, resulted in: D = 16.5 × 10 −9 m 2 /s and η = 0.116 cP. Clearly, even though this model results in the same equilibrium distribution, i.e. the same RDF, as the SPC model, it does not represent the transport properties correctly.
Adding the DPD thermostat, we tuned σ separately for each r c , such that the diffusion became approximately equal to the SPC value. The resulting σ values and the corresponding r c values are listed in Table I. Keeping D constant, we examined the variation of the corresponding viscosities as a function of r c and σ . The results are summarized in Table II. The viscosity is at a minimum for small values of r c and then increases moderately with increasing r c . The minimum value of the UA-DPD viscosity is close to the viscosity of the coarse-grained SPC system: 0.49 cP compared to 0.42 cP. In particular when comparing to the UA model, it is clear that the addition of a carefully tuned DPD thermostat plays an important role in preserving transport properties.
We would like to stress that while our choice of weight function ω(r) was an appropriate starting point, it might be too simple to capture the system in the best possible way. Our results show that the transport properties depend on the functional form of ω(r). In the light of this it is expected, rather than surprising, that the viscosity is not matched exactly. In another study 44 , we have suggested how the force covariance can be used to obtain an estimate of ω(r). A more elaborate theoretical framework based on this idea is currently in progress.   [45] 0.8909 [46]

B. Escape time distribution
In addition to determining the diffusion coefficient and the viscosity, we have also measured the escape probabilities for pairs of particles, that is, how large the probability is for two particles to separate beyond a given distance in a given time interval, given the initial separation. This gives a more detailed view of the local dynamics than the other transport properties. The escape probabilities for SPC water, united atoms with the Berendsen thermostat, and the UA-DPD using ω(r) from Eq. (17) are shown as level plots in Fig. 4. The different levels in these figures represent the probability that two particles, originally separated by the distance r, remain within a separation distance of 0.5 nm at a later time, t, displayed on the y-axis.
Generally, the escape time is decreasing with increasing initial separation r. Note especially that the average escape time is much larger in the region r 0.3 nm than for r > 0.3 nm. This is consistent with the well-known caging effect. Particles repeatedly bounce against their closest neighbors in the fluid, and for a pair of nearby particles to separate beyond  the typical neighbor distance (about 0.28 nm for SPC water) requires either a concerted motion of several particles, or a very high kinetic energy of the particles. Note that this barrier is largely missing from the united atoms system with the Berendsen thermostat; this is the main reason why the diffusion rate is too high compared to the SPC water. In contrast, the UA-DPD thermostat fitted to the diffusion rate and viscosity (Fig. 4C) exhibit escape times similar to those of the SPC model. Hence, matching diffusion and viscosity seems to confer agreement with respect to more general transport properties.

VI. SUMMARY AND DISCUSSION
We have applied a coarse-graining scheme to the SPC water model. The effective dynamics of the center of mass of the water molecules is Langevin-like with a conservative (drift) component. The conservative part of the dynamics is represented by a pairwise potential that can be determined from the radial distribution function. Stochastic and dissipative forces, motivated by the Mori-Zwanzig projection operator formal-ism, are also included as a DPD thermostat. We have demonstrated that the radial dependence of the stochastic and dissipative forces strongly affect the transport properties. The main conclusion is that the diffusion rate and the viscosity of the original system can be approximately matched by the coarsegrained model. This was obtained by tuning the magnitude and functional form of the random and dissipative forces.
We argue that in meso-and microscale simulations of coarse-grained Hamiltonian systems, the stochastic and dissipative components of the interactions should normally be considered as an integrated part of the effective dynamics. For coarse-grained models with only pairwise interactions, a dissipative particle dynamics ansatz is well suited for representing the non-conservative part of the dynamics, since it is the most general pairwise interaction that respects the local conservation or transport of linear and angular momentum, i.e. hydrodynamic behavior. The DPD interactions depend on the inter-particle distance. We argue that this functional dependence should the tuned so that the transport properties are consistent with the microscopic dynamics. More generally we show explicitly in Eq. 10 how the relaxation towards equilibrium depends on the structure of the thermostat. Finally, we point out that all practical molecular dynamics methods, not only united atoms models, are ultimately coarse-grained from some more fundamental level, either from quantum mechanics or by removing fast degrees of freedom (e.g. vibration modes in covalent bonds) from a classical mechanics model. The interaction with these fast degrees of freedom plays an active role when the system equilibrate. The approach used in this paper can therefore possibly also be applied to atomistic molecular dynamics simulations.