1 Introduction

The development of the computational capacities in the last decades has allowed physicists to use numerical simulations to study physical properties at the atomic scale with the help of statistical physics. In particular, molecular dynamics (MD) consists in integrating the equation of motions for the atoms in order to sample probability measures in a high dimensional space^[1-3]. However, traditional microscopic methods suffer from limitations in terms of accessible time and length scales, which drives the development of mesoscopic coarse-grained methods. These mesoscopic models aim at greatly reducing the number of degrees of freedom explicitly described, and thus the computational cost, while retaining some properties absent from more macroscopic models such as hydrodynamics. Smoothed dissipative particle dynamics (SDPD)^[4] belongs to this class of the mesoscopic method. It couples a particle Lagrangian discretization of the Navier-Stokes, smoothed particle hydrodynamics (SPH)^[5-6], and the thermal fluctuations from models like dissipative particle dynamics with energy (DPDE) conservation^[7-8]. It is thus able to deal with hydrodynamics at nanoscale and has been shown to give results consistent with MD for a wide range of resolutions, at equilibrium and for shock waves^[9], or for dynamical properties such as the diffusion coefficient of a colloid in an SDPD bath^[10-11]. SDPD has in particular been used to study colloids^[10-12], polymer suspensions^[13], and fluid mixtures^[14].

One of the main challenges for mesoscopic models incorporating fluctuations is to develop efficient, stable, and parallelizable numerical schemes for the integration of their stochastic dynamics. Most schemes are based on a splitting strategy^[15-16] where the Hamiltonian part is integrated through a velocity Verlet scheme^[17]. A traditional and popular algorithm first proposed for dissipative particle dynamics^[18] and later extended to DPDE^[19] relies on a pairwise treatment of the fluctuation/dissipation part^[20]. The adaptation of this scheme to dynamics preserving various invariants has led to a class of schemes called Shardlow-like splitting algorithms (SSA)^[21]. A major drawback in this strategy is the complexity of its parallelization^[22]. Other schemes have been recently proposed in Ref. [23] to enhance its use in parallel simulations.

All these schemes are, however, hindered by instabilities when internal energies become negative. This especially happens at low temperatures or when small heat capacities are considered, typically for small mesoparticles. It has been proposed to use Monte Carlo principles to sample the invariant measure of DPDE, by resampling the velocities along the lines of centers according to a Maxwell-Boltzmann distribution and redistributing the energy variation into internal energies according to some prescriptions^[24]. However, this approach leads to dynamics which is not consistent with DPDE. It was proposed in Ref. [25] to correct discretization schemes for DPDE by rejecting unlikely or forbidden moves through a Metropolis procedure, which prevents the appearance of negative internal energies and improves the stability of the integration schemes.

There exist relatively few references in the literature about the integration of full SDPD. Most works focus on numerical schemes in the isothermal setting^[26], avoiding the need to preserve the total energy during the simulation. In a previous article^[9], we introduced an adaptation of the Shardlow splitting to SDPD, allowing a good control of the energy conservation. The aim of this work is to provide more details about the possible integration of SDPD in an energy conserving framework and most importantly to increase the stability for small particle sizes by adapting the metropolization procedure described in Ref. [25].

This article is organized as follows. We first present in Section 2 the equations of SDPD as reformulated in Ref. [9]. In Section 3, we recall the Shardlow splitting for SDPD and introduce a Metropolis step to enhance the stability of the algorithm. We evaluate the properties of the Shardlow and Metropolis schemes by means of numerical simulations in Section 4. Our conclusions are gathered in Section 5.

2 SDPD

At the hydrodynamic scale, dynamics of the fluid is governed by the Navier-Stokes equations as below, which reads in their Lagrangian form when the heat conduction is neglected (for the time t≥0 and the position x in a domain Ω⊂³),

(1)

In these equations, the material derivative used in the Lagrangian description is defined as

The unknowns are the density of the fluid ρ(t, x) ∈ ₊, its velocity v(t, x) ∈ ³, its internal energy u(t, x) ∈ , and the stress tensor σ(t, x) ∈ ^3×3,

(2)

where P is the pressure of the fluid, η is the shear viscosity, ζ is the bulk viscosity, and I is the identity matrix.

In the following, we first present the SPH discretization of the Navier-Stokes equations in Subsection 2.1 before introducing the SDPD equations reformulated in terms of internal energies^[9] in Subsection 2.2.

