diff --git a/src/MANYBODY/pair_sw.cpp b/src/MANYBODY/pair_sw.cpp
index 666230362c3f25b90d7f19a5a4d2793dd6fc1bae..da2f5da37e003119d2f77521f883bbb28630bc20 100644
--- a/src/MANYBODY/pair_sw.cpp
+++ b/src/MANYBODY/pair_sw.cpp
@@ -51,6 +51,9 @@ PairSW::PairSW(LAMMPS *lmp) : Pair(lmp)
   params = NULL;
   elem2param = NULL;
   map = NULL;
+
+  maxshort = 10;
+  neighshort = NULL;
 }
 
 /* ----------------------------------------------------------------------
@@ -70,6 +73,7 @@ PairSW::~PairSW()
   if (allocated) {
     memory->destroy(setflag);
     memory->destroy(cutsq);
+    memory->destroy(neighshort);
     delete [] map;
   }
 }
@@ -102,6 +106,8 @@ void PairSW::compute(int eflag, int vflag)
   numneigh = list->numneigh;
   firstneigh = list->firstneigh;
 
+  double fxtmp,fytmp,fztmp;
+
   // loop over full neighbor list of my atoms
 
   for (ii = 0; ii < inum; ii++) {
@@ -111,17 +117,36 @@ void PairSW::compute(int eflag, int vflag)
     xtmp = x[i][0];
     ytmp = x[i][1];
     ztmp = x[i][2];
+    fxtmp = fytmp = fztmp = 0.0;
 
     // two-body interactions, skip half of them
 
     jlist = firstneigh[i];
     jnum = numneigh[i];
+    int 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;
+
+      jtype = map[type[j]];
+      ijparam = elem2param[itype][jtype][jtype];
+      if (rsq >= params[ijparam].cutsq) {
+        continue;
+      } else {
+        neighshort[numshort++] = j;
+        if (numshort >= maxshort) {
+          maxshort += maxshort/2;
+          memory->grow(neighshort,maxshort,"pair:neighshort");
+        }
+      }
+
+      jtag = tag[j];
       if (itag > jtag) {
         if ((itag+jtag) % 2 == 0) continue;
       } else if (itag < jtag) {
@@ -132,21 +157,11 @@ void PairSW::compute(int eflag, int vflag)
         if (x[j][2] == ztmp && x[j][1] == ytmp && x[j][0] < xtmp) continue;
       }
 
-      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;
-
-      ijparam = elem2param[itype][jtype][jtype];
-      if (rsq >= params[ijparam].cutsq) continue;
-
       twobody(&params[ijparam],rsq,fpair,eflag,evdwl);
 
-      f[i][0] += delx*fpair;
-      f[i][1] += dely*fpair;
-      f[i][2] += delz*fpair;
+      fxtmp += delx*fpair;
+      fytmp += dely*fpair;
+      fztmp += delz*fpair;
       f[j][0] -= delx*fpair;
       f[j][1] -= dely*fpair;
       f[j][2] -= delz*fpair;
@@ -155,22 +170,22 @@ void PairSW::compute(int eflag, int vflag)
                            evdwl,0.0,fpair,delx,dely,delz);
     }
 
-    jnumm1 = jnum - 1;
+    jnumm1 = numshort - 1;
 
     for (jj = 0; jj < jnumm1; jj++) {
-      j = jlist[jj];
-      j &= NEIGHMASK;
+      j = neighshort[jj];
       jtype = map[type[j]];
       ijparam = elem2param[itype][jtype][jtype];
       delr1[0] = x[j][0] - xtmp;
       delr1[1] = x[j][1] - ytmp;
       delr1[2] = x[j][2] - ztmp;
       rsq1 = delr1[0]*delr1[0] + delr1[1]*delr1[1] + delr1[2]*delr1[2];
-      if (rsq1 >= params[ijparam].cutsq) continue;
 
-      for (kk = jj+1; kk < jnum; kk++) {
-        k = jlist[kk];
-        k &= NEIGHMASK;
+      double fjxtmp,fjytmp,fjztmp;
+      fjxtmp = fjytmp = fjztmp = 0.0;
+
+      for (kk = jj+1; kk < numshort; kk++) {
+        k = neighshort[kk];
         ktype = map[type[k]];
         ikparam = elem2param[itype][ktype][ktype];
         ijkparam = elem2param[itype][jtype][ktype];
@@ -179,24 +194,29 @@ void PairSW::compute(int eflag, int vflag)
         delr2[1] = x[k][1] - ytmp;
         delr2[2] = x[k][2] - ztmp;
         rsq2 = delr2[0]*delr2[0] + delr2[1]*delr2[1] + delr2[2]*delr2[2];
-        if (rsq2 >= params[ikparam].cutsq) continue;
 
         threebody(&params[ijparam],&params[ikparam],&params[ijkparam],
                   rsq1,rsq2,delr1,delr2,fj,fk,eflag,evdwl);
 
-        f[i][0] -= fj[0] + fk[0];
-        f[i][1] -= fj[1] + fk[1];
-        f[i][2] -= fj[2] + fk[2];
-        f[j][0] += fj[0];
-        f[j][1] += fj[1];
-        f[j][2] += fj[2];
+        fxtmp -= fj[0] + fk[0];
+        fytmp -= fj[1] + fk[1];
+        fztmp -= fj[2] + fk[2];
+        fjxtmp += fj[0];
+        fjytmp += fj[1];
+        fjztmp += fj[2];
         f[k][0] += fk[0];
         f[k][1] += fk[1];
         f[k][2] += fk[2];
 
         if (evflag) ev_tally3(i,j,k,evdwl,0.0,fj,fk,delr1,delr2);
       }
+      f[j][0] += fjxtmp;
+      f[j][1] += fjytmp;
+      f[j][2] += fjztmp;
     }
+    f[i][0] += fxtmp;
+    f[i][1] += fytmp;
+    f[i][2] += fztmp;
   }
 
   if (vflag_fdotr) virial_fdotr_compute();
@@ -211,7 +231,7 @@ void PairSW::allocate()
 
   memory->create(setflag,n+1,n+1,"pair:setflag");
   memory->create(cutsq,n+1,n+1,"pair:cutsq");
-
+  memory->create(neighshort,maxshort,"pair:neighshort");
   map = new int[n+1];
 }
 
diff --git a/src/MANYBODY/pair_sw.h b/src/MANYBODY/pair_sw.h
index ee1fc123ff5bccb70c117aa15a1d20734a0a4772..8d921a26e81a6cdbf20ce0b16af59ac603b5b2e8 100644
--- a/src/MANYBODY/pair_sw.h
+++ b/src/MANYBODY/pair_sw.h
@@ -55,6 +55,8 @@ class PairSW : public Pair {
   int nparams;                  // # of stored parameter sets
   int maxparam;                 // max # of parameter sets
   Param *params;                // parameter set for an I-J-K interaction
+  int maxshort;                 // size of short neighbor list array
+  int *neighshort;              // short neighbor list array
 
   virtual void allocate();
   void read_file(char *);
diff --git a/src/MANYBODY/pair_tersoff.cpp b/src/MANYBODY/pair_tersoff.cpp
index 4481685ec157686382d48c3acf21aa1586d386e9..062eed3f999f01237b7a4702e59a50a87ee27f68 100644
--- a/src/MANYBODY/pair_tersoff.cpp
+++ b/src/MANYBODY/pair_tersoff.cpp
@@ -52,6 +52,9 @@ PairTersoff::PairTersoff(LAMMPS *lmp) : Pair(lmp)
   params = NULL;
   elem2param = NULL;
   map = NULL;
+
+  maxshort = 10;
+  neighshort = NULL;
 }
 
 /* ----------------------------------------------------------------------
@@ -71,6 +74,7 @@ PairTersoff::~PairTersoff()
   if (allocated) {
     memory->destroy(setflag);
     memory->destroy(cutsq);
+    memory->destroy(neighshort);
     delete [] map;
   }
 }
@@ -98,12 +102,15 @@ 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;
   numneigh = list->numneigh;
   firstneigh = list->firstneigh;
 
+  double fxtmp,fytmp,fztmp;
+
   // loop over full neighbor list of my atoms
 
   for (ii = 0; ii < inum; ii++) {
@@ -113,17 +120,32 @@ void PairTersoff::compute(int eflag, int vflag)
     xtmp = x[i][0];
     ytmp = x[i][1];
     ztmp = x[i][2];
+    fxtmp = fytmp = fztmp = 0.0;
 
     // two-body interactions, skip half of them
 
     jlist = firstneigh[i];
     jnum = numneigh[i];
+    int 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 >= maxshort) {
+          maxshort += maxshort/2;
+          memory->grow(neighshort,maxshort,"pair:neighshort");
+        }
+      }
+
+      jtag = tag[j];
       if (itag > jtag) {
         if ((itag+jtag) % 2 == 0) continue;
       } else if (itag < jtag) {
@@ -135,20 +157,14 @@ 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;
+      if (rsq >= params[iparam_ij].cutsq) continue;
 
       repulsive(&params[iparam_ij],rsq,fpair,eflag,evdwl);
 
-      f[i][0] += delx*fpair;
-      f[i][1] += dely*fpair;
-      f[i][2] += delz*fpair;
+      fxtmp += delx*fpair;
+      fytmp += dely*fpair;
+      fztmp += delz*fpair;
       f[j][0] -= delx*fpair;
       f[j][1] -= dely*fpair;
       f[j][2] -= delz*fpair;
@@ -159,10 +175,10 @@ void PairTersoff::compute(int eflag, int vflag)
 
     // three-body interactions
     // skip immediately if I-J is not within cutoff
+    double fjxtmp,fjytmp,fjztmp;
 
-    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];
 
@@ -170,16 +186,16 @@ void PairTersoff::compute(int eflag, int vflag)
       delr1[1] = x[j][1] - ytmp;
       delr1[2] = x[j][2] - ztmp;
       rsq1 = delr1[0]*delr1[0] + delr1[1]*delr1[1] + delr1[2]*delr1[2];
-      if (rsq1 > params[iparam_ij].cutsq) continue;
+      if (rsq1 >= params[iparam_ij].cutsq) continue;
 
       // accumulate bondorder zeta for each i-j interaction via loop over k
 
+      fjxtmp = fjytmp = fjztmp = 0.0;
       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];
 
@@ -187,7 +203,7 @@ void PairTersoff::compute(int eflag, int vflag)
         delr2[1] = x[k][1] - ytmp;
         delr2[2] = x[k][2] - ztmp;
         rsq2 = delr2[0]*delr2[0] + delr2[1]*delr2[1] + delr2[2]*delr2[2];
-        if (rsq2 > params[iparam_ijk].cutsq) continue;
+        if (rsq2 >= params[iparam_ijk].cutsq) continue;
 
         zeta_ij += zeta(&params[iparam_ijk],rsq1,rsq2,delr1,delr2);
       }
@@ -196,22 +212,21 @@ void PairTersoff::compute(int eflag, int vflag)
 
       force_zeta(&params[iparam_ij],rsq1,zeta_ij,fpair,prefactor,eflag,evdwl);
 
-      f[i][0] += delr1[0]*fpair;
-      f[i][1] += delr1[1]*fpair;
-      f[i][2] += delr1[2]*fpair;
-      f[j][0] -= delr1[0]*fpair;
-      f[j][1] -= delr1[1]*fpair;
-      f[j][2] -= delr1[2]*fpair;
+      fxtmp += delr1[0]*fpair;
+      fytmp += delr1[1]*fpair;
+      fztmp += delr1[2]*fpair;
+      fjxtmp -= delr1[0]*fpair;
+      fjytmp -= delr1[1]*fpair;
+      fjztmp -= delr1[2]*fpair;
 
       if (evflag) ev_tally(i,j,nlocal,newton_pair,
                            evdwl,0.0,-fpair,-delr1[0],-delr1[1],-delr1[2]);
 
       // 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];
 
@@ -219,24 +234,30 @@ void PairTersoff::compute(int eflag, int vflag)
         delr2[1] = x[k][1] - ytmp;
         delr2[2] = x[k][2] - ztmp;
         rsq2 = delr2[0]*delr2[0] + delr2[1]*delr2[1] + delr2[2]*delr2[2];
-        if (rsq2 > params[iparam_ijk].cutsq) continue;
+        if (rsq2 >= params[iparam_ijk].cutsq) continue;
 
         attractive(&params[iparam_ijk],prefactor,
                    rsq1,rsq2,delr1,delr2,fi,fj,fk);
 
-        f[i][0] += fi[0];
-        f[i][1] += fi[1];
-        f[i][2] += fi[2];
-        f[j][0] += fj[0];
-        f[j][1] += fj[1];
-        f[j][2] += fj[2];
+        fxtmp += fi[0];
+        fytmp += fi[1];
+        fztmp += fi[2];
+        fjxtmp += fj[0];
+        fjytmp += fj[1];
+        fjztmp += fj[2];
         f[k][0] += fk[0];
         f[k][1] += fk[1];
         f[k][2] += fk[2];
 
         if (vflag_atom) v_tally3(i,j,k,fj,fk,delr1,delr2);
       }
+      f[j][0] += fjxtmp;
+      f[j][1] += fjytmp;
+      f[j][2] += fjztmp;
     }
+    f[i][0] += fxtmp;
+    f[i][1] += fytmp;
+    f[i][2] += fztmp;
   }
 
   if (vflag_fdotr) virial_fdotr_compute();
@@ -251,7 +272,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,maxshort,"pair:neighshort");
   map = new int[n+1];
 }
 
diff --git a/src/MANYBODY/pair_tersoff.h b/src/MANYBODY/pair_tersoff.h
index dd78bb2eb91fa90a24237e438fe8fff342e40524..f44d9f3c5dfa40c854714f17f5005b2f1a45cc78 100644
--- a/src/MANYBODY/pair_tersoff.h
+++ b/src/MANYBODY/pair_tersoff.h
@@ -59,6 +59,8 @@ class PairTersoff : public Pair {
   int nelements;                // # of unique elements
   int nparams;                  // # of stored parameter sets
   int maxparam;                 // max # of parameter sets
+  int maxshort;                 // size of short neighbor list array
+  int *neighshort;              // short neighbor list array
 
   virtual void allocate();
   virtual void read_file(char *);
diff --git a/src/MANYBODY/pair_vashishta.cpp b/src/MANYBODY/pair_vashishta.cpp
index 19f6017907b4767bf6f92fff95672ec861e66894..6c5cc253351d9ea656aab7a17ab15398c788e1f4 100644
--- a/src/MANYBODY/pair_vashishta.cpp
+++ b/src/MANYBODY/pair_vashishta.cpp
@@ -52,6 +52,10 @@ PairVashishta::PairVashishta(LAMMPS *lmp) : Pair(lmp)
   params = NULL;
   elem2param = NULL;
   map = NULL;
+
+  r0max = 0.0;
+  maxshort = 10;
+  neighshort = NULL;
 }
 
 /* ----------------------------------------------------------------------
@@ -69,6 +73,7 @@ PairVashishta::~PairVashishta()
   if (allocated) {
     memory->destroy(setflag);
     memory->destroy(cutsq);
+    memory->destroy(neighshort);
     delete [] map;
   }
 }
@@ -95,12 +100,15 @@ void PairVashishta::compute(int eflag, int vflag)
   int *type = atom->type;
   int nlocal = atom->nlocal;
   int newton_pair = force->newton_pair;
+  const double cutshortsq = r0max*r0max;
 
   inum = list->inum;
   ilist = list->ilist;
   numneigh = list->numneigh;
   firstneigh = list->firstneigh;
 
+  double fxtmp,fytmp,fztmp;
+
   // loop over full neighbor list of my atoms
 
   for (ii = 0; ii < inum; ii++) {
@@ -110,17 +118,32 @@ void PairVashishta::compute(int eflag, int vflag)
     xtmp = x[i][0];
     ytmp = x[i][1];
     ztmp = x[i][2];
+    fxtmp = fytmp = fztmp = 0.0;
 
     // two-body interactions, skip half of them
 
     jlist = firstneigh[i];
     jnum = numneigh[i];
+    int 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 >= maxshort) {
+          maxshort += maxshort/2;
+          memory->grow(neighshort,maxshort,"pair:neighshort");
+        }
+      }
+
+      jtag = tag[j];
       if (itag > jtag) {
         if ((itag+jtag) % 2 == 0) continue;
       } else if (itag < jtag) {
@@ -132,20 +155,14 @@ void PairVashishta::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;
-
       ijparam = elem2param[itype][jtype][jtype];
-      if (rsq > params[ijparam].cutsq) continue;
+      if (rsq >= params[ijparam].cutsq) continue;
 
       twobody(&params[ijparam],rsq,fpair,eflag,evdwl);
 
-      f[i][0] += delx*fpair;
-      f[i][1] += dely*fpair;
-      f[i][2] += delz*fpair;
+      fxtmp += delx*fpair;
+      fytmp += dely*fpair;
+      fztmp += delz*fpair;
       f[j][0] -= delx*fpair;
       f[j][1] -= dely*fpair;
       f[j][2] -= delz*fpair;
@@ -154,11 +171,10 @@ void PairVashishta::compute(int eflag, int vflag)
       			   evdwl,0.0,fpair,delx,dely,delz);
     }
 
-    jnumm1 = jnum - 1;
+    jnumm1 = numshort - 1;
 
     for (jj = 0; jj < jnumm1; jj++) {
-      j = jlist[jj];
-      j &= NEIGHMASK;
+      j = neighshort[jj];
       jtype = map[type[j]];
       ijparam = elem2param[itype][jtype][jtype];
       delr1[0] = x[j][0] - xtmp;
@@ -167,9 +183,11 @@ void PairVashishta::compute(int eflag, int vflag)
       rsq1 = delr1[0]*delr1[0] + delr1[1]*delr1[1] + delr1[2]*delr1[2];
       if (rsq1 >= params[ijparam].cutsq2) continue;
 
-      for (kk = jj+1; kk < jnum; kk++) {
-        k = jlist[kk];
-        k &= NEIGHMASK;
+      double fjxtmp,fjytmp,fjztmp;
+      fjxtmp = fjytmp = fjztmp = 0.0;
+
+      for (kk = jj+1; kk < numshort; kk++) {
+        k = neighshort[kk];
         ktype = map[type[k]];
         ikparam = elem2param[itype][ktype][ktype];
         ijkparam = elem2param[itype][jtype][ktype];
@@ -183,19 +201,25 @@ void PairVashishta::compute(int eflag, int vflag)
         threebody(&params[ijparam],&params[ikparam],&params[ijkparam],
                   rsq1,rsq2,delr1,delr2,fj,fk,eflag,evdwl);
 
-        f[i][0] -= fj[0] + fk[0];
-        f[i][1] -= fj[1] + fk[1];
-        f[i][2] -= fj[2] + fk[2];
-        f[j][0] += fj[0];
-        f[j][1] += fj[1];
-        f[j][2] += fj[2];
+        fxtmp -= fj[0] + fk[0];
+        fytmp -= fj[1] + fk[1];
+        fztmp -= fj[2] + fk[2];
+        fjxtmp += fj[0];
+        fjytmp += fj[1];
+        fjztmp += fj[2];
         f[k][0] += fk[0];
         f[k][1] += fk[1];
         f[k][2] += fk[2];
 
-	if (evflag) ev_tally3(i,j,k,evdwl,0.0,fj,fk,delr1,delr2);
+        if (evflag) ev_tally3(i,j,k,evdwl,0.0,fj,fk,delr1,delr2);
       }
+      f[j][0] += fjxtmp;
+      f[j][1] += fjytmp;
+      f[j][2] += fjztmp;
     }
+    f[i][0] += fxtmp;
+    f[i][1] += fytmp;
+    f[i][2] += fztmp;
   }
 
   if (vflag_fdotr) virial_fdotr_compute();
@@ -210,6 +234,7 @@ void PairVashishta::allocate()
 
   memory->create(setflag,n+1,n+1,"pair:setflag");
   memory->create(cutsq,n+1,n+1,"pair:cutsq");
+  memory->create(neighshort,maxshort,"pair:neighshort");
 
   map = new int[n+1];
 }
@@ -529,13 +554,15 @@ void PairVashishta::setup_params()
     params[m].c0 = params[m].cut*params[m].dvrc - params[m].vrc;
   }
 
-  // set cutmax to max of all params
+  // set cutmax to max of all cutoff params. r0max only for r0
 
   cutmax = 0.0;
+  r0max = 0.0;
   for (m = 0; m < nparams; m++) {
     if (params[m].cut > cutmax) cutmax = params[m].cut;
-    if (params[m].r0 > cutmax) cutmax = params[m].r0;
+    if (params[m].r0 > r0max) r0max = params[m].r0;
   }
+  if (r0max > cutmax) cutmax = r0max;
 }
 
 /* ---------------------------------------------------------------------- */
