Pairwise adaptive thermostats for improved accuracy and stability in dissipative particle dynamics

We examine the formulation and numerical treatment of dissipative particle dynamics (DPD) and momentum-conserving molecular dynamics. We show that it is possible to improve both the accuracy and the stability of DPD by employing a pairwise adaptive Langevin thermostat that precisely matches the dynamical characteristics of DPD simulations (e.g., autocorrelation functions) while automatically correcting thermodynamic averages using a negative feedback loop. In the low friction regime, it is possible to replace DPD by a simpler momentum-conserving variant of the Nos\'{e}--Hoover--Langevin method based on thermostatting only pairwise interactions; we show that this method has an extra order of accuracy for an important class of observables (a superconvergence result), while also allowing larger timesteps than alternatives. All the methods mentioned in the article are easily implemented. Numerical experiments are performed in both equilibrium and nonequilibrium settings; using Lees--Edwards boundary conditions to induce shear flow.


Introduction
Classical molecular dynamics (MD), where the motion of individual atoms is governed by Newton's law in the microcanonical ensemble (where energy, i.e., the Hamiltonian of the system, is conserved), has been widely used in molecular simulations [5,30]. However, the constant energy setting is not relevant to a real-world laboratory environment since energy, as an extensive variable, depends on system size. In typical cases, one replaces the microcanonical ensemble by the canonical one, where temperature is conserved using suitable "thermostat" techniques.
One popular thermostat is Langevin dynamics, whereby each particle is subject to dissipative and collisional interactions with the particles of an artificial "heat bath" and modeled by supplementing the conservative Newtonian equations of motion with balanced damping and stochastic terms in such a way that the desired target system temperature is maintained. However, as pointed in [35], in order to be consistent with hydrodynamics, a particle model should respect Galilean invariance, and, in particular, should conserve momentum, something that Langevin dynamics fails to do. Fundamentally, Langevin dynamics and its overdamped Brownian dynamics limit, are appropriate for modeling systems in or near thermodynamic equilibrium and therefore do not take into account the possibility of an underlying fluid flow, thereby precluding their use in situations where the flow of the soft matter system cannot be predicted beforehand (e.g., when dealing with interfaces or nonuniform flow). Moreover, it has been reported in [22] that, due to the violation of global momentum conservation, Langevin dynamics can lead to nonphysical screening of hydrodynamic interactions with a screening length proportional to the inverse square root of the friction coefficient of the algorithm. In order to control simulation artifacts, one is led to use a very small friction coefficient, effectively reducing Langevin dynamics to Hamiltonian dynamics in the microcanonical ensemble. Therefore, when hydrodynamics is of interest, standard thermostats should be replaced by momentum-conserving thermostats, in particular the so-called dissipative particle dynamics (DPD) method of Hoogerbrugge and Koelman [35].
DPD was first proposed in order to recover the properties of isotropy (and Galilean invariance) that were broken in the so-called lattice-gas automata method [31]. In DPD, each body is regarded as a coarse-grained particle. These particles interact in a soft (and short-ranged) potential, allowing larger integration timesteps than would be possible in MD, while simultaneously decreasing the number of degrees of freedom required. As in Langevin dynamics, a thermostat consisting of well-balanced damping and stochastic terms is applied to each particle. However, unlike in Langevin dynamics, both terms are pairwise and the damping term is based on relative velocities, leading to the conservation of both the angular momentum and the linear momentum. The property of Galilean invariance (i.e., the dependence on the relative velocity) makes DPD a profile-unbiased thermostat (PUT) [25,26] by construction and thus it is an ideal thermostat for nonequilibrium molecular dynamics (NEMD) [75]. The momentum is expected to propagate locally (while global momentum is conserved) and thus the correct hydrodynamics is expected in DPD [75], as demonstrated previously in [23]. Due to the aforementioned properties, DPD has been widely used to recover thermodynamic, dynamical, and rheological properties of complex fluids, with applications in polymer solutions [78], colloidal suspensions [63], multiphase flows [62], and biological systems [52]. DPD has been compared with Langevin dynamics for out-of-equilibrium simulations of polymeric systems in [64], where as expected the correct dynamic fluctuations of the polymers were obtained with the former but not with the latter.
Given its promising prospects from the applications perspective, and its widespread use in large scale simulations, the optimal design of numerical methods for DPD becomes crucially important, in particular the numerical efficiency in practice [8,14,60,83]. Numerous numerical schemes [4,8,16,55,65,72,76] have been proposed in the last two decades following the introduction of DPD, which are intended to reduce nonphysical artifacts (especially in the large stepsize regime) induced by the discretization error. Recently, we have systematically examined the performance (in terms of accuracy, efficiency, and robustness) of a number of the most popular methods in the literature [50].
In addition, we have proposed in [50] an alternative stochastic momentum-conserving thermostat, the pairwise Nosé-Hoover-Langevin (PNHL) thermostat. This method mimics the DPD system in the regime of low friction, however achieving much higher accuracy and computational efficiency. One contribution of the current article is a perturbation analysis showing that averages of observables of a certain (common) form performed using a nonsymmetric splitting of the PNHL system (i.e., the PNHL-N method [50]) have unexpected second order accuracy (as a power of the stepsize), justifying the enhanced performance of PNHL observed in simulations.
The second important contribution of this article is a new pairwise adaptive Langevin (PAdL) thermostat to replace DPD in the regime of moderate or high friction. This method draws on work on adaptive thermostats [38,51,71], by supplementing a DPD type pairwise stochastic perturbation by an auxiliary control law (also pairwise) to maintain the thermodynamic state. The new method fully captures the dynamics of DPD (for example, autocorrelation functions in DPD are precisely reproduced) and thus can be directly applied in the setting of momentum-conserving simulations as a replacement for DPD. We describe a simple splitting-based numerical method for PAdL. While PAdL has similar per-timestep computational cost, the method is shown to generate substantially more accurate approximations to thermodynamic averages at the same stepsize as DPD (as much as an order of magnitude). Moreover, and perhaps more significantly, the stepsize can be increased by around 50% using PAdL, for similar accuracy, resulting in a much more efficient overall simulation method.
Furthermore, we discuss the proper treatment of Lees-Edwards boundary conditions in the DPD setting, an essential part of modeling shear flow.
The rest of the article is organized as follows. In Section 2, we review the formulation of DPD and the momentum-conserving PNHL method, and introduce the newly proposed PAdL thermostat that mimics the dynamics of DPD. We investigate in Section 3 numerical methods for PNHL and PAdL and give results on the order of accuracy of various schemes, in particular showing that the PNHL-N method is second order in its approximation of ergodic (long time) averages of a certain class of observables. Section 4 presents numerical experiments in both equilibrium and nonequilibrium cases, comparing the performance of numerous popular numerical methods in practice. The proper treatment of Lees-Edwards boundary conditions in the context of momentum-conserving thermostats is also discussed in Section 4. Our findings are summarized in Section 5.