2.1 SPH

Hydrodynamics^[5-6] is a Lagrangian discretization of the Navier-Stokes equations (1), where a finite number N of fluid particles play the role of interpolation nodes. These fluid particles are associated with a portion of the fluid of mass m. They are located at the position q_i ∈ Ω and have a momentum p_i ∈³. The internal degrees of freedom are represented by an internal energy ε_i. In general, the energies are bounded below. Upon shifting the minimum of the internal energy, we may consider that the internal energies remain positive (ε_i>0).

2.1.1 Approximation of field variables and their gradients

A key ingredient in the SPH discretization is the use of a particle-based interpolation of the field variables. This leads to an approximation of the field variables by averaging over their values at the particle positions weighted by a smoothing kernel function W. The kernel is generally required to be non-negative, regular, normalized as and with finite support^[27]. We introduce the smoothing length h defined such that W(r)=0 if |r| ≥ h. In the sequel, we use the notation r = |r|. In this work, we rely on a cubic spline^[28], whose expressions read

(3)

The field variables are then approximated as

(4)

where f_i denotes the value of the field f on the particle i.

The approximation of the gradient ∇_xf is obtained by deriving (4), which yields

In order to have more explicit expressions, we introduce the function F such that ∇_r W(r) = -F(|r|)r. In the case of the cubic spline (3), it reads

The gradient approximation can then be rewritten as

To simplify the notation, we define the following quantities for two particles i and j:

We can associate a density ρ_i and a volume to each particle as

(5)

The corresponding approximations of the density gradient evaluated at the particle points read

(6)

2.1.2 Thermodynamic closure

As in the Navier-Stokes equations, an equation of state is required to close the set of equations provided by the SPH discretization. This equation of state relates the entropy S_i of the mesoparticle i with its density ρ_i(q) (as defined by (5)) and its internal energy ε_i through an entropy function,

(7)

The equation of the state can be computed by microscopic simulations or by an analytic expression modeling the material behavior. Using the partial derivatives ∂_ε and ∂_ρ of the entropy with respect to the energy and the density, it is possible to assign to each particle the temperature

the pressure

and the heat capacity at the constant volume

To simplify the notation, we omit in Subsection 2.2 the dependence of T_i, P_i, and C_i on the variables ε_i and q.

2.2 Equations of motion for SDPD

SDPD^[4] is a top-down mesoscopic method relying on the SPH discretization of the Navier-Stokes equations with the addition of thermal fluctuations which are modeled by a stochastic force. In its energy reformulation^[9], SDPD is a set of stochastic differential equations for the following variables: the positions q_i∈Ω⊂³, the momenta p_i∈³, and the energies ε_i∈ for i=1, 2, ..., N.

The dynamics can be split into several elementary dynamics, the first one being conservative dynamics deriving from the pressure part of the stress tensor (2) (see Subsection 2.2.1) and a set of pairwise fluctuation and dissipation dynamics stemming from the viscous terms in (2) coupled with random fluctuations (see Subsection 2.2.2).

2.2.1 Conservative forces

The elementary force between particles i and j arising from the discretization of the pressure gradient in the Navier-Stokes momentum equation reads

(8)

In its original formulation^[4], this conservative dynamics clearly appears as the Hamiltonian dynamics with a potential relating the energy with the positions and the particle entropy S_i. The entropies are then invariants of this subdynamics. In the energy reformulation, the entropies are no longer considered as such. Instead, the focus is on the total energy,

which is preserved by the dynamics. This can be ensured by computing the variation of the particle volume as

leading to the variation of the internal energy given by

This allows us to write the conservative part of the dynamics as

(9)

This dynamics preserves by construction the total momentum and the total energy E(q, p, ε).

2.2.2 Fluctuation and dissipation

In order to give the expression of the viscous and fluctuating part of the dynamics, we define the relative velocity for a pair in particles i and j as

The viscous term in the Navier-Stokes equations (1) is discretized by a pairwise dissipative force, while the thermal fluctuations are modeled by a pairwise stochastic force. In the spirit of DPDE, the pairwise fluctuation and dissipation dynamics for i < j can be chosen as follows:

(10)

