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 = delxdelvx + delydelvy + 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]wdwddotrinv;
fpair += sigma[itype][jtype]wdrandnumdtinvsqrt;
fpair = factor_dpdrinv;
f[i][0] += delxfpair;
f[i][1] += delyfpair;
f[i][2] += delzfpair;
if (newton_pair || j < nlocal) {
f[j][0] -= delxfpair;
f[j][1] -= delyfpair;
f[j][2] -= delzfpair;
…
…
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.