Transport properties controlled by a thermostat: An extended dissipative particle dynamics thermostat

Abstract

We introduce a variation of the dissipative particle dynamics (DPD) thermostat that allows for controlling transport properties of molecular fluids. The standard DPD thermostat acts only on a relative velocity along the interatomic axis. Our extension includes the damping of the perpendicular components of the relative velocity, whilst keeping the advantages of conserving Galilei invariance and within our error bar also hydrodynamics. This leads to a second friction parameter for tuning the transport properties of the system. Numerical simulations of a simple Lennard-Jones fluid and liquid water demonstrate a very sensitive behaviour of the transport properties, e.g., viscosity, on the strength of the new friction parameter. We envisage that the new thermostat will be very useful for the coarse-grained and adaptive resolution simulations of soft matter, where the diffusion constants and viscosities of the coarse-grained models are typically too high/low, respectively, compared to all-atom simulations.

---

I. Introduction

Using the optimal set of degrees of freedom (DOFs) in computer simulations of soft matter guarantees efficiency, accuracy and avoids huge amounts of unnecessary detail, which might even obscure the underlying physics. This very idea is exploited in systematic coarse-graining efforts for modern computational materials science and biophysics problems,^1–3 where atomistic simulations are usually beyond the possibilities of current and near future computers. A similar philosophy of reducing the number of DOFs is also employed by representing clusters of molecules with soft particles when simulating fluids on a mesoscopic scale using dissipative particle dynamics (DPD).^4–8 However, in typical soft matter systems different time- and length-scales are intrinsically interconnected and a multiscale modeling approach is required to tackle such problems in the most efficient way.^9–16

Recently, we have proposed an adaptive resolution scheme (AdResS) that couples the atomistic and coarse-grained levels of detail.^17–22 Due to the reduction in DOFs upon coarse-graining, which eliminates the fluctuating forces associated with those missing molecular DOFs, the coarse-grained molecules typically move faster than the corresponding atomistically resolved ones.^23,24 Although the accelerated dynamics is advantageous in some cases it can turn out to be problematic if one is really interested in dynamics in situations where two levels of resolutions are used within one simulation, as in the case of AdResS. To overcome this problem one can couple different classes of DOFs to the Langevin thermostat with different friction constants.^25–28 We have shown in the example of liquid water that the coarse-grained dynamics can be slowed down by increasing the effective friction in the coarse-grained system.²⁹ However, it is well known that the Langevin thermostat does not reproduce the correct hydrodynamics, i.e., the hydrodynamic interactions are unphysically screened. In order to correctly describe hydrodynamic interactions, one has to resort to the DPD thermostat.³⁰

In the past years DPD has established itself as a useful thermostat for soft matter simulations. The DPD thermostat is known to have several good properties, i.e., the thermostat satisfies Newton's third law by construction and owing to mass, momentum and temperature conservation, hydrodynamics is also correctly reproduced.³¹ As it turns out, however, the DPD thermostat in its standard form is not capable of controlling liquid properties such as viscosity and diffusion constant.³⁰ The aim of this work is to extend the DPD equations in such a way that these quantities can be tuned by changing the parameters of the thermostat. We consider the most general version of DPD⁶ and exploit the terms, which are not used in the standard DPD approach, i.e., the damping of the perpendicular components of the relative particle velocities. This allows us to tune the viscosity of the coarse-grained liquid to match that of the all-atom counterpart while conserving the virtues of the standard DPD thermostat. This should be very useful for both coarse-grained and adaptive resolution simulations of soft matter.

The article is organized as follows: In Section II the standard and new DPD thermostats are presented. The simulation setup is described in Section III. The results of molecular dynamics (MD) simulations of the Lennard-Jones fluid and liquid water are reported in Section IV, followed by the conclusions in Section V.

II. Thermostats

A. Standard dissipative particle dynamics thermostat