diff --git a/src/MANYBODY/pair_vashishta.h b/src/MANYBODY/pair_vashishta.h
index 87077011e630e22f5c83c12004869276f51a8e71..b6d9c574a51cc9e399827e1aa9ee9c3374a0be09 100644
--- a/src/MANYBODY/pair_vashishta.h
+++ b/src/MANYBODY/pair_vashishta.h
@@ -1,4 +1,4 @@
-/* ----------------------------------------------------------------------
+/* -*- c++ -*- ----------------------------------------------------------
    LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
    http://lammps.sandia.gov, Sandia National Laboratories
    Steve Plimpton, sjplimp@sandia.gov
@@ -54,6 +54,9 @@ class PairVashishta : public Pair {
   int nparams;                  // # of stored parameter sets
   int maxparam;                 // max # of parameter sets
   Param *params;                // parameter set for an I-J-K interaction
+  double r0max;                 // largest value of r0
+  int maxshort;                 // size of short neighbor list array
+  int *neighshort;              // short neighbor list array
 
   void allocate();
   void read_file(char *);
diff --git a/src/MANYBODY/pair_vashishta_table.cpp b/src/MANYBODY/pair_vashishta_table.cpp
index 9ba2f258a59691b525898911671ffd06a8a1b1d2..c58f1286d8c741c9f77c2cf22ea3168b0e8fdb9d 100644
--- a/src/MANYBODY/pair_vashishta_table.cpp
+++ b/src/MANYBODY/pair_vashishta_table.cpp
@@ -37,9 +37,6 @@ using namespace LAMMPS_NS;
 
 PairVashishtaTable::PairVashishtaTable(LAMMPS *lmp) : PairVashishta(lmp)
 {
-  neigh3BodyMax = 0;
-  neigh3BodyCount = NULL; 
-  neigh3Body = NULL;
   forceTable = NULL;
   potentialTable = NULL;
 }
@@ -52,8 +49,6 @@ PairVashishtaTable::~PairVashishtaTable()
 {
   memory->destroy(forceTable);
   memory->destroy(potentialTable);
-  memory->destroy(neigh3BodyCount);
-  memory->destroy(neigh3Body);
 }
 
 /* ---------------------------------------------------------------------- */
