From 569136e1604d9637a2e1668dadc05f50d3e3a1c1 Mon Sep 17 00:00:00 2001 From: jtclemm Date: Wed, 6 Jul 2022 14:35:16 -0600 Subject: [PATCH 1/4] Redoing groups in compute contact/atom --- doc/src/compute_contact_atom.rst | 10 +++++- src/GRANULAR/compute_contact_atom.cpp | 52 ++++++++++++++++----------- src/GRANULAR/compute_contact_atom.h | 4 +++ 3 files changed, 44 insertions(+), 22 deletions(-) diff --git a/doc/src/compute_contact_atom.rst b/doc/src/compute_contact_atom.rst index 20dcbfae29..23b5f4639d 100644 --- a/doc/src/compute_contact_atom.rst +++ b/doc/src/compute_contact_atom.rst @@ -8,10 +8,11 @@ Syntax .. parsed-literal:: - compute ID group-ID contact/atom + compute ID group-ID contact/atom group2-ID * ID, group-ID are documented in :doc:`compute ` command * contact/atom = style name of this compute command +* group2-ID = optional argument to restrict which atoms to consider for contacts (see below) Examples """""""" @@ -19,6 +20,7 @@ Examples .. code-block:: LAMMPS compute 1 all contact/atom + compute 1 all contact/atom mygroup Description """"""""""" @@ -45,6 +47,9 @@ overview of LAMMPS output options. The per-atom vector values will be a number >= 0.0, as explained above. +The optional *group2-ID* argument allows to specify from which group atoms +contribute to the coordination number. Default setting is group 'all'. + Restrictions """""""""""" @@ -63,4 +68,7 @@ Related commands Default """"""" +*group2-ID* = all + + none diff --git a/src/GRANULAR/compute_contact_atom.cpp b/src/GRANULAR/compute_contact_atom.cpp index 5a07d14eb8..26018f7bfb 100644 --- a/src/GRANULAR/compute_contact_atom.cpp +++ b/src/GRANULAR/compute_contact_atom.cpp @@ -18,6 +18,7 @@ #include "comm.h" #include "error.h" #include "force.h" +#include "group.h" #include "memory.h" #include "modify.h" #include "neigh_list.h" @@ -34,7 +35,16 @@ ComputeContactAtom::ComputeContactAtom(LAMMPS *lmp, int narg, char **arg) : Compute(lmp, narg, arg), contact(nullptr) { - if (narg != 3) error->all(FLERR,"Illegal compute contact/atom command"); + if (narg != 3 && narg != 4) error->all(FLERR,"Illegal compute contact/atom command"); + + jgroup = group->find("all"); + jgroupbit = group->bitmask[jgroup]; + if (narg == 4) { + group2 = utils::strdup(arg[3]); + jgroup = group->find(group2); + if (jgroup == -1) error->all(FLERR, "Compute contact/atom group2 ID does not exist"); + jgroupbit = group->bitmask[jgroup]; + } peratom_flag = 1; size_peratom_cols = 0; @@ -120,28 +130,28 @@ void ComputeContactAtom::compute_peratom() for (ii = 0; ii < inum; ii++) { i = ilist[ii]; - if (mask[i] & groupbit) { - xtmp = x[i][0]; - ytmp = x[i][1]; - ztmp = x[i][2]; - radi = radius[i]; - jlist = firstneigh[i]; - jnum = numneigh[i]; - for (jj = 0; jj < jnum; jj++) { - j = jlist[jj]; - j &= NEIGHMASK; + xtmp = x[i][0]; + ytmp = x[i][1]; + ztmp = x[i][2]; + radi = radius[i]; + jlist = firstneigh[i]; + jnum = numneigh[i]; - delx = xtmp - x[j][0]; - dely = ytmp - x[j][1]; - delz = ztmp - x[j][2]; - rsq = delx*delx + dely*dely + delz*delz; - radsum = radi + radius[j]; - radsumsq = radsum*radsum; - if (rsq <= radsumsq) { - contact[i] += 1.0; - contact[j] += 1.0; - } + for (jj = 0; jj < jnum; jj++) { + j = jlist[jj]; + j &= NEIGHMASK; + + delx = xtmp - x[j][0]; + dely = ytmp - x[j][1]; + delz = ztmp - x[j][2]; + rsq = delx * delx + dely * dely + delz * delz; + radsum = radi + radius[j]; + radsumsq = radsum * radsum; + if (rsq <= radsumsq) { + // Only tally for atoms in compute group (groupbit) if neighbor is in group2 (jgroupbit) + if ((mask[i] & groupbit) && (mask[j] & jgroupbit)) contact[i] += 1.0; + if ((mask[j] & groupbit) && (mask[i] & jgroupbit)) contact[j] += 1.0; } } } diff --git a/src/GRANULAR/compute_contact_atom.h b/src/GRANULAR/compute_contact_atom.h index f01815e666..04157b5455 100644 --- a/src/GRANULAR/compute_contact_atom.h +++ b/src/GRANULAR/compute_contact_atom.h @@ -37,6 +37,10 @@ class ComputeContactAtom : public Compute { private: int nmax; + + char *group2; + int jgroup, jgroupbit; + class NeighList *list; double *contact; }; From 5496781ecd58b0c59f60bc9d0326128107ae4253 Mon Sep 17 00:00:00 2001 From: jtclemm Date: Thu, 7 Jul 2022 12:44:36 -0600 Subject: [PATCH 2/4] Optimizing group filters --- src/GRANULAR/compute_contact_atom.cpp | 45 ++++++++++++++++----------- 1 file changed, 27 insertions(+), 18 deletions(-) diff --git a/src/GRANULAR/compute_contact_atom.cpp b/src/GRANULAR/compute_contact_atom.cpp index 26018f7bfb..a502630fd4 100644 --- a/src/GRANULAR/compute_contact_atom.cpp +++ b/src/GRANULAR/compute_contact_atom.cpp @@ -125,33 +125,42 @@ void ComputeContactAtom::compute_peratom() int *mask = atom->mask; int nlocal = atom->nlocal; int nall = nlocal + atom->nghost; + bool update_i_flag, update_j_flag; for (i = 0; i < nall; i++) contact[i] = 0.0; for (ii = 0; ii < inum; ii++) { i = ilist[ii]; - xtmp = x[i][0]; - ytmp = x[i][1]; - ztmp = x[i][2]; - radi = radius[i]; - jlist = firstneigh[i]; - jnum = numneigh[i]; + // Only proceed if i is either part of the compute group or will contribute to contacts + if ((mask[i] & groupbit) || (mask[i] & jgroupbit)) { + xtmp = x[i][0]; + ytmp = x[i][1]; + ztmp = x[i][2]; + radi = radius[i]; + jlist = firstneigh[i]; + jnum = numneigh[i]; - for (jj = 0; jj < jnum; jj++) { - j = jlist[jj]; - j &= NEIGHMASK; + for (jj = 0; jj < jnum; jj++) { + j = jlist[jj]; + j &= NEIGHMASK; - delx = xtmp - x[j][0]; - dely = ytmp - x[j][1]; - delz = ztmp - x[j][2]; - rsq = delx * delx + dely * dely + delz * delz; - radsum = radi + radius[j]; - radsumsq = radsum * radsum; - if (rsq <= radsumsq) { // Only tally for atoms in compute group (groupbit) if neighbor is in group2 (jgroupbit) - if ((mask[i] & groupbit) && (mask[j] & jgroupbit)) contact[i] += 1.0; - if ((mask[j] & groupbit) && (mask[i] & jgroupbit)) contact[j] += 1.0; + update_i_flag = (mask[i] & groupbit) && (mask[j] & jgroupbit); + update_j_flag = (mask[j] & groupbit) && (mask[i] & jgroupbit); + + if (update_i_flag || update_j_flag) { + delx = xtmp - x[j][0]; + dely = ytmp - x[j][1]; + delz = ztmp - x[j][2]; + rsq = delx * delx + dely * dely + delz * delz; + radsum = radi + radius[j]; + radsumsq = radsum * radsum; + if (rsq <= radsumsq) { + if (update_i_flag) contact[i] += 1.0; + if (update_j_flag) contact[j] += 1.0; + } + } } } } From 4c7825e5c7ad29dec09a5de51064ed5ea1d6d9bc Mon Sep 17 00:00:00 2001 From: jtclemm Date: Fri, 15 Jul 2022 09:20:07 -0600 Subject: [PATCH 3/4] Inverting if statements --- src/GRANULAR/compute_contact_atom.cpp | 49 +++++++++++++-------------- 1 file changed, 24 insertions(+), 25 deletions(-) diff --git a/src/GRANULAR/compute_contact_atom.cpp b/src/GRANULAR/compute_contact_atom.cpp index a502630fd4..91511f57ec 100644 --- a/src/GRANULAR/compute_contact_atom.cpp +++ b/src/GRANULAR/compute_contact_atom.cpp @@ -133,34 +133,33 @@ void ComputeContactAtom::compute_peratom() i = ilist[ii]; // Only proceed if i is either part of the compute group or will contribute to contacts - if ((mask[i] & groupbit) || (mask[i] & jgroupbit)) { - xtmp = x[i][0]; - ytmp = x[i][1]; - ztmp = x[i][2]; - radi = radius[i]; - jlist = firstneigh[i]; - jnum = numneigh[i]; + if (! (mask[i] & groupbit) && ! (mask[i] & jgroupbit)) continue; - for (jj = 0; jj < jnum; jj++) { - j = jlist[jj]; - j &= NEIGHMASK; + xtmp = x[i][0]; + ytmp = x[i][1]; + ztmp = x[i][2]; + radi = radius[i]; + jlist = firstneigh[i]; + jnum = numneigh[i]; - // Only tally for atoms in compute group (groupbit) if neighbor is in group2 (jgroupbit) - update_i_flag = (mask[i] & groupbit) && (mask[j] & jgroupbit); - update_j_flag = (mask[j] & groupbit) && (mask[i] & jgroupbit); + for (jj = 0; jj < jnum; jj++) { + j = jlist[jj]; + j &= NEIGHMASK; - if (update_i_flag || update_j_flag) { - delx = xtmp - x[j][0]; - dely = ytmp - x[j][1]; - delz = ztmp - x[j][2]; - rsq = delx * delx + dely * dely + delz * delz; - radsum = radi + radius[j]; - radsumsq = radsum * radsum; - if (rsq <= radsumsq) { - if (update_i_flag) contact[i] += 1.0; - if (update_j_flag) contact[j] += 1.0; - } - } + // Only tally for atoms in compute group (groupbit) if neighbor is in group2 (jgroupbit) + update_i_flag = (mask[i] & groupbit) && (mask[j] & jgroupbit); + update_j_flag = (mask[j] & groupbit) && (mask[i] & jgroupbit); + if (! update_i_flag && ! update_j_flag) continue; + + delx = xtmp - x[j][0]; + dely = ytmp - x[j][1]; + delz = ztmp - x[j][2]; + rsq = delx * delx + dely * dely + delz * delz; + radsum = radi + radius[j]; + radsumsq = radsum * radsum; + if (rsq <= radsumsq) { + if (update_i_flag) contact[i] += 1.0; + if (update_j_flag) contact[j] += 1.0; } } } From 8b1734d4c3db10aab43e41576c4bdc834aba1fc5 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Fri, 15 Jul 2022 12:18:15 -0400 Subject: [PATCH 4/4] whitespace --- src/GRANULAR/compute_contact_atom.h | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/GRANULAR/compute_contact_atom.h b/src/GRANULAR/compute_contact_atom.h index 04157b5455..5262c1b29b 100644 --- a/src/GRANULAR/compute_contact_atom.h +++ b/src/GRANULAR/compute_contact_atom.h @@ -37,10 +37,10 @@ class ComputeContactAtom : public Compute { private: int nmax; - + char *group2; - int jgroup, jgroupbit; - + int jgroup, jgroupbit; + class NeighList *list; double *contact; };