Newton's equations of motion are used in microcanonical NVE MD simulations and generate dynamics with constant energy. In order to run simulations in the canonical NVT ensemble the equations of motion have to be modified. In Langevin dynamics two additional forces are introduced, a damping and a random force, whose ratio defines the temperature. The DPD equations of motion (used by the DPD thermostat) are then given by^4,6andwhere denotes the conservative force on the ith particle. The damping and random forces, can be split up in particle pair forces aswhere the dissipative force reads as⁶and the random force is given by

In these equations the relative velocity between the ith and jth particle was introduced, while denotes the unit vector of the interatomic axis . ζ^∥ is the friction constant and σ^∥ the noise strength. w^D(r_ij) and w^R(r_ij) are r-dependent weight functions. These are connected by the fluctuation–dissipation theorem (see eqn (9)). The variable Θ_ij is symmetric in the particle indices (Θ_ij = Θ_ji) and has the following first

〈Θ_ij(t)〉 = 0 (7)

and second moment

〈Θ_ij(t)Θ_kl(t′) = 2(δ_ikδ_jl + δ_ilδ_jk)δ(t−t′). (8)

The fluctuation–dissipation theorem⁶ reads as

(σ^∥)² = k_BTζ^∥ (9)

and

(w^R(r))² = w^D(r). (10)

The above DPD equations conserve the total momentum and reproduce correctly the hydrodynamics interactions in the system. However, previous studies^30,32 have shown that the strength of the friction ζ^∥ does not influence the viscosity in linear order. In order to be able to tune the value of the viscosity of the system while preserving all of the virtues of the standard DPD thermostat presented briefly above we introduce in the next subsection its extension named “Transverse DPD thermostat”.

B. Transverse dissipative particle dynamics thermostat

