Hi:

I’m trying to modify the DPD package (since it seems like one of the two pair_styles that communicates velocities) and use this modified dpd package along with LJ using pair_style hybrid/overlay

Reading the pair_dpd.cpp file brings me to this section

…

…

dot = delx*delvx + dely*delvy + delz*delvz;

wd = 1.0 - r/cut[itype][jtype];

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;
fpair -= gamma[itype][jtype] wdwd*dot

*rinv;*

fpair += sigma[itype][jtype]dtinvsqrt;

fpair += sigma[itype][jtype]

*wd*randnumfpair

*= 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] -= delxfpair;

if (newton_pair || j < nlocal) {

f[j][0] -= delx

f[j][1] -= dely

*fpair;*

f[j][2] -= delzfpair;

f[j][2] -= delz

…

…

which seems like the force on each atom being decremented according to the DPD formulation. I change this appropriately.

The next bit, however, I’m not so sure about. The name evdwl seems to suggest a van-der-waals energy, but I’m really not sure what the variable contains. Could someone please explain? Setting it to zero does not cause lammps to crash, but if I compute the pairwise energy using

compute x all pe/atom

the values it produces seem strange.

Many thanks.

if (eflag) {

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

evdwl *= factor_dpd;

}

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

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

}

Many thanks.