Edinburgh Research Explorer On the numerical treatment of dissipative particle dynamics and related systems

We review and compare numerical methods that simultaneously control temperature while preserving the momentum, a family of particle simulation methods commonly used for the modelling of complex ﬂuids and polymers. The class of methods considered includes dissipative particle dynamics (DPD) as well as extended stochastic-dynamics models incorporating a generalized pairwise thermostat scheme in which stochastic forces are eliminated and the coeﬃcient of dissipation is treated as an additional auxiliary variable subject to a feedback (kinetic energy) control mechanism. In the latter case, we consider the addition of a coupling of the auxiliary variable, as in the Nosé–Hoover– Langevin (NHL) method, with stochastic dynamics to ensure ergodicity, and ﬁnd that the convergence of ensemble averages is substantially improved. To this end, splitting methods are developed and studied in terms of their thermodynamic accuracy, two-point correlation functions, and convergence. In terms of computational eﬃciency as measured by the ratio of thermodynamic accuracy to CPU time, we report signiﬁcant advantages in simulation for the pairwise NHL method compared to popular alternative schemes (up to an 80% improvement), without degradation of convergence rate. The momentum- conserving thermostat technique described here provides a consistent hydrodynamic model in the low-friction regime, but it will also be of use in both equilibrium and nonequilibrium molecular simulation applications owing to its eﬃciency and simple numerical implementation.


Introduction
Stochastic momentum-conserving thermostats, which correctly capture long-ranged hydrodynamic interactions, are increasingly popular tools for simulation of simple and complex fluids [1]. The first important scheme of this type was dissipative particle dynamics (DPD), introduced by Hoogerbrugge and Koelman [2] in 1992 for simulating complex hydrodynamic behavior at a mesoscopic level that is not accessible by conventional molecular dynamics (MD) [3,4]. In DPD, a collection of fluid molecules are grouped at the coarse-grained level and treated as a discrete particle. These particles interact at short range in a soft potential, thereby allowing larger timesteps than would be possible in MD, while simultaneously decreasing the number of degrees of freedom required. DPD thus bridges the gap between microscale (atomistic methods, e.g. molecular dynamics) and macroscale (continuum methods, e.g. Navier-Stokes) models and can be used to recover thermodynamic, dynamical and rheological properties of complex fluids, with applications to colloidal particles [5], polymer molecules [6] and fluid mixtures [7].
A great deal of effort has been devoted to the design of simple, efficient and accurate numerical methods to solve the DPD system due to its promising prospects from the applications perspective. In the example of lipid bilayers, new phenomena arise as the time scale of the system that we are investigating is increased [8]. However, not all algorithms are rigorously founded and may not perform satisfactorily in large scale simulations (see discussions in [9][10][11][12]). All the methods exhibit pronounced artifacts with increasing integration stepsizes due to the discretization error, typically manifest in the form of statistical bias in the calculation of thermodynamic averages. A previous study [13] suggested that, without performing serious checks for each method, the only reliable approach is to use vanishingly small stepsizes, since most of the schemes proposed are formally convergent at some order of accuracy. However, we argue that using very small stepsize significantly limits the time scales accessible for DPD simulation, particularly in large scale simulations.
One of the key questions we tackle in this article is "What is the largest integration stepsize that can be used without damaging both static and dynamical properties?" Answering this question leads to a better understanding of the overall efficiency of each method and the validity of the schemes in the computation-intensive setting (large particle number/longtime interval). Recently, a systematic approach to thermodynamic bias in numerical computations has been used to study the accuracy and efficiency of methods for Langevin dynamics [14][15][16]. The approach suggested is to determine the order of accuracy of a stochastic scheme in relation to its effective invariant distribution (and, thus, with respect to steady state averages computed using numerical trajectories). This technique has led to greatly improved numerical methods for Langevin dynamics, and it is in principle applicable to momentum-conserving thermostat schemes as well. However these results are based on asymptotic expansion and thus are only relevant in the limit of small stepsize (whereas we are interested in the large stepsize threshold). Moreover we find that the analytical computations necessary to perform expansions in the DPD and pairwise NHL cases are in typical cases highly complex; we therefore restrict ourselves in this article to outlining some fundamental and illustrative applications of the theory. In particular for certain symmetric methods, we can demonstrate the even-order approximation of long-time averages. The same conclusion may be reached for additional schemes of a specific structure (related to the Geometric Langevin Algorithms of [14]).
A new stochastic momentum-conserving thermostat is introduced in this article that can be used in place of DPD in the low-friction regime or in nonequilibrium molecular dynamics (NEMD) based on stochastic extension of a scheme in a recent paper by Allen and Schmid [1]. This method is particularly inexpensive to implement and is found to have a very high stepsize stability threshold compared to alternatives.
In typical cases, for understanding the stepsize stability threshold and the performance of different schemes, we are forced to rely on numerical experiment. An excellent survey of the performances of a number of methods was undertaken by Nikunen et al. [11] in 2003. Since then, and despite many additions to the arsenal of methods, such a comparison has been lacking. To this end, we test a number of popular methods [9,17,18] from the DPD literature together with additional methods [19][20][21] that are used in popular software packages. For each method, we examine calculations such as kinetic and configurational temperatures, average potential energy, radial distribution function, velocity autocorrelation function and transverse momentum autocorrelation function, which gives information on the rotational relaxation process.
The outline of the article is as follows. We first review the formulation and various integration schemes of DPD in Section 2. Two extended variable momentum-conserving thermostats are presented in Section 3 with the latter newly proposed in this article. In Section 4, we demonstrate the definitions of error in DPD simulation and summarize the orders of accuracy both theoretically and numerically in a number of methods. Various numerical experiments are carried out in Section 5 to compare all the schemes described in the article. We summarize our findings in Section 6.