where B_ij is a three-dimensional vector of standard Brownian motion, and Γ_ij and Γ_ij are 3×3 symmetric matrices. By construction, (10) preserves the total momentum in the system since dp_i = -dp_j. Furthermore, as in DPDE, the equations for the energy variables are determined to ensure the conservation of the total energy E(q, p, ε). As , Itȏ calculus yields the resulting equations in (10).

We consider friction and fluctuation matrices of the form,

(11)

with the projection matrices and given by

Introducing the coefficients

which are defined from the fluid viscosities η and ζ appearing in the stress tensor (2), we can choose the friction and fluctuations coefficients as

(12)

As shown in Ref. [9], this ensures that measures of the form (with g a given smooth function)

(13)

are left invariant by the elementary dynamics (10). Alternative fluctuation/dissipation relations are possible (such as constant parameters σ), but the relations (12) allow us to retrieve the original SPDP^[4].

2.2.3 Complete equations of motion

By gathering all the terms, the SDPD equations of motion reformulated in the position, momentum, and internal energy variables read

(14)

with Σ_ij and Γ_ij given by (11) and (12). The dynamics (14) preserves the total momentum and the total energy E(q, p, ε) since all the elementary subdynamics can ensure these conservations.

2.3 Reduced units for SDPD

In SDPD, the mass m of the fluid particles allows us to change the resolution of the method. We introduce the particle size , where m₀ is the mass of one microscopic particle (e.g., a molecule). Since we deal with different particle sizes in the following, it is convenient to introduce reduced units for each size K as follows:

(15)

where is the mass unit, is the length unit, is the energy unit, and ρ is the average density of the fluid. With such a set of reduced units, the time unit is

In the following, we select time steps before expressing them in terms of , where K is the particle size used in the simulations. This explains the use of non-round time steps in Section 4.

The smoothing length h_K defining the cut-off radius in (3) also needs to be adapted to the size of the SDPD particles so that the approximations (4) continue to make sense. In order to keep the average number of neighbors roughly constant in the smoothing sum, h_K should be rescaled as

In this work, we take h=2.5, which corresponds to a typical number of from 60 neighbors to 70 neighbors, a commonly accepted number^[28].

3 Integration schemes

In the following, we describe several numerical schemes for the integration of SDPD. They all rely on a splitting strategy^[15-16], where the full dynamics can be divided in simpler elementary dynamics that can be consecutively integrated. Since the conservative part of the dynamics (9) can be viewed as Hamiltonian dynamics, it is natural to resort to a symplectic scheme such as the widely used velocity-Verlet scheme^[17] which ensures a good energy conservation in the long term^[29-30]. This algorithm is briefly described in Subsection 3.1.

There is, however, no definite way to deal with the fluctuation/dissipation part described in Subsection 2.2.2. Due to its close similarities with DPDE, we propose in the following to adapt some schemes devoted to the integration of DPDE to the SDPD setting. One approach to integrate SDPD, described in Ref. [9], is based on the algorithm proposed by Shardlow^[20] for dissipative particle dynamics and its subsequent adaptations to DPDE^[19]. The dynamics can be split into a Hamiltonian part, discretized through a velocity-Verlet algorithm (16), and elementary pairwise fluctuation/dissipation dynamics that can be successively integrated. We first recall in Subsection 3.2 the SSA used in Ref. [9]. While this provides a way to integrate SDPD preserving its invariants (approximately for the energy), it suffers from stability issues especially for small particle sizes, when the internal and kinetic energies are of the same scale. We thus explore methods to improve the stability of these integration algorithms in Subsection 3.3, relying on the ideas developed in Ref. [25] where a Metropolis acceptance-rejection step is included to correct the biases of the numerical discretization of the fluctuation/dissipation part.

3.1 Integrating Hamiltonian part of dynamics

It is convenient to consider the conservative part of the dynamics (9) in its original formulation in the position, momentum, and entropy variables^[4] in order to take advantage of the conservation of the entropy S_i. The internal energies are related to the positions and entropies by an energy function , which allows us to write the Hamiltonian as

The dynamics of (9) can thus be recast in a Hamiltonian form as

The velocity-Verlet scheme^[17] allows us to integrate such dynamics while preserving on average the Hamiltonian H. This corresponds to the following integration scheme:

(16)

3.2 SSA

