diff --git a/src/USER-DPD/fix_shardlow.cpp b/src/USER-DPD/fix_shardlow.cpp
index bf8959fa9f..56597697f7 100644
--- a/src/USER-DPD/fix_shardlow.cpp
+++ b/src/USER-DPD/fix_shardlow.cpp
@@ -148,6 +148,7 @@ void FixShardlow::init()
   int irequest = neighbor->request(this,instance_me);
   neighbor->requests[irequest]->pair = 0;
   neighbor->requests[irequest]->fix  = 1;
+  neighbor->requests[irequest]->ghost= 1;
   neighbor->requests[irequest]->ssa  = 1;
 }
 
@@ -498,7 +499,7 @@ void FixShardlow::ssa_update_dpde(
 
 void FixShardlow::initial_integrate(int vflag)
 {
-  int i,ii,inum;
+  int i,ii,inum,anum;
   int *ilist;
 
   int nlocal = atom->nlocal;
@@ -531,10 +532,12 @@ void FixShardlow::initial_integrate(int vflag)
   v_t0 = (double (*)[3]) memory->smalloc(sizeof(double)*3*nghost, "FixShardlow:v_t0");
 
   inum = list->inum;
+  anum = inum + list->gnum;
   ilist = list->ilist;
 
   dtsqrt = sqrt(update->dt);
 
+  ii = 0;
   //Loop over all 14 directions (8 stages)
   for (airnum = 1; airnum <=8; airnum++){
 
@@ -549,15 +552,16 @@ void FixShardlow::initial_integrate(int vflag)
       }
     }
 
-    // Loop over neighbors of my atoms
-    for (ii = 0; ii < inum; ii++) {
+    // process neighbors in this AIR
+    while (ii < anum) {
       i = ilist[ii];
-      int start = (airnum < 2) ? 0 : list->ndxAIR_ssa[i][airnum - 2];
-      int len = list->ndxAIR_ssa[i][airnum - 1] - start;
+      if (atom->ssaAIR[i] > airnum) break; /* done with curent AIR */
+      int len = list->numneigh[i];
       if (len > 0) {
-        if (useDPDE) ssa_update_dpde(i, &(list->firstneigh[i][start]), len);
-        else ssa_update_dpd(i, &(list->firstneigh[i][start]), len);
+        if (useDPDE) ssa_update_dpde(i, &(list->firstneigh[i][0]), len);
+        else ssa_update_dpd(i, &(list->firstneigh[i][0]), len);
       }
+      ii++;
     }
 
     // Communicate the ghost deltas to the atom owners
diff --git a/src/USER-DPD/nbin_ssa.cpp b/src/USER-DPD/nbin_ssa.cpp
index 73da5e0df3..c2d780bac6 100644
--- a/src/USER-DPD/nbin_ssa.cpp
+++ b/src/USER-DPD/nbin_ssa.cpp
@@ -33,14 +33,13 @@ NBinSSA::NBinSSA(LAMMPS *lmp) : NBinStandard(lmp)
   bins_ssa = NULL;
   maxhead_ssa = 0;
   binhead_ssa = NULL;
-  gbinhead_ssa = NULL;
+  for (int i = 0; i < 9; i++) gairhead_ssa[i] = -1;
 }
 
 NBinSSA::~NBinSSA()
 {
   memory->destroy(bins_ssa);
   memory->destroy(binhead_ssa);
-  memory->destroy(gbinhead_ssa);
 }
 
 /* ----------------------------------------------------------------------
@@ -62,8 +61,11 @@ void NBinSSA::bin_atoms()
 
   last_bin = update->ntimestep;
 
+  for (i = 0; i < 9; i++) {
+    gairhead_ssa[i] = -1;
+  }
+
   for (i = 0; i < mbins; i++) {
-    gbinhead_ssa[i] = -1;
     binhead_ssa[i] = -1;
   }
 
@@ -73,19 +75,19 @@ void NBinSSA::bin_atoms()
     int bitmask = group->bitmask[includegroup];
     int nowned = atom->nlocal; // NOTE: nlocal was set to atom->nfirst above
     for (i = nall-1; i >= nowned; i--) {
-      if (ssaAIR[i] < 2) continue; // skip ghost atoms not in AIR
+      ibin = ssaAIR[i];
+      if (ibin < 2) continue; // skip ghost atoms not in AIR
       if (mask[i] & bitmask) {
-        ibin = coord2bin(x[i]);
-        bins_ssa[i] = gbinhead_ssa[ibin];
-        gbinhead_ssa[ibin] = i;
+        bins_ssa[i] = gairhead_ssa[ibin];
+        gairhead_ssa[ibin] = i;
       }
     }
   } else {
     for (i = nall-1; i >= nlocal; i--) {
-      if (ssaAIR[i] < 2) continue; // skip ghost atoms not in AIR
-      ibin = coord2bin(x[i]);
-      bins_ssa[i] = gbinhead_ssa[ibin];
-      gbinhead_ssa[ibin] = i;
+      ibin = ssaAIR[i];
+      if (ibin < 2) continue; // skip ghost atoms not in AIR
+      bins_ssa[i] = gairhead_ssa[ibin];
+      gairhead_ssa[ibin] = i;
     }
   }
   for (i = nlocal-1; i >= 0; i--) {
@@ -103,10 +105,8 @@ void NBinSSA::bin_atoms_setup(int nall)
 
   if (mbins > maxhead_ssa) {
     maxhead_ssa = mbins;
-    memory->destroy(gbinhead_ssa);
     memory->destroy(binhead_ssa);
     memory->create(binhead_ssa,maxhead_ssa,"binhead_ssa");
-    memory->create(gbinhead_ssa,maxhead_ssa,"gbinhead_ssa");
   }
 
   if (nall > maxbin_ssa) {
@@ -125,7 +125,6 @@ bigint NBinSSA::memory_usage()
   if (maxbin_ssa) bytes += memory->usage(bins_ssa,maxbin_ssa);
   if (maxhead_ssa) {
     bytes += memory->usage(binhead_ssa,maxhead_ssa);
-    bytes += memory->usage(gbinhead_ssa,maxhead_ssa);
   }
   return bytes;
 }
diff --git a/src/USER-DPD/nbin_ssa.h b/src/USER-DPD/nbin_ssa.h
index f0699b3a7a..5a2562d305 100644
--- a/src/USER-DPD/nbin_ssa.h
+++ b/src/USER-DPD/nbin_ssa.h
@@ -32,7 +32,7 @@ class NBinSSA : public NBinStandard {
   int *bins_ssa;             // index of next atom in each bin
   int maxbin_ssa;            // size of bins_ssa array
   int *binhead_ssa;          // index of 1st local atom in each bin
-  int *gbinhead_ssa;         // index of 1st ghost atom in each bin
+  int gairhead_ssa[9];       // index of 1st ghost atom in each AIR
   int maxhead_ssa;           // size of binhead_ssa and gbinhead_ssa arrays
 
   NBinSSA(class LAMMPS *);
diff --git a/src/USER-DPD/npair_half_bin_newton_ssa.cpp b/src/USER-DPD/npair_half_bin_newton_ssa.cpp
index fd67b66e9b..4c9dc95308 100644
--- a/src/USER-DPD/npair_half_bin_newton_ssa.cpp
+++ b/src/USER-DPD/npair_half_bin_newton_ssa.cpp
@@ -32,12 +32,6 @@
 
 using namespace LAMMPS_NS;
 
-// allocate space for static class variable
-// prototype for non-class function
-
-static int *ssaAIRptr;
-static int cmp_ssaAIR(const void *, const void *);
-
 /* ---------------------------------------------------------------------- */
 
 NPairHalfBinNewtonSSA::NPairHalfBinNewtonSSA(LAMMPS *lmp) : NPair(lmp) {}
@@ -64,9 +58,7 @@ void NPairHalfBinNewtonSSA::build(NeighList *list)
   tagint **special = atom->special;
   int **nspecial = atom->nspecial;
   int nlocal = atom->nlocal;
-  int nall = nlocal + atom->nghost;
   if (includegroup) nlocal = atom->nfirst;
-  int *ssaAIR = atom->ssaAIR;
 
   int *molindex = atom->molindex;
   int *molatom = atom->molatom;
@@ -89,16 +81,18 @@ void NPairHalfBinNewtonSSA::build(NeighList *list)
   if (!nb_ssa) error->one(FLERR, "NBin wasn't a NBinSSA object");
   int *bins_ssa = nb_ssa->bins_ssa;
   int *binhead_ssa = nb_ssa->binhead_ssa;
-  int *gbinhead_ssa = nb_ssa->gbinhead_ssa;
+  int *gairhead_ssa = &(nb_ssa->gairhead_ssa[0]);
 
   int inum = 0;
+  int gnum = 0;
+  int xbin,ybin,zbin,xbin2,ybin2,zbin2;
+  int **stencilxyz = ns_ssa->stencilxyz;
 
   ipage->reset();
 
   // loop over owned atoms, storing half of the neighbors
 
   for (i = 0; i < nlocal; i++) {
-    int AIRct[8] = { 0 };
     n = 0;
     neighptr = ipage->vget();
 
@@ -175,51 +169,6 @@ void NPairHalfBinNewtonSSA::build(NeighList *list)
         }
       }
     }