Formulation of DPD
The original version of DPD [2] updated the system in discrete time steps and was later reformulated by Español and Warren in 1995 [22] as a system of stochastic differential equations (SDEs).
Consider N particles with positions q i , momenta p i and masses m i for i = 1, . . . , N evolving in dimension d. The time evolution of the DPD particles, each of which represents a cluster of molecules, is governed by Newton's equations where F i is the total interparticle force acting on particle i due to the presence of the other particles. The force is composed of three pairwise contributions where F C ij , F D ij and F R ij represent conservative, dissipative and random forces, acting on particle i due to the j-th particle, respectively.
The conservative force controlling the thermodynamics of the DPD system is normally [23] chosen as where parameters a ij are symmetric (a ij = a ji ≥ 0) representing the maximum repulsion strength between particle i and particle j, and r c is a certain cutoff radius. The relative positions are denoted by q ij = q i − q j with length q ij = |q ij | and the unit direction from q j to q i by q ij = q ij /q ij . It should be noted that there are many alternative choices for the mesoscopic force field. For example, a Lennard-Jones type potential has been used in combination with DPD thermostat in simulating equilibrium and nonequilibrium molecular dynamics [24]. Potentials may also be obtained by direct coarse-graining of a microscopic system [25,26]. Although we assume the simple force law (3), many of the observations of this article would hold for more general conservative force fields.
The dissipative and random forces behave as a thermostat to keep the temperature constant and are given by where γ and σ represent the dissipative and random strengths respectively, ω D and ω R are position-dependent weight functions, and the relative velocities are denoted by term with the following stochastic property and is chosen independently for each pair of interacting particles at each time step. The canonical ensemble was not preserved in the original formulation of DPD of [2]. This has been corrected by Español and Warren [22] who showed that certain conditions have to be satisfied to guarantee that the DPD system has the same invariant distribution as would be obtained in cases where the dissipative and random forces are put to zero, namely: where k B is the Boltzmann constant and T the equilibrium temperature. This is the fluctuation-dissipation theorem for the DPD thermostat involving dissipative and random forces only. Then it can be easily shown that the canonical ensemble is preserved with invariant distribution defined by the density where β −1 = k B T , Z is the partition function and H(q, p) is the Hamiltonian defined as where the "soft" pair potential energy U (q ij ) corresponding to the conservative force (3) is defined as Although we write the density (8) as an exponential, we note that if the total momentum is conserved, the density should be replaced by where π = (π x , π y , π z ) is the total momentum vector. A similar modification would be needed if the angular momentum were also conserved. It is worth pointing out that the ergodicity of the DPD system has only been demonstrated in the case of high particle density in one dimension by Shardlow and Yan [27]. Due to the fact that the algorithm depends on relative velocities and the interactions between particles are symmetric, both total and angular momenta are conserved, DPD is therefore an isotropic Galilean-invariant thermostat which preserves hydrodynamics [1,28]. If periodic boundary conditions are used, the angular momentum will be destroyed as a conserved quantity, but the symmetric character of the force field has important consequences e.g. for the relaxation of a localized rotational excitation. An early study of constant-temperature molecular dynamics with momentum conservation can be found in [29], but the techniques discussed there are not that relevant to the present work.
One of the two weight functions can be chosen arbitrarily without changing thermodynamic equilibrium. A simple and widely-used choice reads (11) in which case the conservative force (3) can be written in a compact way To make the presentation simpler, a compact form of the SDEs of the DPD system (1) (for particle i) may be used [1] where F C i (q) is the total conservative forces acting on particle i (13) and V i (q, p) and R i (q, p, t) are defined respectively as (15) where dW ij (t) = dW ji (t) are independent increments of a Wiener process with mean zero and variance dt [23].

Numerical integration schemes
Due to the soft repulsive potential (10), the major advantage of the DPD method is that the stepsize used in simulations may be much larger than those of conventional MD simulations with Lennard-Jones internuclear potentials for instance. This feature is crucial especially when a very long simulation time is required. However, large stepsizes may result in errors in computed thermodynamic quantities. Whereas for molecular dynamics, simulations are performed at or near the stability threshold defined by stiff components such as harmonic bonds [16], DPD simulations may be perfectly stable over a wide range of stepsizes for which errors in averages are very large. There have been many attempts to develop accurate and efficient numerical methods that allow larger stepsizes. This is currently an active field of research.
In the early days of DPD, a number of integration schemes were proposed based on the well-recognized velocity Verlet scheme [30] widely used in classical MD simulations [3,4]. Specific examples are the integrator of Groot and Warren [23] (GW), the method of Gibson et al. [31] (GCC) and the DPD velocity-Verlet integrator, which refers to as DPD-VV, of Besold et al. [9] (see more discussions in [10]). Both the GW and GCC integrators incorporate a parameter λ, which has to be chosen carefully for specific model parameters, to reduce unphysical artifacts. Nevertheless, as reported in [9,10], all the integrators mentioned above display pronounced artifacts, especially when the stepsize is large, due to the fact that the velocities and dissipative forces depend on each other implicitly and thus need to be updated in a self-consistent fashion.
In the spirit of the self-consistent leap-frog integrator introduced by Pagonabarraga et al. [32], Besold et al. [9] proposed a self-consistent velocity Verlet scheme, which we label SC-VV, in which, at the end of each iteration step, the velocity is "corrected" based on the newly-calculated dissipative force until the deviation between the instantaneous kinetic temperature and the target temperature is less than a certain value of "tolerance". It should be noted that there is no such "correction" in the DPD-VV scheme but still the dissipative force is recalculated once, using the up-to-date velocities (momenta). A variant of the SC-VV scheme, SC-Th, introduced in [9], couples the original DPD system to an auxiliary Nosé-Hoover thermostat [33][34][35][36] to provide direct kinetic temperature control. Overall, the self-consistent schemes do reduce unphysical artifacts to some extent, however, it is also well-documented [11,12] that they can be substantially slower than standard methods, depending how small the tolerance is. Therefore, computationally expensive self-consistent methods are not the choices present in typical software packages.
Although it has been demonstrated that relatively small stepsizes must be used in the DPD-VV scheme to produce correct static and dynamical properties [10,11], the DPD-VV scheme remains one of the most popular methods for DPD system in software packages due to its efficiency and ease of implementation (particularly in parallel computing). For this reason, the DPD-VV method has been chosen as the "benchmark" to compare with other schemes in this article.
Several novel integration schemes have been proposed over the years, such as the approach by den Otter and Clarke [37], the extended DPD by Cotter and Reich [38], the multiple time step schemes by Jakobsen et al. [39], the Trottersplitting methods by Thalmann and Farago [40] and most recently the algorithm by Goga et al. [41]. In this article, we focus on schemes that have been included in popular software packages (i.e. DPD-VV [9], Peters thermostat [20], Lowe-Andersen thermostat [19] and Nosé-Hoover-Lowe-Andersen thermostat [21] in the DL_MESO [42] package) with two promising splitting methods by Shardlow [17] and De Fabritiis et al. [18], respectively. Particularly, Nikunen et al. [11] in 2003 showed that the performances of Lowe-Andersen thermostat and Shardlow's scheme are superior to those of several other schemes for a number of different observables. Recently, Sharlow-like splitting algorithms have been further applied in DPD with various fixed conditions [43], and its parallel implementation has also been developed in [44].
The full details of the integration steps of each method described in this article are presented in a common language in Appendix B.

