Hi,

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

protected:

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)];

j &= NEIGHMASK;

delx = xtmp - x[j][0];

dely = ytmp - x[j][1];

delz = ztmp - x[j][2];

rsq = delx*delx + dely*dely + 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 = delx*delvx + dely*delvy + 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

force.

------------------------------------------------------

*/*

fpair -= (2well_depth[itype][jtype]*bp/widthsq)

fpair -= (2

*exp(-1*pow(bp,2)/widthsq); // SW

fpair -= gamma[itype][jtype]*wd*wd*dot*rinv;

fpair += sigma[itype][jtype]*wd*randnum*dtinvsqrt;

fpair *= factor_dpd*rinv;

f[i][0] += delx*fpair;
f[i][1] += dely*fpair;

f[i][2] += delz*fpair;

if (newton_pair || j < nlocal) {

f[j][0] -= delx*fpair;
f[j][1] -= dely*fpair;

f[j][2] -= delz*fpair;

}

if (eflag) {

// unshifted eng of conservative term:

// evdwl = -a0[itype][jtype]*r * (1.0-0.5*r/cut[itype][jtype]);

// eng shifted to 0.0 at cutoff

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

evdwl *= factor_dpd;

}

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

evdwl,0.0,fpair,delx,dely,delz);

}

}

}

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;

force->bounds(arg[0],atom->ntypes,ilo,ihi);

force->bounds(arg[1],atom->ntypes,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;

count++;

}

}

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 <http://gnu.org/licenses/gpl.html>

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:

<http://bugs.launchpad.net/gdb-linaro/>…

Reading symbols from /home/steve/verify_mod_lammps/cyberstar/lmp_ubuntu_mod…done.

(gdb) run < water.in

Starting program: /home/steve/verify_mod_lammps/cyberstar/lmp_ubuntu_mod < water.in

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

Steve

pair_dpd_mod.h (2.84 KB)

pair_dpd_mod.cpp (13.1 KB)