LAMMPS Implementation of Constant Energy Dissipative 

Particle Dynamics (DPD-E)

by James P. Larentzos, John K. Brennan, Joshua D. Moore, and 

William D. Mattson

ARL-TR-6863 March 2014

Approved for public release; distribution is unlimited. 
NOTICES

Disclaimers

The findings in this report are not to be construed as an official Department of the Army position unless 

so designated by other authorized documents.

Citation of manufacturer’s or trade names does not constitute an official endorsement or approval of the 

use thereof.

Destroy this report when it is no longer needed. Do not return it to the originator.
Army Research Laboratory

Aberdeen Proving Ground, MD 21005-5069

ARL-TR-6863 March 2014

LAMMPS Implementation of Constant Energy Dissipative 

Particle Dynamics (DPD-E)

James P. Larentzos

Dynamics Research Corporation (DRC)

John K. Brennan, Joshua D. Moore, and William D. Mattson

Weapons and Materials Research Directorate, ARL

Approved for public release; distribution is unlimited. 
ii

REPORT DOCUMENTATION PAGE Form Approved

OMB No. 0704-0188

Public reporting burden for this collection of information is estimated to average 1 hour per response, including the time for reviewing instructions, searching existing data sources, gathering 

and maintaining the data needed, and completing and reviewing the collection information. Send comments regarding this burden estimate or any other aspect of this collection of information, 

including suggestions for reducing the burden, to Department of Defense, Washington Headquarters Services, Directorate for Information Operations and Reports (0704-0188), 1215 Jefferson 

Davis Highway, Suite 1204, Arlington, VA 22202-4302. Respondents should be aware that notwithstanding any other provision of law, no person shall be subject to any penalty for failing to 

comply with a collection of information if it does not display a currently valid OMB control number.

PLEASE DO NOT RETURN YOUR FORM TO THE ABOVE ADDRESS.

1. REPORT DATE (DD-MM-YYYY)

March 2014

2. REPORT TYPE

Final

3. DATES COVERED (From - To)

10 January 2012–13 September 2013

4. TITLE AND SUBTITLE

LAMMPS Implementation of Constant Energy Dissipative Particle Dynamics 

(DPD-E)

5a. CONTRACT NUMBER

5b. GRANT NUMBER

5c. PROGRAM ELEMENT NUMBER

6. AUTHOR(S)

James P. Larentzos,

* John K. Brennan, Joshua D. Moore, and William D. Mattson

5d. PROJECT NUMBER

5e. TASK NUMBER

5f. WORK UNIT NUMBER

7. PERFORMING ORGANIZATION NAME(S) AND ADDRESS(ES)

U.S. Army Research Laboratory

ATTN: RDRL-WML-B

Aberdeen Proving Ground, MD 21005-5069

8. PERFORMING ORGANIZATION

 REPORT NUMBER

ARL-TR-6863

9. SPONSORING/MONITORING AGENCY NAME(S) AND ADDRESS(ES) 10. SPONSOR/MONITOR’S ACRONYM(S)

11. SPONSOR/MONITOR'S REPORT

 NUMBER(S)

12. DISTRIBUTION/AVAILABILITY STATEMENT

Approved for public release; distribution is unlimited.

13. SUPPLEMENTARY NOTES

*

Dynamics Research Corporation (DRC), High Performance Technologies Group at the U.S. Army Research Laboratory, 

Aberdeen Proving Ground, MD 21005

14. ABSTRACT

A general framework is presented for implementing the constant-energy Dissipative Particle Dynamics (DPD-E) method into 

the highly scalable Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) simulation software to efficiently 

model systems under isoenergetic conditions. The current LAMMPS velocity-Verlet (VV) integration scheme for isothermal 

DPD simulations is extended to the isoenergetic case. A description is given of the implementation of the Shardlow-splitting 

algorithm (SSA) to enable longer time steps with comparable accuracy. In addition, a description of example benchmark 

problems is provided, along with a discussion about the tradeoffs between the DPD version of the VV and the SSA integration 

schemes in terms of performance and stability.

15. SUBJECT TERMS

dissipative particle dynamics, Shardlow splitting, parallelization

16. SECURITY CLASSIFICATION OF: 

17. LIMITATION

OF ABSTRACT

UU

18. NUMBER

OF PAGES

44

19a. NAME OF RESPONSIBLE PERSON

John K. Brennan

a. REPORT

Unclassified

b. ABSTRACT

Unclassified

c. THIS PAGE

Unclassified

19b. TELEPHONE NUMBER (Include area code)

(410) 306-0678

Standard Form 298 (Rev. 8/98)

Prescribed by ANSI Std. Z39.18
iii

Contents

List of Figures v

List of Tables v

Acknowledgments vi

1. Introduction 1

2. Methods 2

2.1 Constant Energy Dissipative Particle Dynamics (DPD-E) .............................................2

2.2 Parallel Implementation of DPD-E Into the LAMMPS Framework...............................5

2.2.1 Implementation of the DPD Atom Style .............................................................5

2.2.2 Implementation of the DPD and DPD/Atom Compute Commands....................6

2.2.3 Implementation of DPD-E With the VV Integration Scheme.............................7

2.2.4 Implementation of DPD-E With the VV-SSA Integration Scheme ....................9

3. Example Input Files for VV and VV-SSA Integration 10

4. Results 13

5. Conclusions 13

6. References 14

Appendix A. Code Differences Between the Atom Class Files as Compared to the Native 

LAMMPS Code 15

Appendix B. Summary of LAMMPS Code Differences Between the AtomVecDPD Class 

as Compared to the Native AtomVecAtomic Class 17

Appendix C. Summary of LAMMPS Code Differences Between the PairDPDE Class as 

Compared to the Native PairDPD Class 21

Appendix D. Summary of LAMMPS Code Differences Between the FixDPDE Class as 

Compared to the Native FixNVE Class 25
iv

Appendix E. Summary of the FixEOScv Class 27

Appendix F. Summary of LAMMPS Code Differences Between the PairDPDE 

Conservative Class as Compared to the Native PairDPD Class 29

Appendix G. Summary of LAMMPS Code Differences Between the FixDPDEShardlow 

Class as Compared to the Native FixNVE Class 31

Distribution List 35
v

List of Figures

Figure 1. Example input file for simulating a DPD fluid under isoenergetic conditions with 

the VV integration scheme.......................................................................................................11

Figure 2. Example input file for simulating a DPD fluid under isoenergetic conditions with 

the VV-SSA integration scheme..............................................................................................12

List of Tables

Table 1. Model and system parameters for DPD-E simulations of the DPD fluid model.............13
vi

Acknowledgments

We acknowledge Martin Lísal (Institute of Chemical Process Fundamentals of the ASCR and J. 

E. Purkinje University, Ústí nad Labem, Czech Republic), Timothy Sirk (U.S. Army Research 

Laboratory), and Timothy Mattox (Dynamics Research Corporation) for useful discussions in 

parallelizing the Shardlow splitting algorithm. We also acknowledge the computational resources 

and Productivity Enhancement, Technology Transfer and Training (PETTT) software support 

from the High Performance Computing Modernization office. Joshua D. Moore acknowledges 

support in part by an appointment to the U.S. Army Research Laboratory Postdoctoral 

Fellowship Program administered by the Oak Ridge Associated Universities through a 

cooperative agreement with the U.S. Army Research Laboratory. John K. Brennan and Joshua

Moore acknowledge support in part by the Office of Naval Research and the Department of 