Dissipative particle dynamics and pairwise thermostats
In this section, we review the formulation of DPD and the momentum-conserving PNHL thermostat, followed by the introduction of the PAdL thermostat.

Dissipative particle dynamics (DPD)
The original DPD system was updated in discrete time steps and was later reformulated by Español and Warren [24] as a system of Itō stochastic differential equations (SDEs).
We write the DPD system in a compact (vector) form: where q and p are dN -dimensional vectors, d being the underlying dimensionality of the physical space (typically d = 3) and N being the number of particles, representing positions and momenta of particles, respectively, M is the diagonal mass matrix, −∇U is the conservative force given in terms of a potential energy function U = U (q), γ is the friction coefficient and σ represents the strength of the random forces, W is a vector of S = dN (N − 1)/2 independent Wiener processes (note that, starting from the first element, every d consecutive elements are identical, intended for each interacting pair, and the symmetry of dW ij = dW ji is required to ensure the momentum conservation), and the matrices Γ(q) ∈ R dN ×dN and Σ(q) ∈ R dN ×S satisfy the following relation: which can be thought of as a generalized fluctuation-dissipation relation. The matrix Γ(q) may be written explicitly as where is the weight function defined in the DPD system and E ij =q ijq T ij , wherê q ij = (q i − q j )/r ij is the unit vector with r ij = q i − q j being the distance between two particles, is the d by d projection matrix on particles i and j.
In the standard treatment of DPD, the conservative forces are derived from a sum of pair potentials a typical choice of ϕ being [33] ϕ(r ij ) = where parameters a ij represent the maximum repulsion strengths between interacting pairs, and r c denotes a cutoff radius that is used in order to reduce the computational cost by restricting the number of pairs that need to be involved in the force computation. In more recent studies, the conservative force is frequently obtained by coarse-graining procedures (e.g., [44,53,57]) and may be represented by tabled data or an interpolation thereof [77].
The weight function can be arbitrarily chosen without violating the thermal equilibrium. A simple, popular choice is

Statistical properties of the DPD system in equilibrium
Let H represent the system Hamiltonian and assume the following fluctuation-dissipation relation: where k B is the Boltzmann constant and T the temperature, it is easy to show that the DPD system (1) preserves the momentum-constrained canonical ensemble with density where Z is a suitable normalizing constant (the partition function), β −1 = k B T , and π = (π x , π y , π z ) is the linear momentum vector. Additional constraints should be included if the angular momentum is also conserved. In an open system model, DPD conserves both angular and linear momenta (due to the fact that the interactions between particles depend on relative velocities) thus DPD is an isotropic Galilean-invariant thermostat which also preserves hydrodynamics [4,58,75]. However, it should be pointed out that if periodic boundary conditions are used (even along one direction), the conservation of angular momentum is lost, which is also expected when Lees-Edwards boundary conditions are applied. Due to the use of coarse-grained potentials, and depending on the choice of weight functions, the ergodicity of stochastic thermostats for mesoscale modeling cannot be taken for granted. The ergodicity of the DPD system has been demonstrated only in the case of high particle density in one dimension by Shardlow and Yan [73]. It is also worth noting that the simulation study of Pastorino et al. [64] appears to contradict ergodicity when (i) soft DPD potentials are used and (ii) the number of interactions is limited sufficiently.

