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 <Firi>.
For the virial part, I try to use pairwise interatomic <Fijrij> 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];
ytmp = x[i];
ztmp = x[i];
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];
dely = ytmp - x[j];
delz = ztmp - x[j];
rsq = delxdelx + delydely + 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 = delxdelxfpair;
vir = delydelyfpair;
vir = delzdelzfpair;
The virial part will be summed up from vir to vir 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.