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

Re: [lammps-users] Pressure / Virial calculation


From: Axel Kohlmeyer <akohlmey@...24...>
Date: Tue, 29 Aug 2017 18:31:21 -0400



On Tue, Aug 29, 2017 at 5:53 PM, Xiao Jia <xiaoj@...1722...08...> 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];
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.​

 

Thanks!

------------------------------------------------------------------------------
Check out the vibrant tech community on one of the world's most
engaging tech sites, Slashdot.org! http://sdm.link/slashdot
_______________________________________________
lammps-users mailing list
lammps-users@...6297....sourceforge.net
https://lists.sourceforge.net/lists/listinfo/lammps-users




--
Dr. Axel Kohlmeyer  akohlmey@...24...  http://goo.gl/1wk0
College of Science & Technology, Temple University, Philadelphia PA, USA
International Centre for Theoretical Physics, Trieste. Italy.