Defense High Performance Computing Modernization Program Software Application Institute 

for Multiscale Reactive Modeling of Insensitive Munitions. 
1

1. Introduction

Dissipative Particle Dynamics (DPD) is an attractive method for simulating soft condensed 

matter, including polymers, surfactants, and colloids at the mesoscale. DPD is a particle-based 

mesoscopic method that can be understood as a coarse-graining technique that correctly predicts 

the hydrodynamic nature of a fluid. The general idea behind DPD is that particles interact with 

each other through a set of stochastic differential equations with conservative, dissipative, and 

random forces.

While the original DPD formalism is isothermal (1, 2), it has been extended to the isoenergetic 

(DPD-E) (3, 4) case, which is of practical interest for simulating materials at nonisothermal

conditions. One challenge requiring special consideration is the numerical integration of the 

stochastic equations-of-motion (EOMs). The most commonly applied numerical integration 

scheme for DPD applications originates from the velocity-Verlet (VV) integrator used for 

molecular dynamics applications and is extended to DPD by including the dissipative and 

random forces in the overall force expression (5). This approach is currently used in the Largescale Atomic/Molecular Massively Parallel Simulator (LAMMPS) software package (6) but is 

limited only to the constant-temperature DPD method. An alternative integration scheme using 

the Shardlow-splitting algorithm (7) separates the EOMs into stochastic and deterministic 

integration steps. It is readily extendable to all DPD variants and has been found to be the most 

effective integration scheme available to date, especially when considering DPD simulations 

under isoenergetic conditions to attain the length and time scales necessary to model complex 

systems.

In this work, we present a general framework for implementing the DPD-E method into the 

highly scalable LAMMPS simulation software to efficiently model systems under isoenergetic 

conditions. We extend the current LAMMPS VV integration scheme for isothermal DPD 

simulations to the isoenergetic case. In addition, we describe the implementation of the 

Shardlow-splitting algorithm (SSA) to enable longer time steps with comparable accuracy. 

Finally, a description of example benchmark problems is provided, along with a discussion about 

the tradeoffs between the DPD version of the VV and the SSA integration schemes in terms of 

performance and stability. 
2

2. Methods

2.1 Constant Energy Dissipative Particle Dynamics (DPD-E)

In the DPD-E method (3, 4) particles are defined by a mass mi , position i r , momentum pi , and

particle internal energy ui. The variation of dui is taken as the sum of the two contributions that 

correspond to the mechanical work done on the system mech dui and the heat conduction between 

particles cond dui . For constant-energy and constant-volume conditions, the evolution of DPD 

particles in time t is governed by the following EOM:

( )

q

ij

Rq

ij

Dq

j i i j

ij

cond

i

ij ij

ij

R ij

ij

j i

R

i j

ij

ij

ij

D ij

ij

mech

i

j i

ij

ij

R ij

ij

ij

ij

ij

ij

D ij

ij

C

i ij

i

i

i

u t W

W

r

t

m m

d

t

r

u

W

r

t

r r

t

t

m

d d 1 1 d

d

d 1 1

2

d

2

1 d

d d d d

d d

2

2 2

ω α ω

θ θ

κ

σ ω

ω

σ

γ ω

γ ω σ ω

+ 













 = −















 − ⋅















 − + 













 = ⋅

+ 













 = − ⋅

=

∑

∑

∑

≠

≠

≠

v

r

v

r

r r

v

r

p F

p r

(i = 1,..., N) (1)

where ij i j r = r − r is the separation vector between particle i and particle j , and ij ij r = r .

j

j

i

i

ij m m

p p

v = − is the relative velocity between the particle pair. ij γ and σ ij are the friction 

coefficient and noise amplitude between particle i and particle j , respectively. Wij Wji d = d are 

the increments of the Wiener processes associated with momentum variations. The weight 

functions (r) D ω and (r) R ω vanish for c r ≥ r , where cr is the cutoff radius, and are typically 

chosen as 

ωD(r) = [ωR(r)]2 = �

�1 − r

rc

�

2

, r < rc

0 , r ≥ rc

(2)

Similarly, the weight functions (r) Dq ω and (r) Rq ω are defined in an equivalent manner as 

equation 2. θ i is the particle internal temperature and is related to ui through a mesoparticle 

equation of state. κ ij and αij are the mesoscopic thermal conductivity and the noise amplitude 
3

between particle i and particle j , respectively, and q

ji

q dWij = −dW are the increments of the 

Wiener processes associated with thermal conduction. The conservative force C Fij is given as the 

negative derivative of a coarse-grain potential, CG uij , and is expressed in LAMMPS as

Fij

C = − duij

CG

drij

= AijωR (3)

where Aij is the force constant. Bonet Avalos and Mackie (3) demonstrated that thermodynamic 

consistency requires the following fluctuation-dissipation relations to be satisfied

( ) [ ( )]

( ) [ ( )]

2

B

2

2

B

2

2

2

r r

k

r r

k

Dq Rq

ij ij

D R

ij ij ij

ω ω

α κ

ω ω

σ γ

=

=

=

= Θ

(4)

where the relevant temperature is 















Θ = + −

i j

ij θ θ

1 1

2

1 1 . In the following outline, the EOMs are 

formulated for the VV integration scheme under isoenergetic conditions.

1. For i =1,...,N

a. i i p' ← p

b. i i i

t

p p F

2

∆

← + , where ∑≠

= + +

j i

R

ij

D

ij

C Fi Fij F F

c.

i

i

i i

m

t

p r ← r + ∆

d. ∑≠

∆

+

∆

















← + − j i

q

ij

Rq

ij

Dq

i j

ij

cond

i

cond

i

t t

u u

2 2

1 1

ω α ω ς θ θκ

e.

( )

2

2

1 1

2 2 2

1 2

2 2

t

r

t

m m

t

r

u u

ij ij

ij

R ij

ij

j i

R

i j

ij

ij

ij

D ij

ij

mech

i

mech

i

∆















 − ⋅

∆















 − +

∆

















← − ∑ ⋅ ≠

σ ω ς

ω

σ

γ ω

v

r

v

r

2. Force Calculation: { }

N Fi i=1 , where ∑≠

= + +

j i

R

ij

D

ij

C Fi Fij F F

3. For i =1,...,N

a. i i p' ← p
4

b. i i i

t

p p F

2

∆

← + , where ∑≠

= + +

j i

R

ij

D

ij

C Fi Fij F F

c.

i

i

i i

m

t

p r ← r + ∆

d. ∑≠

∆

+

∆

















← + − j i

q

ij

Rq

ij

Dq

i j

ij

cond

i

cond

i

t t

u u

2 2

1 1

ω α ω ς θ θκ

e.

( )

2

2

1 1

2 2 2

1 2

2 2

t

r

t

m m

t

r

u u

ij ij

ij

R ij

ij

j i

R

i j

ij

ij

ij

D ij

ij

mech

i

mech

i

∆















 − ⋅

∆















 − +

∆

















← − ∑ ⋅ ≠

σ ω ς

ω

σ

γ ω

v

r

v

r

Here, the EOMs are formulated for the VV integration scheme using the Shardlow-splitting 

algorithm (VV-SSA) (8).

1. Stochastic Integration: For all i − j pairs of particles

a. i i p' ← p , j j p' ← p

b.

2 2

t

r r r

t

ij

ij

ij

R

ij

ij

ij

ij

ij

D ij

i i ij

∆