We present here a first possibility for the integration of the fluctuation/dissipation dynamics introduced in Ref. [9] based on existing schemes for DPD^[20] and DPDE^[19]. If we neglect the dependence of Γ and Σ on ε_i, the elementary dynamics (10) on the momenta can be viewed as a standard Ornstein-Uhlenbeck process and solved analytically. We provided in Ref. [9] the corresponding expression for the updated momenta after a time step Δt as

(17)

where G_ijⁿ is a standard three-dimensional Gaussian variable and for

The integration of the momenta with (17) induces a variation of the kinetic energy which is then redistributed symmetrically in the internal energies as suggested in Refs. [19] and ^[31]. This guarantees the exact conservation of the energy during this elementary step. The new internal energies are finally given by

Thermodynamic variables like the temperatures T_i, T_j and heat capacities C_i, C_j are updated with the equation of state using the new internal energies, before turning to another pair of particles.

Let us, however, remark that the pairwise Shardlow-like algorithm is sequential by nature, and its parallelization requires a convoluted method^[22]. Moreover, and maybe more importantly, there is no mechanism preventing the apparition of negative energies during the simulation. This situation happens when the fluctuations are large with respect to the internal energies, typically at low temperatures or when the particle sizes are small (so that their heat capacities are small as well). This leads to stability issues unless very small time steps are used.

3.3 Metropolized integration scheme

To avoid instabilities related to negative internal energies while keeping reasonable time steps, it has been proposed to include a Metropolis step to reject impossible or unlikely moves^[25]. In the following, we show how this procedure can be used for SDPD. First, we reformulate the pairwise dynamics (10) as overdamped Langevin dynamics in the relative velocity v_ij variable only (see Subsection 3.3.1). We then construct proposed moves for the metropolized scheme and compute the corresponding acceptance ratio in Subsection 3.3.2. A simplified version of the metropolized scheme is introduced in Subsection 3.3.3, where the computation of the Metropolis ratio is avoided, and the rejection occurs only to avoid negative internal energies.

3.3.1 Reformulation of fluctuation and dissipation dynamics as overdamped Langevin dynamics

In order to simplify the metropolization of the integration scheme, we show that the elementary fluctuation/dissipation dynamics (10) can be described only in terms of the relative velocity v_ij and formulated as overdamped Langevin dynamics.

Since the dynamics (10) preserves the momentum p_i+p_j, the momenta p_i and p_j can be rewritten as a function of v_ij,

This has already shown how to express the momenta p_i and p_j in terms of v_ij. In addition, the kinetic energy formulated in the relative velocity reads

The conservation of the energy and the fact that dε_i = dε_j provide the expression of the internal energies as a function of the relative velocity,

(18)

Using this relation, the dynamics (10) can in fact be rewritten as effective dynamics on the relative velocity only,

(19)

where Γ_ij and Σ_ij are functions of the relative velocity through (18). We claim that the dynamics (19) can be written more explicitly as overdamped Langevin dynamics under the form of

(20)

with the diffusion matrix

and the potential

where

Let us emphasize that the reformulation (20) is the key element for the Metropolis stabilization. We now check that (20) holds. By definition,

It therefore suffices to check that

We first compute the gradient of the potential ,

Upon application of the matrix M,

The divergence of M with respect to the relative velocity reads

The desired result follows from

In view of (20), it can be immediately deduced that the following measure on the relative velocity, at fixed momenta p_i, p_j and fixed internal energies ε_i, ε_j, is left invariant by the overdamped Langevin dynamics (19):

(21)

3.3.2 Metropolis ratio

We consider (17) as the proposed move. In terms of the relative velocity, it reads

(22)

with

and

The momenta and internal energies can then be updated as

(23)

The internal energies T_i, T_j and heat capacities C_i, C_j are updated accordingly.

In order to decide whether we update the configuration with the proposed move or keep the current one, we first check whether ε_iⁿ⁺¹ and ε_jⁿ⁺¹ are negative, in which case the proposal is rejected. Otherwise, we compute a Metropolis ratio that is an acceptance probability. The probability to accept the proposed move from v to v' is min(1, A_Δt(v, v')) with

where T_Δt is the transition kernel associated with the proposal. In the following, we omit all the dependence on the position qⁿ, which remains constant in this subdynamics, to simplify the notation.

The probability that (22) proposes v' starting from v is given by

with the inverse matrix and the determinant . For the direct move, the transition probability simply reads