DPD velocity-Verlet
For integration stepsize h, the simple DPD-VV integrator [9] reads In fact, the DPD-VV scheme has two differences compared to the standard velocity Verlet method [30]: (1) the forces are not just the conventional conservative forces, but include dissipative and random forces, as well; (2) the dissipative forces have to be updated for a second time at the end of each integration step by using the up-to-date velocities (momenta), F D i (q n+1 , p n+1 ), with the first update taking place right after the positions are updated at each step (see more details in Appendix B). It has been shown that the performance of the DPD-VV method would be significantly improved [9,10] by simply doing the additional update of the dissipative forces in each step, which is actually not time-consuming if one makes use of computation-saving devices such as Verlet neighbor lists [30]. Note that both the GW integrator [23] of Groot and Warren and the modified Verlet method mentioned by Shardlow in [17] do not incorporate the additional update.
It is important to observe that, unlike the standard velocity Verlet method (second order), the DPD-VV scheme should only give first order convergence for the invariant measure (see more details in Section 4) due to the fact that the momentum is not updated in a symmetric manner. However, second order convergence for averages was observed in the numerical experiments in Section 5. Moreover, the term √ h/2 multiplying the random forces, which would be expected to be different in this type of splitting of Langevin dynamics, must be present when random forces are reused in the subsequent step; it ensures that the diffusion coefficient of the particles is independent of the integration timestep (see [23] for further discussion).

Shardlow's splitting method
Splitting techniques were studied by Shardlow [17] based on dividing the vector field of the DPD system into three parts, the first two of which represent the vector field of the Hamiltonian system associated to kinetic and potential energies, and the last term is the remaining Langevin equation (actually Ornstein-Uhlenbeck process with positions fixed) involving the dissipative and random forces. Two integrators, termed there S1 and S2, have been proposed for treating this system in [17]. Only the S1 method will be examined here as it is more efficient than S2. This scheme relies on the method of Brünger, Brooks and Karplus (BBK) [45] to solve the Langevin part, followed by the standard velocity Verlet scheme for the conservative part.
In describing Shardlow's method or another splitting scheme, we use the formal notation of the generator of the diffusion as in, for example, [18,46,40].
We first separate the system of stochastic differential equations for DPD (12) into three pieces, which we label as A, B and O: The generators for each part of the SDE may be written out as follows: Thus, the generator for the DPD system can be written as The flow map (or phase space propagator) of the system may be written in the shorthand notation where the exponential map is here used formally to denote the solution operator. Approximations of F t may then be obtained as products (taken in different arrangements) of exponentials of the various terms of the splitting. For example, the phase space propagation of Shardlow's S1 splitting method [17], termed as DPD-S1, can be written as where h is the stepsize and exp(hL f ) represents the phase space propagator associated to the corresponding vector field f .
In Shardlow's approach, the vector field O is further split into each interacting pair. Therefore, the propagation of the O part is broken down into many terms: where the operators associated to each interacting pair are defined as Each interacting pair preserves the invariant distribution ρ β (8). As a shorthand, we may term the DPD-S1 method OBAB (similarly, the S2 method of Shardlow would be equivalent to BABOBAB in the same language), where, in both cases, it is to be understood that the steplengths associated to various operations are uniform and span the interval h. Thus the B step in OBAB is taken with a step length of h/2, while O with a steplength of h.

DPD-Trotter scheme
The Trotter formula [47] that has been widely used in molecular simulations was investigated and further applied to split the DPD generator (23) in an "optimal" way to reduce artifacts and maintain good temperature control [18,46]. A new scheme, referred to as DPD-Trotter, was introduced but few numerical simulations have been presented and therefore we incorporate it in the comparisons.
In the stochastic DPD-Trotter scheme, the standard DPD system (12) is split into two parts, which are labeled "A" and "S" as indicated below The corresponding operator of part A is exactly the same as in Shardlow's method, while the operator of part S is actually the sum of the operators of B and O defined above As in Shardlow's method, the vector field S is further split into each interacting pair; these pair interactions are exactly solvable (in the sense of distributional fidelity). In fact, for j > i, subtracting dv j from dv i and multiplying q ij on both sides The above equation is an Ornstein-Uhlenbeck process with the exact (in the sense of distributions) solution [48] v ij (t) = and the corresponding momenta can be updated by which defines the propagator e hL S i, j for each interacting pair. Overall, the propagator of the DPD-Trotter scheme can be written as where the momentum part is defined by

Lowe-Andersen thermostat
The Schmidt number, S c , which is the ratio of the kinematic viscosity ν to the diffusion coefficient D, is an important quantity that characterizes the dynamical behavior of fluids. In a typical fluid flow, water for example, momentum can be transported more rapidly than particles, resulting in Schmidt number of the order 10 3 . However, as it was reported in [23], the standard DPD system (12) produces a gas-like Schmidt number of the order 1. Depending on the application, this could be a significant disadvantage for simulating more fluid-like system, although recent work by Fan et al. [49] indicates that the Schmidt number of the standard DPD system can be varied by modifying the weight function and increasing the cutoff radius.
To overcome the issue of low Schmidt number in standard DPD system, instead of using a Langevin thermostat to re-equilibrate the system, Lowe [19] employed a pairwise stochastic momentum-conserving Andersen thermostat [50], in which after updating the positions and momenta due to the conservative forces only by using the standard velocity Verlet method, the momenta are updated, with probability P = Γ h, by reselecting the relative velocities for interacting pairs from the Maxwell-Boltzmann distribution, where R ij (t) are Gaussian random variables with zero mean and unit variance. The parameter Γ can be thought of as the stochastic randomization frequency with upper limit 1/h. Lowe's method is frequently referred to as the Lowe-Andersen (LA) thermostat, which still conserves the momentum and hydrodynamics. The additional Andersen thermostat does not change the distribution of the system [50], therefore the same invariant distribution (8) as in the standard DPD system is maintained. Most important, the LA thermostat is capable of varying the Schmidt number by changing the parameter Γ .
When the probability of further updating the momentum is high (large P ), the viscosity is high and diffusion coefficient is low, resulting in large Schmidt number in the regime of typical fluids. The LA thermostat has been applied in molecular dynamics simulation by Koopman and Lowe [51]. It is worth mentioning that the generator of the LA thermostat does not converge to that of the standard DPD system as h → 0.

