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

# Re: [lammps-users] Pressure / Virial calculation

 From: Axel Kohlmeyer Date: Tue, 29 Aug 2017 18:31:21 -0400

On Tue, Aug 29, 2017 at 5:53 PM, Xiao Jia wrote:
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.

​the morse pair style will request a "half" neighbor list. where i,j pairs (for local atoms, for pairs where one atom is a ghost atom, it depends on the newton setting) are only present once.​  you can see this reported at the beginning of a run (here from lj/cut):

Neighbor list info ...
update every 20 steps, delay 0 steps, check no
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 2.8
ghost atom cutoff = 2.8
binsize = 1.4, bins = 12 12 12
1 neighbor lists, perpetual/occasional/extra = 1 0 0
(1) pair lj/cut, perpetual
attributes: half, newton on
pair build: half/bin/atomonly/newton
stencil: half/bin/3d/newton
bin: standard

​the same thing with a full neighbor list looks like this:

Neighbor list info ...
update every 20 steps, delay 0 steps, check no
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 2.8
ghost atom cutoff = 2.8
binsize = 1.4, bins = 12 12 12
1 neighbor lists, perpetual/occasional/extra = 1 0 0
(1) pair lj/cut/full, perpetual
attributes: full, newton off
pair build: full/bin/atomonly
stencil: full/bin/3d
bin: standard

in that case, the virial tally is indeed computed differently to avoid the double counting:​

if (evflag) ev_tally(i,j,nlocal,/* newton_pair */ 1,
0.5*evdwl,0.0,0.5*fpair,delx,dely,delz);

​axel.​

Thanks!

------------------------------------------------------------------------------
Check out the vibrant tech community on one of the world's most