Dear LAMMPS users,

I am working on the pressure calculation of nickel with morse potential.

I am considering two terms in the mathematical expression, including thermal part (nkT) and the virial part <Fi*ri>.
For the virial part, I try to use pairwise interatomic <Fij*rij> to calculate the results.

I tried to borrow the code to compute the interatomic force within pair_morse.cpp for this calculation, which is

inum = force->pair->list->inum;

ilist = force->pair->list->ilist;

numneigh = force->pair->list->numneigh;

firstneigh = force->pair->list->firstneigh;

for (ii = 0; ii < inum; ii++)

{

i = ilist[ii];

xtmp = x[i][0];

ytmp = x[i][1];

ztmp = x[i][2];

itype = type[i];

jlist = firstneigh[i];

jnum = numneigh[i];

for (jj = 0; jj < jnum; jj++)

{

j = jlist[jj];

factor_lj = 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);

dr = r - r0[itype][jtype];

dexp = exp(-alpha[itype][jtype] * dr);

fpair = factor_lj * morse1[itype][jtype] * (dexp*dexp - dexp) / r;

vir[0] = delx*delx*fpair;

vir[1] = dely*dely*fpair;

vir[2] = delz*delz*fpair;

}

}

}

The virial part will be summed up from vir[0] to vir[3] and averaged.

I can understand that for every atom i, the neighbor list is stored in j, so that in this manner every pair of interatomic interaction can be calculated.

But I am not sure about how this neighbor list is generated.

Actually if j can be i, and i becomes j later on, there will be double counting of virial between i and j.

This introduces my question that whether this calculation double counts the virial, and do I need to add the 0.5 factor in front of every atom virial.

Thanks!