Pairwise Nosé-Hoover-Langevin (PNHL) thermostat
We have recently proposed the PNHL thermostat [50] in an attempt to improve the accuracy and stability in DPD simulations. Comparing to the standard formulation of DPD, in PNHL, the stochastic term was completely removed and the constant friction coefficient in the damping term was replaced by an additional dynamical variable that was driven by the difference between the instantaneous temperature based on relative velocities and the target temperature. Moreover, a Langevin thermostat was applied to the auxiliary variable in order to enhance the ergodicity. As in DPD, the PNHL thermostat conserves the momentum and is Galilean-invariant, thus correct hydrodynamics is expected to be preserved.
The equations of motion of the PNHL thermostat are given by where ξ is an auxiliary dynamical variable and G(q, p) denotes the accumulated deviation of the instantaneous temperature away from the target temperature where µ is a coupling parameter which is referred to as the "thermal mass", coefficient constantsγ andσ satisfy the following fluctuation-dissipation relation: and W = W(t) is a standard Wiener process. The special caseγ =σ = 0 of the PNHL thermostat reduces the system to the pairwise Nosé-Hoover (PNH) thermostat by Allen and Schmid [4]. However, as mentioned above, ergodicity is not expected in some coarse-grained models, even with the addition of stochastic forces, thus, in particular, the PNH thermostat is likely to fail for many choices of weight functions and potentials. The inclusion of noise through the auxiliary variable has been rigorously shown to restore ergodicity to the system, albeit only in the case of a nonpairwise scheme [49].
The PNHL thermostat (10) preserves the canonical ensemble with a modified invariant distribution (comparing to ρ β (9)) in the standard DPD system We conjecture that the PNHL system is ergodic for this distribution if the weight functions have a sufficiently large support. Two different splitting methods have been proposed in [50] for the PNHL thermostat, with the first being of a symmetric manner, labeled as PNHL-S, and the other nonsymmetric, PNHL-N (see details in Appendix A). Both PNHL methods have been compared to various popular schemes for a number of physical quantities in [50], and it turns out that both methods (especially the PNHL-N method) achieve significant improvements in terms of accuracy, robustness, and numerical efficiency over alternatives. Both mathematical theory and numerical experiment with methods for Langevin dynamics (see, e.g., [48]) have repeatedly shown the efficiency advantage of symmetric numerical methods, thus the numerical results reported in [50] were, until now, a curious anomaly. In the next section of this article, we finally demonstrate that PNHL-N has a superconvergence property (an extra order of accuracy) for averages of an important class of observables, thus explaining its numerical performance in practice.
The dynamical properties of the PNHL formulation correspond to those of the standard DPD system only in the low friction regime. In what follows, we propose a new momentumconserving thermostat in order to have full control of the dynamics.

Thermodynamically-corrected and hydrodynamics-preserving pairwise adaptive Langevin (PAdL) thermostat
In order to maximize numerical efficiency (especially for large scale simulations), a large timestep is always preferred to discretize the system of interest. However, as already mentioned, a large timestep can result in pronounced nonphysical artifacts. A recent article [74] interprets the effect of the discretization error in Langevin dynamics as a means of driving the system away from a desired invariant distribution-excess energy is pumped into the system in the form of "shadow work" to prevent it from maintaining the desired temperature.
Furthermore, if the forces are computed "on-the-fly" (see, e.g., [54]), they are likely subject to substantial errors, which would also effectively "heat" the system up (see more discussions in [51]). As an illustration, the hybrid quantum mechanics/molecular mechanics (QM/MM) method introduces localized heating due to the force-mixing at the boundary of the coupled QM and MM regions [59]. Moreover, the use of external fields in nonequilibrium models causes viscous heating (i.e., the energy pumped into the system leads to a temperature rise under steady perturbation), thus, in those cases, proper thermostats are required to drain the excess energy in order to maintain the correct system temperature [7]. One way to regulate excess heat in mechanical systems is via the use of negative feedback loop controls. Nosé-Hoover [36,61] is one such feedback control system, but the observation of [38] is that feedback loop controls can be introduced in tandem with a Langevin thermostat (their so-called "adaptive" Langevin dynamics), benefiting from the strong ergodicity properties of Langevin dynamics together with an automatic regulation of the kinetic energy. These methods were further explored in [51,71]. Notably, in [51] it was shown that an adaptive Langevin device could be used to dissipate excess heat due to noisy forces, while providing statistical convergence properties very similar to those of Langevin dynamics.
It is natural to consider extending the adaptive thermostat idea to DPD, thus hopefully providing correction for thermodynamic observables while mimicking, to a large extent, the properties of DPD. To this end, we propose here a momentum-conserving PAdL thermostat, whose equations of motion are given by where σ is a constant coefficient as in DPD and the function G is given in (11). An additional Langevin thermostat could also be added on the auxiliary variable ξ as in PNHL, but it is not necessary here because of the presence of the additional stochastic term directly contacting the physical variables.
It can be demonstrated, modifying the argument of [38], that the PAdL system preserves the canonical ensemble with a modified invariant distributioñ whereγ can be thought of as the "effective friction" and the following fluctuation-dissipation relation is satisfied as in DPD: The invariant distribution (15) implies that the auxiliary variable ξ is Gaussian distributed with meanγ and variance β −1 µ −1 . The auxiliary variable will fluctuate around its mean value during simulation and we can tune the value of the effective friction in order to recover the dynamics of DPD in a wide range of friction regimes. Therefore, we can think of the PAdL thermostat as the standard DPD system with an adaptive friction coefficient. It can be seen that the PAdL thermostat inherits essential properties of DPD (such as Galilean invariance and momentum conservation) required for consistent hydrodynamics. In the large thermal mass limit (i.e., µ → ∞), the PAdL thermostat effectively reduces to the standard DPD formulation (1).

Numerical methods for pairwise thermostats
A great deal of effort has been devoted to developing accurate and efficient numerical methods to solve DPD and related systems. We have compared a number of popular schemes in [50].
In what follows we briefly review the numerical methods for the PNHL thermostat that will be included for further investigation in this article, followed by the derivation of the numerical scheme for the newly proposed PAdL thermostat. We also discuss the accuracy of equilibrium averages for those methods, in particular, we show the unexpected second order convergence of the PNHL-N method for certain observables.

