add option to restrict coordination number by group

doc/src/compute_coord_atom.txt

@@ -15,8 +15,9 @@ compute ID group-ID coord/atom cstyle args ... :pre
 ID, group-ID are documented in "compute"_compute.html command :ulb,l
 coord/atom = style name of this compute command :l
 cstyle = {cutoff} or {orientorder} :l
-  {cutoff} args = cutoff typeN
+  {cutoff} args = cutoff \[group group2-ID\] typeN
     cutoff = distance within which to count coordination neighbors (distance units)
+    group {group2-ID} = select group-ID to restrict which atoms to cound for coordination number (optional)
     typeN = atom type for Nth coordination count (see asterisk form below)
   {orientorder} args = orientorderID threshold
     orientorderID = ID of an orientorder/atom compute
@@ -28,6 +29,7 @@ cstyle = {cutoff} or {orientorder} :l
 compute 1 all coord/atom cutoff 2.0
 compute 1 all coord/atom cutoff 6.0 1 2
 compute 1 all coord/atom cutoff 6.0 2*4 5*8 *
+compute 1 solute coord/atom cutoff 2.0 group solvent
 compute 1 all coord/atom orientorder 2 0.5 :pre
 
 [Description:]
@@ -38,9 +40,14 @@ meaning of the resulting value depend on the {cstyle} keyword used.
 
 The {cutoff} cstyle calculates one or more traditional coordination
 numbers for each atom.  A coordination number is defined as the number
-of neighbor atoms with specified atom type(s) that are within the
-specified cutoff distance from the central atom.  Atoms not in the
-specified group are included in the coordination number tally.
+of neighbor atoms with specified atom type(s), and optionally within
+the specified group, that are within the specified cutoff distance from
+the central atom. The compute group selects only the central atoms; all
+neighboring atoms, unless selected by type, type range, or group option,
+are included in the coordinations number tally.
+
+The optional {group} keyword allows to specify from which group atoms
+contribute to the coordination number. Default setting is group 'all'. 
 
 The {typeN} keywords allow specification of which atom types
 contribute to each coordination number.  One coordination number is
@@ -122,7 +129,9 @@ explained above.
 "compute cluster/atom"_compute_cluster_atom.html
 "compute orientorder/atom"_compute_orientorder_atom.html
 
-[Default:] none
+[Default:]
+
+group = all
 
 :line
 

src/compute_coord_atom.cpp

@@ -25,6 +25,7 @@
 #include "force.h"
 #include "pair.h"
 #include "comm.h"
+#include "group.h"
 #include "memory.h"
 #include "error.h"
 
@@ -37,10 +38,12 @@ using namespace LAMMPS_NS;
 ComputeCoordAtom::ComputeCoordAtom(LAMMPS *lmp, int narg, char **arg) :
   Compute(lmp, narg, arg),
   typelo(NULL), typehi(NULL), cvec(NULL), carray(NULL),
-  id_orientorder(NULL), normv(NULL)
+  group2(NULL), id_orientorder(NULL), normv(NULL)
 {
   if (narg < 5) error->all(FLERR,"Illegal compute coord/atom command");
 
+  jgroup = group->find("all");
+  jgroupbit = group->bitmask[jgroup];
   cstyle = NONE;
 
   if (strcmp(arg[3],"cutoff") == 0) {
@@ -48,18 +51,29 @@ ComputeCoordAtom::ComputeCoordAtom(LAMMPS *lmp, int narg, char **arg) :
     double cutoff = force->numeric(FLERR,arg[4]);
     cutsq = cutoff*cutoff;
 
-    ncol = narg-5 + 1;
+    int iarg = 5;
+    if ((narg > 6) && (strcmp(arg[5],"group") == 0)) {
+      int len = strlen(arg[6])+1;
+      group2 = new char[len];
+      strcpy(group2,arg[6]);
+      iarg += 2;
+      jgroup = group->find(group2);
+      if (jgroup == -1)
+        error->all(FLERR,"Compute coord/atom group2 ID does not exist");
+      jgroupbit = group->bitmask[jgroup];
+    }
+
+    ncol = narg-iarg + 1;
     int ntypes = atom->ntypes;
     typelo = new int[ncol];
     typehi = new int[ncol];
 
-    if (narg == 5) {
+    if (narg == iarg) {
       ncol = 1;
       typelo[0] = 1;
       typehi[0] = ntypes;
     } else {
       ncol = 0;
-      int iarg = 5;
       while (iarg < narg) {
         force->bounds(FLERR,arg[iarg],ntypes,typelo[ncol],typehi[ncol]);
         if (typelo[ncol] > typehi[ncol])
@@ -106,6 +120,7 @@ ComputeCoordAtom::ComputeCoordAtom(LAMMPS *lmp, int narg, char **arg) :
 
 ComputeCoordAtom::~ComputeCoordAtom()
 {
+  delete [] group2;
   delete [] typelo;
   delete [] typehi;
   memory->destroy(cvec);
@@ -229,13 +244,15 @@ void ComputeCoordAtom::compute_peratom()
             j = jlist[jj];
             j &= NEIGHMASK;
 
-            jtype = type[j];
-            delx = xtmp - x[j][0];
-            dely = ytmp - x[j][1];
-            delz = ztmp - x[j][2];
-            rsq = delx*delx + dely*dely + delz*delz;
-            if (rsq < cutsq && jtype >= typelo[0] && jtype <= typehi[0])
-              n++;
+            if (mask[j] & jgroupbit) {
+              jtype = type[j];
+              delx = xtmp - x[j][0];
+              dely = ytmp - x[j][1];
+              delz = ztmp - x[j][2];
+              rsq = delx*delx + dely*dely + delz*delz;
+              if (rsq < cutsq && jtype >= typelo[0] && jtype <= typehi[0])
+                n++;
+            }
           }
 
           cvec[i] = n;

src/compute_coord_atom.h

@@ -45,6 +45,9 @@ class ComputeCoordAtom : public Compute {
   double *cvec;
   double **carray;
 
+  char *group2;
+  int jgroup,jgroupbit;
+
   class ComputeOrientOrderAtom *c_orientorder;
   char *id_orientorder;
   double threshold;