git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@8562 f3b2605a-c512-4ea7-a41b-209d697bcdaa

src/compute_coord_atom.cpp

@@ -34,23 +34,48 @@ using namespace LAMMPS_NS;
 ComputeCoordAtom::ComputeCoordAtom(LAMMPS *lmp, int narg, char **arg) :
   Compute(lmp, narg, arg)
 {
-  if (narg != 4) error->all(FLERR,"Illegal compute coord/atom command");
+  if (narg < 4) error->all(FLERR,"Illegal compute coord/atom command");
 
   double cutoff = atof(arg[3]);
   cutsq = cutoff*cutoff;
 
+  ncol = narg-4 + 1;
+  int ntypes = atom->ntypes;
+  typelo = new int[ncol];
+  typehi = new int[ncol];
+
+  if (narg == 4) {
+    ncol = 1;
+    typelo[0] = 1;
+    typehi[0] = ntypes;
+  } else {
+    ncol = 0;
+    int iarg = 4;
+    while (iarg < narg) {
+      force->bounds(arg[iarg],ntypes,typelo[ncol],typehi[ncol]);
+      if (typelo[ncol] > typehi[ncol])
+        error->all(FLERR,"Illegal compute coord/atom command");
+      ncol++;
+    }
+  }
+
   peratom_flag = 1;
-  size_peratom_cols = 0;
+  if (ncol == 1) size_peratom_cols = 0;
+  else size_peratom_cols = ncol;
 
   nmax = 0;
-  coordination = NULL;
+  cvec = NULL;
+  carray = NULL;
 }
 
 /* ---------------------------------------------------------------------- */
 
 ComputeCoordAtom::~ComputeCoordAtom()
 {
-  memory->destroy(coordination);
+  delete [] typelo;
+  delete [] typehi;
+  memory->destroy(cvec);
+  memory->destroy(carray);
 }
 
 /* ---------------------------------------------------------------------- */
@@ -90,19 +115,27 @@ void ComputeCoordAtom::init_list(int id, NeighList *ptr)
 
 void ComputeCoordAtom::compute_peratom()
 {
-  int i,j,ii,jj,inum,jnum,n;
+  int i,j,m,ii,jj,inum,jnum,jtype,n;
   double xtmp,ytmp,ztmp,delx,dely,delz,rsq;
   int *ilist,*jlist,*numneigh,**firstneigh;
+  double *count;
 
   invoked_peratom = update->ntimestep;
 
   // grow coordination array if necessary
 
   if (atom->nlocal > nmax) {
-    memory->destroy(coordination);
+    if (ncol == 1) {
-    nmax = atom->nmax;
+      memory->destroy(cvec);
-    memory->create(coordination,nmax,"coord/atom:coordination");
+      nmax = atom->nmax;
-    vector_atom = coordination;
+      memory->create(cvec,nmax,"coord/atom:cvec");
+      vector_atom = cvec;
+    } else {
+      memory->destroy(carray);
+      nmax = atom->nmax;
+      memory->create(carray,nmax,ncol,"coord/atom:carray");
+      array_atom = carray;
+    }
   }
 
   // invoke full neighbor list (will copy or build if necessary)
@@ -114,35 +147,71 @@ void ComputeCoordAtom::compute_peratom()
   numneigh = list->numneigh;
   firstneigh = list->firstneigh;
 
-  // compute coordination number for each atom in group
+  // compute coordination number(s) for each atom in group
   // use full neighbor list to count atoms less than cutoff
 
   double **x = atom->x;
+  int *type = atom->type;
   int *mask = atom->mask;
 
-  for (ii = 0; ii < inum; ii++) {
+  if (ncol == 1) {
-    i = ilist[ii];
+    for (ii = 0; ii < inum; ii++) {
-    if (mask[i] & groupbit) {
+      i = ilist[ii];
-      xtmp = x[i][0];
+      if (mask[i] & groupbit) {
-      ytmp = x[i][1];
+        xtmp = x[i][0];
-      ztmp = x[i][2];
+        ytmp = x[i][1];
-      jlist = firstneigh[i];
+        ztmp = x[i][2];
-      jnum = numneigh[i];
+        jlist = firstneigh[i];
+        jnum = numneigh[i];
 
-      n = 0;
+        n = 0;
-      for (jj = 0; jj < jnum; jj++) {
+        for (jj = 0; jj < jnum; jj++) {
-        j = jlist[jj];
+          j = jlist[jj];
-        j &= NEIGHMASK;
+          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++;
+        }
+        
+        cvec[i] = n;
+      } else cvec[i] = 0.0;
+    }
 
-        delx = xtmp - x[j][0];
+  } else {
-        dely = ytmp - x[j][1];
+    for (ii = 0; ii < inum; ii++) {
-        delz = ztmp - x[j][2];
+      i = ilist[ii];
-        rsq = delx*delx + dely*dely + delz*delz;
+      count = carray[i];
-        if (rsq < cutsq) n++;
+      for (m = 0; m < ncol; m++) count[m] = 0.0;
+
+      if (mask[i] & groupbit) {
+        xtmp = x[i][0];
+        ytmp = x[i][1];
+        ztmp = x[i][2];
+        jlist = firstneigh[i];
+        jnum = numneigh[i];
+
+
+        for (jj = 0; jj < jnum; jj++) {
+          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) {
+            for (m = 0; m < ncol; m++)
+              if (jtype >= typelo[m] && jtype <= typehi[m])
+                count[m] += 1.0;
+          }
+        }
       }
-
+    }
-      coordination[i] = n;
-    } else coordination[i] = 0.0;
   }
 }
 
@@ -152,6 +221,6 @@ void ComputeCoordAtom::compute_peratom()
 
 double ComputeCoordAtom::memory_usage()
 {
-  double bytes = nmax * sizeof(double);
+  double bytes = ncol*nmax * sizeof(double);
   return bytes;
 }

src/compute_coord_atom.h

@@ -34,10 +34,13 @@ class ComputeCoordAtom : public Compute {
   double memory_usage();
 
  private:
-  int nmax;
+  int nmax,ncol;
   double cutsq;
   class NeighList *list;
-  double *coordination;
+
+  int *typelo,*typehi;
+  double *cvec;
+  double **carray;
 };
 
 }