Merge pull request #1896 from pmla/user-ptm-group-fix

Apply compute group correctly for compute ptm/atom

doc/src/compute_ptm_atom.rst

@@ -9,12 +9,13 @@ Syntax
 
 .. parsed-literal::
 
-   compute ID group-ID ptm/atom structures threshold
+   compute ID group-ID ptm/atom structures threshold group2-ID
 
 * ID, group-ID are documented in :doc:`compute <compute>` command
 * ptm/atom = style name of this compute command
 * structures = structure types to search for
 * threshold = lattice distortion threshold (RMSD)
+* group2-ID determines which group is used for neighbor selection (optional, default "all")
 
 Examples
 """"""""
@@ -22,8 +23,8 @@ Examples
 
 .. parsed-literal::
 
-   compute 1 all ptm/atom default 0.1
-   compute 1 all ptm/atom fcc-hcp-dcub-dhex 0.15
+   compute 1 all ptm/atom default 0.1 all
+   compute 1 all ptm/atom fcc-hcp-dcub-dhex 0.15 all
    compute 1 all ptm/atom all 0
 
 Description
@@ -82,7 +83,9 @@ The neighbor list needed to compute this quantity is constructed each
 time the calculation is performed (e.g. each time a snapshot of atoms
 is dumped).  Thus it can be inefficient to compute/dump this quantity
 too frequently or to have multiple compute/dump commands, each with a
-*ptm/atom* style.
+*ptm/atom* style. By default the compute processes **all** neighbors
+unless the optional *group2-ID* argument is given, then only members
+of that group are considered as neighbors.
 
 **Output info:**
 
@@ -106,6 +109,7 @@ The interatomic distance is computed from the scale factor in the RMSD equation.
 The (qw,qx,qy,qz) parameters represent the orientation of the local structure
 in quaternion form.  The reference coordinates for each template (from which the
 orientation is determined) can be found in the *ptm\_constants.h* file in the PTM source directory.
+For atoms that are not within the compute group-ID, all values are set to zero.
 
 Restrictions
 """"""""""""

src/USER-PTM/compute_ptm_atom.cpp

@@ -33,6 +33,7 @@ under
 #include "neigh_request.h"
 #include "neighbor.h"
 #include "update.h"
+#include "group.h"
 
 #include "ptm_functions.h"
 
@@ -61,7 +62,7 @@ static const char cite_user_ptm_package[] =
 
 ComputePTMAtom::ComputePTMAtom(LAMMPS *lmp, int narg, char **arg)
     : Compute(lmp, narg, arg), list(NULL), output(NULL) {
-  if (narg != 5)
+  if (narg < 5 || narg > 6)
     error->all(FLERR, "Illegal compute ptm/atom command");
 
   char *structures = arg[3];
@@ -122,6 +123,14 @@ ComputePTMAtom::ComputePTMAtom(LAMMPS *lmp, int narg, char **arg)
   if (rmsd_threshold == 0)
     rmsd_threshold = INFINITY;
 
+  char* group_name = (char *)"all";
+  if (narg > 5) {
+    group_name = arg[5];
+  }
+  int igroup = group->find(group_name);
+  if (igroup == -1) error->all(FLERR,"Could not find fix group ID");
+  group2bit = group->bitmask[igroup];
+
   peratom_flag = 1;
   size_peratom_cols = NUM_COLUMNS;
   create_attribute = 1;
@@ -168,6 +177,8 @@ typedef struct
   int **firstneigh;
   int *ilist;
   int nlocal;
+  int *mask;
+  int group2bit;
 
 } ptmnbrdata_t;
 
@@ -184,6 +195,8 @@ static bool sorthelper_compare(ptmnbr_t const &a, ptmnbr_t const &b) {
 static int get_neighbours(void* vdata, size_t central_index, size_t atom_index, int num, size_t* nbr_indices, int32_t* numbers, double (*nbr_pos)[3])
 {
   ptmnbrdata_t* data = (ptmnbrdata_t*)vdata;
+  int *mask = data->mask;
+  int group2bit = data->group2bit;
 
   double **x = data->x;
   double *pos = x[atom_index];
@@ -203,6 +216,9 @@ static int get_neighbours(void* vdata, size_t central_index, size_t atom_index,
 
   for (int jj = 0; jj < jnum; jj++) {
     int j = jlist[jj];
+    if (!(mask[j] & group2bit))
+      continue;
+
     j &= NEIGHMASK;
     if (j == atom_index)
       continue;
@@ -265,7 +281,11 @@ void ComputePTMAtom::compute_peratom() {
 
   double **x = atom->x;
   int *mask = atom->mask;
-  ptmnbrdata_t nbrlist = {x, numneigh, firstneigh, ilist, atom->nlocal};
+  ptmnbrdata_t nbrlist = {x, numneigh, firstneigh, ilist, atom->nlocal, mask, group2bit};
+
+  // zero output
+
+  memset(output,0,nmax*NUM_COLUMNS*sizeof(double));
 
   for (int ii = 0; ii < inum; ii++) {
 

src/USER-PTM/compute_ptm_atom.h

@@ -39,6 +39,7 @@ class ComputePTMAtom : public Compute {
   double rmsd_threshold;
   class NeighList *list;
   double **output;
+  int group2bit;
 };
 
 }