# Stable and accurate schemes for Smoothed Dissipative Particle Dynamics

###### Abstract

Abstract Smoothed Dissipative Particle Dynamics (SDPD) is a mesoscopic particle method which allows to select the level of resolution at which a fluid is simulated.
The numerical integration of its equations of motion still suffers from the lack of numerical schemes satisfying all the desired properties such as energy conservation and stability.
The similarities between SDPD and Dissipative Particle Dynamics with Energy conservation (DPDE), which is another coarse-grained model, enable the adaptation of recent numerical schemes developed for DPDE to the SDPD setting.
In this article, we introduce a Metropolis step in the integration of the fluctuation/dissipation part of SDPD to improve its stability.

Key words Smoothed Dissipative Particle Dynamics, Numerical integration, Metropolis algorithm

2010 Mathematics Subject Classification 82-08, 65C30, 82C80

## 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, 2, 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 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 conservation (DPDE) [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 a 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 [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 small temperatures or when small heat capacities are considered, typically for small mesoparticles. It has been proposed to use Monte Carlo principles 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 prescription [24]. This approach leads however to a dynamics which is not consistent with DPDE. It was proposed in [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 the full SDPD dynamics. Most work 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 [25].

This article is organized as follows. We first present in Section 2 the equations of SDPD as reformulated in [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 Smoothed Dissipative Particle Dynamics

At the hydrodynamic scale, the dynamics of the fluid is governed by the Navier-Stokes equations (1), which read in their Lagrangian form when the heat conduction is neglected (for times and positions 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, its velocity, its internal energy and the stress tensor:

(2) |

where is the pressure of the fluid, the shear viscosity and the bulk viscosity.

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

### 2.1 Smoothed Particle Hydrodynamics

Smoothed Particle Hydrodynamics [5, 6] is a Lagrangian discretization of the Navier-Stokes equations (1) on a finite number of fluid particles playing the role of interpolation nodes. These fluid particles are associated with a portion of fluid of mass . They are located at positions and have a momentum . The internal degrees of freedom are represented by an internal energy . In general, the energies are bounded below. Upon shifting the minimum of the internal energy, we may consider that the internal energies remain positive ().

#### 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 . The kernel is generally required to be non-negative, regular, normalized as and with finite support [27]. We introduce the smoothing length defined such that if . In the sequel, we use the notation . In this work, we rely on a cubic spline [28], whose expression reads

(3) |

The field variables are then approximated as

(4) |

where denotes the value of the field on the particle .

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

In order to have more explicit expressions, we introduce the function such that . In the case of the cubic spline (3), it reads

The gradient approximation can then be rewritten as

In order to simplify the notation, we define the following quantities for two particles and :

We can associate a density and 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 of the mesoparticle with its density (as defined by (5)) and its internal energy through an entropy function

(7) |

The equation of state can be computed by microscopic simulations or by an analytic expression modeling the material behavior. It is then possible to assign to each particle a temperature

pressure

and heat capacity at constant volume

To simplify the notation, we omit in Sections 2.2 the dependence of , and on the variables and .

### 2.2 Equations of motion for SDPD

Smoothed Dissipative Particle Dynamics [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 , the momenta and the energies for .

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

#### 2.2.1 Conservative forces

The elementary force between particles and 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 a Hamiltonian dynamics with a potential relating the energy with the positions and the particle entropy . 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 preserve by construction the total momentum and the total energy .

#### 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 of particles and 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 is chosen of the following form:

(10) |

where is a -dimensional vector of standard Brownian motions and , are symmetric matrices. by construction, (10) preserves the total momentum in the system since Furthermore, as in DPDE, the equations for the energy variables are determined to ensure the conservation of the total energy . 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

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

(12) | ||||

As shown is [9], this ensures that measures of the form (with 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 to retrieve the original SPDP [4].

#### 2.2.3 Complete equations of motion

### 2.3 Reduced units for SDPD

In SDPD, the mass of the fluid particles allows us to change the resolution of the method. We introduce the particle size , where 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 :

(15) | ||||

where is the mass unit, the length unit, the energy unit and 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 , with the particle size used in the simulations. This explains the use of non round time steps in Section 4.

The smoothing length 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, should be rescaled as

In this work, we take , which corresponds to a typical number of 60-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 is divided in simpler elementary dynamics that are consecutively integrated. Since the conservative part of the dynamics (9) can be viewed as a 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 Section 3.1.

There is however no definite way to deal with the fluctuation/dissipation part described in Section 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 [9], is based on the algorithm proposed by Shardlow [20] for DPD and its subsequent adaptations to DPDE [19]. The dynamics is split into a Hamiltonian part, discretized through a Velocity-Verlet algorithm (16), and elementary pairwise fluctuation/dissipation dynamics that are successively integrated. We first recall in Section 3.2 the Shardlow-like splitting scheme (SSA) used in [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 energy are of the same scale. We thus explore methods to improve the stability of these integration algorithms in Section 3.3, relying on the ideas developed in [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 the Hamiltonian part of the 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 entropies . The internal energies are related to the positions and entropies by an energy function , which allows us to write the Hamiltonian as

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

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

(16) |

### 3.2 Shardlow-like Splitting Algorithm

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

(17) |

where is a standard 3-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 [31, 19]. 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 , and heat capacities , 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 temperature or when the particle sizes are small (so that their heat capacity are small as well). This leads to stability issues unless very small timesteps 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 an overdamped Langevin dynamics in the relative velocity variable only, see Section 3.3.1. We then construct proposed moves for the Metropolized scheme and compute the corresponding acceptance ratio in Section 3.3.2. A simplified version of the Metropolized scheme is introduced in Section 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 the fluctuation and dissipation dynamics as an 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 and formulated as an overdamped Langevin dynamics.

Since the dynamics (10) preserve the momentum , the momenta and can be rewritten as a function of as:

and

This already shows how to express the momenta and in terms of . In addition, the kinetic energy formulated in the relative velocity reads

The conservation of the energy and the fact that provides the expression of the internal energies as a function of the relative velocity as

(18) |

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

(19) |

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

(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 ,

The divergence of with respect to the relative velocity reads

The desired result follows from

#### 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