@@ -78,25 +73,14 @@ void PairVashishtaTable::compute(int eflag, int vflag)
   int *type = atom->type;
   int nlocal = atom->nlocal;
   int newton_pair = force->newton_pair;
+  const double cutshortsq = r0max*r0max;
 
   inum = list->inum;
   ilist = list->ilist;
   numneigh = list->numneigh;
   firstneigh = list->firstneigh;
 
-  // reallocate 3-body neighbor list if necessary
-  // NOTE: using 1000 is inefficient
-  //       could make this a LAMMPS paged neighbor list
-
-  if (nlocal > neigh3BodyMax) {
-    neigh3BodyMax = atom->nmax;
-    memory->destroy(neigh3BodyCount);
-    memory->destroy(neigh3Body);
-    memory->create(neigh3BodyCount,neigh3BodyMax,
-                   "pair:vashishta:neigh3BodyCount");
-    memory->create(neigh3Body,neigh3BodyMax,1000,
-                   "pair:vashishta:neigh3Body");
-  }
+  double fxtmp,fytmp,fztmp;
 
   // loop over full neighbor list of my atoms
 
@@ -107,36 +91,32 @@ void PairVashishtaTable::compute(int eflag, int vflag)
     xtmp = x[i][0];
     ytmp = x[i][1];
     ztmp = x[i][2];