Peters thermostat
Based on various numerical studies on the DPD system, all the numerical methods based on discretization of the equation of motion are dependent on the stepsize. In order to reduce the dependence, Peters generalized the Lowe-Andersen (LA) thermostat and presented another approach to re-equilibrate the system [20]. Following the update of the conservative part by using the standard velocity Verlet scheme, the momenta of all interacting pairs will be further updated (not in random order as in the original paper, which does not have much effect on the results but reduces computational cost) as follows where R ij (t) are the standard Gaussian random variables as in Lowe-Andersen thermostat, and the coefficients γ ij and σ ij satisfies the following condition which reduces to the fluctuation-dissipation theorem (7) in standard DPD in the limit of h → 0. Two possible choices of the coefficients was presented in the original paper [20], but only "Scheme II", which has less restriction on the choice of stepsize h, was chosen to compare with other methods in this article. In the h → 0 limit, it can be easily shown that the generator of the Peters thermostat converges to that of the standard DPD system and therefore conserves the canonical ensemble. Unfortunately, numerical simulations in the original paper showed that the thermostat still exhibits significant deviation both in static (kinetic temperature) and dynamical (diffusion coefficient) quantities in standard model settings (model B in the language of [11]) when the timestep is above 0.05.

Nosé-Hoover-Lowe-Andersen thermostat
Recently, Stoyanov and Groot combined the Lowe-Andersen (LA) thermostat with a Nosé-Hoover-like thermostat to construct a local Galilean invariant stochastic momentum-conserving thermostat, Nosé-Hoover-Lowe-Andersen (NHLA) thermostat [21], to achieve direct kinetic temperature control. Unlike the LA thermostat, a modified version of the velocity Verlet scheme is used to update the positions and velocities at the start of each integration step which is followed by the calculation of the updated conservative forces. After that, the fraction (1 − P ) of interacting pairs that do not have their relative velocities stochastically reselected are thermalized by a deterministic method instead. For each such pair, the dissipative force is calculated where α is a coupling parameter chosen as 0.9/(ρh) in this article such that, overall, the dissipative force defined above is the same as that in the original paper [21], and ρ is the particle density. The dissipative force of each particle is updated Then, after the further update of the velocities the momenta are corrected by where T 0 is the desired temperature and the momentary kinetic temperature T k is calculated from the relative velocities at the time of calculating the conservative force to enhance computational efficiency and save memory space (this is slightly different from the approach in the original paper which uses the updated Verlet neighbor lists but the stored velocities from the previous integration step) where ω(q ij ) is a smearing function defined as Finally, the momenta are further updated as in LA thermostat (Eqs. (33)- (34)).
The factor (1 −T k /T 0 ) in (44) acts like the "friction coefficient" to tune the kinetic temperature of the system. It is actually not a dynamical variable as in the standard Nosé-Hoover thermostat, instead is more closely related to the Berendsen thermostat [52]. As reported in [21], the NHLA thermostat maintains an order of magnitude improvement in kinetic temperature and can also vary the Schmidt number by several orders of magnitude as in the LA thermostat. However, with large stepsizes that maintain good kinetic temperature control (1% relative error) in the NHLA thermostat, substantial errors in configurational temperature were reported [53], which indicates that the system temperature was not sampled correctly. It is worthy of mention that a slightly modified integration strategy was used in [53], which does not have much effect on the results as discussed in the original paper [21]. Moreover, the generator of the NHLA thermostat does not converge to that of the standard DPD system as h → 0, and, it has not been shown that the NHLA thermostat preserves the canonical ensemble.

Pairwise Nosé-Hoover thermostat
In all the standard DPD methods (Section 2.2) and alternative methods (Section 2.3) in DPD simulations, a random number has to be generated for each interacting pair, which can be very time-consuming depending on the particle density and cutoff radius, and, becomes trickier when parallel computing techniques (multiple processors for domain-decomposed cells) are used: it requires additional, even substantial, effort to communicate interacting particles in different cells [54,55]. Based on the Nosé-Hoover-Lowe-Andersen (NHLA) thermostat by Stoyanov and Groot [21], Allen and Schmid [1] presented a new thermostat of the Nosé-Hoover type, in which stochastic terms were totally removed and the constant friction coefficient was replaced by a dynamical variable that was driven by the difference between the instantaneous kinetic temperature and the target temperature. The so-called pairwise Nosé-Hoover (PNH) thermostat conserves the momentum and is also Galilean-invariant, therefore correct hydrodynamics are still expected to be generated and it can be used in DPD simulations. Moreover, one may find it useful in nonequilibrium molecular dynamics (NEMD) to reduce unphysical behaviors (see more discussions in [1]).
The equation of motion of the PNH thermostat (for particle i) is given by where ξ is the dynamical variable and G(q, p) is the instantaneous accumulated deviation of the kinetic temperature away from the target temperature [1] where μ is a coupling parameter which is referred to as the "thermal mass". Canonical ensemble is still preserved with modified invariant distribution than ρ β (8) in standard DPD system A nonsymmetric integration algorithm (see Appendix B) was applied in the original paper [1] to solve the system.