-    AIRct[0] = n;
-
-    // loop over AIR ghost atoms in all bins in "full" stencil
-    // Note: the non-AIR ghost atoms have already been filtered out
-    // That is a significant time savings because of the "full" stencil
-    // Note2: only non-pure locals can have ghosts as neighbors
-
-    if (ssaAIR[i] == 1) for (k = 0; k < nstencil_full; k++) {
-      for (j = gbinhead_ssa[ibin+stencil[k]]; j >= 0;
-           j = bins_ssa[j]) {
-
-        jtype = type[j];
-        if (exclude && exclusion(i,j,itype,jtype,mask,molecule)) continue;
-
-        delx = xtmp - x[j][0];
-        dely = ytmp - x[j][1];
-        delz = ztmp - x[j][2];
-        rsq = delx*delx + dely*dely + delz*delz;
-
-        if (rsq <= cutneighsq[itype][jtype]) {
-          if (molecular) {
-            if (!moltemplate)
-              which = find_special(special[i],nspecial[i],tag[j]);
-            else if (imol >= 0)
-              which = find_special(onemols[imol]->special[iatom],
-                                   onemols[imol]->nspecial[iatom],
-                                   tag[j]-tagprev);
-            else which = 0;
-            if (which == 0) {
-              neighptr[n++] = j;
-              ++(AIRct[ssaAIR[j] - 1]);
-            } else if (domain->minimum_image_check(delx,dely,delz)) {
-              neighptr[n++] = j;
-              ++(AIRct[ssaAIR[j] - 1]);
-            } else if (which > 0) {
-              neighptr[n++] = j ^ (which << SBBITS);
-              ++(AIRct[ssaAIR[j] - 1]);
-            }
-          } else {
-            neighptr[n++] = j;
-            ++(AIRct[ssaAIR[j] - 1]);
-          }
-        }
-      }
-    }
 
     ilist[inum++] = i;
     firstneigh[i] = neighptr;