+ 













 ⋅ ∆

← − r r

v

r

p p γ ω σ ω ς

c.

2 2

t

r r r

t

ij

ij

ij

R

ij

ij

ij

ij

ij

D ij

j j ij

∆ − 













 ⋅ ∆

← +

r r

v

r

p p γ ω σ ω ς

d.

j

j

i

i

ij m m

p p

v ← −

e.

2 2

2

1 2

t

r

t

r r r

t

t

ij

ij

ij

R

ij

ij

ij

ij

R

ij

ij

ij

ij

ij

ij

ij

D

ij

ij

D

ij

i i

∆

+

















+ ∆ 













 ⋅

+ ∆

∆

← − r r r

v

r

p p σ ω ς σ ω ς

μ

γ ω

μ

γ ω

f.

2 2

2

1 2

t

r

t

r r r

t

t

ij

ij

ij

R

ij

ij

ij

ij

R

ij

ij

ij

ij

ij

ij

ij

D

ij

ij

D

ij

j j

∆ −

















+ ∆ 













 ⋅

+ ∆

∆

← +

r r r

v

r

p p σ ω ς σ ω ς

μ

γ ω

μ

γ ω where

i j

ij m m

1 1

μ = +

g. u u t t q

ij

Rq

ij

Dq

i j

ij

cond

i

cond

i ∆ + ∆ 















← + − ω α ω ς θ θ

κ

1 1

h. u u t t q

ij

Rq

ij

Dq

i j

ij

cond

j

cond

j ∆ − ∆ 















← − − ω α ω ς θ θκ

1 1
5

i.















 ⋅ − ⋅ − ⋅ + ⋅ ← −

j

j j

i

i i

j

j j

i

mech i i

i

mech

i

m

' '

m

' '

m m

u u

2 2 2 2 2

1 p p p p p p p p

j. 













 ⋅ − ⋅ − ⋅ + ⋅ ← −

j

j j

i

i i

j

j j

i

mech i i

j

mech

j m

' '

m

' '

m m

u u

2 2 2 2 2

1 p p p p p p p p

2. Deterministic Integration #1: For i =1,...,N

a.

C

i i i

t

p p F

2

∆

← +

b.

i

i

i i

m

t

p r ← r + ∆

3. Conservative Force Calculation: { }

N

i

C Fi =1 , where ∑≠

= j i

C

ij

C Fi F

4. Deterministic Integration #2: For i =1,...,N

C

i i i

t

p p F

2

∆

← +

2.2 Parallel Implementation of DPD-E Into the LAMMPS Framework

A detailed description of the theoretical foundations and parallelization of the DPD-E method

using the VV and VV-SSA integration schemes can be found elsewhere (8, 9). In this report, the 

relevant source code modifications that were required to implement the VV and VV-SSA 

integration schemes for DPD-E are presented. A summary of the code modifications is provided 

in appendices A–G. The source code is current with the 27 January 2014 version of LAMMPS 

and is available upon written request to the authors.

To implement DPD-E into LAMMPS, a new atom style was created to handle the DPD particle 

attributes required for DPD-E simulations. In addition, new pair styles were created to compute 

the DPD forces, and new fixes were created to integrate the DPD-E equations of motion through 

a VV and VV-SSA integration scheme. Compute functions were created to monitor the DPD 

particle attributes as a simulation progresses. Each new feature to LAMMPS is described in 

detail in the following sections.

2.2.1 Implementation of the DPD Atom Style

For all isoenergetic DPD calculations, a new atom style dpd is required to compute and 

communicate between processors the DPD particle attributes and per-atom arrays. DPD particles 

are specified in the LAMMPS input file via the atom_style command:

atom_style dpd

Selection of the dpd atom style requires the DPD internal temperature to be specified in the 

corresponding LAMMPS data file according to the following format:

• Column 1: particle id
6

• Column 2: particle type

• Column 3: Internal temperature (θi) of particle i

• Columns 4–6: Cartesian coordinates (x, y, z) of particle i

The implementation of the dpd atom style to LAMMPS requires a new AtomVecDPD class to be 

created in order to compute and communicate the DPD particle internal temperature (dpdTheta), 

equation of state flag (eos), conductive energy (uCond), mechanical energy (uMech), as well as 

the differences in the DPD particle’s conductive energy (duCond) and mechanical energy 

(duMech) between two subsequent integration steps. The AtomVecDPD atom style class derives 

from the AtomVec parent class and is similar to and modeled after the existing AtomVecAtomic 

(atom style atomic) class. 

Implementation into LAMMPS requires modification of the existing atom class (LAMMPS 

source files: atom.h and atom.cpp) and creation of the AtomVecDPD class (LAMMPS source 

files: atom_vec_dpd.h and atom_vec_dpd.cpp). The atom class code is modified to define the 

necessary pointers to the DPD attributes. The differences between the modified and the native 

LAMMPS codes are shown in appendix A. Next, the new AtomVecDPD child class was created 

by copying the existing AtomVecAtomic class (LAMMPS source files: atom_vec_atomic.h and 

atom_vec_atomic.cpp), then adjusting the data file format, the data structure sizes, and the

communication buffers to handle the additional per-particle arrays. A summary of the differences 

between the AtomVecDPD and AtomVecAtomic classes is given in appendix B.

2.2.2 Implementation of the DPD and DPD/Atom Compute Commands

For all isoenergetic DPD calculations, it is useful to monitor all internal properties of the 

particles on both a system and per-particle basis. The ComputeDPD and ComputeDPDatom 

classes were created to monitor these properties over the course of the DPD-E simulations. 

The ComputeDPD class (LAMMPS source files: compute_dpd.h and compute_dpd.cpp) 

computes the total particle internal conductive and mechanical energies by summing the perparticle energies. In addition, the particle internal temperature of the system is computed through 

a harmonic average of the per-particle internal temperatures, defined as:

θavg

−1 = 1

N� 1θi

N

i=1

(5)

where N is the number of particles in the system. The ComputeDPD class is accessed through the 

LAMMPS input file via the compute command:

compute 1 all dpd

and returns a vector of size 5 that contains the following particle internal properties:
7

• The total conductive energy of the system

• The total mechanical energy of the system

• The sum of the conductive and mechanical energy of the system

• The harmonic averaged internal temperature of the system

• The number of particles in the group

The ComputeDpdAtom class (LAMMPS source files: compute_dpd_atom.h and 

compute_dpd_atom.cpp) accesses the per-particle internal energies and internal temperature. The 

compute is specified through the LAMMPS input file via the compute command

compute 1 all dpd/atom

and enables access to the following particle properties:

• The per-particle conductive energy

• The per-particle mechanical energy

• The per-particle internal temperature 

2.2.3 Implementation of DPD-E With the VV Integration Scheme

Implementation of the isoenergetic DPD using the VV numerical integration scheme is similar to 

the existing isothermal DPD LAMMPS implementation. The major difference is the calculation 

of the particle internal energies within the dpde pair style compute command and the integration 

of the internal energies within the fix dpde command. The new pair style and fix commands are 

presented together since the classes interact with one another to perform all stages of the VV 

integration. 

2.2.3.1 Implementation of the DPDE Pair Style

Implementation of constant energy DPD with the VV integration scheme requires a new pair 

style to be defined, which computes the conservative, dissipative, and random forces as well as 

the change in particle internal conductive energy and mechanical energy from one integration 

step to the next. The new PairDPDE class (LAMMPS source files: pair_dpde.h and 

pair_dpde.cpp) is specified through the LAMMPS input file via the pair_style command