while for the reverse move,

Use (21) with the reference taken at the iteration n,

Finally, the acceptance ratio is given by

(24)

with

Starting from a configuration (p_iⁿ, p_jⁿ, ε_iⁿ, ε_jⁿ), the overall algorithm (exact Metropolis scheme (EMS)) to integrate the fluctuation/dissipation for a pair of particles (i, j) is organized as follows:

(ⅰ) Compute a proposed move for v_ijⁿ⁺¹ with (22).

(ⅱ) If the following energy bound does not hold:

(25)

the move will be rejected, that is, (p_iⁿ⁺¹, p_jⁿ⁺¹, ε_iⁿ⁺¹, ε_jⁿ⁺¹)=(p_iⁿ, p_jⁿ, ε_iⁿ, ε_jⁿ).

If the bound is satisfied, the algorithm continues.

(ⅲ) Compute the acceptance ratio with (24).

(ⅳ) Draw U_ijⁿ ~ [0, 1] and compare it with A_Δt(v_ijⁿ, v_ijⁿ⁺¹). If U_ijⁿ > A_Δt(v_ijⁿ, v_ijⁿ⁺¹), the move will be rejected, and (p_iⁿ⁺¹, p_jⁿ⁺¹, ε_iⁿ⁺¹, ε_jⁿ⁺¹)=(p_iⁿ, p_jⁿ, ε_iⁿ, ε_jⁿ).

Otherwise, it is accepted, and the momenta and internal energies are updated with (23), along with the internal temperatures T_i, T_j and heat capacities C_i, C_j.

3.3.3 Approximate metropolized scheme (AMS)

Since the computation of the Metropolis ratio may be cumbersome in practical simulations, we propose a simplified and approximate scheme where we only reject moves that cause internal energies to become negative. It avoids the need to actually compute the Metropolis acceptance ratio.

As for the complete metropolized scheme, we use the expression (22) as the proposed evolution for the relative velocities. We then check whether the updated internal energies remain positive and reject the moves that do not satisfy this property. The current configuration at the time is then used as the new configuration and counted as usual in the averages. Otherwise the move is accepted, and the velocities and internal energies, along with the internal temperatures and heat capacities, are updated accordingly. When no stability issues, i.e., negative internal energies, appear, the AMS is equivalent to SSA.

4 Numerical results

In the following, we test the accuracy of our schemes for the ideal gas equation of state given by

(26)

The interest of this model is that the marginal distribution for the internal energies ε_i has an analytic expression,

(27)

where is the heat capacity in the equation of state (26), and Γ is the gamma function. This distribution is plotted in Fig. 1 for various particle sizes. As the size K decreases, very small internal energies become more likely, and stability issues arise.

Our simulations have been carried out on a three-dimensional system of 1 000 particles initialized on a simple cubic lattice with an initial temperature T =1 000 K. The internal energies are chosen so that T_i(ε_i, ρ_i(q)) = T with the density ρ_i(q) evaluated from the initial distribution of the positions; while the velocities are distributed along the Boltzmann distribution. We let the system equilibrate during the time τ_therm = 50 to obtain an equilibrated initial configuration. The shear viscosity is set to η=2×10^-3 Pa·s, and the bulk viscosity ζ is neglected.

4.1 Integrating fluctuation/dissipation dynamics

We first investigate the properties of the integration schemes for the fluctuation/dissipation part only and do not couple the SSA and the metropolized schemes with velocity Verlet. While the SSA is quite stable for large particles (K>10), for which time steps as large as Δt = 5 can be used with no occurrence of a negative internal energy during the simulation time τ_sim, stability issues arise for smaller particles. At K=5, we need a time step Δt < 0.025 to avoid the appearance of negative internal energies with the SSA. As a comparison, the stability limit for the Verlet scheme at K=5 is Δt = 0.8. As the particle size decreases further, it becomes impossible to run simulations, and no admissible time step has been found for K=2. With the rejection of moves provoking negative energies, the metropolized schemes are stable at any time step for every particle sizes.