@@ -227,34 +176,74 @@ void NPairHalfBinNewtonSSA::build(NeighList *list)
     ipage->vgot(n);
     if (ipage->status())
       error->one(FLERR,"Neighbor list overflow, boost neigh_modify one");
-
-    // sort the ghosts in the neighbor list by their ssaAIR number
-
-    ssaAIRptr = atom->ssaAIR;
-    qsort(&(neighptr[AIRct[0]]), n - AIRct[0], sizeof(int), cmp_ssaAIR);
-
-    // do a prefix sum on the counts to turn them into indexes
-
-    list->ndxAIR_ssa[i][0] = AIRct[0];
-    for (int ndx = 1; ndx < 8; ++ndx) {
-      list->ndxAIR_ssa[i][ndx] = AIRct[ndx] + list->ndxAIR_ssa[i][ndx - 1];
-    }
   }
 
   list->inum = inum;
+
+  // loop over AIR ghost atoms, storing their local neighbors
+  // since these are ghosts, must check if stencil bin is out of bounds
+  for (int airnum = 2; airnum <= 8; airnum++) {
+    for (i = gairhead_ssa[airnum]; i >= 0; i = bins_ssa[i]) {
+      n = 0;
+      neighptr = ipage->vget();
+
+      itype = type[i];
+      xtmp = x[i][0];
+      ytmp = x[i][1];
+      ztmp = x[i][2];
+      if (moltemplate) {
+        imol = molindex[i];
+        iatom = molatom[i];
+        tagprev = tag[i] - iatom - 1;
+      }
+
+      ibin = coord2bin(x[i],xbin,ybin,zbin);
+
+      // loop over AIR ghost atoms in all bins in "full" stencil
+      // Note: the non-AIR ghost atoms have already been filtered out
+      for (k = 0; k < nstencil_full; k++) {
+        xbin2 = xbin + stencilxyz[k][0];
+        ybin2 = ybin + stencilxyz[k][1];
+        zbin2 = zbin + stencilxyz[k][2];
+        // since we only care about ghost to local neighbors, these "bounds" could be inset
+        if (xbin2 < 0 || xbin2 >= mbinx ||
+            ybin2 < 0 || ybin2 >= mbiny ||
+            zbin2 < 0 || zbin2 >= mbinz) continue;
+        for (j = binhead_ssa[ibin+stencil[k]]; j >= 0; j = bins_ssa[j]) {
+
+          jtype = type[j];
+          if (exclude && exclusion(i,j,itype,jtype,mask,molecule)) continue;
+
+          delx = xtmp - x[j][0];
+          dely = ytmp - x[j][1];
+          delz = ztmp - x[j][2];
+          rsq = delx*delx + dely*dely + delz*delz;
+
+          if (rsq <= cutneighsq[itype][jtype]) {
+            if (molecular) {
+              if (!moltemplate)
+                which = find_special(special[i],nspecial[i],tag[j]);
+              else if (imol >= 0)
+                which = find_special(onemols[imol]->special[iatom],
+                                     onemols[imol]->nspecial[iatom],
+                                     tag[j]-tagprev);
+              else which = 0;
+              if (which == 0) neighptr[n++] = j;
+              else if (domain->minimum_image_check(delx,dely,delz))
+                neighptr[n++] = j;
+              else if (which > 0) neighptr[n++] = j ^ (which << SBBITS);
+            } else neighptr[n++] = j;
+          }
+        }
+      }
+
+      if (n > 0) ilist[inum + (gnum++)] = i;
+      firstneigh[i] = neighptr;
+      numneigh[i] = n;
+      ipage->vgot(n);
+      if (ipage->status())
+        error->one(FLERR,"Neighbor (ghost) list overflow, boost neigh_modify one");
+    }
+  }
+  list->gnum = gnum;
 }
