LAMMPS WWW Site - LAMMPS Documentation - LAMMPS Mailing List Archives
[lammps-users] Pressure / Virial calculation

# [lammps-users] Pressure / Virial calculation

 From: Xiao Jia Date: Tue, 29 Aug 2017 17:53:27 -0400

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

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!