-
-    // reset the 3-body neighbor list
-
-    neigh3BodyCount[i] = 0;
+    fxtmp = fytmp = fztmp = 0.0;
 
     // two-body interactions, skip half of them
 
     jlist = firstneigh[i];
     jnum = numneigh[i];
+    int numshort = 0;
 
     for (jj = 0; jj < jnum; jj++) {
       j = jlist[jj];
       j &= NEIGHMASK;
-      jtag = tag[j];
-
-      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;
-      ijparam = elem2param[itype][jtype][jtype];
 
-      if (rsq <= params[ijparam].cutsq2) {
-        neigh3Body[i][neigh3BodyCount[i]] = j;
-        neigh3BodyCount[i]++;
+      if (rsq < cutshortsq) {
+        neighshort[numshort++] = j;
+        if (numshort >= maxshort) {
+          maxshort += maxshort/2;
+          memory->grow(neighshort,maxshort,"pair:neighshort");
+        }
       }
 
-      if (rsq > params[ijparam].cutsq) continue;
-
+      jtag = tag[j];
       if (itag > jtag) {
         if ((itag+jtag) % 2 == 0) continue;
       } else if (itag < jtag) {
@@ -147,11 +127,15 @@ void PairVashishtaTable::compute(int eflag, int vflag)
         if (x[j][2] == ztmp && x[j][1] == ytmp && x[j][0] < xtmp) continue;
       }
 
+      jtype = map[type[j]];
+      ijparam = elem2param[itype][jtype][jtype];
+      if (rsq >= params[ijparam].cutsq) continue;
+
       twobody_table(params[ijparam],rsq,fpair,eflag,evdwl);
 
-      f[i][0] += delx*fpair;
-      f[i][1] += dely*fpair;
-      f[i][2] += delz*fpair;
+      fxtmp += delx*fpair;
+      fytmp += dely*fpair;
+      fztmp += delz*fpair;
       f[j][0] -= delx*fpair;
       f[j][1] -= dely*fpair;
       f[j][2] -= delz*fpair;
@@ -160,13 +144,10 @@ void PairVashishtaTable::compute(int eflag, int vflag)
       			   evdwl,0.0,fpair,delx,dely,delz);
     }
 
