[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

From: |
Xiao Jia <xiaoj@...1508...> |

Date: |
Tue, 29 Aug 2017 19:11:36 -0400 |

Dear Axel,

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!

On Tue, Aug 29, 2017 at 6:31 PM, Axel Kohlmeyer <akohlmey@...43...4...> wrote:

On Tue, Aug 29, 2017 at 5:53 PM, Xiao Jia <xiaoj@...1508...> 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 isinum = 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 nomax neighbors/atom: 2000, page size: 100000master list distance cutoff = 2.8ghost atom cutoff = 2.8binsize = 1.4, bins = 12 12 121 neighbor lists, perpetual/occasional/extra = 1 0 0(1) pair lj/cut, perpetualattributes: half, newton onpair build: half/bin/atomonly/newtonstencil: half/bin/3d/newtonbin: standardthe same thing with a full neighbor list looks like this:Neighbor list info ...update every 20 steps, delay 0 steps, check nomax neighbors/atom: 2000, page size: 100000master list distance cutoff = 2.8ghost atom cutoff = 2.8binsize = 1.4, bins = 12 12 121 neighbor lists, perpetual/occasional/extra = 1 0 0(1) pair lj/cut/full, perpetualattributes: full, newton offpair build: full/bin/atomonlystencil: full/bin/3dbin: standardin 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@...655....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.

**Follow-Ups**:**Re: [lammps-users] Pressure / Virial calculation***From:*Axel Kohlmeyer <akohlmey@...24...>

**References**:**[lammps-users] Pressure / Virial calculation***From:*Xiao Jia <xiaoj@...1508...>

**Re: [lammps-users] Pressure / Virial calculation***From:*Axel Kohlmeyer <akohlmey@...24...>

- Prev by Date:
**Re: [lammps-users] Segmentation fault when using Kolmogorov/Crespi/z with REBO** - Next by Date:
**Re: [lammps-users] Pressure / Virial calculation** - Previous by thread:
**Re: [lammps-users] Pressure / Virial calculation** - Next by thread:
**Re: [lammps-users] Pressure / Virial calculation** - Index(es):