From 0659a1fc5a1c5c5efb7d8a0212ad5f5d8aa18cee Mon Sep 17 00:00:00 2001
From: sjplimp <sjplimp@f3b2605a-c512-4ea7-a41b-209d697bcdaa>
Date: Fri, 22 Jul 2016 22:10:50 +0000
Subject: [PATCH] git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@15351
 f3b2605a-c512-4ea7-a41b-209d697bcdaa

---
 src/compute_pair_local.cpp | 40 ++++++++++++++++++++++++++++++++++----
 src/compute_pair_local.h   |  3 +--
 2 files changed, 37 insertions(+), 6 deletions(-)

diff --git a/src/compute_pair_local.cpp b/src/compute_pair_local.cpp
index ad404b3f5e..13dccc5787 100644
--- a/src/compute_pair_local.cpp
+++ b/src/compute_pair_local.cpp
@@ -31,6 +31,7 @@ using namespace LAMMPS_NS;
 #define DELTA 10000
 
 enum{DIST,ENG,FORCE,FX,FY,FZ,PN};
+enum{TYPE,RADIUS};
 
 /* ---------------------------------------------------------------------- */
 
@@ -48,7 +49,8 @@ ComputePairLocal::ComputePairLocal(LAMMPS *lmp, int narg, char **arg) :
   pindex = new int[nvalues];
 
   nvalues = 0;
-  for (int iarg = 3; iarg < narg; iarg++) {
+  int iarg = 3;
+  while (iarg < narg) {
     if (strcmp(arg[iarg],"dist") == 0) pstyle[nvalues++] = DIST;
     else if (strcmp(arg[iarg],"eng") == 0) pstyle[nvalues++] = ENG;
     else if (strcmp(arg[iarg],"force") == 0) pstyle[nvalues++] = FORCE;
@@ -61,9 +63,32 @@ ComputePairLocal::ComputePairLocal(LAMMPS *lmp, int narg, char **arg) :
                              "Invalid keyword in compute pair/local command");
       pstyle[nvalues] = PN;
       pindex[nvalues++] = n-1;
-    } else error->all(FLERR,"Invalid keyword in compute pair/local command");
+
+    } else break;
+
+    iarg++;
   }
 
+  // optional args
+
+  cutstyle = TYPE;
+
+  while (iarg < narg) {
+    if (strcmp(arg[iarg],"cutoff") == 0) {
+      if (iarg+2 > narg) 
+        error->all(FLERR,"Illegal compute pair/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 pair/local command");
+      iarg += 2;
+    } else error->all(FLERR,"Illegal compute pair/local command");
+  }
+
+  // error check
+
+  if (cutstyle == RADIUS && !atom->radius_flag)
+    error->all(FLERR,"Compute pair/local requires atom attribute radius");
+
   // set singleflag if need to call pair->single()
 
   singleflag = 0;
@@ -140,11 +165,12 @@ int ComputePairLocal::compute_pairs(int flag)
   int i,j,m,n,ii,jj,inum,jnum,itype,jtype;
   tagint itag,jtag;
   double xtmp,ytmp,ztmp,delx,dely,delz;
-  double rsq,eng,fpair,factor_coul,factor_lj;
+  double rsq,radsum,eng,fpair,factor_coul,factor_lj;
   int *ilist,*jlist,*numneigh,**firstneigh;
   double *ptr;
 
   double **x = atom->x;
+  double *radius = atom->radius;
   tagint *tag = atom->tag;
   int *type = atom->type;
   int *mask = atom->mask;
@@ -216,7 +242,13 @@ int ComputePairLocal::compute_pairs(int flag)
       delz = ztmp - x[j][2];
       rsq = delx*delx + dely*dely + delz*delz;
       jtype = type[j];
-      if (rsq >= cutsq[itype][jtype]) continue;
+
+      if (cutstyle == TYPE) {
+        if (rsq >= cutsq[itype][jtype]) continue;
+      } else {
+        radsum = radius[i] + radius[j];
+        if (rsq >= radsum*radsum) continue;
+      }
 
       if (flag) {
         if (singleflag)
diff --git a/src/compute_pair_local.h b/src/compute_pair_local.h
index 663c3f931c..d6c2a12f22 100644
--- a/src/compute_pair_local.h
+++ b/src/compute_pair_local.h
@@ -34,8 +34,7 @@ class ComputePairLocal : public Compute {
   double memory_usage();
 
  private:
-  int nvalues;
-  int ncount;
+  int nvalues,ncount,cutstyle;
 
   int *pstyle;              // style of each requested output
   int *pindex;              // for pI, index of the output (0 to M-1)