-    jlist = neigh3Body[i];
-    jnum = neigh3BodyCount[i];
-    jnumm1 = jnum - 1;
+    jnumm1 = numshort - 1;
 
     for (jj = 0; jj < jnumm1; jj++) {
-      j = jlist[jj];
-      j &= NEIGHMASK;
+      j = neighshort[jj];
       jtype = map[type[j]];
       ijparam = elem2param[itype][jtype][jtype];
       delr1[0] = x[j][0] - xtmp;
@@ -175,9 +156,11 @@ void PairVashishtaTable::compute(int eflag, int vflag)
       rsq1 = delr1[0]*delr1[0] + delr1[1]*delr1[1] + delr1[2]*delr1[2];
       if (rsq1 >= params[ijparam].cutsq2) continue;
 
-      for (kk = jj+1; kk < jnum; kk++) {
-        k = jlist[kk];
-        k &= NEIGHMASK;
+      double fjxtmp,fjytmp,fjztmp;
+      fjxtmp = fjytmp = fjztmp = 0.0;
+
+      for (kk = jj+1; kk < numshort; kk++) {
+        k = neighshort[kk];
         ktype = map[type[k]];
         ikparam = elem2param[itype][ktype][ktype];
         ijkparam = elem2param[itype][jtype][ktype];
@@ -191,19 +174,25 @@ void PairVashishtaTable::compute(int eflag, int vflag)
         threebody(&params[ijparam],&params[ikparam],&params[ijkparam],
                   rsq1,rsq2,delr1,delr2,fj,fk,eflag,evdwl);
 
-        f[i][0] -= fj[0] + fk[0];
-        f[i][1] -= fj[1] + fk[1];
-        f[i][2] -= fj[2] + fk[2];
-        f[j][0] += fj[0];
-        f[j][1] += fj[1];
-        f[j][2] += fj[2];
+        fxtmp -= fj[0] + fk[0];
+        fytmp -= fj[1] + fk[1];
+        fztmp -= fj[2] + fk[2];
+        fjxtmp += fj[0];
+        fjytmp += fj[1];
+        fjztmp += fj[2];
         f[k][0] += fk[0];
         f[k][1] += fk[1];
         f[k][2] += fk[2];
 
         if (evflag) ev_tally3(i,j,k,evdwl,0.0,fj,fk,delr1,delr2);
       }
+      f[j][0] += fjxtmp;
+      f[j][1] += fjytmp;
+      f[j][2] += fjztmp;
     }
+    f[i][0] += fxtmp;
+    f[i][1] += fytmp;
+    f[i][2] += fztmp;
   }
 
   if (vflag_fdotr) virial_fdotr_compute();
