# Pressure / Virial calculation

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 <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 = 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[0] = delxdelxfpair;
vir[1] = delydelyfpair;
vir[2] = delzdelzfpair;

}
}
}

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!

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.​

Dear Axel,

Thanks a lot for your detailed explanation.Based on your explanation, I think I can keep the current code for the virial part (no need to half that) for the calculation of pressure.
Could I have another question is that I want to calculate the pressure distribution within the system, not the overall pressure as reported from .log file in lammps.
What I am doing is keeping going calcualtion from i,j pairs, and get the virial value in three dimensions.
When location of atom i is within a certain smaller bin, I will count the virial from i,j pair into the overal virial in this bin.
My question is this method processed properly? and Do I need to require atom j also within the bin?
Additionally, I guess the atoms in ilist are from local atoms, and atoms in jlist can be from both local and ghost atoms,
do I need to apply the condition that (if (j<nlocal), something like this) to count in the virial from a certain pair?

Thanks!

Dear Axel,

Thanks a lot for your detailed explanation.
Based on your explanation, I think I can keep the current code for the
virial part (no need to half that) for the calculation of pressure.
Could I have another question is that I want to calculate the pressure
distribution within the system, not the overall pressure as reported from
.log file in lammps.
What I am doing is keeping going calcualtion from i,j pairs, and get the
virial value in three dimensions.
When location of atom i is within a certain smaller bin, I will count the
virial from i,j pair into the overal virial in this bin.
My question is this method processed properly? and Do I need to require
atom j also within the bin?
Additionally, I guess the atoms in ilist are from local atoms, and atoms
in jlist can be from both local and ghost atoms,
do I need to apply the condition that (if (j<nlocal), something like this)
to count in the virial from a certain pair?

i think ​this is beyond the scope of this mailing list and - more
importantly - beyond my comfort level to give specific ​advice, since i
don't have much experience with this. there are several ways to tally
stress into a spatial grid, specifically if you also want to include time
averaging.
the simplest approach i can imagine would be to compute stress/atom and
then set up a 3d (or 2d or 1d) binning using compute chunk/atom and sum
those with fix ave/chunk. since you know the volume of each bin, you can
then compute the pressure distribution.
other approaches are coded in the AtC package. that is worth taking a look
as well.

i also recall seeing some code and a poster at the LAMMPS workshop about 4
years ago with more elaborate spatial stress/pressure distribution