pair_style dpde <kappa flag> <cutoff> <random number seed>

pair_coeff i j ADPD sigmaij kappa <cutoff>

An energy-independent or energy-dependent (10) kappa model is specified in the pair_style 

command by setting the kappa_flag to 0 or 1, respectively. In the energy-independent model 

(kappa_flag = 0), κij is explicitly given as a pair coefficient in the LAMMPS input file. In the 

energy-dependent model (kappa_flag = 1), κij is given by the equation 
8

κij = κ0

kB

�

ui+uj

2 �

2

(6)

where the κ0 pair coefficient is specified in the LAMMPS input file, ui is the total internal energy 

of particle i, and kB is Boltzmann’s constant.

The PairDPDE pair style class derives from the Pair parent class and is very similar to and 

modeled after the existing PairDPD pair style class, except that it contains additional per-atom 

arrays that update the DPD particle internal temperature (dpdTheta), conductive energy (uCond), 

mechanical energy (uMech), heat capacity (cv), and density (rho), as well as the differences in 

the DPD particle’s conductive energy (duCond) and mechanical energy (duMech) between two 

subsequent integration steps. The main differences between the native PairDPD and modified 

PairDPDE classes are summarized in appendix C.

2.2.3.2 Implementation of DPDE Fix Commands

Implementation of constant energy DPD with the VV integration scheme also requires a new fix 

class to be defined, which accounts for the integration of the conductive and mechanical energy

and applying a mesoparticle equation of state. The new fix is defined as the FixDPDE class

(LAMMPS source files: fix_dpde.h and fix_dpde.cpp) and is specified in the LAMMPS input file 

via the fix command

fix 1 all dpde

The FixDPDE is very similar to and modeled after the existing FixNVE class, except that it 

contains the integration of the conductive and mechanical energy. The main differences between 

the native FixNVE and modified FixDPDE classes are summarized in appendix D.

Once the updated conductive and mechanical energy is computed, the DPD particle internal 

temperatures are updated through the mesoparticle equation of state that relates the total particle 

internal particle energy to the internal temperature. The mesoparticle equation of state must be 

specified with a separate fix command. Currently, only one choice of the mesoparticle equation 

of state has been implemented. It relates the total internal energy to the internal temperature 

through the heat capacity according to the relation

ui = CV,iθi (7)

The new fix is defined as the FixEOScv class (LAMMPS source files: fix_eos_cv.h and 

fix_eos_cv.cpp) and is specified in the LAMMPS input file via the fix command:

fix 1 all eos/cv <cv>

where cv is the value of the heat capacity. The FixEOScv class is summarized in appendix E.
9

2.2.4 Implementation of DPD-E With the VV-SSA Integration Scheme

Implementation of the isoenergetic DPD using the VV-SSA integration scheme requires a 

splitting of the stochastic and deterministic EOMs. The conservative force is computed within 

the dpde/conservative pair style and is integrated deterministically through the VV integration 

scheme with the fix dpde/shardlow command. The random and dissipative forces are computed 

and integrated through the SSA within the fix dpde/shardlow prior to the deterministic 

integration of the conservative force. The new pair style and fix commands are presented 

together since the classes interact with one another to perform all stages of the VV-SSA 

integration. Additional implementation details of the SSA can be found elsewhere (9).

2.2.4.1 Implementation of the DPDE/Conservative Pair Style

Implementation of constant energy DPD using the VV-SSA integration scheme requires a new 

pair style that computes the conservative force to be defined. The new PairDPDEConservative 

class (LAMMPS source files: pair_dpde_conservative.h and pair_dpde_conservative.cpp) is 

specified in the LAMMPS input file via the pair_style command

pair_style dpde/conservative <kappa flag> <cutoff> <random number seed>

pair_coeff i j ADPD sigmaij kappa <cutoff>

An energy-independent or energy-dependent (10) kappa model is specified in the pair_style 

command by setting the kappa_flag to 0 or 1, respectively. In the energy-independent model 

(kappa_flag = 0), κij is explicitly given as a pair coefficient in the LAMMPS input file. In the 

energy-dependent model (kappa_flag = 1), κij is related to the κ0 pair coefficient as given in 

equation 6.

The PairDPDEConservative pair style class derives from the Pair parent class and is very similar 

to and modeled after the existing PairDPD pair style class. The main difference is that random 

and dissipative computations are removed, while the required data structures to compute the 

DPD particle attributes are included. The main differences between the PairDPD and 

PairDPDEConservative classes are summarized in section 2.2.4.2 and are explicitly shown in 

appendix F.

2.2.4.2 Implementation of DPDE/Shardlow Fix Command

Implementation of constant energy DPD using the VV-SSA integration also requires a new fix 

class to be defined, which accounts for the stochastic integration of the random and dissipative 

forces, the deterministic integration of the conservative force, and the integration of the 

conductive and mechanical energy for constant energy simulations. The new fix is defined as the 

FixDPDEShardlow class (LAMMPS source files: pair_dpde_shardlow.h and pair_dpde_

shardlow.cpp) and is specified in the LAMMPS input file via the fix command:

fix 1 all dpde/shardlow
10

The FixDPDEShardlow class is very similar to the existing FixNVE class, except that it contains

an additional function to integrate the stochastic equations of motion, the conductive energy, and 

the mechanical energy. These integrations are handled via a function called stochastic_integrate 

within the FixDPDEShardlow class. The DPD particle internal temperatures are updated through 

the mesoparticle equation of state that relates the total particle internal particle energy to the 

internal temperature. The mesoparticle equation of state must be specified with a separate fix 

command. Currently, only one choice of the mesoparticle equation of state has been 

implemented and can be specified in the LAMMPS input file via the fix command

fix 1 all eos/cv <cv>

The main differences between the native FixNVE and modified FixDPDEShardlow classes are 

summarized in appendix G.

3. Example Input Files for VV and VV-SSA Integration

Example input files running a DPD-E simulation in LAMMPS with the VV and VV-SSA 

integration schemes are provided in figures 1 and 2, respectively. 
11

Figure 1. Example input file for simulating a DPD fluid under isoenergetic conditions with the VV integration 

scheme. 

# Input File for DPD fluid under isoenergetic conditions using the VV integration scheme

boundary p p p

units metal # ev, ps

atom_style dpd

read_data initial.conf.DPDfluid

communicate single vel yes

mass 1 125.9

pair_style dpde 1 8.60 234324

pair_coeff 1 1 0.0752 0.0223 4.55E-05 8.60

neighbor 2.0 bin

neigh_modify every 1 delay 0 check no once no

# Time in ps for metal units

timestep 0.001

compute dpdU all dpd

variable totEnergy equal pe+ke+c_dpdU[3]

thermo 1

thermo_style custom step temp pe ke c_dpdU[1] c_dpdU[2] c_dpdU[3] v_totEnergy c_dpdU[4]

thermo_modify format float %15.10E

fix 1 all dpde

fix 2 all eos/cv 0.00517041

run 100
12

Figure 2. Example input file for simulating a DPD fluid under isoenergetic conditions with the VV-SSA 

integration scheme.

# Input File for DPD fluid under isoenergetic conditions using the VV-SSA integration scheme

boundary p p p

units metal # ev, ps

atom_style dpd

read_data initial.conf.DPDfluid

communicate single vel yes

mass 1 125.9

pair_style dpde/conservative 1 8.60 234324

