Issues with implementing a modified dpd pair style


I am trying to implement a modified version of the DPD pair style to include a Gaussian well component in the conservative force (as described in Magee and Siperstein’s article “Formation of Ordered Mesoporous Materials under Slow Aggregation Conditions”). The modified conservative force requires three additional pairwise variables to be defined: a well width parameter, a well center parameter, and a well depth parameter.

In an effort to implement this simple force update, I copied the available source code for the pair_dpd.cpp/h files and basically defined my three new parameters similarly to the other DPD variables:

From pair_dpd_mod.h:
// Code snippet

double cut_global,temperature;
int seed;
double **cut;
double **a0,**gamma;
double **sigma;
double **well_width; // SW
double **well_center; // SW
double **well_depth; // SW

class RanMars *random;

void allocate();

// End code snippet

Then I added the appropriate lines for the constructor and destructor, and added in my force update shown below in the compute member function:

From pair_dpd_mod.cpp
// Code snippet

void PairDPDmod::compute(int eflag, int vflag) // SW
int i,j,ii,jj,inum,jnum,itype,jtype;
double xtmp,ytmp,ztmp,delx,dely,delz,evdwl,fpair;
double vxtmp,vytmp,vztmp,delvx,delvy,delvz;
double rsq,r,rinv,dot,wd,randnum,factor_dpd,bp,widthsq; // SW
int *ilist,*jlist,*numneigh,**firstneigh;

evdwl = 0.0;
if (eflag || vflag) ev_setup(eflag,vflag);
else evflag = vflag_fdotr = 0;

double **x = atom->x;
double **v = atom->v;
double **f = atom->f;
int *type = atom->type;
int nlocal = atom->nlocal;
double *special_lj = force->special_lj;
int newton_pair = force->newton_pair;
double dtinvsqrt = 1.0/sqrt(update->dt);

inum = list->inum;
ilist = list->ilist;
numneigh = list->numneigh;
firstneigh = list->firstneigh;

// loop over neighbors of my atoms