When they are not coupled to the velocity Verlet scheme, the SSA and metropolized schemes preserve exactly the energy by construction. We can, however, compare the bias in the distributions of internal energies for the different schemes. Figure 2 shows the distributions of internal energy for the exact and approximate metropolized approaches using the ideal gas equation of state (26) with K=5 obtained with the simulation time τ_sim=20 000. The results are compared with the analytic distribution (27). In practice, the distribution ν_Δt obtained from the numerical simulations is approximated using histograms computed on 50 configurations extracted at regular time intervals from the simulations. This ensures a constant number of sampling points for all time steps. The histograms consist in N_bins=150 bins uniformly distributed between ε_min = 0.1k_BT and ε_max = 20k_BT. A more quantitative measurement of the bias is to evaluate the quadratic error with respect to the theoretical distribution,

(28)

Due to the Metropolis procedure, the only source of error for the exact metropolized scheme (EMS) is of statistical nature. This is not guaranteed for the AMS but no systematic bias is apparent up to Δt = 5. The agreement with the theoretical distribution is a bit deteriorated when the full acceptance ratio is not computed, and the AMS displays a 30% larger error compared with the EMS.

We also observe the scaling of the rejection rate with the time step in Fig. 3. The metropolized scheme displays rejection rates between 0.1% and 0.2% for time steps between Δt = 1 and Δt = 5. For our system size (1 000 particles), it means that, on average, several fluctuation/dissipation interactions are rejected each time step. Most of these rejections are due to the Metropolis ratio and not to the appearance of negative internal energies which accounts for approximately one rejection every thousand. When we only reject forbidden moves that would cause negative energies, the rejection rate of the AMS is between 0.005% and 0.01%, which is about 20 times smaller than the overall rejection rate of the EMS. This is however two orders of magnitude larger than the occurrence of negative internal energies with the EMS. By rejecting authorized but unlikely moves (leading to small energies for instance), the EMS is less prompt to the appearance of negative internal energies. A linear fit in log scale shows that the total rejection rates for both the EMS and AMS are Δt^0.42. The rejection of negative energies with the EMS roughly follows the same scaling with Δt^0.5.

4.2 Integrating full SDPD

We now turn to the numerical integration of the full dynamics and study the behavior of the full schemes coupling the velocity Verlet scheme for the integration of the conservative part of the dynamics and either SSA or its metropolized versions for the fluctuation/dissipation dynamics.

To evaluate the energy conservation and the scheme stability, we run n_sim=10 independent simulation with K=5 during the time τ_sim=1 000 for each time step Δt and average the results. The schemes obtained by superimposing (16) and either SSA or the metropolized scheme lead to linear energy drifts which have already been observed in DPDE^{[21, 23]} and in SDPD^[9]. This is illustrated in Fig. 4 where the total energy with respect to the time is displayed for different time steps in the case of the metropolized scheme. We characterize the energy drift by fitting the time evolution of the energy on a linear function and plot the resulting slope in Fig. 5 for the SSA, Metropolis and its approximate version. We observe a similar energy drift for all methods. As in the integration of the fluctuation/dissipation only, the SSA is limited to small time steps to prevent stability issues to arise while the metropolized schemes greatly increase the admissible time steps. However, the time steps for which the dynamics is stable are much smaller with the full dynamics than those reported in Subsection 4.1. Although the conservative interactions are bounded in SDPD unlike the DPDE simulations in Ref. [25], they still induce a stringent stability limit on the time steps. This observation leads us to consider the multiple time step (MTS) implementations where the fluctuation/dissipation is integrated with a larger time step. We introduce the time step Δt_VV used to integrate the conservative part with the velocity Verlet scheme and Δt_FD = θΔt_VV used for the discretization of the fluctuation/dissipation with a metropolized scheme (EMS or AMS). We test this approach with θ= 5 and θ= 10. The algorithm then reads as follows:

(ⅰ) θ consecutive steps of velocity Verlet with Δt =Δt_VV.

(ⅱ) One step of EMS or AMS with Δt = θΔt_VV.

We plot in Fig. 6 the slope of the energy drift compared with their single time step (STS) version (with the time step Δt = Δt_VV). For both EMS and AMS, the energy drift rate is smaller for the MTS approach when the time step is large enough. Moreover, the reduction of the energy drift is enhanced for larger θ with a division by 6 of the rate for θ= 10 and Δt = 0.392 when it is only halved for θ= 5.

