LAMMPS WWW Site - LAMMPS Documentation - LAMMPS Mailing List Archives
[lammps-users] Pressure / Virial calculation
[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

[lammps-users] Pressure / Virial calculation


From: Xiao Jia <xiaoj@...1508...>
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];
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!