Thanks! I am trying that. I thought to share the piece of code I modified which is pasted below. I checked each term by term energy and force and it is doing correct.

What I found is that dmax sets maximum distance a particle can move in the line minimization. If a particle needs to move larger than the given dmax value in order to find minimum in line-minimization, it can not, rather it returns with a higher value. I could check this as I have my own minimization code and comparing results of both from lammps and my code. I need lammps to work because I need to run for very large system size. Anyway thanks.

Regards,

Moumita

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

randnum = random->gaussian();

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

f[i][1] += delyfpair;

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

printf(“i=%d %f %f %f\n”, i, f[i][0],f[i][1],f[i][2]);

if (newton_pair || j < nlocal) {

f[j][0] -= delxfpair;

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

f[j][2] -= delzfpair;

printf(“j=%d %f %f %f\n”, j, f[j][0],f[j][1],f[j][2]);

}

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 = wd*wd;

//evdwl *= factor_dpd;

}

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

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

}

}

}