-
-/* ----------------------------------------------------------------------
-   comparison function invoked by qsort()
-   accesses static class member ssaAIRptr, set before call to qsort()
-------------------------------------------------------------------------- */
-
-static int cmp_ssaAIR(const void *iptr, const void *jptr)
-{
-  int i = NEIGHMASK & *((int *) iptr);
-  int j = NEIGHMASK & *((int *) jptr);
-  if (ssaAIRptr[i] < ssaAIRptr[j]) return -1;
-  if (ssaAIRptr[i] > ssaAIRptr[j]) return 1;
-  return 0;
-}
-
diff --git a/src/USER-DPD/npair_half_bin_newton_ssa.h b/src/USER-DPD/npair_half_bin_newton_ssa.h
index 13347b33b0..c9ccbc4bd9 100644
--- a/src/USER-DPD/npair_half_bin_newton_ssa.h
+++ b/src/USER-DPD/npair_half_bin_newton_ssa.h
@@ -15,7 +15,7 @@
 
 NPairStyle(half/bin/newton/ssa,
            NPairHalfBinNewtonSSA,
-           NP_HALF | NP_BIN | NP_NEWTON | NP_ORTHO | NP_SSA)
+           NP_HALF | NP_BIN | NP_NEWTON | NP_ORTHO | NP_SSA | NP_GHOST)
 
 #else
 
