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

src/compute_property_local.cpp

@@ -27,6 +27,7 @@
 using namespace LAMMPS_NS;
 
 enum{NONE,NEIGH,PAIR,BOND,ANGLE,DIHEDRAL,IMPROPER};
+enum{TYPE,RADIUS};
 
 #define DELTA 10000
 
@@ -39,15 +40,14 @@ ComputePropertyLocal::ComputePropertyLocal(LAMMPS *lmp, int narg, char **arg) :
 
   local_flag = 1;
   nvalues = narg - 3;
-  if (nvalues == 1) size_local_cols = 0;
-  else size_local_cols = nvalues;
-
   pack_choice = new FnPtrPack[nvalues];
 
   kindflag = NONE;
 
   int i;
-  for (int iarg = 3; iarg < narg; iarg++) {
+  nvalues = 0;
+  int iarg = 3;
+  while (iarg < narg) {
     i = iarg-3;
 
     if (strcmp(arg[iarg],"natom1") == 0) {
@@ -206,8 +206,28 @@ ComputePropertyLocal::ComputePropertyLocal(LAMMPS *lmp, int narg, char **arg) :
                    "Compute property/local cannot use these inputs together");
       kindflag = IMPROPER;
 
-    } else error->all(FLERR,
+    } else break;
-                      "Invalid keyword in compute property/local command");
+
+    nvalues++;
+    iarg++;
+  }
+
+  if (nvalues == 1) size_local_cols = 0;
+  else size_local_cols = nvalues;
+
+  // optional args
+
+  cutstyle = TYPE;
+
+  while (iarg < narg) {
+    if (strcmp(arg[iarg],"cutoff") == 0) {
+      if (iarg+2 > narg) 
+        error->all(FLERR,"Illegal compute property/local command");
+      if (strcmp(arg[iarg+1],"type") == 0) cutstyle = TYPE;
+      else if (strcmp(arg[iarg+1],"radius") == 0) cutstyle = RADIUS;
+      else error->all(FLERR,"Illegal compute property/local command");
+      iarg += 2;
+    } else error->all(FLERR,"Illegal compute property/local command");
   }
 
   // error check
@@ -229,6 +249,8 @@ ComputePropertyLocal::ComputePropertyLocal(LAMMPS *lmp, int narg, char **arg) :
   if (kindflag == IMPROPER && atom->avec->impropers_allow == 0)
     error->all(FLERR,
                "Compute property/local for property that isn't allocated");
+  if (cutstyle == RADIUS && !atom->radius_flag)
+    error->all(FLERR,"Compute property/local requires atom attribute radius");
 
   nmax = 0;
   vector = NULL;
@@ -336,10 +358,11 @@ int ComputePropertyLocal::count_pairs(int allflag, int forceflag)
 {
   int i,j,m,ii,jj,inum,jnum,itype,jtype;
   tagint itag,jtag;
-  double xtmp,ytmp,ztmp,delx,dely,delz,rsq;
+  double xtmp,ytmp,ztmp,delx,dely,delz,rsq,radsum;
   int *ilist,*jlist,*numneigh,**firstneigh;
 
   double **x = atom->x;
+  double *radius = atom->radius;
   tagint *tag = atom->tag;
   int *type = atom->type;
   int *mask = atom->mask;
@@ -404,7 +427,14 @@ int ComputePropertyLocal::count_pairs(int allflag, int forceflag)
       delz = ztmp - x[j][2];
       rsq = delx*delx + dely*dely + delz*delz;
       jtype = type[j];
-      if (forceflag && rsq >= cutsq[itype][jtype]) continue;
+      if (forceflag) {
+        if (cutstyle == TYPE) {
+          if (rsq >= cutsq[itype][jtype]) continue;
+        } else {
+          radsum = radius[i] + radius[j];
+          if (rsq >= radsum*radsum) continue;
+        }
+      }
 
       if (allflag) {
         indices[m][0] = i;

src/compute_property_local.h

@@ -34,7 +34,7 @@ class ComputePropertyLocal : public Compute {
   double memory_usage();
 
  private:
-  int nvalues,kindflag;
+  int nvalues,kindflag,cutstyle;
 
   int nmax;
   double *vector;