Pressure / Virial calculation

Xiao_Jia · August 29, 2017, 9:53pm

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][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 = 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!

akohlmey · August 29, 2017, 10:31pm

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.

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.

Xiao_Jia · August 29, 2017, 11:11pm

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!

akohlmey · August 30, 2017, 12:56am

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
computation for molecular systems. you should look for papers published by
Takenobu Nakamura. there are likely more publications on that kind of
topic. i think comparing your plans with what other people have done, that
have spent more time on this than somebody like me, is going to be much
more helpful and informative.

axel.

Xiao_Jia · August 30, 2017, 1:32am

Dear Axel,

Thanks for your kind advice.I will follow your recommended materials for more information about this problem.

Thanks!