diff --git a/src/USER-DPD/npair_halffull_newton_ssa.cpp b/src/USER-DPD/npair_halffull_newton_ssa.cpp
index 2c9de3e50f..d0be1685b6 100644
--- a/src/USER-DPD/npair_halffull_newton_ssa.cpp
+++ b/src/USER-DPD/npair_halffull_newton_ssa.cpp
@@ -64,6 +64,10 @@ void NPairHalffullNewtonSSA::build(NeighList *list)
   int inum_full = list->listfull->inum;
 
   int inum = 0;
+
+  error->one(FLERR,"NPairHalffullNewtonSSA not yet implemented for ghosts with neighbors.");
+  return;
+
   ipage->reset();
 
   // loop over parent full list
diff --git a/src/USER-DPD/nstencil_half_bin_2d_newton_ssa.cpp b/src/USER-DPD/nstencil_half_bin_2d_newton_ssa.cpp
index df379a109a..254339bffc 100644
--- a/src/USER-DPD/nstencil_half_bin_2d_newton_ssa.cpp
+++ b/src/USER-DPD/nstencil_half_bin_2d_newton_ssa.cpp
@@ -46,8 +46,12 @@ void NStencilHalfBin2dNewtonSSA::create()
   for (j = 0; j <= sy; j++)
     for (i = -sx; i <= sx; i++)
       if (j > 0 || (j == 0 && i > 0))
-        if (bin_distance(i,j,0) < cutneighmaxsq)
+        if (bin_distance(i,j,0) < cutneighmaxsq) {
+          stencilxyz[pos][0] = i;
+          stencilxyz[pos][1] = j;
+          stencilxyz[pos][2] = 0;
           stencil[pos++] = j*mbinx + i;
+        }
 
   nstencil_half = pos; // record where normal half stencil ends
 
@@ -56,8 +60,12 @@ void NStencilHalfBin2dNewtonSSA::create()
   for (j = -sy; j <= 0; j++)
     for (i = -sx; i <= sx; i++) {
       if (j == 0 && i > 0) continue;
-      if (bin_distance(i,j,0) < cutneighmaxsq)
+      if (bin_distance(i,j,0) < cutneighmaxsq) {
+        stencilxyz[pos][0] = i;
+        stencilxyz[pos][1] = j;
+        stencilxyz[pos][2] = 0;
         stencil[pos++] = j*mbinx + i;
+      }
     }
 
   nstencil = pos; // record where full stencil ends
diff --git a/src/USER-DPD/nstencil_half_bin_2d_newton_ssa.h b/src/USER-DPD/nstencil_half_bin_2d_newton_ssa.h
index 30901bb3e2..1d5cc3f6b2 100644
--- a/src/USER-DPD/nstencil_half_bin_2d_newton_ssa.h
+++ b/src/USER-DPD/nstencil_half_bin_2d_newton_ssa.h
@@ -15,7 +15,7 @@
 
 NStencilStyle(half/bin/2d/newton/ssa,
               NStencilHalfBin2dNewtonSSA,
-              NS_HALF | NS_BIN | NS_2D | NS_NEWTON | NS_SSA | NS_ORTHO)
+              NS_HALF | NS_BIN | NS_2D | NS_NEWTON | NS_SSA | NS_ORTHO | NS_GHOST)
 
 #else
 