Numerical Treatment of PNHL
Since we are unable to solve the PNHL system (10) "exactly", we instead split the vector field of the system into pieces (A, B, C, D, and O) in such a way that we can solve each subsystem exactly (see more details in [50]). We can write down the generators associated with each piece, respectively, as the sum of which is the generator of the PNHL system: The flow map (or phase space propagator) of the system is given by where the exponential map represents the solution operator. Various approximations of F t can be obtained as products (taken in different arrangements) of exponentials of the splitting terms, however, it turns out in a number of studies [45-47, 50, 51] that different splittings and/or combinations give dramatically different performance in practice. Two different splitting methods of PNHL have been proposed in [50], termed PNHL-S and PNHL-N, respectively: and The detailed integration steps of both methods can be found in Appendix A. It is important to note that the steplengths associated with various operations are uniform and span the interval h. Thus the O step in either of the two methods is taken with a steplength of h, while the other pieces with a steplength of h/2. Note also that these two integrators differ only from the order of integrating the last two pieces. However, an additional force calculation is required (using the updated positions) in the PNHL-N scheme just before updating the last B piece at the end of each integration step. The additional force calculation could be costly, however, it was found to be offset by a great increase in accuracy and usable steplength [50].

Numerical Treatment of PAdL
As in PNHL schemes, we separate the vector field of the PAdL system (14) into pieces, which we label as A, B, O, and D, respectively: Note that the generators of pieces A, B, and D here are exactly the same as defined in PNHL (Eqs. (18), (19), and (21), respectively), while the generator for the remaining piece is given by where : denotes the Frobenius product for matrices, i.e., A : B = Tr(A T B). Note that the generator O here is different from that of PNHL, although they both represent Ornstein-Uhlenbeck processes. Overall, the generator for the PAdL system can be written as It is straightforward to note that pieces A, B, and D can be solved exactly as in PNHL schemes. In fact, each interacting pair in O is also exactly solvable (in the sense of distributional fidelity) [70]. For each interacting pair, subtracting dv j , where v j is the velocity of particle j (with corresponding mass m j ), from dv i and multiplying byq ij on both sides yields where m ij = m i m j /(m i +m j ) is the "reduced mass" and v ij =q ij ·v ij . With ξ being a positive constant, the above equation is a standard Ornstein-Uhlenbeck process with the exact (in the sense of distributions) solution [41] v whereτ = ξω D /m ij , v ij (0) are the initial relative velocities, and R ij are normally distributed variables with zero mean and unit variance. As demonstrated in [51], the exact solution above (31) is still valid for ξ < 0. When ξ = 0, one can simply replace (1 − e −2τ t )/(2ξm ij ) by its well-defined asymptotic limit, in which case (31) becomes Then the velocity increments can be obtained as and subsequently the corresponding momenta can be updated by which defines the propagator e hL O i,j for each interacting pair.
Various splitting methods could be constructed, we propose in this article a symmetric PAdL method (alternatively, the "ABODOBA" method), whose propagator can be written as where the O part associated with interacting pairs is given by The detailed integration steps can be found in Appendix A.

Accuracy of equilibrium averages
We have demonstrated in [50] that a symmetric splitting method gives at least second order convergence to the invariant measure, and similar observations would hold for the symmetric methods described above, PNHL-S (25) and the PAdL scheme (36). Superconvergence properties (i.e., fourth order convergence) have also been proven in both Langevin dynamics [45] and adaptive Langevin thermostat [51] for configurational sampling. Nonsymmetric splitting methods of geometric Langevin algorithms [11] type can exhibit high order ergodic approximations, but first order convergence is generally expected for other nonsymmetric splitting methods. This raises the obvious question of why we observed second order accuracy using the PNHL-N method (26) in simulations performed in [50]. The framework of long-time Talay-Tubaro expansion has been widely used in analyzing the accuracy of ergodic averages (i.e., averages with respect to the invariant measure) in stochastic numerical methods [1, 2, 17, 45-48, 51, 79]. Therefore, in this section, we adopt the procedures to verify the second order convergence of the PNHL-N method for certain observables.
For a splitting method described by L = L α +L β +· · ·+L ζ , its associated effective operator L † with stepsize h is given by which can be computed using the Baker-Campbell-Hausdorff (BCH) expansion and can thus be viewed as a perturbation of the exact Fokker-Planck operator L † : for some perturbation operators L † i . We also define the invariant distributionρ associated with the numerical method as an approximation of the target invariant distribution ρ β : for some correction functions f i satisfying f i = 0, where · denotes the average with respect to the target invariant distribution. Thus, substitutingL † andρ into the stationary Fokker-Planck equation gives Since the exact Fokker-Planck operator preserves the target invariant distribution, i.e., L † ρ β = 0, we obtain by equating first order terms in h.
We are able to find the perturbation operator L † 1 by using the BCH expansion for any particular integration scheme. Then we can calculate its action on ρ β . However, it is generally very hard to solve the above PDE (43) in order to obtain the leading correction function f 1 (see examples in Langevin dynamics [45]). Note that for symmetric splitting methods, including the PAdL scheme (36), f 1 is always equal to zero, thereby ensuring second order convergence to the invariant measure [51].
Consider now the PNHL-N method (26), which can be written as where and By using the BCH expansion, we obtain where the notation [A, B] = AB − BA denotes the commutator of operators A and B, and subsequentlyL † Thus the leading perturbation operator of the PNHL-N scheme is whose action on the invariant distribution of the PNHL system (13) reads (assuming M = I for simplicity) Recall the Fokker-Planck operator of the PNHL system: Although in this case the right-hand side of the PDE (43) is relatively simple, it is still very challenging to solve the PDE in order to obtain the corresponding leading correction function f 1,PNHL−N . However, the additional variable ξ is normally distributed with mean zero and variance β −1 µ −1 . Thus, the variance of ξ will be small if the thermal mass µ is large. Therefore, following [51], we consider projecting the Fokker-Planck equation and its solution by integrating with respect to the Gaussian distribution of ξ in the ergodic limit. That is, we apply the projection operator [32] Pν(q, p, ξ) := Ω ξρ β (q, p, ξ)ν(q, p, ξ) dξ where ν is an arbitrary function, to the PDE (43). Effectively, this leads to the reduced where the operatorĽ † , acting on functions of q and p only, is just the operator L † reduced by the action of the projection. In fact, now L † PNHL is simply reduced tǒ while the right-hand side is unchanged due to the fact that ξ is not present there. Finally, we can solve the reduced PDE to obtain the leading correction function: which leads to the following proposition: Proposition 1. For sufficiently large thermal mass, the PNHL-N method (26) exhibits second order convergence for thermodynamic averages of certain observables, i.e., in the form of φ(q, p) = p 2k ϑ(q), where k is an integer and ϑ(q) can be constant.
The class of observables includes kinetic temperature and observables that only depend on the configurations. In other words, for those observables, assuming the asymptotic expansion holds, the computed average (in the large thermal mass limit) reads which is of order two. This is fully consistent with what we have observed numerically for a number of observables in [50] and which we further verify in the following section.