Pairwise Nosé-Hoover-Langevin thermostat
In order to enhance the ergodicity and have a better temperature control, we have reformulated the Pairwise Nosé-Hoover (PNH) thermostat to form a pairwise Nosé-Hoover-Langevin (PNHL) thermostat by adding a Langevin thermostat to the additional variable ξ in such a way that the invariant distribution (49) is not violated. As in the PNH thermostat, the PNHL thermostat has the potential of being useful in NEMD, but we focus on the application of DPD in this article.
The SDEs of the PNHL system (for particle i) is given by where coefficient constants γ * and σ * satisfy the fluctuation-dissipation theorem in standard Langevin dynamics and W = W(t) is a standard Wiener process.
The vector field of the PNHL system can be split into five pieces below such that each piece can be solved "exactly", Note that the operators of A and B are exactly the same as defined in (20) and (21), respectively. We can also write down the operators for the remaining pieces as The generator O here is slightly different from that in (22) which involves pairwise terms in DPD. Overall, the generator for the PNHL system can be written as There are a variety of approaches to splitting this system. For example, we could use the same technique as in DPD-Trotter scheme (Section 2.2.3) to solve part C, but without the conservative and stochastic terms. Also, the O part may be solved exactly using the analytical weak solution of the Ornstein-Uhlenbeck process [48] ξ(t) = e −γ * t ξ(0) + k B T 1 − e −2γ * t /μR(t), (54) where ξ(0) is the initial value of the additional variable and R(t) are uncorrelated independent standard normal random variables.
Interestingly as noted in the setting of Langevin dynamics [15,16], integrating those different splitting pieces in different orders gives dramatically different performances in terms of kinetic temperature control and other configurational quantities. We present here two approaches to integrate the PNHL system: first in a symmetric manner, termed PNHL-S, and the other nonsymmetric, termed PNHL-N. The propagators of the two schemes (see more details in Appendix B) may be defined as It is important to mention that the only difference between these two integrators is the order of integrating the last two pieces. In particular, an additional force calculation is needed in the PNHL-N scheme just before updating the last B piece at the end of each integration step. In experiments, the high per-timestep cost of PNHL-N was found to be offset by a great increase in accuracy and usable steplength. Detailed numerical comparisons will be presented in Section 5.

Error in DPD simulations
To our knowledge, little attention has been paid to the mathematical analysis of the DPD system, or more generally stochastic momentum-conserving thermostats. Since the spectral properties of the Fokker-Planck operators in the case of DPD is not available, a rigorous study of the order of convergence of numerical methods in this context has been lacking. Because of the inclusion of stochastic terms, it is not reliable to depend directly on intuition regarding the error behavior of deterministic schemes [56,57]. Here we perform a few first steps toward the analysis of stochastic DPD integrators by extending a framework for investigating the perturbation of long-time average computed using numerical methods in Langevin dynamics proposed recently by Leimkuhler and Matthews [15,16].
In DPD simulations, the error in averages related to the evolving distribution is generally of interest, i.e. the weak error (finite-time error in the weak sense) for nonequilibrium and dynamical properties, and, long-time error (error in the invariant distribution obtained as t → ∞) for thermodynamics. We next give definitions of these two errors following [58].
For the weak error, consider a finite time interval [0, τ ] with τ = nh. The probability measure associated to certain system is described by a probability density ρ(z, t) which evolves according to the Fokker-Planck equation ∂ρ ∂t where L † is the adjoint of the generator of the system (the Fokker-Planck operator of the standard DPD system is given in [22]). Assuming ergodicity, the solution ρ(z, t) evolves from an initial probability distribution ρ(z, 0) to the steady state (invariant distribution) ρ(z, ∞) = ρ β . For a smooth and bounded test function φ of a suitable class, the average of φ at time τ may be defined bȳ The discretization scheme can also be viewed as giving rise to an evolving sequence of probability distributions ρ 1 , ρ 2 , · · ·. With stepsize h, the average at time τ = nh is defined aŝ φ(τ , h) = φ(z)ρ n (z)dz. (56) Thus, we could define the weak error as the difference between (55) and (56) where the coefficient K depends on the time interval and p is the order of the weak error. To be more precise, K also depends on the initial distribution ρ(z, 0) as well as the particular observable φ. The asymptotic (τ → ∞) behavior [59] of K can be used to describe the performance of the numerical method for computing averages with respect to the invariant distribution. Hence, the long-time error in averages can be written as