diff --git a/src/USER-DPD/nstencil_half_bin_3d_newton_ssa.cpp b/src/USER-DPD/nstencil_half_bin_3d_newton_ssa.cpp
index 76c9931ab2..1e2c18c66a 100644
--- a/src/USER-DPD/nstencil_half_bin_3d_newton_ssa.cpp
+++ b/src/USER-DPD/nstencil_half_bin_3d_newton_ssa.cpp
@@ -47,8 +47,12 @@ void NStencilHalfBin3dNewtonSSA::create()
     for (j = -sy; j <= sy; j++)
       for (i = -sx; i <= sx; i++)
         if (k > 0 || j > 0 || (j == 0 && i > 0))
-          if (bin_distance(i,j,k) < cutneighmaxsq)
+          if (bin_distance(i,j,k) < cutneighmaxsq) {
+            stencilxyz[pos][0] = i;
+            stencilxyz[pos][1] = j;
+            stencilxyz[pos][2] = k;
             stencil[pos++] = k*mbiny*mbinx + j*mbinx + i;
+          }
 
   nstencil_half = pos; // record where normal half stencil ends
 
@@ -57,8 +61,12 @@ void NStencilHalfBin3dNewtonSSA::create()
   for (k = -sz; k < 0; k++)
     for (j = -sy; j <= sy; j++)
       for (i = -sx; i <= sx; i++)
-        if (bin_distance(i,j,k) < cutneighmaxsq)
+        if (bin_distance(i,j,k) < cutneighmaxsq) {
+          stencilxyz[pos][0] = i;
+          stencilxyz[pos][1] = j;
+          stencilxyz[pos][2] = k;
           stencil[pos++] = k*mbiny*mbinx + j*mbinx + i;
+        }
 
   // For k==0, make sure to skip already included bins
 
@@ -66,8 +74,12 @@ void NStencilHalfBin3dNewtonSSA::create()
   for (j = -sy; j <= 0; j++)
     for (i = -sx; i <= sx; i++) {
       if (j == 0 && i > 0) continue;
-      if (bin_distance(i,j,k) < cutneighmaxsq)
+      if (bin_distance(i,j,k) < cutneighmaxsq) {
+        stencilxyz[pos][0] = i;
+        stencilxyz[pos][1] = j;
+        stencilxyz[pos][2] = k;
         stencil[pos++] = k*mbiny*mbinx + j*mbinx + i;
+      }
     }
 
   nstencil = pos; // record where full stencil ends
diff --git a/src/USER-DPD/nstencil_half_bin_3d_newton_ssa.h b/src/USER-DPD/nstencil_half_bin_3d_newton_ssa.h
index 7765b256d3..450a696e46 100644
--- a/src/USER-DPD/nstencil_half_bin_3d_newton_ssa.h
+++ b/src/USER-DPD/nstencil_half_bin_3d_newton_ssa.h
@@ -15,7 +15,7 @@
 
 NStencilStyle(half/bin/3d/newton/ssa,
               NStencilHalfBin3dNewtonSSA,
-              NS_HALF | NS_BIN | NS_3D | NS_NEWTON | NS_SSA | NS_ORTHO)
+              NS_HALF | NS_BIN | NS_3D | NS_NEWTON | NS_SSA | NS_ORTHO | NS_GHOST)
 
 #else
 
diff --git a/src/USER-DPD/nstencil_ssa.h b/src/USER-DPD/nstencil_ssa.h
index 9fcd19ee26..e6dfce60f4 100644
--- a/src/USER-DPD/nstencil_ssa.h
+++ b/src/USER-DPD/nstencil_ssa.h
@@ -20,7 +20,7 @@ namespace LAMMPS_NS {
 
 class NStencilSSA : public NStencil {
  public:
-  NStencilSSA(class LAMMPS *lmp) : NStencil(lmp) { }
+  NStencilSSA(class LAMMPS *lmp) : NStencil(lmp) { xyzflag = 1; }
   ~NStencilSSA() {}
   virtual void create() = 0;