We measure the time per iteration and per particle for the different Metropolis schemes with Δt=2.24×10^-2 or Δt_VV=2.24×10^-2 and gather them in Table 1. For the MTS algorithms, the number of iteration is the number of Verlet steps (which is thus the same with that in the STS case). The integration of the fluctuation/dissipation dynamics with the SSA represents a quarter of the total computational time. The integration of the fluctuation/dissipation part with the EMS is about twice as long since we need to compute the reverse move and estimate the Metropolis ratio. This results in an overall increase by 30% of the total simulation time. However, much larger time steps can be chosen with the EMS while the SSA suffers from stringent stability limitations. There is almost no overhead when resorting to the AMS which also greatly improves the stability and is as good as the EMS in terms of energy conservation. With the MTS strategy, the time needed for the fluctuation/dissipation is greatly reduced, as expected, by a factor θ.

4.3 Simulation of nonequilibrium systems

While all the previous simulations were carried out in an equilibrium situation, the metro-polization procedure we propose is also suited for nonequilibrium settings. In particular, SDPD has been applied to model shock waves^[9]. Stability is a crucial issue for these phenomena since they involve dramatic changes in the thermodynamic states of the material. We thus illustrate the enhanced stability of the metropolized schemes in nonequilibrium situations by simulating a shock wave. The system is initialized as previously mentioned but with N=23 400 particles organized on 10 ×10 ×234 lattices with periodic boundary conditions in the x-and y-directions. In the z-direction, two walls are located at each end of the system and formed of "virtual" SDPD particles as described in Refs. [9] and ^[12]. These virtual particles interact with the real SDPD particles through the conservative forces (8) and a repulsive Lennard-Jones potential that ensures the impermeability of the walls. After the system equilibration during τ_therm = 50, the lower wall is given a constant velocity v_P = 1 661 m/s in the z-direction.

To obtain the profiles of physical properties at the given time, the simulation box is divided into n_sl=100 slices regularly distributed along the z-axis in which the physical properties are averaged. Since a shock wave is a stationary process in the reference frame of the shock front, we can average profiles over the time after shifting the position of the shock front to z=0. We plot in Fig. 7 the density profile for the metropolized schemes (EMS and AMS) with Δt = 0.045. This choice is governed by the piston velocity v_P and ensures that the piston does not move by more than 20% of the characteristic distance between two particles. This avoids instabilities in the conservative part of the dynamics. The Rankine-Hugoniot (RH) relations make use of the conservation of mass, momentum, and energy to predict the thermodynamic properties in the shocked state from the initial thermodynamic conditions and the velocity of the particles in the shocked region. The properties estimated from the SDPD simulations, namely, the velocity of the shock wave v_S, the density ρ_S, the pressure P_S, and the temperature T_S in the shocked state, all agree very well with the theoretical predictions as seen in Table 2. Let us point out that this simulation would not have been possible with the SSA since negative internal energies appear very early in the simulation even for time steps as small as Δt = 10^-4. The metropolization procedure can thus be particularly useful in nonequilibrium simulations where stability issues are aggravated.

5 Conclusions

In this work, we have introduced a Metropolis procedure for the integration of the fluctuation/dissipation part in SDPD. This adaptation of the metropolized schemes for DPDE has led to a significant increase in the stability of the dynamics for small particle sizes. This allows us to carry out simulations that traditional schemes could not achieve due to the very stringent limitation on the time step they must respect. It appears that an approximate version of the Metropolis step where only the negative internal energies are rejected is actually enough to ensure the stability of the algorithm and does not display a larger bias or energy drift than its exact version.

With the addition of the Metropolis step, the integration of the fluctuation/dissipation is stable for very large time steps, and the limit on the admissible time steps emerges from the conservative part. An MTS approach has been tested where a smaller time step was used for velocity Verlet and a larger one for the metropolized SSA scheme. This resulted in a similar energy drift at a reduced computational cost.

The relevance of the metropolization in nonequilibrium situations has been illustrated by the simulation of a shock wave. While traditional schemes such as the SSA fail to perform the simulation for small particles, due to the aggravated stability issues, both EMS and AMS have allowed us to recover the correct physical properties of the shock wave.

Our Metropolis schemes still suffer from the difficult parallelization of the SSA on which they are based. It would be most beneficial to adapt a stabilization procedure based on rejecting moves leading to negative energies to parallel schemes^[22-23] in order to deal with larger systems.

Acknowledgements We thank J. B. MAILLET for helpful discussion.