add short neighbor list support to pair style tersoff

src/MANYBODY/pair_tersoff.cpp

@@ -52,6 +52,10 @@ PairTersoff::PairTersoff(LAMMPS *lmp) : Pair(lmp)
   params = NULL;
   elem2param = NULL;
   map = NULL;
+
+  sizeshort = 10;
+  numshort = 0;
+  neighshort = NULL;
 }
 
 /* ----------------------------------------------------------------------
@@ -71,6 +75,7 @@ PairTersoff::~PairTersoff()
   if (allocated) {
     memory->destroy(setflag);
     memory->destroy(cutsq);
+    memory->destroy(neighshort);
     delete [] map;
   }
 }
@@ -98,6 +103,7 @@ void PairTersoff::compute(int eflag, int vflag)
   int *type = atom->type;
   int nlocal = atom->nlocal;
   int newton_pair = force->newton_pair;
+  const double cutshortsq = cutmax*cutmax;
 
   inum = list->inum;
   ilist = list->ilist;
@@ -118,12 +124,26 @@ void PairTersoff::compute(int eflag, int vflag)
 
     jlist = firstneigh[i];
     jnum = numneigh[i];
+    numshort = 0;
 
     for (jj = 0; jj < jnum; jj++) {
       j = jlist[jj];
       j &= NEIGHMASK;
-      jtag = tag[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 < cutshortsq) {
+        neighshort[numshort++] = j;
+        if (numshort > sizeshort) {
+          sizeshort += sizeshort/2;
+          memory->grow(neighshort,sizeshort,"pair:neighshort");
+        }
+      }
+
+      jtag = tag[j];
       if (itag > jtag) {
         if ((itag+jtag) % 2 == 0) continue;
       } else if (itag < jtag) {
@@ -135,12 +155,6 @@ void PairTersoff::compute(int eflag, int vflag)
       }
 
       jtype = map[type[j]];
-
-      delx = xtmp - x[j][0];
-      dely = ytmp - x[j][1];
-      delz = ztmp - x[j][2];
-      rsq = delx*delx + dely*dely + delz*delz;
-
       iparam_ij = elem2param[itype][jtype][jtype];
       if (rsq > params[iparam_ij].cutsq) continue;
 
@@ -160,9 +174,8 @@ void PairTersoff::compute(int eflag, int vflag)
     // three-body interactions
     // skip immediately if I-J is not within cutoff
 
-    for (jj = 0; jj < jnum; jj++) {
-      j = jlist[jj];
-      j &= NEIGHMASK;
+    for (jj = 0; jj < numshort; jj++) {
+      j = neighshort[jj];
       jtype = map[type[j]];
       iparam_ij = elem2param[itype][jtype][jtype];
 
@@ -176,10 +189,9 @@ void PairTersoff::compute(int eflag, int vflag)
 
       zeta_ij = 0.0;
 
-      for (kk = 0; kk < jnum; kk++) {
+      for (kk = 0; kk < numshort; kk++) {
         if (jj == kk) continue;
-        k = jlist[kk];
-        k &= NEIGHMASK;
+        k = neighshort[kk];
         ktype = map[type[k]];
         iparam_ijk = elem2param[itype][jtype][ktype];
 
@@ -208,10 +220,9 @@ void PairTersoff::compute(int eflag, int vflag)
 
       // attractive term via loop over k
 
-      for (kk = 0; kk < jnum; kk++) {
+      for (kk = 0; kk < numshort; kk++) {
         if (jj == kk) continue;
-        k = jlist[kk];
-        k &= NEIGHMASK;
+        k = neighshort[kk];
         ktype = map[type[k]];
         iparam_ijk = elem2param[itype][jtype][ktype];
 
@@ -251,7 +262,7 @@ void PairTersoff::allocate()
 
   memory->create(setflag,n+1,n+1,"pair:setflag");
   memory->create(cutsq,n+1,n+1,"pair:cutsq");
-
+  memory->create(neighshort,sizeshort,"pair:neighshort");
   map = new int[n+1];
 }
 

src/MANYBODY/pair_tersoff.h

@@ -59,6 +59,9 @@ class PairTersoff : public Pair {
   int nelements;                // # of unique elements
   int nparams;                  // # of stored parameter sets
   int maxparam;                 // max # of parameter sets
+  int sizeshort;                // size of short neighbor list array
+  int numshort;                 // current count of neighbors in short list
+  int *neighshort;              // short neighbor list array
 
   virtual void allocate();
   virtual void read_file(char *);