We generalize eqn (5) and eqn (6) asandwhere ζ and σ are the friction constant and the noise strength of the generalized thermostat, respectively (see the text below). is a projection operatorwhich is symmetric in the particle indices . The scalar noise (see eqn (8)) is replaced by a noise vector which is antisymmetric in the particle indices due to the symmetry of the projection operator and the antisymmetry of the pair force (Newton's third law). The corresponding Fokker–Planck operator is a sum of two parts: the deterministic part and the generalized DPD part

The equilibrium condition then yields the dissipation–fluctuation theorem in the same form as given by eqn (9) and eqn (10).

For the case where we choose the projector along the interatomic axis between particle i and j

we retain the standard DPD thermostat.

Alternatively, one can project on the plane perpendicular to the interatomic axis

The space defined by the projector (eqn (18)) is orthogonal to the case of the standard DPD and introduces an extension of the DPD thermostat, i.e., the Transverse DPD thermostat. Note that owing to this orthogonality the new thermostat can be used in combination with the standard one. This enables us to adjust at the same time two friction constants ζ^∥ and ζ^⊥ for the standard and Transverse DPD thermostats, respectively. Galilei invariance remains valid by construction while hydrodynamics is conserved within our error bar.³³

Our basic assumption is that in contrast to the standard DPD the viscosity is very sensitive to the damping perpendicular to the interatomic axis. This damping mimics the shear of those DOFs that were integrated out in the coarse-graining procedure. In a system with two particles the stochastic forces of the Transverse DPD thermostat act in the same direction as the shear forces. The mean force acting on a particle in the sheared system with more than two particles is hence a sum of two contributions: a force coming from the Transverse DPD thermostat and another one originating from the shearing of the probe. Therefore, the shear viscosity in a simulation with the Transverse DPD thermostat is always higher than with the standard one. In the Green–Kubo picture this additional viscosity arises from the projected velocity–velocity autocorrelation function, which is derived by the Mori–Zwanzig formalism.^35–37 The exact derivation is beyond the scope of this paper and will be presented elsewhere. Here, we shall demonstrate this by the results of our numerical experiments presented in the results section.

III. Simulation setup

All simulations of the Lennard-Jones (LJ) liquid and liquid water are performed using the ESPResSo package.³⁸

A. Lennard-Jones fluid

We use the repulsive Weeks–Chandler–Andersen potentialwith the cutoff at r_c = 2^1/6σ, σ and ε being the standard LJ parameters of length and energy.

We chose as the weight function (see eqn (10)) for both thermostats a step function

The simulations are carried out with a system consisting of 1000, 2000 and 4000 LJ particles at a temperature T = 1.2ε/k_B and density ρ = N_part/V = 1/(1.05σ)³ in a cubic box with periodic boundary conditions.

B. Liquid water

All-atom water NVT simulations at ambient conditions are performed using the rigid TIP3P water model.³⁹ The electrostatics is described by the reaction field (RF) method, in which all molecules outside a spherical cavity of a molecular based cutoff radius R_c = 9Å are treated as a dielectric continuum with a dielectric constant ε_RF = 80.^40–42 Typically, all-atom water simulations are carried out using global thermostats, e.g., Berendsen,⁴³ Nosé–Hoover,^44,45 Nosé–Hoover chains⁴⁶ thermostats, that dissipate the energy uniformly in the system. Local thermostats, e.g., Langevin, DPD,³⁰ Andersen,⁴⁷ Lowe–Andersen,^48–50 Nosé–Hoover–Lowe–Andersen⁵¹ thermostats, that dissipate energy on a spatially localized scale are usually used in coarse-grained simulations. Here, in order to reproduce the hydrodynamics correctly we employed the DPD thermostat³⁰ with the friction constant ζ^∥ = 0.038 ps⁻¹, cutoff radius R_c, and the weight function given by eqn (20). The constant ζ^∥ is much smaller than the intrinsic friction coefficient ξ of the TIP3P water system, i.e., ξ = 288.6 ps⁻¹, so that the stochastic dynamics yields the correct dynamics.^29,52,53 For the coarse-grained water simulations with the Transverse DPD thermostat we have used the single-site water model from Ref. 22, which reproduces essential thermodynamics and structural properties, e.g., the pressure, density, and radial distribution functions, of the all-atom rigid TIP3P water model at standard conditions. Other simulation details are the same as given in Ref. 22.

IV. Results

A. Lennard-Jones fluid

First, we checked in an equilibrium simulation the dependency of pressure and temperature on the strength of the friction. We set the reference temperature to 1.2 ε/k_B and measured the instantaneous temperature defined aswhere E_kin and N_part are the kinetic energy and the number of particles of the system, respectively. The relative deviation between the measured and reference temperature was smaller then 1.2% for all strengths of friction and all combination of thermostats. The mean pressure at that temperature turned out to be 9.8 ± 0.2 ε/σ³, which is in a perfect agreement with the results of previous studies.⁵⁴

Next, we studied the dependency of the liquid transport properties, i.e., the diffusion constant and shear viscosity, on the friction constants ζ^∥ and ζ^⊥ for the standard and Transverse DPD thermostats, respectively.

1. Diffusion constant. The diffusion constant was computed from the particle displacements using the Einstein relation

A small influence on this constant from the standard DPD thermostat³⁰ is expected, but a considerable one from the new Transverse DPD thermostat. Former results³⁰ for the standard DPD thermostat were confirmed. The value of the diffusion constant approaches the equilibrium value (D = 0.08 σ²/τ) for vanishing friction of the Transverse DPD thermostat (see Fig. 1). The diffusion constant D is very sensitive on the friction ζ^⊥ for the Transverse DPD thermostat. By changing ζ^⊥ it is therefore possible to tune the diffusion constant.

Fig. 1 Diffusion constant (4000 LJ particles) as a function of the friction measured in equilibrium for different thermostats. In the case of the combined DPD thermostat only the strength of the friction parameter ζ^⊥ was varied while the friction for the standard one was held constant at the value ζ^∥ = 1.0 ετσ⁻². The errors were obtained by averaging over several runs and Jackknife analysis.⁵⁵

2. Shear viscosity. The viscosity was measured in nonequilibrium molecular dynamics (NEMD) simulation by shearing the system with a constant shear rate in the z-direction³⁰where u_x is a fluid streaming velocity in x direction. The viscosity can then be determined by the following simple formulawhere F is the mean force (momentum transfer per time unit) and L_x and L_y are the box lengths in x and y directions, respectively. The apparent shear viscosity η measured in NEMD simulation approaches the equilibrium viscosity with decreasing shear rate.

We found nearly no dependency on the strength of the friction ζ^∥ for the standard DPD thermostat, (see also Ref. 30). In contrast the friction ζ^⊥ for the Transverse DPD thermostat gives a very sensitive means of controlling the viscosity (see Fig. 2). In the case of the combination of both thermostats the shear viscosity is mostly controlled by the Transverse DPD thermostat. In the limit of a vanishing shear rate a value of 2.45 ± 0.07 ετ/σ³ was extrapolated for the standard DPD thermostat, which matches former results.⁵⁴ This extrapolation has to be done due to short characteristic timescale of the LJ system. Additionally, we also checked that the equilibrium correlation of the pressure tensor is in accordance with the Green–Kubo picture: all (non-auto) off-diagonal – (off-)diagonal elements are uncorrelated, (60 out of 81 elements are zero).

Fig. 2 Shear viscosity measured with the NEMD algorithm for different thermostats (4000 LJ particles) and the different shear rates (0.1 τ⁻¹ and 0.01 τ⁻¹ shown), which then are extrapolated to vanishing shear rate. The errors are obtained by Jackknife analysis.⁵⁵

For higher values of ζ^⊥ the apparent viscosity becomes increasingly dependent on the shear rate. This is due to the fact that the dynamics is controlled in this regime by the thermostat forces (that are linear in ζ^⊥) and hence we end up measuring the “viscosity of the thermostat”.

B. Tuning the dynamics of water

Having shown that the new thermostat enables us to tune the diffusion constant and viscosity of a simple liquid we apply it to an important physical example, i.e., liquid water at ambient conditions.

We first check that the structural properties do not depend on the thermostat and also that we obtain consistency between the coarse-grained and atomistic simulations. The center-of-mass radial distribution functions of the all-atom and coarse-grained system using different values of ζ^⊥ match within the line thickness (see Fig. 3).

Fig. 3 The center-of-mass radial distribution functions of the all-atom (ζ^∥ = 0.038 ps⁻¹, ζ^⊥ = 0) and several coarse-grained simulations (ζ^∥ = 0 and ζ^⊥₁ = 0.5 ps⁻¹, ζ^⊥₂ = 0.75 ps⁻¹ and ζ^⊥₃ = 1.0 ps⁻¹).

There is an intrinsic timescale difference in the diffusive dynamics of the coarse-grained water system because of the reduced number of DOFs, i.e., the self-diffusion constant for the coarse-grained water model is approximately 2 times larger then the corresponding all-atom one using the Langevin thermostat with the same background friction in both cases.²² We used ζ^⊥ = 0.75 ps⁻¹ for the Transverse DPD thermostat (ζ^∥ = 0) to match the diffusion constant of the coarse-grained water model to the corresponding value D = 3.2 × 10⁻⁵ cm² s⁻¹ obtained from the all-atom simulation with the standard DPD thermostat (see Fig. 4). With a friction strength of ζ^⊥ = 0.6 ps⁻¹, which is as desired very close to the above value for matching the diffusion constants, we were also able to match the viscosity to the desired value η = (0.5 ± 0.1) × 10⁻³ Pa s for the TIP3P water model (from our atomistic simulation with ζ^∥ = 0.038 ps⁻¹). In this case the characteristic times for the atomistic model are much longer than the time scale of the shearing (1/[small gamma, Greek, dot above]). Therefore even with a shear rate of [small gamma, Greek, dot above] = 0.01 τ⁻¹, where τ is 1.579 ps, we are in the no shear limit and hence no extrapolation is required. The obtained diffusion constants and viscosities are in good agreement with the published data.⁵⁶

Fig. 4 The mean square displacements over time plot of the all-atom (ζ^∥ = 0.038 ps⁻¹, ζ^⊥ = 0) and several coarse-grained simulations (ζ^∥ = 0 and ζ^⊥₁ = 0.5 ps⁻¹, ζ^⊥₂ = 0.75 ps⁻¹ and ζ^⊥₃ = 1.0 ps⁻¹).

Thus, employing the new thermostat one can reproduce both the structure and the dynamics of the all-atom liquid water with the single-site coarse-grained water model. This is essential for synchronizing the timescales of the all-atom and coarse-grained regimes in the adaptive resolution MD simulations.²²

V. Conclusions

In this paper we introduced an extension of the dissipative particle dynamics thermostat that allows for controlling the transport properties of molecular liquids, e.g., water, while preserving the hydrodynamics.

The presented Galilean invariant local thermostat can be used in coarse-grained simulations to tune the diffusion constant and viscosity of the system to the desired values. This opens up the possibility of reproducing the atomistic dynamics by coarse-grained simulations, as it is required, for example, in recently introduced adaptive resolution simulations.

Acknowledgements

We are grateful to R. Delgado-Buscalioni and B. Dünweg for discussions at early stage of this work. We also thank J. Kirkpatrick for critical reading of the manuscript.

References

1. V. A. Harmandaris, N. P. Adhikari, N. F. A. van der Vegt and K. Kremer, Macromolecules, 2006, 39, 6708–6719.
2. B. Hess, S. Leon, N. F. A. van der Vegt and K. Kremer, Soft Matter, 2006, 2, 409–414.
3. N. F. A. van der Vegt, C. Peter and K. Kremer, Structure-based coarse- and fine-graining in soft matter simulations, in Coarse-graining of Condensed Phase and Biomolecular Systems, ed. G. A. Voth, Chapman and Hall/CRC Press, Taylor and Francis Group, 2007.
4. P. J. Hoogerbrugge and J. M. V. A. Koelman, Europhys. Lett., 1992, 19, 155–160.
5. J. M. V. A. Koelman and P. J. Hoogerbrugge, Europhys. Lett., 1993, 21, 363–368.
6. P. Español and P. Warren, Europhys. Lett., 1995, 30, 191–196.
7. R. D. Groot and P. B. Warren, J. Chem. Phys., 1997, 107, 4423–4435.
8. P. Español, Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top., 1998, 57, 2930–2948.
9. L. Delle Site, C. F. Abrams, A. Alavi and K. Kremer, Phys. Rev. Lett., 2002, 89, 156103.
10. R. Delgado-Buscalioni and P. V. Coveney, Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top., 2003, 67, 046704.
11. S. Barsky, R. Delgado-Buscalioni and P. V. Coveney, J. Chem. Phys., 2004, 121, 2403–2411.
12. R. Delgado-Buscalioni, E. G. Flekkoy and P. V. Coveney, Europhys. Lett., 2005, 69, 959–965.
13. M. Neri, C. Anselmi, M. Cascella, A. Maritan and P. Carloni, Phys. Rev. Lett., 2005, 95, 218102.
14. E. Villa, A. Balaeff and K. Schulten, Proc. Natl. Acad. Sci. U. S. A., 2005, 102, 6783–6788.
15. G. D. Fabritiis, R. Delgado-Buscalioni and P. V. Coveney, Phys. Rev. Lett., 2006, 97, 134501.
16. Q. Shi, S. Izvekov and G. A. Voth, J. Phys. Chem. B, 2006, 110, 15045–15048.
17. M. Praprotnik, L. Delle Site and K. Kremer, J. Chem. Phys., 2005, 123, 224106.
18. M. Praprotnik, L. Delle Site and K. Kremer, Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top., 2006, 73, 066701.
19. M. Praprotnik, K. Kremer and L. Delle Site, Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top., 2007, 75, 017701.
20. M. Praprotnik, K. Kremer and L. Delle Site, J. Phys. A: Math. Theor., 2007, 40, F281–F288.
21. M. Praprotnik, L. Delle Site and K. Kremer, J. Chem. Phys., 2007, 126, 134902.
22. M. Praprotnik, S. Matysiak, L. Delle Site, K. Kremer and C. Clementi, J. Phys.: Condens. Matter, 2007, 19, 292201.
23. W. Tschöp, K. Kremer, J. Batoulis, T. Bürger and O. Hahn, Acta Polym., 1998, 49, 61–74.
24. S. Izvekov and G. A. Voth, J. Chem. Phys., 2006, 125, 151101.
25. B. Carmeli and A. Nitzan, Chem. Phys. Lett., 1983, 102, 517–522.
26. M. J. Moix and R. Hernandez, J. Chem. Phys., 2005, 122, 114111.
27. P. Grigolini, J. Chem. Phys., 1988, 89, 4300–4308.
28. J. B. Straus and G. A. Voth, J. Chem. Phys., 1992, 96, 5460–5470.
29. S. Matysiak, C. Clementi, M. Praprotnik, K. Kremer and L. Delle Site, ArXiv:0710.0048v1, 2007.
30. T. Soddemann, B. Dünweg and K. Kremer, Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top., 2003, 68, 046702.
31. P. Español, Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top., 1995, 52, 1734–1742.
32. C. A. Marsh, G. Backx and M. H. Ernst, Europhys. Lett., 1997, 38, 411–415.
33. Note that due to the Transverse DPD thermostat, forces between molecules are not central and therefore the local angular momentum is not conserved. Nevertheless, simulation results show that these angular momentum fluctuations cancel out on average and hence the total angular momentum of the system is almost a conserved quantity. Moreover, after the submission of this manuscript it has been shown that only in special circumstances, i.e., when the boundary conditions on walls are given by forces, fluids with different viscosities are in contact, or finite-sized objects rotate in fluids, the angular momentum conservation plays a significant role.³⁴ For the adaptive resolution simulations, the present thermostat can be used to match the viscosities to avoid the above mentioned problems.
34. I. O. Götze, H. Noguchi and G. Gompper, ArXiv:0709.1600v1, 2007.
35. R. Zwanzig, J. Chem. Phys., 1960, 33, 1338–1341.
36. R. Zwanzig, Phys. Rev., 1961, 124, 983–992.
37. D. Forster, Hydrodynamic Fluctuation, Broken Symmetry, and Correlation Functions, Benjamin, Reading, MA, 1975.
38. H. J. Limbach, A. Arnold, B. A. Mann and C. Holm, Comput. Phys. Commun., 2006, 174, 704–727 http://www.espresso.mpg.de .
39. W. L. Jorgensen, J. Chandrasekhar, J. D. Madura, R. W. Impey and M. L. Klein, J. Chem. Phys., 1983, 79, 926–935.
40. M. Neumann, Mol. Phys., 1983, 50, 841–841.
41. M. Neumann, J. Chem. Phys., 1985, 82, 5663–5672.
42. I. G. Tironi, R. Sperb, P. E. Smith and W. F. van Gunsteren, J. Chem. Phys., 1995, 102, 5451–5459.
43. H. Berendsen, J. Postma, W. V. Gunsteren, A. D. Nola and J. Haak, J. Chem. Phys., 1984, 81, 3684–3690.
44. S. Nosé, J. Chem. Phys., 1984, 81, 511–519.
45. W. G. Hoover, Phys. Rev. A, 1985, 31, 1695–1697.
46. G. J. Martyna, M. L. Klein and M. Tuckermann, J. Chem. Phys., 1992, 97, 2635–2643.
47. H. C. Andersen, J. Chem. Phys., 1980, 72, 2384–2393.
48. C. P. Lowe, Europhys. Lett., 1999, 47, 145–151.
49. E. A. J. F. Peters, Europhys. Lett., 2004, 66, 311–318.
50. E. A. Koopman and C. P. Lowe, J. Chem. Phys., 2006, 124, 204103.
51. S. D. Stoyanov and R. D. Groot, J. Chem. Phys., 2005, 122, 114112.
52. K. Kremer and G. S. Grest, J. Chem. Phys., 1990, 92, 5057–5086.
53. C. Pastorino, T. Kreer, M. Müller and K. Binder, Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top., 2007, 76, 026706.
54. B. Dünweg and K. Kremer, J. Chem. Phys., 1993, 99, 6983–6997.
55. E. Efron, The Jackknife, the Bootstrap and other Resampling Plans, SIAM, Philadelphia, 1982.
56. Y. Wu, H. L. Tepper and G. A. Voth, J. Chem. Phys., 2006, 124, 024503.

---