pair_coeff 1 1 0.0752 0.0223 4.55E-05 8.60

neighbor 2.0 bin

neigh_modify every 1 delay 0 check no once no

# Time in ps for metal units

timestep 0.001

compute dpdU all dpd

variable totEnergy equal pe+ke+c_dpdU[3]

thermo 1

thermo_style custom step temp pe ke c_dpdU[1] c_dpdU[2] c_dpdU[3] v_totEnergy c_dpdU[4]

thermo_modify format float %15.10E

fix 1 all dpde/shardlow

fix 2 all eos/cv 0.00517041

run 100

The example input files contain the commands to simulate a 10,125 particle system in a cubic 

box with volume 129.0 Å3 at a temperature of 300 K under isoenergetic conditions with a time 

step of 0.001 ps. The eos/cv equation of state is specified along with an energy-dependent kappa 

model. The neighbor lists are reconstructed after every time step, and the thermodynamics are 

output every 0.001 ps. The total particle internal energies and internal temperature are included 

in the thermodynamics via the compute dpd command. The DPD-E parameters are summarized 

in table 1.
13

Table 1. Model and system parameters for DPD-E simulations 

of the DPD fluid model.

System Property Real Units

ρDPD 4.72 × 10-3 Å3

TSP = Tinit 300 K

mass 125.9 g mol-1

ADPD 7.52 × 10-2 eV Å-1

rcut 8.6 Å

σ 2.23 × 10-2 eV ps1/2 Å-1

κ0 4.55 × 10-5 ps-1

CV 5.17 × 10-3 eV K-1

4. Results

The example benchmarks described in section 5 are compared in terms of performance and 

accuracy in Lísal et al. (8), where DPD fluid calculations were performed under isoenergetic 

conditions with 10,125 particles at time steps ranging from 0.001 to 0.4 ps for a total of 1 ns of 

simulation time. Lísal et al. (8) shows that the VV-SSA integrator is considerably more stable 

than the VV integrator, where comparable accuracy is achieved with time steps that are 10–100 

times longer than the traditional VV integrator. However, the VV-SSA implementation is shown 

to be more computationally expensive than the serial VV implementation on a per time step basis 

and carries additional parallelization overhead. Analysis of the trade-off between computation 

performance and time-step size stability showed that the VV-SSA integration scheme can 

decrease the overall time-to-solution by a factor of 10–100 for a DPD-E simulation, thereby 

justifying practical and regular implementation of the approach. 

5. Conclusions

In this report we have provided the user with the following information and resources:

• A detailed description of the implementation of the DPD-E method using the VV and 

VV-SSA integration schemes into LAMMPS. 

• Documentation that highlights the major modifications to the source code. 

• Example input scripts to run the software. 

• A benchmark study to give an indication of the performance and accuracy of DPD-E with 

the VV and VV-SSA integration schemes. 
14

6. References

1. Hoogerbrugge, P. J.; Koelman, J. M. V. A. Simulating Microscopic Hydrodynamic 

Phenomena With Dissipative Particle Dynamics. Europhys. Lett. 1992, 19, 155.

2. Koelman, J. M. V. A.; Hoogerbrugge, P. J. Dynamic Simulation of Hard Sphere Suspensions 

Under Steady Shear. Europhys. Lett. 1993, 21, 363.

3. Bonet Avalos, J.; Mackie, A. D. Dissipative Particle Dynamics With Energy Conservation. 

Europhys. Lett. 1997, 40 (2), 141.

4. Español, P. Dissipative Particle Dynamics With Energy Conservation. Europhys. Lett. 1997,

40 (6), 631.

5. Groot, R. D.; Warren, P. B. Dissipative Particle Dynamics: Bridging the Gap Between 

Atomistic and Mesoscopic Simulation. J. Chem. Phys. 1997, 107, 4423.

6. Plimpton, S. Fast Parallel Algorithms for Short-Range Molecular Dynamics. J. Comp. Phys.

1995, 117, 1–19.

7. Shardlow, T. Splitting for Dissipative Particle Dynamics. SIAM J. Sci. Comput. 2003, 24,

1267.

8. Lísal, M.; Brennan, J. K.; Avalos, J. B. Dissipative Particle Dynamics at Isothermal, 

Isobaric, Isoenergetic, and Isoenthalpic Conditions Using Shardlow-Like Splitting 

Algorithms. J. Chem. Phys. 2011, 135, 204105. 

9. Larentzos, J. P.; Brennan, J. K.; Moore, J. D.; Lísal, M.; Mattson, W. D. Parallel 

Implementation of Isothermal and Isoenergetic Dissipative Particle Dynamics Using 

Shardlow-Like Splitting Algorithms. Comput. Phys. Commun., in press.

10. Hernando, R. Kinetic Theory of Dissipative Particle Dynamics Models. Ph.D Thesis, 

Universidad Nacional de Educaciόn a Distancia, 2002.
15

Appendix A. Code Differences Between the Atom Class Files as Compared to 

the Native LAMMPS Code 

 

This appendix appears in its original form, without editorial change.
16

diff --git a/src/atom.h b/src/atom.h

index cc044f1..1d122ec 100644

---a/src/atom.h

+++ b/src/atom.h