for (ii = 0; ii < inum; ii++) {
i = ilist[ii];
xtmp = x[i][0];
ytmp = x[i][1];
ztmp = x[i][2];
vxtmp = v[i][0];
vytmp = v[i][1];

vztmp = v[i][2];
itype = type[i];
jlist = firstneigh[i];
jnum = numneigh[i];

for (jj = 0; jj < jnum; jj++) {
j = jlist[jj];
factor_dpd = special_lj[sbmask(j)];

delx = xtmp - x[j][0];
dely = ytmp - x[j][1];
delz = ztmp - x[j][2];
rsq = delxdelx + delydely + delz*delz;
jtype = type[j];

if (rsq < cutsq[itype][jtype]) {
r = sqrt(rsq);
if (r < EPSILON) continue; // r can be 0.0 in DPD systems
rinv = 1.0/r;
delvx = vxtmp - v[j][0];
delvy = vytmp - v[j][1];
delvz = vztmp - v[j][2];
dot = delxdelvx + delydelvy + delz*delvz;
wd = 1.0 - r/cut[itype][jtype];
bp = r - well_center[itype][jtype]; // SW added line

widthsq = pow(well_width[itype][jtype],2); // SW added line
randnum = random->gaussian();

// conservative force = a0 * wd
// drag force = -gamma * wd^2 * (delx dot delv) / r
// random force = sigma * wd * rnd * dtinvsqrt;

fpair = a0[itype][jtype]wd;
SW: The next update to fpair establishes the
gaussian well component of the conservative
fpair -= (2
well_depth[itype][jtype]*bp/widthsq)exp(-1pow(bp,2)/widthsq); // SW

fpair -= gamma[itype][jtype]wdwddotrinv;
fpair += sigma[itype][jtype]wdrandnum*dtinvsqrt;
fpair = factor_dpdrinv;

f[i][0] += delxfpair;
f[i][1] += dely
f[i][2] += delz*fpair;
if (newton_pair || j < nlocal) {

f[j][0] -= delxfpair;
f[j][1] -= dely
f[j][2] -= delz*fpair;

if (eflag) {
// unshifted eng of conservative term:
// evdwl = -a0[itype][jtype]r * (1.0-0.5r/cut[itype][jtype]);
// eng shifted to 0.0 at cutoff
evdwl = 0.5*a0[itype][jtype]cut[itype][jtype] * wdwd;
evdwl *= factor_dpd;

if (evflag) ev_tally(i,j,nlocal,newton_pair,

if (vflag_fdotr) virial_fdotr_compute();

// End code snippet

In an effort to verify that the modified code behaves the same as the original code, I wanted to compare data produced when the well parameters are all set to 0.0. My well parameters were passed similarly to the program from the input file as the other DPD parameters as shown below:

From pair_dpd_mod.cpp
// Code Snippet

void PairDPDmod::coeff(int narg, char **arg) // SW
if (narg < 4 || narg > 8) error->all(FLERR,“Incorrect args for pair coefficients”); // SW 8 replaced 5
if (!allocated) allocate();

int ilo,ihi,jlo,jhi;

double a0_one = force->numeric(arg[2]);
double gamma_one = force->numeric(arg[3]);
double cut_one = force->numeric(arg[4]); // SW now cut-off MUST be defined

double well_width_one = force->numeric(arg[5]); // SW
double well_center_one = force->numeric(arg[6]); // SW
double well_depth_one = force->numeric(arg[7]); // SW

int count = 0;
for (int i = ilo; i <= ihi; i++) {
for (int j = MAX(jlo,i); j <= jhi; j++) {
a0[i][j] = a0_one;
gamma[i][j] = gamma_one;
cut[i][j] = cut_one;

well_width[i][j] = well_width_one; // SW
well_center[i][j] = well_center_one; // SW
well_depth[i][j] = well_depth_one; // SW

setflag[i][j] = 1;

if (count == 0) error->all(FLERR,“Incorrect args for pair coefficients”);


// End code snippet

The issue that I am having is that (running with a fix NVT) the initial pressure cannot be calculated and a number of the other components fail to be updated and appear in the output as -nan.
Furthermore, I was getting a seg fault which gdb pointed out to be due to my fix rdf line in my input file. It is confusing though because using the original dpd pair style does not produce this seg fault.

From a gdb session:

[email protected]…3899… ~/verify_mod_lammps/cyberstar » gdb lmp_ubuntu_mod
GNU gdb (Ubuntu/Linaro 7.3-0ubuntu2) 7.3-2011.08
Copyright © 2011 Free Software Foundation, Inc.
License GPLv3+: GNU GPL version 3 or later <>
This is free software: you are free to change and redistribute it.
There is NO WARRANTY, to the extent permitted by law. Type “show copying”
and “show warranty” for details.
This GDB was configured as “x86_64-linux-gnu”.
For bug reporting instructions, please see:
Reading symbols from /home/steve/verify_mod_lammps/cyberstar/lmp_ubuntu_mod…done.
(gdb) run <
Starting program: /home/steve/verify_mod_lammps/cyberstar/lmp_ubuntu_mod <
[Thread debugging using libthread_db enabled]
LAMMPS (19 Dec 2011)
Reading data file …
orthogonal box = (0 0 0) to (20 20 20)
1 by 1 by 1 MPI processor grid
24000 atoms
Setting up run …
Memory usage per processor = 51.5927 Mbytes
Step Temp PotEng KinEng TotEng Press Pxx Pyy Pzz Volume Lx Ly Lz
0 0 7.8431557 0 7.8431557 -nan -nan -nan -nan 8000 20 20 20
1 -nan 0 -nan -nan -nan -nan -nan -nan 8000 20 20 20
2 -nan 0 -nan -nan -nan -nan -nan -nan 8000 20 20 20
3 -nan 0 -nan -nan -nan -nan -nan -nan 8000 20 20 20
4 -nan 0 -nan -nan -nan -nan -nan -nan 8000 20 20 20
5 -nan 0 -nan -nan -nan -nan -nan -nan 8000 20 20 20


95 -nan 0 -nan -nan -nan -nan -nan -nan 8000 20 20 20
96 -nan 0 -nan -nan -nan -nan -nan -nan 8000 20 20 20
97 -nan 0 -nan -nan -nan -nan -nan -nan 8000 20 20 20
98 -nan 0 -nan -nan -nan -nan -nan -nan 8000 20 20 20
99 -nan 0 -nan -nan -nan -nan -nan -nan 8000 20 20 20

Program received signal SIGSEGV, Segmentation fault.
LAMMPS_NS::ComputeRDF::compute_array (this=0xa56a20) at compute_rdf.cpp:254
254 compute_rdf.cpp: No such file or directory.
in compute_rdf.cpp
(gdb) quit

Commenting out the fix rdf line in the input file solved the seg fault issue, but I believe the issue lies somewhere in my modified pair style code. Moreover, without the fix rdf line, the same behavior that led to the -nan output still occured.

I was just wondering if you had any input on what may be stopping the system variables from updating and why my code is producing the seg fault and why the original file is not.

I am new to the forum and I am unsure if I can send attachments, but I will try to send my two modified source files regardless.

Any feedback would be greatly appreciated.

Thank you,


pair_dpd_mod.h (2.84 KB)

pair_dpd_mod.cpp (13.1 KB)

If you debug with valgrind, you will be able to see where
you use a bad or undefined value for the first time.

The new 2d arrays you defined (wellwidth, etc) need
to be defined for all possible I,J type pairs. That is done in init_one(),
to set J,I = I,J. I presume you did that for your variables?

When you have it working, send us the files for the new pair style and
we can release
them in the main distro.


I was looking into the DPD code as well to see how difficult it would be to implement a couple of different weighting function. The reason for doing something like that can be found here for example:

Would there be any interest to get something like this working? I will probably take this own for my own work but I have to sort out c++ first (sometimes it looks like ancient Egyptian Hieroglyphs to me)


I was looking into the DPD code as well to see how difficult it would be to
implement a couple of different weighting function. The reason for doing
something like that can be found here for example:

Would there be any interest to get something like this working? I will

probably take this own for my own work but I have to sort out c++ first
(sometimes it looks like ancient Egyptian Hieroglyphs to me)

have you seen a program in APL. that *is* written with hieroglyphs. :wink:


Hi Steve,

I did debug with valgrind but I am still having some trouble figuring out what’s going on.

Before I spend more time trying to debug it, I realized that I didn’t update the “single” member function in the pair style. In pair_dpd.cpp it appears that single is calculating a force “fforce” that looks similar to the conservative force and another parameter “phi”. I am not quite sure what is going on here in the code, what single is used for and if it is necessary at all for the general operation of the code. I was just wondering if you could help provide some insight into what this function is doing.

Could this be responsible for why my system pressures and energies are being printed as “-nan” when I test the modified pair style? Could this furthermore be responsible for why the program has a segmentation fault when I use compute RDF?


  • Steve

pair::single() is used by only a couple commands, like compute group/group
to calculate pairwise interactions on a one by one basis, rather than
for the entire set of atoms. Eventually your new pair style should
have a valid single() method, but it's unlikely you are even invoking it
in your script. A print statement added to single() would verify that.