≥1(2)
According to the Baker-Campbell-Hausdorff (BCH) formula, a nonsymmetric splitting method generally gives only first order convergence for the invariant measure (long-time error), whereas second order is anticipated in symmetric splittings in the asymptotic limit of small stepsize [60]. In our numerical studies, we have verified second order convergence for a number of nonsymmetric methods (Shardlow's splitting method, Lowe-Andersen thermostat and Peters thermostat) which are similar to the family of Geometric Langevin Algorithms (GLA) following [14] (see more details in Appendix A). We compute the long-time error of various observables, including kinetic and configurational temperatures and average potential energy, and demonstrate the convergence order for each method in Table 1. The results shown in Table 1 are based on the numerical experiments in Section 5. All the symmetric methods show second order as well as those three nonsymmetric ones (DPD-S1, LA and Peters). Some other nonsymmetric methods, which are not of GLA type, exhibit second order convergence in calculated quantities; this observation remains to be demonstrated rigorously. It is entirely possible that the superconvergence observed in these special cases is related to the form of the observable we have used in our simulation test.

Numerical experiments
To investigate the performance of all the numerical methods described in this article, we perform systematic numerical experiments in this section.

Simulation details
Tests have been carried out by using the standard parameter set [23] that is commonly used in algorithms tests [10,11,17,20,21,46]. A system of N = 500 identical particles (m i = m = 1) was simulated in a cubic box (length L = 5) with periodic boundary conditions (particle density ρ = 4). The cut-off radius is r c = 1 and k B T = 1. The potential repulsion parameter a ij was set to 25, while dissipative strength γ was chosen as 4.5, resulting the value of random strength σ to be 3. It is worthy of mention that we did investigate the influence of other values of the friction coefficient, such as γ = 0.5 and γ = 40.5, but not all the results will be presented unless necessary. The initial positions of the particles were IID (independent, identically distributed) uniformly distributed over the box, while the initial momenta were IID normally distributed with mean zero and variance k B T . Verlet neighbor lists [30] were used in all the simulations as long as possible.
In particular, the thermal mass μ in PNH and PNHL thermostats has to be chosen with care. For PNH thermostat, we used μ = 200 to maintain as good stability as possible for the integration scheme, while μ = 10 and γ * = 4.5 were used in PNHL thermostat. Since the focus of this article is on DPD simulations, the stochastic randomization frequency Γ in LA and NHLA thermostats was set to be 0.44 as in [11,13] so that similar translational diffusion properties of the fluid were obtained.

Static
The kinetic temperature T k appears to be the most popular quantity to be tested in DPD literature, which is defined as But in practice, the kinetic temperature is not that important, instead those configurational quantities play more crucial roles. Recent studies [61,62] in simulations showed that the system temperature can be measured from static snapshots of its constituents' instantaneous configurations rather than their momenta. We test the configurational temperature T c [63], which can be defined as where the angle brackets denote the averages, ∇ i U and ∇ 2 i U are respectively the gradient and Laplacian of the potential energy U with respect to the position of particle i. To test the correctness of codes and/or algorithms in simulations, both kinetic and configurational temperatures can be used to calculate the system temperature (in principle they should both be equal to the target). However, it turned out that the configurational temperature is more reliable [61,62] since it can rapidly ) Log-log plot of the relative error in kinetic temperature against stepsize by using various numerical methods. The system was simulated for 1000 reduced time units but only last 80% data were collected to calculate the static quantity to make sure the system was well equilibrated. The stepsizes tested began at h = 0.05 and were increased incrementally by 15% until all methods either were above 100% relative error or became unstable. and accurately track changes in system temperature even when the system is not in global thermodynamic equilibrium. It becomes more crucial when it comes to experimental studies of soft condensed matter systems [64,65] most notably due to their applicability to overdamped systems whose instantaneous momenta may not be accessible. It was den Otter and Clarke that first investigated the deviations of the kinetic and configurational temperatures from the system temperature in DPD system [37]. Since then, little attention has been paid to the configurational temperature in DPD simulations until Allen recently argued that the configurational temperature should be added to the list of diagnostic tests applied to DPD simulations [53]. In addition, we also calculate the average potential energy U and the radial distribution function (RDF) g(r) [3,4], both of which are very important configurational quantities in simulations.

Dynamical
To have a deeper understanding of how the physical system evolves, it is not enough to just evaluate the static quantities described above. It is essential to measure and compare the relevant dynamical properties. In general, in simulation, one relies on various Green-Kubo formulas [66] to calculate various transport coefficients; in this article we restrict ourselves to two special cases.
The velocity autocorrelation function (VAF) is an important measure of dynamical fidelity that numerical methods should be able to reproduce, particularly if they are to be used for nonequilibrium (transient-phase simulation). The VAF, which characterizes the translational diffusion of the system, is defined as where v i (0) is the initial velocity picked up after the system is well equilibrated. By integrating the VAF (Green-Kubo relation [67,68]), we can compute the translational diffusion coefficient, which can also be obtained by using Einstein's relation [69] giving the diffusion coefficient as the slope of the mean square displacement as a function of time t.
To investigate the rotational relaxation process of the system, we measure the transverse momentum autocorrelation function (TMAF) [70,71], which is related to the shear viscosity, η, in the hydrodynamic limit, i.e. small wavenumber k and large time t, and is defined by where ρ is the particle density, m is the particle mass and p x (k, t) represents the transverse component of the momentum, where p j,x denotes the x-component of the momentum of particle j and similarly q j,z represents the z-component of its position (see more details in [70]). Note that, i here is the imaginary unit in complex number.

Results
The kinetic temperature control of various methods was tested and shown in Fig. 1. According to the black dashed second order line in the figure, all the methods seem to have second order convergence to the invariant measure, which verifies the  error analysis results on nonsymmetric DPD-S1, LA and Peters thermostats in Section 4, but is a bit surprising for DPD-VV, NHLA, PNH and PNHL-N methods that were also based on nonsymmetric splittings.
The performance of standard DPD methods (DPD-VV, DPD-S1 and DPD-Trotter) and the Peters thermostat that converges to the standard DPD system in the limit of h → 0, are almost indistinguishable with the tendency to quickly blow up after the stepsize reaches h = 0.1. Both LA and NHLA thermostats show similar behaviors and maximal stepsizes that can be used are limited around h = 0.1. The latter displays a better accuracy when stepsize is large due to the additional Nosé-Hoover-like thermostat (44). The PNH thermostat illustrates the worst kinetic temperature accuracy and we can also see the clear low stability threshold for the PNH thermostat due to the lack of ergodicity. Most surprisingly and interestingly, the PNHL-S and PNHL-N methods, based on different splitting orders on the same pairwise Nosé-Hoover-Langevin (PNHL) thermostat (50), show dramatically different kinetic temperature control: the symmetric PNHL-S method is worse than most of the methods, whereas the nonsymmetric PNHL-N method outstandingly maintains almost two orders of magnitude improvement on the accuracy of kinetic temperature than all the other methods. Configurational quantities, configurational temperature and average potential energy, were compared in Fig. 2. Again, all the methods seem to have second order convergence to the invariant measure except the PNH thermostat, which displays a stability threshold of h = 0.05. Most of the methods, including standard DPD methods, LA and Peters thermostats, are indistinguishable and cross the 100% barrier in configurational temperature and 10% barrier in average potential energy respectively around h = 0.1. As in the case of kinetic temperature, the performance of the NHLA thermostat on those configurational quantities is slightly better than the majority due to the additional thermostat. Unlike the kinetic temperature case, the PNHL-S method does have very good accuracy in configurational quantities with almost one order of magnitude improvement than the majority. Incredibly, the PNHL-N method manages more than one order of magnitude accuracy enhancement in configurational temperature and almost two orders of magnitude in average potential energy. Fig. 3 compares the radial distribution function (RDF) that characterizes the structure of the fluids. The DPD-S1 and PNHL-N methods were used for the standard DPD and PNHL systems respectively. Given that very small stepsizes were used, all the methods display exactly the same RDF, which indicates that different systems maintain the same structure of the fluids without the impacts of discretization errors. Expectedly, discretization errors start to destroy the structure  methods with three different values of friction coefficient were calculated by using DPD-S1 method to compare with other methods. 100 different runs were averaged to reduce the sampling errors after the system was well equilibrated. of the fluids with increasing stepsizes as shown in Fig. 4. Standard DPD methods and Peters thermostat exhibit similar behaviors with the RDFs starting to show artifacts around h = 0.09 and being heavily destroyed around h = 0.13. The LA and PNH thermostats again show lower stability in maintaining the correct structure, blowing up around h = 0.11 and 0.07, respectively. The performance of PNHL-S and NHLA methods are slightly better than the majority of the schemes, while the PNHL-N method only starts to show pronounced artifacts above h = 0.25, which is remarkably more than two times larger than the stepsize usable in the standard DPD methods.
The velocity autocorrelation function (VAF) of various numerical methods were calculated in Fig. 5 to compare with standard DPD methods with three different values of friction coefficient. The DPD-S1 method was used to calculate the VAF of standard DPD methods since there is no difference with other two if very small stepsizes were used to reduce the effects of discretization errors. Similarly, the uncorrupted dynamics of PNHL-S and PNHL-N methods should be exactly the same, and the latter was used. The VAF of the PNH and PNHL systems are indistinguishable and consistent with standard DPD methods in the regime of small friction coefficient γ = 0.5. This is not surprising since the average of the dynamical variable ξ , which is controlled by an additional thermostat and can be thought of as the "dynamical friction", tends to zero.
As we expected, the Peters thermostat shares the same VAF as standard DPD methods in the regime of the commonly-used friction coefficient γ = 4.5 if the same value of γ was chosen. The stochastic randomization frequency Γ = 0.44 was used in LA and NHLA thermostats to maintain similar translational diffusion properties as in standard DPD system with γ = 4.5. The VAF of standard DPD system with large friction coefficient γ = 40.5 is also shown in the figure, which indicates that the larger the friction is the faster the VAF goes down (the system loses memory faster).
The effects of the discretization error on each method were also investigated in Fig. 6. The results are largely consistent with those observed for the configurational temperature. The only surprise is that the PNHL-S method allows a maximum stepsize which is similar to those of the thermostats considered (and much lower than the useful stepsize for PNHL-N).
Moreover, we observe that the PNH thermostat begins to display non-physical artifacts, at a stepsize of just h = 0.05.
Likewise the LA thermostat has a lower stability threshold. Among the various schemes, PNHL-N is again by far the most stable scheme, exhibiting only a mild deviation from the reference VAF at h = 0.17.
One of the most important features of DPD simulations is the correct handling of rotational relaxation which will be important for resolving correct vortical motion and long-ranged interactions. We therefore investigate the computation of transverse momentum autocorrelation function (TMAF, see Section 5.2.2) for each scheme. The results here are in many ways similar to those obtained for the velocity autocorrelation function (Fig. 5), thus only the results of DPD and PNHL methods are shown in Fig. 7: the PNHL method is expectedly consistent with small friction (γ = 0.5) in standard DPD simulation, and only a minor difference is observed between the PNHL method and standard DPD with a moderate value of friction (γ = 4.5). One may notice that the TMAF shown in Fig. 7 is not as smooth as the VAF in Fig. 5 even by averaging 1000 times more different runs. We emphasize here that collective (N-particle) correlations, such as TMAF and stress autocorrelation function [12], fluctuate rapidly and thus are always determined with poorer statistics than single-particle ones, such as VAF.
We emphasize the following facts regarding our numerical tests: • A detailed statistical analysis of the results presented has not been incorporated in this article, due to the extensive computational requirements of doing so. An early study [17] has already suggested that highly reliable estimation of the kinetic temperature can be typically obtained by various methods in the stepsize regime of our interests, i.e. h ≥ 0.05.
In terms of convergence rate of thermodynamic properties to distribution, all the methods perform similarly in practice (see Figs. 1-2). PNHL methods with stepsize h = 0.05. DPD-S1 and PNHL-N were used to solve the DPD and PNHL systems respectively. The ratio of the curves on the right panel is proportional to the corresponding shear viscosity. The wavenumber k was chosen as 2π /L and 100,000 different runs were averaged to reduce the sampling errors after the system was well equilibrated.

Table 2
Comparisons of the computational efficiency of the various numerical methods. "Critical stepsize" is the stepsize beyond which the numerical method starts to show pronounced artifacts, while "maximal stepsize" is the stepsize stability threshold. The "numerical efficiency" of each method was scaled to that of the benchmark DPD-VV method.  • Based on the results for the VAF and TMAF, we have observed very good agreement of the dynamical properties between PNHL and DPD with relatively small friction coefficient, particularly γ = 0.5. Although the PNHL method may not be able to fully recover the hydrodynamics as DPD, we have already seen that PNHL offers substantially improved stability in simulations. Moreover, the PNHL method, as a general momentum-conserving thermostat, has a valid motivation in a broader context than just comparison with DPD; in particular it will be useful in NEMD and in other cases where fluid dynamics per se is not at issue.

Computational efficiency
The computational efficiency of the various methods was tested with simulation details in Section 5.1. All the tests were run on an HP Z600 Workstation with 15.7 GB RAM. As shown in Table 2, we calculated the CPU time (milliseconds) taken with the use of Verlet neighbor lists for the integration of a single time step of h = 0.05 (averaged over 10 000 consecutive time steps). Note that only the time for the main integration step (without calculating any physical quantity) was counted.
In order to quantitatively compare the overall performance of each method, we define a quantity, the "numerical efficiency", that measures the amount of simulation time accessible per unit of computational work, i.e.
Numerical Efficiency = Critical Stepsize

CPU Time Per Step
, where the "critical stepsize" is defined as the stepsize beyond which pronounced artifacts become apparent in some observable thus rendering the simulation unusable. For practical purposes, we use a threshold of 10% relative error in configurational temperature in determining the critical stepsize. Our reasoning in choosing the configurational temperature as the quantity for determination of the critical stepsize is that (a) it is a sensitive observable and difficult to control in simulation, (b) accuracy of configurational temperature seems to us to be intuitively to be a core requirement for canonical simulation methods, (c) good control of the configurational temperature appears to imply good results in other stationary computations. The choice of 10% as the bar for accuracy is clearly arbitrary. If a smaller error threshold were used, the results may change slightly, but the ordering of the methods in terms of efficiency would remain essentially the same.
We also show the "maximal stepsize" in Table 2, which is the stepsize at which the numerical method is either above 100% relative error in configurational temperature or unstable. Only the PNHL-N method needs to calculate the force twice in each integration step, which is the reason why an almost doubled time was needed than other methods. Finally, the computed "numerical efficiency" of each method was scaled to that of the DPD-VV method since it is still the most popular method in DPD simulations and we use that method as the benchmark in the comparisons.
As can be seen from the table, those standard DPD methods and the Peters thermostat have comparable "numerical efficiency". The LA and PNH thermostats are slightly better than the commonly-used DPD-VV method with around 15% and 27% improvements, respectively. Both NHLA and PNHL-S methods maintain about 50% enhancement. The improvement of the "numerical efficiency" of the PNHL-N is remarkably 87% in comparison to the DPD-VV method. Similarly, we may also compute the enhancements of the PNHL-N method based on kinetic temperature and average potential energy, using the same critical stepsize h = 0.05 as in configurational temperature: 43% for the former and 120% for the latter.

Conclusions
We have reviewed a number of numerical methods that are widely used in DPD simulations and have also proposed a new stochastic momentum-conserving thermostat, pairwise Nosé-Hoover-Langevin (PNHL) thermostat. Two favorable splitting methods of the PNHL thermostat were introduced and compared with existing methods in the computation of various physical quantities (both static and dynamical).
We have observed that, for the PNHL thermostat proposed here, the PNHL-N method based on a nonsymmetric arrangement of the terms of a splitting, gives an enormous stability benefit. The PNHL-N method needs to calculate the force twice in each integration step, which is computationally costly in the model setting used in this article; nevertheless when the computational overhead is costed carefully the PNHL-N method outperforms the alternatives. To measure the practical performance of numerical methods in DPD simulations quantitatively, we have defined the "numerical efficiency", based on which we have reported substantially improvements of both methods of the newly proposed thermostat, with the symmetric PNHL-S method 46% more efficient than the commonly-used DPD-VV method and the nonsymmetric PNHL-N method incredibly 87% better than the benchmark method in DPD simulations. It should be noted that, if the force calculation is not that expensive in other model settings, the gain of using the PNHL-N method could be further exploited.
Based on the numerical experiments of the velocity and transverse momentum autocorrelation functions which characterize the translational and rotational diffusions of the system respectively, DPD and PNHL give rather similar dynamical properties in practice. Although PNHL as formulated is not based on a hydrodynamic interaction model [72,73], we have seen that it is an effective replacement for DPD in the low-friction regime. Moreover, we point out that the projection along interacting particle pairs in PNHL could be replaced by alternatives to achieve further control of transport properties. The method is also potentially useful more broadly in molecular simulation applications, whenever momentum conservation is at issue.
The only difference between PNH and PNHL is the additional Langevin thermostat acting on the dynamical friction, however they perform very differently. This is probably because the cutoffs used in PNH and the simple potentials there provided insufficient internal mechanisms for the system to achieve an ergodic sampling of the canonical distribution. In molecular dynamics with steep pair potentials (e.g. Lennard-Jones), the ergodic properties develop more naturally and the "L" in PNHL can be redundant in some instances. It is also worth mentioning that both PNH and PNHL are not able to vary the Schmidt number since the average of the dynamical friction tends to zero, whereas the Schmidt number can be tuned in other schemes.
We have also investigated the order of convergence of the long-time averages to the invariant measure for a couple of methods described in this article. By extending the framework recently introduced in Langevin dynamics, we can infer (and verify using numerics) the second order convergence for those nonsymmetric GLA-like methods (DPD-S1, LA and Peters thermostats). However, rigorous investigation on other nonsymmetric methods (DPD-VV, NHLA, PNH and PNHL-N methods) that surprisingly obtained second order convergence remains to be established. Overall, we claim here that PNHL thermostat indeed can be used (and may be preferred in some typical cases) as an alternative to low-friction DPD simulations with substantially improved computational efficiency and no degradation of convergence rate.

Appendix A. Order of long-time error
To investigate the order of convergence of the invariant distribution, we define the operator L † associated to propagation under the numerical method with stepsize h for an SDE. It can be thought of as a perturbation of the exact Fokker-Planck for some perturbation operators L † i .
We also view the invariant distribution ρ associated to the numerical method as a perturbation of the target canonical distribution ρ β for some correction function f i satisfying f i = 0. Substituting L † and ρ into the stationary Fokker-Planck equation Since the exact operator preserves the canonical distribution, i.e. L † Exact ρ β = 0, we obtain where [X, Y ] = X Y − Y X is the commutator. Therefore, a symmetric splitting (such as DPD-Trotter) guarantees at least second order convergence automatically whereas a nonsymmetric one could only typically gives first order. In what follows we show that the nonsymmetric schemes widely used in DPD simulations provide second order convergence due to cancellations in the error expansions.
In the example of DPD-S1, which can be termed as OBAB and clearly is not symmetric, each interacting pair preserves Thus, by applying the BCH expansion on the operator of the DPD-S1 scheme (24), we obtain L † DPD-S1 = L † Hence, which gives the only solution of the PDE (A.3) to the DPD-S1 scheme Given that higher order perturbations in (A.2) are not generally equal to zero, we have shown that the nonsymmetric DPD-S1 scheme has second order convergence to its invariant distribution. Similarly, we can also demonstrate that both the Lowe-Andersen and Peters thermostats (in the fashion of BABO) maintain second order convergence to the invariant distribution. Since all the three methods involve a second-order symplectic Verlet method for the deterministic part and other actions only on the Ornstein-Uhlenbeck (OU) process, we may refer to them as generalized Geometric Langevin Algorithms of order two (GLA-2) [14] in the context of stochastic momentum-conserving thermostats. It should be emphasized that the OU process in stochastic momentum-conserving thermostats is pairwise and thus different from the standard setting of Langevin dynamics. Nevertheless, it has been demonstrated that second order is still achieved by the combination of a symplectic method for the deterministic part and an exact solve for the OU process. More discussions regarding the long-time accuracy of the GLA-type methods in Langevin dynamics can be found in [14,60,74].

Appendix B. Integration schemes
We list detailed integration steps for each method described in the article here. Verlet neighbor lists [30] are used throughout each method as long as possible. , p) and F R i (q) are conservative, dissipative and random forces, respectively, in standard DPD system. Note that, at the end of each integration step, the dissipative forces F D i (q n+1 , p n+1 ) are further updated by using the up-to-date velocities (momenta).

B.3. DPD-Trotter scheme: DPD-Trotter
For each interacting pair within cutoff radius (q ij < r c ) p n+1/2 i = p n i + hm ij v ij q n , p n q n ij /2, p n+1/2 j = p n j − hm ij v ij q n , p n q n ij /2, where τ = γ ω D /m ij , m ij = m i m j /(m i + m j ) and R ij (t) are normally distributed variables with zero mean and unit variance.

B.5. Peters scheme II: Peters
For each particle i p n+1/3 i = p n i + hF C i q n /2, For each interacting pair within cutoff radius (q ij < r c ) where γ and σ are the same dissipation and random coefficients respectively as in standard DPD system.
where α is is a coupling parameter chosen as 0.9/(ρh), and ρ is the particle density.
where T k is the momentary kinetic temperature and T 0 is the desired temperature. For the remaining ( P ) fraction interacting pairs within cutoff radius (q ij < r c )