@@ -50,6 +50,10 @@ class Atom : protected Pointers {

 int *type,*mask;

 imageint *image;

 double **x,**v,**f;

+ double *uCond, *uMech;

+ double *duCond, *duMech;

+ double *dpdTheta;

+ int *eos;

 tagint *molecule;

 int *molindex,*molatom;

diff --git a/src/atom.cpp b/src/atom.cpp

index 86ae09b..888aebe 100644

---a/src/atom.cpp

+++ b/src/atom.cpp

@@ -72,6 +72,10 @@ Atom::Atom(LAMMPS *lmp) : Pointers(lmp)

 type = mask = NULL;

 image = NULL;

 x = v = f = NULL;

+ uCond = uMech = NULL;

+ duCond = duMech = NULL;

+ dpdTheta = NULL;

+ eos = NULL;

 molecule = NULL;

 molindex = molatom = NULL;
17

Appendix B. Summary of LAMMPS Code Differences Between the 

AtomVecDPD Class as Compared to the Native AtomVecAtomic Class 

 

This appendix appears in its original form, without editorial change.
18

• AtomVecDPD Constructor: 

o Initialized uCond, uMech, duCond, duMech and eos pointers to NULL

o Updated comm_x_only and comm._f_only to 1 to allow for forward communication of rho, cv and dpdTheta

o Updated data buffer size definitions for the forward and border communications (size_forward=9, size_border=12) to 

account for cv, dpdTheta, eos, uCond and uMech

o Updated the number of columns in the data file (size_data_atom = 6). The data file contains the following six columns:

 Column 1: particle id

 Column 2: particle type

 Column 3: Internal temperature (θi

) of particlei

 Columns 4-6: Cartesian coordinates (x, y, z) of particlei

o Updated the column the x data is located (xcol_data = 4)

o Set the rho_flag to be on (atom->rho_flag = 1)

• void AtomVecDPD::grow(int n)

o Added memory allocation for rho, drho, cv, dpdTheta, eos, uCond, uMech, duCond and duMech

• void AtomVecDPD::grow_reset()

o Added pointers for rho, drho, cv, dpdTheta, eos, uCond, uMech, duCond and duMech

• void AtomVecDPD::copy(int i, int j, int delflag) 

o Copies particle i rho, drho, cv, dpdTheta, eos, uCond, uMechto particle j

• int AtomVecDPD::pack_comm(int n, int *list, double *buf, int pbc_flag, int *pbc)

o Packs rho, cv, dpdTheta, eos, uCond and uMech in the forward communication buffer

• int AtomVecDPD::pack_comm_vel(int n, int *list, double *buf, int pbc_flag, int *pbc)

o Packs rho, cv, dpdTheta, eos, uCond and uMech in the forward communication buffer, along with velocities

• void AtomVecDPD::unpack_comm(int n, int first, double *buf)

o Unpacks rho, cv, dpdTheta, eos, uCond and uMech from forward communication buffer

• void AtomVecDPD::unpack_comm_vel(int n, int first, double *buf)

o Unpacks rho, cv, dpdTheta, eos, uCond and uMech from forward communication buffer, along with velocities

• int AtomVecDPD::pack_border(int n, int *list, double *buf, int pbc_flag, int *pbc)

o Packs rho, cv, dpdTheta, eos, uCond and uMech in the border communication buffer. 

• int AtomVecDPD::pack_border_vel(int n, int *list, double *buf, int pbc_flag, int *pbc)

o Packs rho, cv, dpdTheta, eos, uCond and uMech in the border communication buffer, along with velocities

• void AtomVecDPD::pack_comm_hybrid(int n, int *list, double *buf)

o Packs rho, cv, dpdTheta, eos, uCond and uMech for the hybrid atom style communication buffer

• void AtomVecDPD::pack_border_hybrid(int n, int *list, double *buf)

o Packs rho, cv, dpdTheta, eos, uCond and uMech in the hybrid atom style border communication buffer

• void AtomVecDPD::unpack_border(int n, int first, double *buf)

o Unpacks rho, cv, dpdTheta, eos, uCond and uMech from the border communication buffer

• void AtomVecDPD::unpack_border_vel(int n, int first, double *buf)

o Unpacks rho, cv, dpdTheta, eos, uCond and uMech from the border communication buffer, along with velocities

• void AtomVecDPD::unpack_comm_hybrid(int n, int first, double *buf)

o Unpacks rho, cv, dpdTheta, eos, uCond and uMech from the hybrid atom style communication buffer

• void AtomVecDPD::unpack_border_hybrid(int n, int first, double *buf)

o Unpacks rho, cv, dpdTheta, eos, uCond and uMech from the hybrid atom style border communication buffer

• int AtomVecDPD::pack_exchange(int i, double *buf)

o Packs rho, cv, dpdTheta, eos, uCond and uMech in the exchange communication buffer. 

• int AtomVecDPD::unpack_exchange(double *buf)
19

o Unpacks rho, cv, dpdTheta, eos, uCond and uMech from exchange communication buffer. 

• int AtomVecDPD::size_restart()

o Adjusted the size of the restart data to include rho, cv, dpdTheta, eos, uCond and uMech

 Changed from 11 to 17

• int AtomVecDPD::pack_restart(int i, double *buf)

o Packs rho, cv, dpdTheta, eos, uCond and uMech in the restart buffer

• int AtomVecDPD::unpack_restart(double *buf)

o Unpacks rho, cv, dpdTheta, eos, uCond and uMech from the restart buffer

• void AtomVecDPD::create_atom(int itype, double *coord)

o Creates a new particle and sets default rho, drho, cv, dpdTheta, eos, uCond, uMech, duCond and duMech values 

• void AtomVecDPD::data_atom(double *coord, tagint imagetmp, char **values)

o Read internal temperature (dpdTheta) from data file; initialize rho, cv, uCond[i] and uMech[i] to zeroand eos to -1

• int AtomVecDPD::data_atom_hybrid(int nlocal, char **values)

o Read internal temperature (dpdTheta) from hydrid atom style, data file 

• void AtomVecDPD::pack_data(double **buf)

o Packs internal temperature (dpdTheta) into the data file buffer

• int AtomVecDPD::pack_data_hybrid(int i, double *buf)

o Packs internal temperature (dpdTheta) into the hydrid atom style, data file buffer

• void AtomVecDPD::write_data(FILE *fp, int n, double **buf)

o Writes internal temperature (dpdTheta) to the data file

• int AtomVecDPD::write_data_hybrid(FILE *fp, double *buf)

o Writes internal temperature (dpdTheta) to the hybrid atom info to data file

• bigint AtomVecDPD::memory_usage()

o Added memory checks for rho, drho, cv, dpdTheta, eos, uCond, uMech, duCond and duMech
20

INTENTIONALLY LEFT BLANK.
21

Appendix C. Summary of LAMMPS Code Differences Between the PairDPDE 

Class as Compared to the Native PairDPD Class 

 

This appendix appears in its original form, without editorial change.
22

• PairDPDE constructor: PairDPDE::PairDPDE(LAMMPS *lmp) : Pair(lmp)

o Set the size of the pair forward communication to 3(duCondi, duMechi, dpdThetai)

o Set the size of the pair reverse communication to 2 (duCondi

, duMechi)

• PairDPDE destructor: PairDPDE::~PairDPDE()

o Deallocate kappa array; remove gamma deallocation

• void PairDPDE::compute(int eflag, int vflag)

o Initialize duCond and duMech by setting the arrays to 0.0 followed by a forward communication 

o Define wR and wD: ωD(r) = [ωR(r)]2 = �

�1 − r

rc

�

2

, r < rc

0 , r ≥ rc

where rc is the cutoff radius for the pair of interacting particles.

NOTE: In PairDPD class, wd definition differs

o Compute theta_ij: Θij

−1 = 1

2 � 1θi

+ 1

θj

�

o Compute gamma_ij: γij = σij

2

2kBθij

o Compute the conservative, dissipative and random forces

R

ij

D

ij

C Fij = Fij + F + F

Fij

C = ADPDωR(r) rij

rij

( )

ij

ij

ij

ij

ij

ij

D

ij

D

ij r r

r

r

v

r

F 













 = −γ ω ⋅

Fij

R = σijωR(r)ςij

rij

rij

1

√Δt

where ςijj is a random number sampling a Gaussian distribution and Δt is the timestep used in the simulation

o Compute kappa_ij: κijor κij = κ0

kB

�

ui+uj

2 �

2

o Compute alpha_ij: αij = �2kBκij

o Compute mu_ij: μij = � 1

mi

− 1

mj

�

o Compute duMechi and duMechj

dui

mech = 1

2 �γijωD �

rij

rij

∙ vij�

2

− σij

2

2 μijωD − σijωR �

rij

rij

∙ vij� ςij

2

√Δt�

duj

mec h = duimech

o Compute duCondi and duCondj

dui

cond = κij � 1

θi

− 1

θj

� ωDq + αijωRqςij

q 1

√Δt

 

duj

cond = −duicond

o Compute the interaction energy between the particle pair

 evdwl = 0.5*a0[itype][jtype]*cut[itype][jtype] * wd;

o Reverse communication of duCond and duMech

• void PairDPDE::allocate()

o Allocate memory for kappa array; Remove gamma allocation

• void PairDPDE::settings(int narg, char **arg)
23

o Read 3 arguments associated with pair style (kappa_flag, cut_global and seed); Remove temperature

o Ensure kappa_flag is set to either 0 or 1

• void PairDPDE::coeff(int narg, char **arg)

o Reads in 3-4 three pair coefficients

 ADPD 

 Sigma (σij); removed gamma

 Kappa_ij (κij) or Kappa0 (κ0)

 cutoff (optional)

• void PairDPDE::init_style()

o Ensures that the Newton pair is set to “on” for DPD calculations

• double PairDPDE::init_one(int i, int j)

o Initializes the pair coefficients for all particle pairs; Add kappa, sigma and remove gamma

• void PairDPDE::write_restart(FILE *fp)

o Updates the pair coefficient information to the restart file; Add kappa, sigma and remove gamma

• void PairDPDE::read_restart(FILE *fp)

o Read pair coefficients from restart file and broadcast across processors; Add kappa, sigma and remove gamma

• void PairDPDE::write_restart_settings(FILE *fp)

o Removed the temperature variable 

• void PairDPDE::read_restart_settings(FILE *fp)

o Removed temperature variable

• void PairDPD::write_data(FILE *fp)

o Removed gamma; Added kappa and sigma to list of variable to write to the data file

• void PairDPD::write_data_all(FILE *fp)

o Removed gamma; Added kappa and sigma to list of variable to write to the data file

• double PairDPDE::single(int i, int j, int itype, int jtype, double rsq, double factor_coul, double factor_dpd, double &fforce)

o Updated definition of wr and wd

• int PairDPDE::pack_comm(int n, int *list, double *buf, int pbc_flag, int *pbc)

o Packs the buffer for forward communication(duCond, duMech and dpdTheta)

• void PairDPDE::unpack_comm(int n, int first, double *buf)

o Unpacks the buffer for forward communication (duCond, duMech and dpdTheta)

• int PairDPDE::pack_reverse_comm(int n, int first, double *buf)

o Packs the buffer for reverse communication(duCond and duMech)

• void PairDPDE::unpack_reverse_comm(int n, int *list, double *buf)

o Unpacks the buffer for reverse communication (duCond and duMech)
24

INTENTIONALLY LEFT BLANK.
25

Appendix D. Summary of LAMMPS Code Differences Between the FixDPDE 

Class as Compared to the Native FixNVE Class 

 

This appendix appears in its original form, without editorial change.
26

• FixDPDE constructor

o Ensure 3 keywords are given in the fix dpde command

• int FixDPDE::setmask()

o Removed all reference to RESPA

• void FixDPDE::init()

o Ensure that an equationof state is specified in another fix; otherwise, print error message and stop

o Remove all references to RESPA

• void FixDPDE::initial_integrate(int vflag)

o Integrate uCond and uMech by the relations:

uCondi = uCondi + 1⁄2 * dt * duCondi

uMechi = uMechi + 1⁄2 * dt * duMechi

• void FixDPDE::final_integrate()

o Integrated uCond and uMech by the relations:

uCondi = uCondi + 1⁄2 * dt * duCondi

uMechi = uMechi + 1⁄2 * dt * duMechi
27

Appendix E. Summary of the FixEOScv Class 

 

This appendix appears in its original form, without editorial change.
28

• FixEOScv constructor 

o Ensure 4 keywords are specified in the fix eos/cv command

o Read the heat capacity, cv

• int FixEOScv::setmask()

o Set the mask to POST_INTEGRATE and END_OF_STEP

• void FixEOScv::init()

o Initialize the cv array

o Initialize the eos array

o Initialize uCond and uMech to be 1⁄2 * cv[i] * dpdTheta[i]

• void FixEOScv::post_integrate()

o Apply the EOS to compute the DPD particle internal temperature. The eos/cv defines the following relation:

θi = (uCondi + uMechi

) / Cv,i

where θi is the internal particle temperature (dpdTheta). 

• void FixEOScv::end_of_step()

o Apply the EOS to compute the DPD particle internal temperature. The eos/cv defines the following relation:

θi = (uCondi + uMechi

) / Cv,i

 where θi is the internal particle temperature (dpdTheta). 
29

Appendix F. Summary of LAMMPS Code Differences Between the PairDPDE

Conservative Class as Compared to the Native PairDPD Class 

 

This appendix appears in its original form, without editorial change.
30

• PairDPDEConservative destructor: PairDPDE::~PairDPDE()

o Deallocate kappa array; remove gamma deallocation

• void PairDPDEConservative::compute(int eflag, int vflag)

o Removed dissipative and random force calculations

o Define wR and wD: ωD(r) = [ωR(r)]2 = �

�1 − r

rc

�

2

, r < rc

0 , r ≥ rc

where rc is the cutoff radius for the pair of interacting particles.

NOTE: In PairDPD class, wd definition differs

o Compute the conservative force: Fij

C = ADPDωR(r)

o Compute the interaction energy between the particle pair

 evdwl = 0.5*a0[itype][jtype]*cut[itype][jtype] * wd;

• void PairDPDEConservative::allocate()

o Allocate memory for kappa array; Remove gamma allocation

• void PairDPDEConservative::settings(int narg, char **arg)

o Read 3 pair style arguments (kappa_flag, cut_global and seed); Remove the temperature argument

o Ensure kappa_flag is set to either 0 or 1

• void PairDPDEConservative::coeff(int narg, char **arg)

o Reads in 3-4 three pair coefficients

 ADPD 

 Sigma (σij); removed gamma

 Kappa_ij (κij) or Kappa0 (κ0) 

 cutoff (optional)

• void PairDPDEConservative::init_style()

o Ensures that ghost velocities are stored and Newton pair is set to “on” for DPD calculations

• double PairDPDEConservative::init_one(int i, int j)

o Initializes the pair coefficients for all particle pairs; Add kappa, sigmaand remove gamma

• void PairDPDEConservative::write_restart(FILE *fp)

o Updates the pair coefficient information to the restart file; Add kappa, sigma and remove gamma

• void PairDPDEConservative::read_restart(FILE *fp)

o Reads pair coefficients from restart file and broadcast across processors; Add kappa, sigma and remove gamma

• void PairDPDEConservative::write_restart_settings(FILE *fp)

o Remove the temperature variable 

• void PairDPDEConservative::read_restart_settings(FILE *fp)

o Remove temperature variable

• void PairDPDEConservative::write_data(FILE *fp)

o Removed gamma; Added kappa and sigma to list of variable to write to the data file

• void PairDPDEConservative::write_data_all(FILE *fp)

o Removed gamma; Added kappa and sigma to list of variable to write to the data file

• double PairDPDEConservative::single(int i, int j, int itype, int jtype, double rsq, double factor_coul, double factor_dpd, double &fforce)

o changed definition of wd and wr to be consistent with manuscript
31

Appendix G. Summary of LAMMPS Code Differences Between the 

FixDPDEShardlow Class as Compared to the Native FixNVE Class 

 

This appendix appears in its original form, without editorial change.
32

• FixDPDEShardlow Constructor

o Ensure 3 keywords are given in the fix dpde/shardlow command

o Set the size of the fix forward communication to 10 (dvSSA[i][0-2], v[i][0-2], duCond[i], duMech[i], uCond[i], uMech[i])

o Set the size of the fix reverse communication to 5 (dvSSA[i][0-2],duCond[i], duMech[i])

• int FixDPDEShardlow::setmask()

o Remove all references to RESPA

• void FixDPDEShardlow::init()

o Ensure that an equation of state is specified in another fix; otherwise, print error message and stop

o Set the pairDPDE pointer to access the pairDPDEConservative class data structures

• void FixDPDEShardlow::initial_integrate(int vflag)

o Added a shardlow_integrate() function prior to the velocity integration

o Integrate the velocity and position at a half time step (Same as FixNVE)

• void FixDPDEShardlow::final_integrate()

o Integrate the velocities at a half time step (Same as FixNVE)

• void FixDPDEShardlow::shardlow_integrate(): Member function to perform the stochastic integration of velocities

o Compute the length of the bounding box dimensions (bbx, bby, bbz)

 Ensure all bounding box dimensions are larger than interaction cutoff radius; Print error if not satisfied

o Allocate memory for velocity (dvSSA[j][0-2])

o Initialize (dvSSA[j][0-2], duCond[j], duMech[j]) arrays and forward communicate to ensure all processors are initialized

o Allocate memory for the particle pair active interaction regions (jbin)

 Count the number of pairs (count), then allocate active interaction region bins, jbin[count]

 Assign each pair to an active interaction region depending on the coordinates of the pair. One atom will be in 

central box; assign the AIR based on the 2nd atom ,which may be in the central box or a neighboring box

o Loop over the eight active interaction regions (Stages)

 If the pair lies in the current active interaction region, then ...

• Compute the pair separation distance, rij, and check if it is less than the interaction cutoff

o Store the current velocities (vx0i, vy0i, vz0i, vx0j, vy0j, vz0j)

o Compute the velocity difference between i and j: 

j

j

i

i

ij m m

p p

v ← − = vi - vj

o Compute the dot product, dot = rij • vij

o Compute wR and wD: ωD(r) = [ωR(r)]2 = �

�1 − r

rc

�

2

, r < rc

0 , r ≥ rc

 where rc is the cutoff radius for the pair of interacting particles

o Apply EOS to compute the current internal temperature.

o Compute theta_ij: Θij

−1 = 1

2 � 1θi

+ 1

θj

�

o Compute gamma_ij: γij = σij

2

2kBθij

o Obtain a random number ςij from a Gaussian distribution 

o Compute the momentum (velocity) update

2 2

t

r r r

t

ij

ij

ij

R

ij

ij

ij

ij

ij

D ij

i i ij

∆

+ 













 ⋅ ∆

← − r rv

r

p p γ ω σ ω ς
33

2 2

t

r r r

t

ij

ij

ij

R

ij

ij

ij

ij

ij

D ij

j j ij

∆ − 













 ⋅ ∆

← +

r r

v

r

p p γ ω σ ω ς

o Re-Compute the velocity difference between i and j: 

j

j

i

i

ij m m

p p

v ← − = vi - vj

o Re-compute the dot product: dot = rij • vij

o Compute mu_ij: μij = � 1

mi

− 1

mj

�

o Compute the momentum (velocity) update

2 2

2

1 2

t

r

t

r r r

t

t

ij

ij

ij

R

ij

ij

ij

ij

R

ij

ij

ij

ij

ij

ij

ij

D

ij

ij

D

ij

i i

∆

+

















+ ∆ 













 ⋅

+ ∆

∆

← − r r rv

r

p p σ ω ς σ ω ς

μ

γ ω

μ

γ ω

2 2

2

1 2

t

r

t

r r r

t

t

ij

ij

ij

R

ij

ij

ij

ij

R

ij

ij

ij

ij

ij

ij

ij

D

ij

ij

D

ij

j j

∆ −

















+ ∆ 













 ⋅

+ ∆

∆

← +

r r r

v

r

p p σ ω ς σ ω ς

μ

γ ω

μ

γ ω

o Compute a new random number, ςij

q

o Compute kappa_ij: κij = κ0

kB

�

ui+uj

2 �

2

 (for T-dependent model)

o Compute alpha_ij: αij = �2kBκij

o Compute the conductive energy update

u u t t q

ij

Rq

ij

Dq

i j

ij

cond

i

cond

i ∆ + ∆ 















← + − ω α ω ς θ θ

κ

1 1

u u t t q

ij

Rq

ij

Dq

i j

ij

cond

j

cond

j ∆ − ∆ 















← − − ω α ω ς θ θ

κ

1 1

o Compute the mechanical energy update















 ⋅ − ⋅ − ⋅ + ⋅ ← −

j

j j

i

i i

j

j j

i

mech i i

i

mech

i

m

' '

m

' '

m m

u u

2 2 2 2 2

1 p p p p p p p p















 ⋅ − ⋅ − ⋅ + ⋅ ← −

j

j j

i

i i

j

j j

i

mech i i

j

mech

j m

' '

m

' '

m m

u u

2 2 2 2 2

1 p p p p p p p p

• Reverse communicate the velocity updates. 

• Update the velocities and energies (v[j][0-2], uCond[j], uMech[j])

• Reinitialize the arrays to zero (dvSSA[j][0-2], duCond[j], duMech[j])

• Forward communicate the velocity/energy updates

 Delete the memory allocations for jbin and dvSSA

• int FixDPDEShardlow::sort_bin(double rx, double ry, double rz)

o Compare the particle position to the reference points of the bounding box

 Determine whether the point lies inside the bounding box, or outside the bounding box 

 Return a number to indicate the active interaction region

• 0 = Inside bounding box

• 1 = Top

• 2 = Right
34

• 3 = Top-right and Bottom-right

• 4 = Back 

• 5 = Top-back and Bottom-back

• 6 = Left-back and Right-back 

• 7 = Back Corners

• int FixDPDEShardlow::pack_comm(int n, int *list, double *buf, int pbc_flag, int *pbc)

o Packs the buffer for forward communication(dvSSA, v, duCond, duMech, uCond, uMech)

• void FixDPDEShardlow::unpack_comm(int n, int first, double *buf)

o Unpacks the buffer for forward communication (dvSSA, v, duCond, duMech, uCond, uMech)

• int FixDPDEShardlow::pack_reverse_comm(int n, int first, double *buf)

o Packs the buffer for reversecommunication(dvSSA, duCond, duMech)

• void FixDPDEShardlow::unpack_reverse_comm(int n, int *list, double *buf)

o Unpacks the buffer for reversecommunication (dvSSA, duCond, duMech)
NO. OF NO. OF

COPIES ORGANIZATION COPIES ORGANIZATION

35

1 DEFENSE TECHNICAL

(PDF) INFORMATION CTR

DTIC OCA

1 DIRECTOR

(PDF) US ARMY RESEARCH LAB

IMAL HRA

1 DIRECTOR

(PDF) US ARMY RESEARCH LAB

RDRL CIO LL

1 GOVT PRINTG OFC

(PDF) A MALHOTRA

 1 U OF MD-COLLEGE PARK

 (PDF) DEPT MECH ENG

P CHUNG

32 DIR USARL

(PDF) RDRL CIH C

L MUNDAY

RDRL CIH S

J LARENTZOS

RDRL WM

P BAKER

B FORCH

J MCCAULEY

RDRL WML B

I BATYREV

J BRENNAN

S BUNTE

E BYRD

S IZVEKOV

W MATTSON

J MOORE

J MORRIS

R PESCE-RODRIGUEZ

B RICE

R SAUSA

N TRIVEDI

N WEINGARTEN

RDRL WML C

S AUBERT

K MCNESBEY

T PIEHLER

RDRL WML D

R BEYER

RDRL WML H

J NEWILL

RDRL WMM

J ZABINSKI

RDRL WMP C

R BECKER

J CLAYTON

M GREENFIELD

RDRL WMM G

J ANDZELM

T SIRK

Y SLIOZBERG

RDRL WMP G

S KUKUCK

B HOMAN
36

INTENTIONALLY LEFT BLANK.