diff --git a/src/.gitignore b/src/.gitignore
index fe23bc1f55446fc6873390ef1fca32af96d94d8d..a9dc1c8fa14ee901ca24ebf2aea7581fcf6e29e9 100644
--- a/src/.gitignore
+++ b/src/.gitignore
@@ -547,6 +547,7 @@
 /fix_shake.h
 /fix_shardlow.cpp
 /fix_shardlow.h
+/random_external_state.h
 /fix_smd.cpp
 /fix_smd.h
 /fix_species.cpp
diff --git a/src/KOKKOS/fix_shardlow_kokkos.cpp b/src/KOKKOS/fix_shardlow_kokkos.cpp
index cc1bd6bede9fabe94f270464cd1ead9558ecad8e..571f4880230caefaaed8c2561f16f0681807ac54 100644
--- a/src/KOKKOS/fix_shardlow_kokkos.cpp
+++ b/src/KOKKOS/fix_shardlow_kokkos.cpp
@@ -61,6 +61,7 @@
 
 using namespace LAMMPS_NS;
 using namespace FixConst;
+using namespace random_external_state;
 
 #define EPSILON 1.0e-10
 #define EPSILON_SQUARED ((EPSILON) * (EPSILON))
@@ -89,10 +90,6 @@ FixShardlowKokkos<DeviceType>::FixShardlowKokkos(LAMMPS *lmp, int narg, char **a
 //   if(k_pairDPDE){
     comm_forward = 3;
     comm_reverse = 5;
-    maxRNG = 0;
-#ifdef DPD_USE_RAN_MARS
-    pp_random = NULL;
-#endif
 //   } else {
 //     comm_forward = 3;
 //     comm_reverse = 3;
@@ -121,13 +118,6 @@ template<class DeviceType>
 FixShardlowKokkos<DeviceType>::~FixShardlowKokkos()
 {
   ghostmax = 0;
-#ifdef DPD_USE_RAN_MARS
-  if (pp_random) {
-    for (int i = 1; i < maxRNG; ++i) delete pp_random[i];
-    delete[] pp_random;
-    pp_random = NULL;
-  }
-#endif
 }
 
 /* ---------------------------------------------------------------------- */
@@ -279,11 +269,7 @@ void FixShardlowKokkos<DeviceType>::ssa_update_dpd(
   int start_ii, int count, int id
 )
 {
-#ifdef DPD_USE_RAN_MARS
-  class RanMars *pRNG = pp_random[id];
-#else
-  rand_type rand_gen = rand_pool.get_state(id);
-#endif
+  es_RNG_t RNGstate = d_rand_state(id);
 
   int ct = count;
   int ii = start_ii;
@@ -345,12 +331,7 @@ void FixShardlowKokkos<DeviceType>::ssa_update_dpd(
         double halfsigma_ij = STACKPARAMS?m_params[itype][jtype].halfsigma:params(itype,jtype).halfsigma;
         double halfgamma_ij = halfsigma_ij*halfsigma_ij*boltz_inv*theta_ij_inv;
 
-        double sigmaRand = halfsigma_ij*wr*dtsqrt*ftm2v *
-#ifdef DPD_USE_RAN_MARS
-            pRNG->gaussian();
-#else
-            rand_gen.normal();
-#endif
+        double sigmaRand = halfsigma_ij*wr*dtsqrt*ftm2v * es_normal(RNGstate);
 
         const double mass_j = masses(massPerI ? j : jtype);
         double massinv_j = 1.0 / mass_j;
@@ -412,9 +393,7 @@ void FixShardlowKokkos<DeviceType>::ssa_update_dpd(
     v(i, 2) = vzi;
   }
 
-#ifndef DPD_USE_RAN_MARS
-  rand_pool.free_state(rand_gen);
-#endif
+  d_rand_state(id) = RNGstate;
 }
 #endif
 
@@ -431,11 +410,7 @@ void FixShardlowKokkos<DeviceType>::ssa_update_dpde(
   int start_ii, int count, int id
 ) const
 {
-#ifdef DPD_USE_RAN_MARS
-  class RanMars *pRNG = pp_random[id];
-#else
-  rand_type rand_gen = rand_pool.get_state(id);
-#endif
+  es_RNG_t RNGstate = d_rand_state(id);
 
   int ct = count;
   int ii = start_ii;
@@ -506,12 +481,7 @@ void FixShardlowKokkos<DeviceType>::ssa_update_dpde(
         double halfsigma_ij = STACKPARAMS?m_params[itype][jtype].halfsigma:params(itype,jtype).halfsigma;
         double halfgamma_ij = halfsigma_ij*halfsigma_ij*boltz_inv*theta_ij_inv;
 
-        double sigmaRand = halfsigma_ij*wr*dtsqrt*ftm2v *
-#ifdef DPD_USE_RAN_MARS
-            pRNG->gaussian();
-#else
-            rand_gen.normal();
-#endif
+        double sigmaRand = halfsigma_ij*wr*dtsqrt*ftm2v * es_normal(RNGstate);
 
         const double mass_j = masses(massPerI ? j : jtype);
         double mass_ij_div_neg4_ftm2v = mass_j*mass_i_div_neg4_ftm2v;
@@ -520,12 +490,7 @@ void FixShardlowKokkos<DeviceType>::ssa_update_dpde(
         // Compute uCond
         double kappa_ij = STACKPARAMS?m_params[itype][jtype].kappa:params(itype,jtype).kappa;
         double alpha_ij = STACKPARAMS?m_params[itype][jtype].alpha:params(itype,jtype).alpha;
-        double del_uCond = alpha_ij*wr*dtsqrt *
-#ifdef DPD_USE_RAN_MARS
-            pRNG->gaussian();
-#else
-            rand_gen.normal();
-#endif
+        double del_uCond = alpha_ij*wr*dtsqrt * es_normal(RNGstate);
 
         del_uCond += kappa_ij*(theta_i_inv - theta_j_inv)*wdt;
         uCond[j] -= del_uCond;
@@ -601,9 +566,7 @@ void FixShardlowKokkos<DeviceType>::ssa_update_dpde(
     ii++;
   }
 
-#ifndef DPD_USE_RAN_MARS
-  rand_pool.free_state(rand_gen);
-#endif
+  d_rand_state(id) = RNGstate;
 }
 
 
@@ -644,20 +607,16 @@ void FixShardlowKokkos<DeviceType>::initial_integrate(int vflag)
     maxWorkItemCt = (int) ssa_gitemLoc.dimension_1();
   }
   if (maxWorkItemCt > maxRNG) {
-#ifdef DPD_USE_RAN_MARS
-    if (pp_random) {
-      for (int i = 1; i < maxRNG; ++i) delete pp_random[i];
-      delete[] pp_random;
-      pp_random = NULL;
-    }
-    pp_random = new RanMars*[maxWorkItemCt];
-    for (int i = 1; i < maxWorkItemCt; ++i) {
-      pp_random[i] = new RanMars(lmp, k_pairDPDE->seed + comm->me + comm->nprocs*i);
+    es_RNG_t serial_rand_state;
+    es_init(serial_rand_state, pairDPDE->seed + comm->me);
+
+    d_rand_state = es_RNGs_type("Kokkos::fix_shardlow::rand_state",maxWorkItemCt);
+    typename es_RNGs_type::HostMirror h_rand_state = create_mirror_view(d_rand_state);
+    for (int i = 0; i < maxWorkItemCt; ++i) {
+      es_genNextParallelState(serial_rand_state, h_rand_state(i));
     }
-    pp_random[0] = k_pairDPDE->random;
-#else
-    rand_pool.init(k_pairDPDE->seed + comm->me, maxWorkItemCt);
-#endif
+    deep_copy(d_rand_state,h_rand_state);
+
     maxRNG = maxWorkItemCt;
   }
 
diff --git a/src/KOKKOS/fix_shardlow_kokkos.h b/src/KOKKOS/fix_shardlow_kokkos.h
index 70dccf2e2d6b9faeba8bae406d26c986d0fac3d5..4e87de6910d7724a0bc1a736cb15b2b5d56c1049 100644
--- a/src/KOKKOS/fix_shardlow_kokkos.h
+++ b/src/KOKKOS/fix_shardlow_kokkos.h
@@ -93,17 +93,6 @@ class FixShardlowKokkos : public FixShardlow {
 #endif
   PairDPDfdtEnergyKokkos<DeviceType> *k_pairDPDE;
 
-  int maxRNG;
-#ifdef DPD_USE_RAN_MARS
-  class RanMars **pp_random;
-#elif defined(DPD_USE_Random_XorShift1024)
-  Kokkos::Random_XorShift1024_Pool<DeviceType> rand_pool;
-  typedef typename Kokkos::Random_XorShift1024_Pool<DeviceType>::generator_type rand_type;
-#else
-  Kokkos::Random_XorShift64_Pool<DeviceType> rand_pool;
-  typedef typename Kokkos::Random_XorShift64_Pool<DeviceType>::generator_type rand_type;
-#endif
-
   Kokkos::DualView<params_ssa**,Kokkos::LayoutRight,DeviceType> k_params;
   typename Kokkos::DualView<params_ssa**,
     Kokkos::LayoutRight,DeviceType>::t_dev_const_um params;
@@ -127,6 +116,10 @@ class FixShardlowKokkos : public FixShardlow {
   typename AT::t_float_1d_randomread masses;
   typename AT::t_efloat_1d dpdTheta;
 
+  // Storage for the es_RNG state variables
+  typedef Kokkos::View<random_external_state::es_RNG_t*,DeviceType> es_RNGs_type;
+  es_RNGs_type d_rand_state;
+
   double dtsqrt; // = sqrt(update->dt);
   int ghostmax;
   int nlocal, nghost;
diff --git a/src/USER-DPD/fix_shardlow.cpp b/src/USER-DPD/fix_shardlow.cpp
index cec53ab15fcec09a05cadd48f9b8c2f1b4ebab7e..39b7ba22246065e6616cac705c48e87fc794a02a 100644
--- a/src/USER-DPD/fix_shardlow.cpp
+++ b/src/USER-DPD/fix_shardlow.cpp
@@ -48,7 +48,6 @@
 #include "neighbor.h"
 #include "neigh_list.h"
 #include "neigh_request.h"
-#include "random_mars.h"
 #include "memory.h"
 #include "domain.h"
 #include "modify.h"
@@ -60,6 +59,7 @@
 
 using namespace LAMMPS_NS;
 using namespace FixConst;
+using namespace random_external_state;
 
 #define EPSILON 1.0e-10
 #define EPSILON_SQUARED ((EPSILON) * (EPSILON))
@@ -87,6 +87,7 @@ static const char cite_fix_shardlow[] =
 
 FixShardlow::FixShardlow(LAMMPS *lmp, int narg, char **arg) :
   Fix(lmp, narg, arg), pairDPD(NULL), pairDPDE(NULL), v_t0(NULL)
+  ,rand_state(NULL)
 {
   if (lmp->citeme) lmp->citeme->add(cite_fix_shardlow);
 
@@ -99,6 +100,7 @@ FixShardlow::FixShardlow(LAMMPS *lmp, int narg, char **arg) :
   if (pairDPDE == NULL)
     pairDPDE = (PairDPDfdtEnergy *) force->pair_match("dpd/fdt/energy/kk",1);
 
+  maxRNG = 0;
   if(pairDPDE){
     comm_forward = 3;
     comm_reverse = 5;
@@ -116,6 +118,8 @@ FixShardlow::FixShardlow(LAMMPS *lmp, int narg, char **arg) :
 
 FixShardlow::~FixShardlow()
 {
+  memory->destroy(rand_state);
+  maxRNG = 0;
 }
 
 /* ---------------------------------------------------------------------- */
@@ -173,12 +177,12 @@ void FixShardlow::setup(int vflag)
    NOTE: only implemented for orthogonal boxes, not triclinic
 ------------------------------------------------------------------------- */
 void FixShardlow::ssa_update_dpd(
-  int i,
-  int *jlist,
-  int jlen
+  int start_ii,
+  int count,
+  int id
 )
 {
-  class RanMars *pRNG;
+  es_RNG_t RNGstate = rand_state[id];
   double **x = atom->x;
   double **v = atom->v;
   double *rmass = atom->rmass;
@@ -192,6 +196,16 @@ void FixShardlow::ssa_update_dpd(
 
   const double dt     = update->dt;
 
+  int ct = count;
+  int ii = start_ii;
+
+while (ct-- > 0) {
+  const int i = list->ilist[ii];
+  const int *jlist = list->firstneigh[ii];
+  const int jlen = list->numneigh[ii];
+  ii++;
+  if (jlen <= 0) continue;
+
   const double xtmp = x[i][0];
   const double ytmp = x[i][1];
   const double ztmp = x[i][2];
@@ -203,7 +217,6 @@ void FixShardlow::ssa_update_dpd(
 
   int itype = type[i];
 
-  pRNG = pairDPD->random;
   cut2_i = pairDPD->cutsq[itype];
   cut_i  = pairDPD->cut[itype];
   sigma_i = pairDPD->sigma[itype];
@@ -254,7 +267,7 @@ void FixShardlow::ssa_update_dpd(
       double halfsigma_ij = 0.5*sigma_i[jtype];
       double halfgamma_ij = halfsigma_ij*halfsigma_ij*boltz_inv*theta_ij_inv;
 
-      double sigmaRand = halfsigma_ij*wr*dtsqrt*ftm2v * pRNG->gaussian();
+      double sigmaRand = halfsigma_ij*wr*dtsqrt*ftm2v * es_normal(RNGstate);
 
       double mass_j = (rmass) ? rmass[j] : mass[jtype];
       double massinv_j = 1.0 / mass_j;
@@ -316,6 +329,9 @@ void FixShardlow::ssa_update_dpd(
   v[i][2] = vzi;
 }
 
+  rand_state[id] = RNGstate;
+}
+
 /* ----------------------------------------------------------------------
    Perform the stochastic integration and Shardlow update for constant energy
    Allow for both per-type and per-atom mass
@@ -323,12 +339,12 @@ void FixShardlow::ssa_update_dpd(
    NOTE: only implemented for orthogonal boxes, not triclinic
 ------------------------------------------------------------------------- */
 void FixShardlow::ssa_update_dpde(
-  int i,
-  int *jlist,
-  int jlen
+  int start_ii,
+  int count,
+  int id
 )
 {
-  class RanMars *pRNG;
+  es_RNG_t RNGstate = rand_state[id];
   double **x = atom->x;
   double **v = atom->v;
   double *rmass = atom->rmass;
@@ -346,6 +362,16 @@ void FixShardlow::ssa_update_dpde(
 
   const double dt     = update->dt;
 
+  int ct = count;
+  int ii = start_ii;
+
+while (ct-- > 0) {
+  const int i = list->ilist[ii];
+  const int *jlist = list->firstneigh[ii];
+  const int jlen = list->numneigh[ii];
+  ii++;
+  if (jlen <= 0) continue;
+
   const double xtmp = x[i][0];
   const double ytmp = x[i][1];
   const double ztmp = x[i][2];
@@ -359,7 +385,6 @@ void FixShardlow::ssa_update_dpde(
   double uCond_i = uCond[i];
   int itype = type[i];
 
-  pRNG = pairDPDE->random;
   cut2_i = pairDPDE->cutsq[itype];
   cut_i  = pairDPDE->cut[itype];
   sigma_i = pairDPDE->sigma[itype];
@@ -415,7 +440,7 @@ void FixShardlow::ssa_update_dpde(
       double halfsigma_ij = 0.5*sigma_i[jtype];
       double halfgamma_ij = halfsigma_ij*halfsigma_ij*boltz_inv*theta_ij_inv;
 
-      double sigmaRand = halfsigma_ij*wr*dtsqrt*ftm2v * pRNG->gaussian();
+      double sigmaRand = halfsigma_ij*wr*dtsqrt*ftm2v * es_normal(RNGstate);
 
       double mass_j = (rmass) ? rmass[j] : mass[jtype];
       double mass_ij_div_neg4_ftm2v = mass_j*mass_i_div_neg4_ftm2v;
@@ -424,7 +449,7 @@ void FixShardlow::ssa_update_dpde(
       // Compute uCond
       double kappa_ij = kappa_i[jtype];
       double alpha_ij = sqrt(boltz2*kappa_ij);
-      double del_uCond = alpha_ij*wr*dtsqrt * pRNG->gaussian();
+      double del_uCond = alpha_ij*wr*dtsqrt * es_normal(RNGstate);
 
       del_uCond += kappa_ij*(theta_i_inv - theta_j_inv)*wdt;
       uCond[j] -= del_uCond;
@@ -499,6 +524,9 @@ void FixShardlow::ssa_update_dpde(
   uCond[i] = uCond_i;
 }
 
+  rand_state[id] = RNGstate;
+}
+
 void FixShardlow::initial_integrate(int vflag)
 {
   int i,ii,inum;
@@ -507,7 +535,6 @@ void FixShardlow::initial_integrate(int vflag)
   int nlocal = atom->nlocal;
   int nghost = atom->nghost;
 
-  int airnum;
   const bool useDPDE = (pairDPDE != NULL);
 
   // NOTE: this logic is specific to orthogonal boxes, not triclinic
@@ -530,6 +557,31 @@ void FixShardlow::initial_integrate(int vflag)
     error->one(FLERR, msg);
   }
 
+  NPairHalfBinNewtonSSA *np_ssa = dynamic_cast<NPairHalfBinNewtonSSA*>(list->np);
+  if (!np_ssa) error->one(FLERR, "NPair wasn't a NPairHalfBinNewtonSSA object");
+  int ssa_phaseCt = np_ssa->ssa_phaseCt;
+  int *ssa_phaseLen = np_ssa->ssa_phaseLen;
+  int **ssa_itemLoc = np_ssa->ssa_itemLoc;
+  int **ssa_itemLen = np_ssa->ssa_itemLen;
+  int ssa_gphaseCt = np_ssa->ssa_gphaseCt;
+  int *ssa_gphaseLen = np_ssa->ssa_gphaseLen;
+  int **ssa_gitemLoc = np_ssa->ssa_gitemLoc;
+  int **ssa_gitemLen = np_ssa->ssa_gitemLen;
+
+  int maxWorkItemCt = np_ssa->ssa_maxPhaseLen;
+  if (maxWorkItemCt > maxRNG) {
+    uint64_t my_seed = comm->me + (useDPDE ? pairDPDE->seed : pairDPD->seed);
+    es_RNG_t serial_rand_state;
+    es_init(serial_rand_state, my_seed);
+
+    memory->grow(rand_state, maxWorkItemCt, "FixShardlow:rand_state");
+    for (int i = 0; i < maxWorkItemCt; ++i) {
+      es_genNextParallelState(serial_rand_state, rand_state[i]);
+    }
+
+    maxRNG = maxWorkItemCt;
+  }
+
 #ifdef DEBUG_SSA_PAIR_CT
   for (int i = 0; i < 2; ++i)
     for (int j = 0; j < 3; ++j)
@@ -545,13 +597,6 @@ void FixShardlow::initial_integrate(int vflag)
 
   dtsqrt = sqrt(update->dt);
 
-  NPairHalfBinNewtonSSA *np_ssa = dynamic_cast<NPairHalfBinNewtonSSA*>(list->np);
-  if (!np_ssa) error->one(FLERR, "NPair wasn't a NPairHalfBinNewtonSSA object");
-  int ssa_phaseCt = np_ssa->ssa_phaseCt;
-  int *ssa_phaseLen = np_ssa->ssa_phaseLen;
-  int **ssa_itemLoc = np_ssa->ssa_itemLoc;
-  int **ssa_itemLen = np_ssa->ssa_itemLen;
-
   // process neighbors in the local AIR
   for (int workPhase = 0; workPhase < ssa_phaseCt; ++workPhase) {
     int workItemCt = ssa_phaseLen[workPhase];
@@ -559,22 +604,14 @@ void FixShardlow::initial_integrate(int vflag)
     for (int workItem = 0; workItem < workItemCt; ++workItem) {
       int ct = ssa_itemLen[workPhase][workItem];
       ii = ssa_itemLoc[workPhase][workItem];
-
-      while (ct-- > 0) {
-        int len = list->numneigh[ii];
-        if (len > 0) {
-          if (useDPDE) ssa_update_dpde(ilist[ii], list->firstneigh[ii], len);
-          else ssa_update_dpd(ilist[ii], list->firstneigh[ii], len);
-        }
-        ii++;
-      }
+      if (useDPDE) ssa_update_dpde(ii, ct, workItem);
+      else ssa_update_dpd(ii, ct, workItem);
     }
   }
 
-  ii = inum;
   //Loop over all 13 outward directions (7 stages)
-  for (airnum = 1; airnum <=7; airnum++){
-    int ct = list->AIRct_ssa[airnum];
+  for (int workPhase = 0; workPhase < ssa_gphaseCt; ++workPhase) {
+    int workItemCt = ssa_gphaseLen[workPhase];
 
     // Communicate the updated velocities to all nodes
     comm->forward_comm_fix(this);
@@ -585,12 +622,11 @@ void FixShardlow::initial_integrate(int vflag)
       memset(&(atom->uMech[nlocal]), 0, sizeof(double)*nghost);
     }
 
-    // process neighbors in this AIR
-    while (ct-- > 0) {
-      int len = list->numneigh[ii];
-      if (useDPDE) ssa_update_dpde(ilist[ii], list->firstneigh[ii], len);
-      else ssa_update_dpd(ilist[ii], list->firstneigh[ii], len);
-      ii++;
+    for (int workItem = 0; workItem < workItemCt; ++workItem) {
+      int ct = ssa_gitemLen[workPhase][workItem];
+      ii = ssa_gitemLoc[workPhase][workItem];
+      if (useDPDE) ssa_update_dpde(ii, ct, workItem);
+      else ssa_update_dpd(ii, ct, workItem);
     }
 
     // Communicate the ghost deltas to the atom owners
@@ -699,6 +735,7 @@ double FixShardlow::memory_usage()
 {
   double bytes = 0.0;
   bytes += sizeof(double)*3*atom->nghost; // v_t0[]
+  bytes += sizeof(*rand_state)*maxRNG; // rand_state[]
   return bytes;
 }
 
diff --git a/src/USER-DPD/fix_shardlow.h b/src/USER-DPD/fix_shardlow.h
index e8e5f484a0c7b1540cce855680d4459ffeeaad6c..21f7569a233de12032311d80b0144ef7fa8acf1e 100644
--- a/src/USER-DPD/fix_shardlow.h
+++ b/src/USER-DPD/fix_shardlow.h
@@ -21,6 +21,8 @@ FixStyle(shardlow,FixShardlow)
 #define LMP_FIX_SHARDLOW_H
 
 #include "fix.h"
+#include "random_external_state.h"
+#include <math.h>
 
 namespace LAMMPS_NS {
 
@@ -52,12 +54,14 @@ class FixShardlow : public Fix {
   class PairDPDfdt *pairDPD;
   class PairDPDfdtEnergy *pairDPDE;
   double (*v_t0)[3];
+  int maxRNG;
 
  private:
+  random_external_state::es_RNG_t *rand_state;
   double dtsqrt; // = sqrt(update->dt);
 
-  void ssa_update_dpd(int, int *, int);  // Constant Temperature
-  void ssa_update_dpde(int, int *, int); // Constant Energy
+  void ssa_update_dpd(int, int, int);  // Constant Temperature
+  void ssa_update_dpde(int, int, int); // Constant Energy
 
 };
 
diff --git a/src/USER-DPD/npair_half_bin_newton_ssa.cpp b/src/USER-DPD/npair_half_bin_newton_ssa.cpp
index a6479d4c4f9c2d225aae248fe82f08c6a7c42fd9..ce50f7603adb0a757bbf9e03d0f31830e0272bd7 100644
--- a/src/USER-DPD/npair_half_bin_newton_ssa.cpp
+++ b/src/USER-DPD/npair_half_bin_newton_ssa.cpp
@@ -42,6 +42,10 @@ NPairHalfBinNewtonSSA::NPairHalfBinNewtonSSA(LAMMPS *lmp) : NPair(lmp)
   ssa_phaseLen = NULL;
   ssa_itemLoc = NULL;
   ssa_itemLen = NULL;
+  ssa_gphaseCt = 7;
+  memory->create(ssa_gphaseLen,ssa_gphaseCt,"NPairHalfBinNewtonSSA:ssa_gphaseLen");
+  memory->create(ssa_gitemLoc,ssa_gphaseCt,1,"NPairHalfBinNewtonSSA:ssa_gitemLoc");
+  memory->create(ssa_gitemLen,ssa_gphaseCt,1,"NPairHalfBinNewtonSSA:ssa_gitemLen");
 }
 
 /* ---------------------------------------------------------------------- */
@@ -54,6 +58,10 @@ NPairHalfBinNewtonSSA::~NPairHalfBinNewtonSSA()
   memory->destroy(ssa_phaseLen);
   memory->destroy(ssa_itemLoc);
   memory->destroy(ssa_itemLen);
+  ssa_gphaseCt = 0;
+  memory->destroy(ssa_gphaseLen);
+  memory->destroy(ssa_gitemLoc);
+  memory->destroy(ssa_gitemLen);
 }
 
 /* ----------------------------------------------------------------------
@@ -236,13 +244,14 @@ void NPairHalfBinNewtonSSA::build(NeighList *list)
 
   if (ssa_phaseCt != workPhase) error->one(FLERR,"ssa_phaseCt was wrong");
 
-  list->AIRct_ssa[0] = list->inum = inum;
+  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 = 1; airnum <= 7; airnum++) {
+  for (workPhase = 0; workPhase < ssa_gphaseCt; workPhase++) {
     int locAIRct = 0;
-    for (i = gairhead_ssa[airnum]; i >= 0; i = bins[i]) {
+    ssa_gitemLoc[workPhase][0] = inum + gnum; // record where workItem starts in ilist
+    for (i = gairhead_ssa[workPhase+1]; i >= 0; i = bins[i]) {
       n = 0;
       neighptr = ipage->vget();
 
@@ -305,7 +314,8 @@ void NPairHalfBinNewtonSSA::build(NeighList *list)
       if (ipage->status())
         error->one(FLERR,"Neighbor (ghost) list overflow, boost neigh_modify one");
     }
-    list->AIRct_ssa[airnum] = locAIRct;
+    ssa_gitemLen[workPhase][0] = locAIRct;
+    ssa_gphaseLen[workPhase] = 1;
   }
   list->gnum = gnum;
 }
diff --git a/src/USER-DPD/npair_half_bin_newton_ssa.h b/src/USER-DPD/npair_half_bin_newton_ssa.h
index ea292316ca740608346ccc8b2c8be5fa09b51370..584d87e3adf3aaeecf876665d88ebf92c8aacab3 100644
--- a/src/USER-DPD/npair_half_bin_newton_ssa.h
+++ b/src/USER-DPD/npair_half_bin_newton_ssa.h
@@ -33,13 +33,18 @@ class NPairHalfBinNewtonSSA : public NPair {
   int *ssa_phaseLen;
   int **ssa_itemLoc;
   int **ssa_itemLen;
+  int ssa_gphaseCt;
+  int *ssa_gphaseLen;
+  int **ssa_gitemLoc;
+  int **ssa_gitemLen;
+
+  int ssa_maxPhaseLen;
 
   NPairHalfBinNewtonSSA(class LAMMPS *);
   ~NPairHalfBinNewtonSSA();
   void build(class NeighList *);
  private:
   int ssa_maxPhaseCt;
-  int ssa_maxPhaseLen;
 };
 
 }
diff --git a/src/USER-DPD/pair_dpd_fdt.h b/src/USER-DPD/pair_dpd_fdt.h
index 5c20f2fc8f3b8fb02efd02a214f0e76c34e9e04d..84b09b0fa2bba261219f2c7d9f5a96905c0209c5 100644
--- a/src/USER-DPD/pair_dpd_fdt.h
+++ b/src/USER-DPD/pair_dpd_fdt.h
@@ -44,11 +44,11 @@ class PairDPDfdt : public Pair {
   double **sigma;
   double temperature;
 
+  int seed;
   class RanMars *random;
 
  protected:
   double cut_global;
-  int seed;
   bool splitFDT_flag;
   bool a0_is_zero;
 
diff --git a/src/USER-DPD/random_external_state.h b/src/USER-DPD/random_external_state.h
new file mode 100644
index 0000000000000000000000000000000000000000..d97d5a17ce7c5242a98263e18e8b05ddf3abb9cb
--- /dev/null
+++ b/src/USER-DPD/random_external_state.h
@@ -0,0 +1,179 @@
+/* -*- c++ -*- ----------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   http://lammps.sandia.gov, Sandia National Laboratories
+   Steve Plimpton, sjplimp@sandia.gov
+
+   Copyright (2003) Sandia Corporation.  Under the terms of Contract
+   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+   certain rights in this software.  This software is distributed under
+   the GNU General Public License.
+
+   See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+/*
+   This random_external_state.h file was derrived from the Kokkos
+   file algorithms/src/Kokkos_Random.hpp and adapted to work
+   without Kokkos installed, as well as being converted to a form
+   that has no internal state.  All RNG state information is kept
+   outside this "class", and is passed in by reference by the caller.
+ */
+/*
+//@HEADER
+// ************************************************************************
+//
+//                        Kokkos v. 2.0
+//              Copyright (2014) Sandia Corporation
+//
+// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
+// the U.S. Government retains certain rights in this software.
+//
+// Redistribution and use in source and binary forms, with or without
+// modification, are permitted provided that the following conditions are
+// met:
+//
+// 1. Redistributions of source code must retain the above copyright
+// notice, this list of conditions and the following disclaimer.
+//
+// 2. Redistributions in binary form must reproduce the above copyright
+// notice, this list of conditions and the following disclaimer in the
+// documentation and/or other materials provided with the distribution.
+//
+// 3. Neither the name of the Corporation nor the names of the
+// contributors may be used to endorse or promote products derived from
+// this software without specific prior written permission.
+//
+// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
+// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
+// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
+// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
+// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
+// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
+// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
+// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
+// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
+// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+//
+// Questions? Contact  H. Carter Edwards (hcedwar@sandia.gov)
+//
+// ************************************************************************
+//@HEADER
+*/
+
+#ifndef LMP_RANDOM_EXTERNALSTATE_H
+#define LMP_RANDOM_EXTERNALSTATE_H
+
+#include <math.h>
+#include "accelerator_kokkos.h"
+
+#ifdef LMP_KOKKOS
+#define RND_INLINE KOKKOS_INLINE_FUNCTION
+#else
+#define RND_INLINE inline
+#endif
+
+
+/// \file random_external_state.h
+/// \brief Pseudorandom number generators
+///
+/// These generators are based on Vigna, Sebastiano (2014). "An
+/// experimental exploration of Marsaglia's xorshift generators,
+/// scrambled."  See: http://arxiv.org/abs/1402.6246
+
+// A replacement for the Kokkos Random_XorShift64 class that uses
+// an external state variable, instead of a class member variable.
+namespace random_external_state {
+  typedef uint64_t es_RNG_t;
+
+  enum {MAX_URAND = 0xffffffffU};
+  enum {MAX_URAND64 = 0xffffffffffffffffULL-1};
+
+  RND_INLINE
+  uint32_t es_urand(es_RNG_t &state_) {
+    state_ ^= state_ >> 12;
+    state_ ^= state_ << 25;
+    state_ ^= state_ >> 27;
+
+    es_RNG_t tmp = state_ * 2685821657736338717ULL;
+    tmp = tmp>>16;
+    return static_cast<uint32_t>(tmp&MAX_URAND);
+  }
+
+  RND_INLINE
+  uint64_t es_urand64(es_RNG_t &state_) {
+    state_ ^= state_ >> 12;
+    state_ ^= state_ << 25;
+    state_ ^= state_ >> 27;
+    return (state_ * 2685821657736338717ULL) - 1;
+  }
+
+  RND_INLINE
+  int es_rand(es_RNG_t &state_) {
+    return static_cast<int>(es_urand(state_)/2);
+  }
+
+  RND_INLINE
+  double es_drand(es_RNG_t &state_) {
+    return 1.0 * es_urand64(state_)/MAX_URAND64;
+  }
+
+  //Marsaglia polar method for drawing a standard normal distributed random number
+  RND_INLINE
+  double es_normal(es_RNG_t &state_) {
+    double S, U;
+    do {
+      U = 2.0*es_drand(state_) - 1.0;
+      const double V = 2.0*es_drand(state_) - 1.0;
+      S = U*U+V*V;
+    } while ((S >= 1.0) || (S == 0.0));
+    return U*sqrt(-2.0*log(S)/S);
+  }
+
+  RND_INLINE
+  double es_normalPair(es_RNG_t &state_, double &second) {
+    double S, U, V;
+    do {
+      U = 2.0*es_drand(state_) - 1.0;
+      V = 2.0*es_drand(state_) - 1.0;
+      S = U*U+V*V;
+    } while ((S >= 1.0) || (S == 0.0));
+    const double fac = sqrt(-2.0*log(S)/S);
+    second = V*fac;
+    return U*fac;
+  }
+
+  // Use es_init() to init a serial RNG, that is then
+  // used to generate the initial state of your k parallel
+  // RNGs with k calls to genNextParallelState()
+  RND_INLINE
+  void es_init(es_RNG_t &serial_state, uint64_t seed) {
+    if(seed==0) seed = uint64_t(1318319);
+    serial_state = seed;
+    for(int i = 0; i < 17; i++) es_rand(serial_state);
+  }
+
+  // Call genNextParallelState() once for each RNG to generate
+  // the initial state for that RNG.
+  RND_INLINE
+  void es_genNextParallelState(es_RNG_t &serial_state, es_RNG_t &new_state) {
+    int n1 = es_rand(serial_state);
+    int n2 = es_rand(serial_state);
+    int n3 = es_rand(serial_state);
+    int n4 = es_rand(serial_state);
+    new_state = ((((static_cast<es_RNG_t>(n1)) & 0xffff)<<00) |
+                 (((static_cast<es_RNG_t>(n2)) & 0xffff)<<16) |
+                 (((static_cast<es_RNG_t>(n3)) & 0xffff)<<32) |
+                 (((static_cast<es_RNG_t>(n4)) & 0xffff)<<48));
+  }
+}
+
+#endif
+
+/* ERROR/WARNING messages:
+
+E: Invalid seed for Marsaglia random # generator
+
+The initial seed for this random number generator must be a positive
+integer less than or equal to 900 million.
+
+*/
diff --git a/src/neigh_list.cpp b/src/neigh_list.cpp
index 934b9f7d9b0409825efc3a51d269df88ba46a15b..93f4b13bf2c3359db2d0a22037e0218b960fe9c4 100644
--- a/src/neigh_list.cpp
+++ b/src/neigh_list.cpp
@@ -88,7 +88,6 @@ NeighList::NeighList(LAMMPS *lmp) : Pointers(lmp)
 
   // USER-DPD package
 
-  for (int i = 0; i < 8; i++) AIRct_ssa[i] = 0;
   np = NULL;
 }
 
diff --git a/src/neigh_list.h b/src/neigh_list.h
index d633ba839e0fc286856f54148a57f64f6559fc0c..755a1bf13415c4cdd7d428f3000c2039047ed081 100644
--- a/src/neigh_list.h
+++ b/src/neigh_list.h
@@ -92,7 +92,6 @@ class NeighList : protected Pointers {
 
   // USER-DPD package and Shardlow Splitting Algorithm (SSA) support
 
-  int AIRct_ssa[8]; // count of how many atoms in each AIR
   class NPair *np;           // ptr to NPair instance I depend on
 
   // methods