diff --git a/src/MANYBODY/pair_vashishta_table.h b/src/MANYBODY/pair_vashishta_table.h
index 8f96829142a0bba12cbc032a020cbbd493ee0e10..a45cac5ae1c719407c3faeaf0885f5f2aa7ef1b4 100644
--- a/src/MANYBODY/pair_vashishta_table.h
+++ b/src/MANYBODY/pair_vashishta_table.h
@@ -39,11 +39,6 @@ class PairVashishtaTable : public PairVashishta {
   double ***forceTable;         // table of forces per element pair
   double ***potentialTable;     // table of potential energies
 
-  int neigh3BodyMax;            // max size of short neighborlist
-  int *neigh3BodyCount;         // # of neighbors in short range 
-                                // 3 particle forces neighbor list
-  int **neigh3Body;             // neighlist for short range 3 particle forces
-
   void twobody_table(const Param &, double, double &, int, double &);
   void setup_params();
   void create_tables();
diff --git a/src/USER-OMP/pair_sw_omp.cpp b/src/USER-OMP/pair_sw_omp.cpp
index ba7eeacae1df8be1c050834610e83e4c05a10b5a..9a63b4f4394846b09988a22ab2f55e00cb05d049 100644
--- a/src/USER-OMP/pair_sw_omp.cpp
+++ b/src/USER-OMP/pair_sw_omp.cpp
@@ -17,6 +17,7 @@
 #include "atom.h"
 #include "comm.h"
 #include "force.h"
+#include "memory.h"
 #include "neighbor.h"
 #include "neigh_list.h"
 
@@ -71,13 +72,13 @@ void PairSWOMP::compute(int eflag, int vflag)
 template <int EVFLAG, int EFLAG>
 void PairSWOMP::eval(int iifrom, int iito, ThrData * const thr)
 {
-  int i,j,k,ii,jj,kk,jnum,jnumm1;
+  int i,j,k,ii,jj,kk,jnum,jnumm1,maxshort_thr;
   tagint itag,jtag;
   int itype,jtype,ktype,ijparam,ikparam,ijkparam;
   double xtmp,ytmp,ztmp,delx,dely,delz,evdwl,fpair;
   double rsq,rsq1,rsq2;
   double delr1[3],delr2[3],fj[3],fk[3];
-  int *ilist,*jlist,*numneigh,**firstneigh;
+  int *ilist,*jlist,*numneigh,**firstneigh,*neighshort_thr;
 
   evdwl = 0.0;
 
@@ -90,6 +91,8 @@ void PairSWOMP::eval(int iifrom, int iito, ThrData * const thr)
   ilist = list->ilist;
   numneigh = list->numneigh;
   firstneigh = list->firstneigh;
+  maxshort_thr = maxshort;
+  memory->create(neighshort_thr,maxshort_thr,"pair_thr:neighshort_thr");
 
   double fxtmp,fytmp,fztmp;
 
@@ -109,12 +112,30 @@ void PairSWOMP::eval(int iifrom, int iito, ThrData * const thr)
 
     jlist = firstneigh[i];
     jnum = numneigh[i];
+    int numshort = 0;
 
     for (jj = 0; jj < jnum; jj++) {
       j = jlist[jj];
       j &= NEIGHMASK;
-      jtag = tag[j];
 
+      delx = xtmp - x[j].x;
+      dely = ytmp - x[j].y;
+      delz = ztmp - x[j].z;
+      rsq = delx*delx + dely*dely + delz*delz;
+
+      jtype = map[type[j]];
+      ijparam = elem2param[itype][jtype][jtype];
+      if (rsq >= params[ijparam].cutsq) {
+        continue;
+      } else {
+        neighshort_thr[numshort++] = j;
+        if (numshort >= maxshort_thr) {
+          maxshort_thr += maxshort_thr/2;
+          memory->grow(neighshort_thr,maxshort_thr,"pair:neighshort_thr");
+        }
+      }
+
+      jtag = tag[j];
       if (itag > jtag) {
         if ((itag+jtag) % 2 == 0) continue;
       } else if (itag < jtag) {
@@ -125,16 +146,6 @@ void PairSWOMP::eval(int iifrom, int iito, ThrData * const thr)
         if (x[j].z == ztmp && x[j].y == ytmp && x[j].x < xtmp) continue;
       }
 
-      jtype = map[type[j]];
-
-      delx = xtmp - x[j].x;
-      dely = ytmp - x[j].y;
-      delz = ztmp - x[j].z;
-      rsq = delx*delx + dely*dely + delz*delz;
-
-      ijparam = elem2param[itype][jtype][jtype];
-      if (rsq >= params[ijparam].cutsq) continue;
-
       twobody(&params[ijparam],rsq,fpair,EFLAG,evdwl);
 
       fxtmp += delx*fpair;
@@ -148,25 +159,22 @@ void PairSWOMP::eval(int iifrom, int iito, ThrData * const thr)
                                evdwl,0.0,fpair,delx,dely,delz,thr);
     }
 
-    jnumm1 = jnum - 1;
+    jnumm1 = numshort - 1;
 
     for (jj = 0; jj < jnumm1; jj++) {
-      j = jlist[jj];
-      j &= NEIGHMASK;
+      j = neighshort_thr[jj];
       jtype = map[type[j]];
       ijparam = elem2param[itype][jtype][jtype];
       delr1[0] = x[j].x - xtmp;
       delr1[1] = x[j].y - ytmp;
       delr1[2] = x[j].z - ztmp;
       rsq1 = delr1[0]*delr1[0] + delr1[1]*delr1[1] + delr1[2]*delr1[2];
-      if (rsq1 >= params[ijparam].cutsq) continue;
 
       double fjxtmp,fjytmp,fjztmp;
       fjxtmp = fjytmp = fjztmp = 0.0;
 
-      for (kk = jj+1; kk < jnum; kk++) {
-        k = jlist[kk];
-        k &= NEIGHMASK;
+      for (kk = jj+1; kk < numshort; kk++) {
+        k = neighshort_thr[kk];
         ktype = map[type[k]];
         ikparam = elem2param[itype][ktype][ktype];
         ijkparam = elem2param[itype][jtype][ktype];
@@ -175,7 +183,6 @@ void PairSWOMP::eval(int iifrom, int iito, ThrData * const thr)
         delr2[1] = x[k].y - ytmp;
         delr2[2] = x[k].z - ztmp;
         rsq2 = delr2[0]*delr2[0] + delr2[1]*delr2[1] + delr2[2]*delr2[2];
-        if (rsq2 >= params[ikparam].cutsq) continue;
 
         threebody(&params[ijparam],&params[ikparam],&params[ijkparam],
                   rsq1,rsq2,delr1,delr2,fj,fk,EFLAG,evdwl);
@@ -200,6 +207,7 @@ void PairSWOMP::eval(int iifrom, int iito, ThrData * const thr)
     f[i].y += fytmp;
     f[i].z += fztmp;
   }
+  memory->destroy(neighshort_thr);
 }
 
 /* ---------------------------------------------------------------------- */
diff --git a/src/USER-OMP/pair_tersoff_omp.cpp b/src/USER-OMP/pair_tersoff_omp.cpp
index 71d00108dfdec0d0d0750f27538c68f90bf728fa..db249b4ba907cd3afb78333141b2137a084504ff 100644
--- a/src/USER-OMP/pair_tersoff_omp.cpp
+++ b/src/USER-OMP/pair_tersoff_omp.cpp
@@ -17,6 +17,7 @@
 #include "atom.h"
 #include "comm.h"
 #include "force.h"
+#include "memory.h"
 #include "neighbor.h"
 #include "neigh_list.h"
 
@@ -73,14 +74,14 @@ void PairTersoffOMP::compute(int eflag, int vflag)
 template <int EVFLAG, int EFLAG, int VFLAG_ATOM>
 void PairTersoffOMP::eval(int iifrom, int iito, ThrData * const thr)
 {
-  int i,j,k,ii,jj,kk,jnum;
+  int i,j,k,ii,jj,kk,jnum,maxshort_thr;
   tagint itag,jtag;
   int itype,jtype,ktype,iparam_ij,iparam_ijk;
   double xtmp,ytmp,ztmp,delx,dely,delz,evdwl,fpair;
   double rsq,rsq1,rsq2;
   double delr1[3],delr2[3],fi[3],fj[3],fk[3];
   double zeta_ij,prefactor;
-  int *ilist,*jlist,*numneigh,**firstneigh;
+  int *ilist,*jlist,*numneigh,**firstneigh,*neighshort_thr;
 
   evdwl = 0.0;
 
@@ -89,10 +90,13 @@ void PairTersoffOMP::eval(int iifrom, int iito, ThrData * const thr)
   const tagint * _noalias const tag = atom->tag;
   const int * _noalias const type = atom->type;
   const int nlocal = atom->nlocal;
+  const double cutshortsq = cutmax*cutmax;
 
   ilist = list->ilist;
   numneigh = list->numneigh;
   firstneigh = list->firstneigh;
+  maxshort_thr = maxshort;
+  memory->create(neighshort_thr,maxshort_thr,"pair_thr:neighshort_thr");
 
   double fxtmp,fytmp,fztmp;
 
@@ -112,12 +116,26 @@ void PairTersoffOMP::eval(int iifrom, int iito, ThrData * const thr)
 
     jlist = firstneigh[i];
     jnum = numneigh[i];
+    int numshort = 0;
 
     for (jj = 0; jj < jnum; jj++) {
       j = jlist[jj];
       j &= NEIGHMASK;
-      jtag = tag[j];
 
+      delx = xtmp - x[j].x;
+      dely = ytmp - x[j].y;
+      delz = ztmp - x[j].z;
+      rsq = delx*delx + dely*dely + delz*delz;
+
+      if (rsq < cutshortsq) {
+        neighshort_thr[numshort++] = j;
+        if (numshort >= maxshort_thr) {
+          maxshort_thr += maxshort_thr/2;
+          memory->grow(neighshort_thr,maxshort_thr,"pair_thr:neighshort_thr");
+        }
+      }
+
+      jtag = tag[j];
       if (itag > jtag) {
         if ((itag+jtag) % 2 == 0) continue;
       } else if (itag < jtag) {
@@ -129,14 +147,8 @@ void PairTersoffOMP::eval(int iifrom, int iito, ThrData * const thr)
       }
 
       jtype = map[type[j]];
-
-      delx = xtmp - x[j].x;
-      dely = ytmp - x[j].y;
-      delz = ztmp - x[j].z;
-      rsq = delx*delx + dely*dely + delz*delz;
-
       iparam_ij = elem2param[itype][jtype][jtype];
-      if (rsq > params[iparam_ij].cutsq) continue;
+      if (rsq >= params[iparam_ij].cutsq) continue;
 
       repulsive(&params[iparam_ij],rsq,fpair,EFLAG,evdwl);
 
@@ -155,9 +167,8 @@ void PairTersoffOMP::eval(int iifrom, int iito, ThrData * const thr)
     // skip immediately if I-J is not within cutoff
     double fjxtmp,fjytmp,fjztmp;
 
-    for (jj = 0; jj < jnum; jj++) {
-      j = jlist[jj];
-      j &= NEIGHMASK;
+    for (jj = 0; jj < numshort; jj++) {
+      j = neighshort_thr[jj];
       jtype = map[type[j]];
       iparam_ij = elem2param[itype][jtype][jtype];
 
@@ -165,17 +176,16 @@ void PairTersoffOMP::eval(int iifrom, int iito, ThrData * const thr)
       delr1[1] = x[j].y - ytmp;
       delr1[2] = x[j].z - ztmp;
       rsq1 = delr1[0]*delr1[0] + delr1[1]*delr1[1] + delr1[2]*delr1[2];
-      if (rsq1 > params[iparam_ij].cutsq) continue;
+      if (rsq1 >= params[iparam_ij].cutsq) continue;
 
       // accumulate bondorder zeta for each i-j interaction via loop over k
 
       fjxtmp = fjytmp = fjztmp = 0.0;
       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_thr[kk];
         ktype = map[type[k]];
         iparam_ijk = elem2param[itype][jtype][ktype];
 
@@ -183,7 +193,7 @@ void PairTersoffOMP::eval(int iifrom, int iito, ThrData * const thr)
         delr2[1] = x[k].y - ytmp;
         delr2[2] = x[k].z - ztmp;
         rsq2 = delr2[0]*delr2[0] + delr2[1]*delr2[1] + delr2[2]*delr2[2];
-        if (rsq2 > params[iparam_ijk].cutsq) continue;
+        if (rsq2 >= params[iparam_ijk].cutsq) continue;
 
         zeta_ij += zeta(&params[iparam_ijk],rsq1,rsq2,delr1,delr2);
       }
@@ -204,10 +214,9 @@ void PairTersoffOMP::eval(int iifrom, int iito, ThrData * const thr)
 
       // 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_thr[kk];
         ktype = map[type[k]];
         iparam_ijk = elem2param[itype][jtype][ktype];
 
@@ -215,7 +224,7 @@ void PairTersoffOMP::eval(int iifrom, int iito, ThrData * const thr)
         delr2[1] = x[k].y - ytmp;
         delr2[2] = x[k].z - ztmp;
         rsq2 = delr2[0]*delr2[0] + delr2[1]*delr2[1] + delr2[2]*delr2[2];
-        if (rsq2 > params[iparam_ijk].cutsq) continue;
+        if (rsq2 >= params[iparam_ijk].cutsq) continue;
 
         attractive(&params[iparam_ijk],prefactor,
                    rsq1,rsq2,delr1,delr2,fi,fj,fk);
@@ -240,6 +249,7 @@ void PairTersoffOMP::eval(int iifrom, int iito, ThrData * const thr)
     f[i].y += fytmp;
     f[i].z += fztmp;
   }
+  memory->destroy(neighshort_thr);
 }
 
 /* ---------------------------------------------------------------------- */
diff --git a/src/USER-OMP/pair_vashishta_omp.cpp b/src/USER-OMP/pair_vashishta_omp.cpp
index 679eceab503669113a532401e355702968ad7bb3..c432738f6239e5079e2298b3055efdf5be0be217 100644
--- a/src/USER-OMP/pair_vashishta_omp.cpp
+++ b/src/USER-OMP/pair_vashishta_omp.cpp
@@ -17,6 +17,7 @@
 #include "atom.h"
 #include "comm.h"
 #include "force.h"
+#include "memory.h"
 #include "neighbor.h"
 #include "neigh_list.h"
 
@@ -71,13 +72,13 @@ void PairVashishtaOMP::compute(int eflag, int vflag)
 template <int EVFLAG, int EFLAG>
 void PairVashishtaOMP::eval(int iifrom, int iito, ThrData * const thr)
 {
-  int i,j,k,ii,jj,kk,jnum,jnumm1;
+  int i,j,k,ii,jj,kk,jnum,jnumm1,maxshort_thr;
   tagint itag,jtag;
   int itype,jtype,ktype,ijparam,ikparam,ijkparam;
   double xtmp,ytmp,ztmp,delx,dely,delz,evdwl,fpair;
   double rsq,rsq1,rsq2;
   double delr1[3],delr2[3],fj[3],fk[3];
-  int *ilist,*jlist,*numneigh,**firstneigh;
+  int *ilist,*jlist,*numneigh,**firstneigh,*neighshort_thr;
 
   evdwl = 0.0;
 
@@ -86,10 +87,13 @@ void PairVashishtaOMP::eval(int iifrom, int iito, ThrData * const thr)
   const tagint * _noalias const tag = atom->tag;
   const int * _noalias const type = atom->type;
   const int nlocal = atom->nlocal;
+  const double cutshortsq = r0max*r0max;
 
   ilist = list->ilist;
   numneigh = list->numneigh;
   firstneigh = list->firstneigh;
+  maxshort_thr = maxshort;
+  memory->create(neighshort_thr,maxshort_thr,"pair_thr:neighshort_thr");
 
   double fxtmp,fytmp,fztmp;
 
@@ -109,12 +113,26 @@ void PairVashishtaOMP::eval(int iifrom, int iito, ThrData * const thr)
 
     jlist = firstneigh[i];
     jnum = numneigh[i];
+    int numshort = 0;
 
     for (jj = 0; jj < jnum; jj++) {
       j = jlist[jj];
       j &= NEIGHMASK;
-      jtag = tag[j];
 
+      delx = xtmp - x[j].x;
+      dely = ytmp - x[j].y;
+      delz = ztmp - x[j].z;
+      rsq = delx*delx + dely*dely + delz*delz;
+
+      if (rsq < cutshortsq) {
+        neighshort_thr[numshort++] = j;
+        if (numshort >= maxshort_thr) {
+          maxshort_thr += maxshort_thr/2;
+          memory->grow(neighshort_thr,maxshort_thr,"pair_thr:neighshort_thr");
+        }
+      }
+
+      jtag = tag[j];
       if (itag > jtag) {
         if ((itag+jtag) % 2 == 0) continue;
       } else if (itag < jtag) {
@@ -126,12 +144,6 @@ void PairVashishtaOMP::eval(int iifrom, int iito, ThrData * const thr)
       }
 
       jtype = map[type[j]];
-
-      delx = xtmp - x[j].x;
-      dely = ytmp - x[j].y;
-      delz = ztmp - x[j].z;
-      rsq = delx*delx + dely*dely + delz*delz;
-
       ijparam = elem2param[itype][jtype][jtype];
       if (rsq >= params[ijparam].cutsq) continue;
 
@@ -148,11 +160,10 @@ void PairVashishtaOMP::eval(int iifrom, int iito, ThrData * const thr)
                                evdwl,0.0,fpair,delx,dely,delz,thr);
     }
 
-    jnumm1 = jnum - 1;
+    jnumm1 = numshort - 1;
 
     for (jj = 0; jj < jnumm1; jj++) {
-      j = jlist[jj];
-      j &= NEIGHMASK;
+      j = neighshort_thr[jj];
       jtype = map[type[j]];
       ijparam = elem2param[itype][jtype][jtype];
       delr1[0] = x[j].x - xtmp;
@@ -164,9 +175,8 @@ void PairVashishtaOMP::eval(int iifrom, int iito, ThrData * const thr)
       double fjxtmp,fjytmp,fjztmp;
       fjxtmp = fjytmp = fjztmp = 0.0;
 
-      for (kk = jj+1; kk < jnum; kk++) {
-        k = jlist[kk];
-        k &= NEIGHMASK;
+      for (kk = jj+1; kk < numshort; kk++) {
+        k = neighshort_thr[kk];
         ktype = map[type[k]];
         ikparam = elem2param[itype][ktype][ktype];
         ijkparam = elem2param[itype][jtype][ktype];
@@ -200,6 +210,7 @@ void PairVashishtaOMP::eval(int iifrom, int iito, ThrData * const thr)
     f[i].y += fytmp;
     f[i].z += fztmp;
   }
+  memory->destroy(neighshort_thr);
 }
 
 /* ---------------------------------------------------------------------- */
diff --git a/src/USER-OMP/pair_vashishta_table_omp.cpp b/src/USER-OMP/pair_vashishta_table_omp.cpp
index b5769438c4709789f7bebec4416dd6eebbde04b4..d2298dca09d1a6bcec326c7eb4638ecd8a2038c2 100644
--- a/src/USER-OMP/pair_vashishta_table_omp.cpp
+++ b/src/USER-OMP/pair_vashishta_table_omp.cpp
@@ -41,20 +41,6 @@ void PairVashishtaTableOMP::compute(int eflag, int vflag)
     ev_setup(eflag,vflag);
   } else evflag = vflag_fdotr = 0;
 
-  // reallocate 3-body neighbor list if necessary
-  // NOTE: using 1000 is inefficient
-  //       could make this a LAMMPS paged neighbor list
-
-  if (atom->nlocal > neigh3BodyMax) {
-    neigh3BodyMax = atom->nmax;
-    memory->destroy(neigh3BodyCount);
-    memory->destroy(neigh3Body);
-    memory->create(neigh3BodyCount,neigh3BodyMax,
-                   "pair:vashishta:neigh3BodyCount");
-    memory->create(neigh3Body,neigh3BodyMax,1000,
-                   "pair:vashishta:neigh3Body");
-  }
-
   const int nall = atom->nlocal + atom->nghost;
   const int nthreads = comm->nthreads;
   const int inum = list->inum;
@@ -86,13 +72,13 @@ void PairVashishtaTableOMP::compute(int eflag, int vflag)
 template <int EVFLAG, int EFLAG>
 void PairVashishtaTableOMP::eval(int iifrom, int iito, ThrData * const thr)
 {
-  int i,j,k,ii,jj,kk,jnum,jnumm1;
+  int i,j,k,ii,jj,kk,jnum,jnumm1,maxshort_thr;
   tagint itag,jtag;
   int itype,jtype,ktype,ijparam,ikparam,ijkparam;
   double xtmp,ytmp,ztmp,delx,dely,delz,evdwl,fpair;
   double rsq,rsq1,rsq2;
   double delr1[3],delr2[3],fj[3],fk[3];
-  int *ilist,*jlist,*numneigh,**firstneigh;
+  int *ilist,*jlist,*numneigh,**firstneigh,*neighshort_thr;
 
   evdwl = 0.0;
 
@@ -101,13 +87,16 @@ void PairVashishtaTableOMP::eval(int iifrom, int iito, ThrData * const thr)
   const tagint * _noalias const tag = atom->tag;
   const int * _noalias const type = atom->type;
   const int nlocal = atom->nlocal;
+  const double cutshortsq = r0max*r0max;
 
   ilist = list->ilist;
   numneigh = list->numneigh;
   firstneigh = list->firstneigh;
+  maxshort_thr = maxshort;
+  memory->create(neighshort_thr,maxshort_thr,"pair_thr:neighshort_thr");
 
   double fxtmp,fytmp,fztmp;
-  
+
   // loop over full neighbor list of my atoms
 
   for (ii = iifrom; ii < iito; ++ii) {
@@ -120,35 +109,30 @@ void PairVashishtaTableOMP::eval(int iifrom, int iito, ThrData * const thr)
     ztmp = x[i].z;
     fxtmp = fytmp = fztmp = 0.0;
 
-    // reset the 3-body neighbor list
-
-    neigh3BodyCount[i] = 0;
-
     // two-body interactions, skip half of them
 
     jlist = firstneigh[i];
     jnum = numneigh[i];
+    int numshort = 0;
 
     for (jj = 0; jj < jnum; jj++) {
       j = jlist[jj];
       j &= NEIGHMASK;
-      jtag = tag[j];
-
-      jtype = map[type[j]];
 
       delx = xtmp - x[j].x;
       dely = ytmp - x[j].y;
       delz = ztmp - x[j].z;
       rsq = delx*delx + dely*dely + delz*delz;
-      ijparam = elem2param[itype][jtype][jtype];
 
-      if (rsq <= params[ijparam].cutsq2) {
-        neigh3Body[i][neigh3BodyCount[i]] = j;
-        neigh3BodyCount[i]++;
+      if (rsq < cutshortsq) {
+        neighshort_thr[numshort++] = j;
+        if (numshort >= maxshort_thr) {
+          maxshort_thr += maxshort_thr/2;
+          memory->grow(neighshort_thr,maxshort_thr,"pair_thr:neighshort_thr");
+        }
       }
 
-      if (rsq >= params[ijparam].cutsq) continue;
-
+      jtag = tag[j];
       if (itag > jtag) {
         if ((itag+jtag) % 2 == 0) continue;
       } else if (itag < jtag) {
@@ -159,6 +143,10 @@ void PairVashishtaTableOMP::eval(int iifrom, int iito, ThrData * const thr)
         if (x[j].z == ztmp && x[j].y == ytmp && x[j].x < xtmp) continue;
       }
 
+      jtype = map[type[j]];
+      ijparam = elem2param[itype][jtype][jtype];
+      if (rsq >= params[ijparam].cutsq) continue;
+
       twobody_table(params[ijparam],rsq,fpair,EFLAG,evdwl);
 
       fxtmp += delx*fpair;
@@ -172,13 +160,10 @@ void PairVashishtaTableOMP::eval(int iifrom, int iito, ThrData * const thr)
                                evdwl,0.0,fpair,delx,dely,delz,thr);
     }
 
-    jlist = neigh3Body[i];
-    jnum = neigh3BodyCount[i];
-    jnumm1 = jnum - 1;
+    jnumm1 = numshort - 1;
 
     for (jj = 0; jj < jnumm1; jj++) {
-      j = jlist[jj];
-      j &= NEIGHMASK;
+      j = neighshort_thr[jj];
       jtype = map[type[j]];
       ijparam = elem2param[itype][jtype][jtype];
       delr1[0] = x[j].x - xtmp;
@@ -190,9 +175,8 @@ void PairVashishtaTableOMP::eval(int iifrom, int iito, ThrData * const thr)
       double fjxtmp,fjytmp,fjztmp;
       fjxtmp = fjytmp = fjztmp = 0.0;
 
-      for (kk = jj+1; kk < jnum; kk++) {
-        k = jlist[kk];
-        k &= NEIGHMASK;
+      for (kk = jj+1; kk < numshort; kk++) {
+        k = neighshort_thr[kk];
         ktype = map[type[k]];
         ikparam = elem2param[itype][ktype][ktype];
         ijkparam = elem2param[itype][jtype][ktype];
@@ -226,6 +210,7 @@ void PairVashishtaTableOMP::eval(int iifrom, int iito, ThrData * const thr)
     f[i].y += fytmp;
     f[i].z += fztmp;
   }
+  memory->destroy(neighshort_thr);
 }
 
 /* ---------------------------------------------------------------------- */