Numerical experiments
In this section, we conduct various numerical experiments, in both equilibrium and nonequilibrium regimes, to compare the newly proposed PAdL method with a number of alternative momentum-conserving schemes described in [50].

Simulation details
We adopted the simulation details used in [50], a standard parameter set commonly used in algorithms tests [33,60,65,70,72,76,83]. Specifically, a system of N = 500 identical particles with unity mass was simulated in a cubic box (length L = 5) with periodic boundary conditions, unless otherwise stated. Particle density ρ d = 4 was used with cutoff radius r c = 1 and k B T = 1. The initial positions of the particles were independent and identically distributed (i.i.d.) with a uniform distribution over the box, while the initial momenta were i.i.d. normal random variables with mean zero and variance k B T . The potential repulsion parameters a ij were set to be 25, while a wide range of (effective) friction coefficients (0.5, 4.5, and 40.5) were used. Verlet neighbor lists [84] were used in each method. In our simulations, the thermal mass in PAdL was chosen to be the same as that of PNHL (whereγ = 4.5 was used), i.e., µ = 10. When comparing different formulations, we have to make sure that similar translational diffusion properties of the fluid were obtained. For the PAdL thermostat, we can always tune the value of σ so that the same (effective) friction coefficient as in DPD was obtained, i.e.,γ = γ.  average potential energy (right) against stepsize by using various numerical methods with (effective) friction coefficient γ = 4.5. The system was simulated for 1000 reduced time units but only the last 80% of the data were collected to calculate the static quantity to make sure the system was well equilibrated. Ten different runs were averaged to further reduce the sampling errors. The stepsizes tested began at h = 0.05 and were increased incrementally by 15% until all methods either started to show significant relative error (100% in configurational temperature or 10% in average potential energy) or became unstable.
We have observed in [50] that standard DPD methods (including the DPD velocity Verlet method (DPD-VV) [8], Shardlow's S1 splitting method (DPD-S1) [72], and the DPD Trotter scheme (DPD-Trotter) [70]) and the Peters thermostat [65] give almost indistinguishable performance in all the quantities that we have tested. Therefore, the DPD-S1 method was used to represent the standard DPD formulation, unless otherwise stated. The PNH thermostat [4] was not included for further comparison because of its stability issue at relatively small stepsize, which is well documented in [50]. Since the dynamics of the PNHL thermostat is consistent with that of DPD in the low friction regime, both PNHL methods were compared to alternatives only in the case of γ = 0.5.
As in [50], we measured the "numerical efficiency", defined as the ratio of the "critical stepsize" and CPU time per step, of each method and then scaled it to the benchmark DPD-VV method, which is widely used in popular software packages. The CPU time (in milliseconds) for the main integration step (without calculating any physical quantity) is the time taken (on an HP Z600 Workstation with 15.7 GB RAM) with the use of Verlet neighbor lists for the integration of a single time step of h = 0.05 (averaged over 10000 consecutive time steps). The critical stepsize was chosen as the stepsize corresponding to 10% relative error in computed configurational temperature [3,12,69,81], an observable function of positions whose average in the canonical ensemble is precisely the target temperature:  Table 1: Comparisons of the computational efficiency of the various numerical methods in the moderate (effective) friction regime of γ = 4.5. "Critical stepsize" is the stepsize beyond which the numerical method starts to show pronounced artifacts (10% relative error in computed configurational temperature), while "maximal stepsize" is the stepsize stability threshold above which the method is unstable. The "numerical efficiency" of each method was scaled to that of the benchmark DPD-VV method. The efficiency figures above quantify the computational work in terms of the number of force computations and correctly take into account the fact that all the methods require one force evaluation per timestep with the sole exception of PNHL-N (included only in the low friction case, see Table 2) which requires two force evaluations.
where ∇ i U and ∇ 2 i U denote the gradient and Laplacian of the potential energy U with respect to the position of particle i, respectively (see more discussions in [50]). We mention that this criterion was not precisely adopted in a recent article by Farago and Grønbech-Jensen [27], where they determined the critical stepsize based on the accuracy control of average potential energy but compared their results to those obtained from configurational temperature in [50]. We would also like to point out that the Nosé-Hoover device employed in PAdL and PNHL is straightforward to implement. In MD simulations, the choice of the thermal mass µ can be a technical impediment to using Nosé-Hoover controls, but we have observed in [49] and also in the experiments of this article that the addition of stochastic noise changes the nature of the coupling parameter and the results obtained are relatively insensitive to its selection.

Equilibrium
Configurational quantities, such as configurational temperature and average potential energy, were compared in Figure 1 with (effective) friction coefficient γ = 4.5. According to the dashed order line, we can see that all the methods tested exhibit second order convergence to the invariant measure for both quantities. More specifically, DPD and the LA thermostat show rather similar behavior, while the NHLA thermostat is slightly better than those two. Quite remarkably, the newly proposed PAdL method (36) achieves one order of magnitude improvement over DPD in terms of numerical accuracy for a fixed stepsize. For certain accuracy (i.e., a fixed relative error), the PAdL method can use doubled stepsize, thus substantially improving the "numerical efficiency" defined in [50] (100% efficiency improvement = double the performance). Our observations were confirmed in Table 1, which shows that the PAdL method indeed has a more than 130% improvement in numerical efficiency over the DPD method. The results on the configurational temperature and average potential energy are rather similar, therefore in what follows we present only configurational temperature results.
We also explore in Figure 2

Method
Critical  system with a fixed value ofγ = 4.5. All the methods compared clearly display second order convergence to the invariant measure, with ABODOBA and BADODAB achieving one order of magnitude improvement in accuracy compared to the other methods. However, it should be noted that in ABODOBA, we needed to update the Verlet neighbor lists only once at each step since the update of A at the end was directly followed by another update of A at the beginning of next step, whereas one additional update of Verlet neighbor lists was required in BADODAB. In our simulation code, the cost of such an update was essentially the same as one force calculation, thus the ABODOBA method (36) was used throughout the current article for the PAdL system. This again illustrates the importance of optimal design of numerical methods. Figure 3 (left) compares the configurational temperature control of various methods in both low and high friction regimes. In the low friction regime, where both PNHL methods were included, again all the methods exhibit second order convergence to the invariant measure.   The NHLA, PNHL-S, and PAdL methods are rather similar to each other, all of which are superior to both DPD and LA methods. The PNHL-N method achieves more than one order of magnitude improvement in numerical accuracy over the DPD method. Although the PNHL-N method requires two force calculations at each step, it still achieves a more than 80% improvement as shown in Table 2. The table also reveals that the PAdL method has an almost 50% improvement in performance compared to the DPD method.

Method
In the high friction regime, the behavior of those methods are rather different from that of other regimes. As shown in Figure 3 (right), surprisingly the most popular DPD-VV method is slightly worse than other standard DPD methods. Both LA and NHLA are indistinguishable from the DPD method. Superconvergence property (i.e., fourth order convergence to the invariant measure) was not observed for the PAdL method in this regime. Nevertheless, the PAdL method still obtains a dramatic improvement over all other schemes. Table 3 shows that the PAdL method has a more than 190% improvement in overall numerical efficiency over the benchmark DPD method.
The control of the dynamical properties of the PAdL method was also tested and plotted in Figure 4. In particular, we compared two important quantities: the velocity autocorrelation The DPD-S1 method was used for DPD with a small stepsize of h = 0.01, while h = 0.05 was used for PAdL. 100 and 100,000 different runs were averaged in the cases of VAF and TMAF (the wavenumber was chosen as 2π/L), respectively, to reduce the sampling errors after the system was well equilibrated.
function (VAF) and the transverse momentum autocorrelation function (TMAF) [34,80], which characterize the translational and rotational diffusions of the system, respectively. The integral of the VAF is related to the diffusion constant, while the logarithm of the TMAF is proportional to the shear viscosity in the hydrodynamic limit (see more discussions in [50]). Unlike the PNHL formulation, which corresponds to the low friction regime, the PAdL system is able to capture the dynamics of DPD in a wide range of friction coefficients-the relaxations of the VAF and the TMAF of both formulations are indeed indistinguishable (only visible with the help of the on-center symbols).

Nonequilibrium
It is well known that nonequilibrium methods, where a steady state is maintained under external perturbations (either stationary fluxes or spatial gradients of some quantities), are more efficient means than equilibrium autocorrelation functions for extracting transport coefficients (e.g., rheological properties such as shear stress and shear viscosity) from fluid dynamics simulations (see more discussions in [5]). Thus, there has been a rapidly growing interest in NEMD [25]. For instance, planar Couette flow, where a simple and steady shear flow is generated, is commonly employed as a numerical "viscometer" in particle-based methods to obtain transport coefficients [5,25]. Furthermore, due to the fact that it is relevant to many real-life phenomena as well as for its simplicity, planar Couette flow has been widely adopted in laboratory experiments. In order to measure the shear viscosity in lab experiments, a linear profile is often imposed at a fixed shear rate and then the resultant shear stress can be measured. However, in computer simulations, simple periodic boundary conditions are unable to maintain a steady linear velocity profile, resulting in problems at the boundaries of the simulation domain. There have been early attempts to generate momentum or energy flows in MD simulations where particles are made to interact with momentum or energy reservoirs (e.g., a velocity profile can be obtained by fixing the average velocity in the extremal slabs of a fluid) [6,15,39,82]. However, these methods are not compatible with periodic boundary conditions, and thus lead to surface effects. Alternatively and more appealingly, one can apply Lees-Edwards boundary conditions (LEBC) [43] to retain periodicity but alter the position and velocity of the periodic images. In this case, a simple shear flow (with a constant shear rate) is generated, which allows the investigation of the dependence of the viscosity on the shear rate [29]. LEBC has been extensively studied in DPD and related systems to study rheological behavior in colloidal suspensions [9,10], polymeric systems [28,64,75], as well as multiphase systems [62] (see more discussions on boundary conditions in DPD in [66][67][68]).
Before we analyze the numerical results obtained by simulating various methods under LEBC (i.e., the system is perturbed by a simple shear flow), we first briefly review LEBC and then discuss two important issues in NEMD: (1) the practical implementation of LEBC in DPD and related momentum-conserving systems, where forces are dependent on relative velocities; (2) the practical measurement of system temperature in NEMD.

Lees-Edwards boundary conditions (LEBC)
In order to generate a simple shear flow in NEMD, the periodic boundary conditions (PBC) have to be modified. A common way to achieve that is to apply the well-known LEBC [43]. In LEBC, the primary cubic box (with lengths L x , L y , and L z ) remains centered at the origin, however, a uniform shear velocity profile is expected [25] u =γye x , where e x is the unit vector in the x-direction andγ is the shear rate defined aṡ LEBC is also called the "sliding brick" boundary conditions. It is important to note that LEBC is applied only in the x-direction, while the other directions (y and z) remain with PBC. Special attention has to be paid in LEBC when a particle is crossing the boundary in the y-direction. In this case, one of the images of the crossing particle will enter through the opposite face, but with both position and velocity modified in a proper way because of the streaming velocity (58).
The periodic boundary crossing is now handled as follows [5]: where N L is the "rounded" number of layers (boxes) between the current position of particle i in the y-direction and the origin, and ∆q x is the displacement of the upper layer during the The minimum image convention should now proceed as follows [5]: Note that when the shear rate is zero, i.e.,γ = 0, LEBC reduces to PBC.

LEBC in pairwise thermostats
A recent article of Chatterjee [13] claimed that, owing to the dependence of inter-particle forces on the relative velocities of the particles, it is problematic to directly apply LEBC to DPD, especially in the high friction regime. As shown in Figure 5 (left), where exactly the same setting as in the original article [13] was used, as the friction increases, the velocity profile starts to (significantly) deviate away from the target linear profile. A simple remedy is to switch off the interactions of dissipative and random forces (i.e., the DPD thermostat) if one particle is within interacting range of an image of other particle near the boundaries where adjacent layers have different streaming velocities (i.e., the y-direction), as proposed in [13]. However, the finding of [13] directly contradicts the principle of LEBC, which is translationally invariant and is intended to overcome the surface effects. In fact, as pointed in [25], in no way can the particles actually sense the boundaries of any given box since the system is spatially homogeneous. Furthermore, our numerical experiments, which are in perfect agreement with theoretically expected behavior as shown in the right panel of Figure 5 even in the high friction regime, suggest that LEBC might have been incorrectly implemented in [13]. One possibility is that when calculating the relative velocity between one particle and an image of another, which is in a layer with different streaming velocity from its origin, the original velocity, rather than the "modified velocity" due to the different streaming velocities in different layers, was used as the velocity of the image particle. By neglecting the necessary modification, we obtained the left panel of Figure 5, while if the velocity of the image particle was properly modified the expected linear velocity profile was recovered in Figure 5 (right) using otherwise exactly the same setting.
Overall, it should be emphasized that if one particle is interacting with the image of another under certain conditions, the relative velocity (in the x-direction) between them should be modified as follows: where function "fabs(·)" returns the absolute value of the argument. Note that: (1) q y ij in (73) has to be evaluated either before the minimum image convention (67)- (71) or by other proper ways to determine the actual distance between two interacting particles; (2) (74) is the relative ("absolute") velocity between the two particles and one has to take into account the effects of the streaming velocity as indicated.
We suspect that the necessary modification (72)- (74) was not correctly implemented in [13], resulting in the nonphysical behavior as shown in Figure 5 (left). It is not surprising that by switching off the DPD thermostat on interactions that cross certain boundaries would recover the expected linear velocity profile as shown in [13], since it directly avoids the situation described in (72)- (74) where special attention has to be paid. Overall, the "workaround" does not provide any physical explanation, and could affect the underlying dynamics of the system, implying that it should be abandoned.

Temperature in nonequilibrium molecular dynamics
Another interesting question in NEMD is that "What is the most appropriate way to measure the system temperature?" Since the streaming velocity should be subtracted from the particle velocity before calculating the kinetic temperature, the normal definition should be modified as follows [25]: where N d is the number of degrees of freedom of the system and u denotes the corresponding streaming velocity at the location of particle i (58). In other words, if u = 0, (75) reduces to the standard form.
If we can assume that the velocity profile is linear, as in uniform shear flow, we can just calculate and subsequently subtract it. However, as pointed in [25,26], at higher shear rates and/or Reynolds numbers (i.e., the ratio of the inertial and viscous forces), the assumption of a linear streaming velocity profile is extremely dubious, even though Lees-Edwards boundary conditions are used. This issue was addressed in [25,26]

Method
Critical  proposed. The PUT allows the simulation itself define the local streaming velocity (see more details in [25,26,85]). However, this is still not completely satisfactory since PUT assumes that the streaming velocity profile is stationary in time, whereas the profile could vary in time.
In DPD and related systems, temperature calculations can be based on relative velocities (e.g., the NHLA thermostat [76]), which do not rely on the relationship between the absolute particle velocity and an underlying streaming velocity. As discussed in [50], it is more desirable, especially in NEMD, to define the temperature solely based on the configurations, which leads to the configurational temperature defined in (57) (see applications in [18,19]). In addition, thermostats based on the configurational temperature have been widely used in NEMD with shear flows [20,21,56]. Therefore, the configurational temperature formulation (57) is used in our numerical experiments. Figure 6 compares the configurational temperature control of various systems described under LEBC with different shear rates. As can be shown from the figure, when the shear rate is relatively small (γ = 0.2, which is larger than that of Figure 5), the behavior is largely similar to that of Figure 1 (left): all the methods appear to show second order convergence to the invariant measure, and the newly proposed PAdL method achieves one order of magnitude improvement in numerical accuracy over both DPD and LA methods, both of which are slightly worse than the NHLA method. The overall numerical efficiency was compared in Table 4. Again the PAdL method is by far the most efficient method of all, which has an about 130% improvement over the benchmark DPD-VV method.

Results
When the shear rate is relatively high (γ = 2), as shown in Figure 6 (right), all the methods appear to lose the clear second order convergence previously observed. The LA thermostat, displaying large relative error even when the stepsize is relatively small, appears to be most vulnerable to the high shear rate, with the DPD method being slightly better. While exhibiting some unexpected behavior, both NHLA and PAdL have better numerical accuracy control than the other two. It can be seen that, as the stepsize increases, the relative error of the PAdL method starts to decrease before growing as would be expected. We believe this unexpected decrease at the beginning is due to sampling errors, which is more likely to be observed in the high accuracy regime, since increasing the system size (i.e., number of particles) would resolve this issue. Nevertheless, the PAdL method consistently achieves an order of magnitude improvement in numerical accuracy over the DPD method.

Conclusions
The most significant contributions of this article are as follows: (i) the introduction of a new negative-feedback-loop-controlled formulation of DPD, i.e., the PAdL thermostat, that preserves the dynamical behavior of DPD while enhancing the accuracy of averages, and (ii) an analysis explaining the second order convergence of the nonsymmetric PNHL method [50] for a broad class of observables. We have, additionally, clarified the treatment of Lees-Edwards boundary conditions for modeling shear flow. PAdL (and PNHL, in the low friction regime) effectively double the performance of DPD, are easy to implement, and rest on a solid foundation of theory and previous numerical experience in the Langevin dynamics context.
PAdL and PNHL use auxiliary dynamical variables (like the "global demon" of Nosé-Hoover dynamics), whose interaction with the physical variables is controlled by a coupling parameter. In the deterministic case, the choice of the coupling coefficient (or "thermal mass") can be an added complication. However, our previous experience with the Nosé-Hoover-Langevin method (as well as in the numerical experiments of this article) has demonstrated that, once stochastic noise is present, the character of the parameter dramatically changes, so that the method properties are robust to a much wider range of values.
The key idea that is exploited here in the formulation of PAdL is that the thermostat can correct for numerical errors. This builds on the observation of [74] that the discretization error can be interpreted as irreversible work that induces effective heating. What is less obvious is that the negative feedback loop controller used here, which introduces additional dynamics that also must be discretized, does not distort the invariant distribution, but we surmise that this is a consequence of the simple form of this equation and its simple discretization which maintains the form of a discrete control law.
The usefulness of the PAdL scheme will be felt most strongly in simulations that have a strong hydrodynamic character, e.g., nonequilibrium systems undergoing shear flow. In other cases, for example, models in the solid or gas states, one would expect that Langevin dynamics or the PNHL method described here will be of more value. A secondary benefit of the PAdL thermostat which we have not explored in this article, is its ability to correct thermodynamics automatically for errors through the computation of the conservative force.
DPD and PAdL would likely have similar ergodic properties which depend on the choice of weight functions and potentials. It will be of interest to explore this issue further in future studies. These choices will also influence the accuracy and stability of the numerical methods and thus the potential efficiency improvements available. Exploration is needed of the numerical behavior when the potential energies arise from tabled or interpolated data.
It should be noted that the idea of "transverse" DPD [40,53] (i.e., including the transverse component of dissipative and random forces) could easily be incorporated into PAdL. The parallelization of the proposed method of PAdL is not completely trivial, but is similar to the task of parallelizing Shardlow-like schemes. The problem has recently been addressed by Larentzos et al. [42].
For each interacting pair within cutoff radius (r ij < r c ),

Pairwise adaptive Langevin thermostat: PAdL
For each particle i, q n+1/2 i = q n i + hm −1 i p n i /2 , p n+1/4 i = p n i + hF C i (q n+1/2 )/2 . For each interacting pair within cutoff radius (r ij < r c ), whereτ = ξω D (r ij )/m ij and R ij are normally distributed variables with zero mean and unit variance; else: ∆v ij = σ ω R (r ij )/m ij h/2R ij .