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 = 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
force.
------------------------------------------------------/
fpair -= (2well_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] += delyfpair;
f[i][2] += delz*fpair;
if (newton_pair || j < nlocal) {
f[j][0] -= delxfpair;
f[j][1] -= delyfpair;
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,
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:
steve@…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)
