diff --git a/src/atom.cpp b/src/atom.cpp
index ac88c4308b..e96e5e2050 100644
--- a/src/atom.cpp
+++ b/src/atom.cpp
@@ -93,6 +93,12 @@ Atom::Atom(LAMMPS *lmp) : Pointers(lmp)
   rho = drho = e = de = cv = NULL;
   vest = NULL;
 
+  // USER-DPD
+
+  uCond = uMech = uChem = uCG = uCGnew = NULL;
+  duCond = duMech = duChem = NULL;
+  dpdTheta = NULL;
+
   // USER-SMD
 
   contact_radius = NULL;
@@ -157,6 +163,7 @@ Atom::Atom(LAMMPS *lmp) : Pointers(lmp)
   cs_flag = csforce_flag = vforce_flag = etag_flag = 0;
 
   rho_flag = e_flag = cv_flag = vest_flag = 0;
+  dpd_flag = 0;
 
   // USER-SMD
 
@@ -273,6 +280,16 @@ Atom::~Atom()
   memory->destroy(eff_plastic_strain_rate);
   memory->destroy(damage);
 
+  memory->destroy(dpdTheta);
+  memory->destroy(uCond);
+  memory->destroy(uMech);
+  memory->destroy(uChem);
+  memory->destroy(uCG);
+  memory->destroy(uCGnew);
+  memory->destroy(duCond);
+  memory->destroy(duMech);
+  memory->destroy(duChem);
+
   memory->destroy(nspecial);
   memory->destroy(special);
 
@@ -2105,6 +2122,8 @@ void *Atom::extract(char *name)
     return (void *) eff_plastic_strain_rate;
   if (strcmp(name, "damage") == 0) return (void *) damage;
 
+  if (strcmp(name,"dpdTheta") == 0) return (void *) dpdTheta;
+
   return NULL;
 }
 
diff --git a/src/atom.h b/src/atom.h
index c0b69c9520..39611adf9e 100644
--- a/src/atom.h
+++ b/src/atom.h
@@ -85,6 +85,12 @@ class Atom : protected Pointers {
   double *eff_plastic_strain_rate;
   double *damage;
 
+  // USER-DPD package
+
+  double *uCond, *uMech, *uChem, *uCGnew, *uCG;
+  double *duCond, *duMech, *duChem;
+  double *dpdTheta;
+
   // molecular info
 
   int **nspecial;               // 0,1,2 = cummulative # of 1-2,1-3,1-4 neighs
@@ -133,6 +139,7 @@ class Atom : protected Pointers {
   int vfrac_flag,spin_flag,eradius_flag,ervel_flag,erforce_flag;
   int cs_flag,csforce_flag,vforce_flag,ervelforce_flag,etag_flag;
   int rho_flag,e_flag,cv_flag,vest_flag;
+  int dpd_flag;
 
   // USER-SMD package
 
diff --git a/src/library.cpp b/src/library.cpp
index 5801cedc5b..44734f95fc 100644
--- a/src/library.cpp
+++ b/src/library.cpp
@@ -150,6 +150,8 @@ void *lammps_extract_global(void *ptr, char *name)
   if (strcmp(name,"yz") == 0) return (void *) &lmp->domain->yz;
   if (strcmp(name,"natoms") == 0) return (void *) &lmp->atom->natoms;
   if (strcmp(name,"nlocal") == 0) return (void *) &lmp->atom->nlocal;
+  if (strcmp(name,"ntimestep") == 0) return (void *) &lmp->update->ntimestep;
+
   return NULL;
 }
 
diff --git a/src/neigh_request.cpp b/src/neigh_request.cpp
index 7c9cd448c7..9a7152a1bf 100644
--- a/src/neigh_request.cpp
+++ b/src/neigh_request.cpp
@@ -48,6 +48,7 @@ NeighRequest::NeighRequest(LAMMPS *lmp) : Pointers(lmp)
   // default is no CUDA neighbor list build
   // default is no multi-threaded neighbor list build
   // default is no Kokkos neighbor list build
+  // default is no Shardlow Splitting Algorithm (SSA) neighbor list build
 
   occasional = 0;
   newton = 0;
@@ -58,6 +59,7 @@ NeighRequest::NeighRequest(LAMMPS *lmp) : Pointers(lmp)
   omp = 0;
   intel = 0;
   kokkos_host = kokkos_device = 0;
+  ssa = 0;
 
   // default is no copy or skip
 
@@ -133,6 +135,7 @@ int NeighRequest::identical(NeighRequest *other)
   if (cudable != other->cudable) same = 0;
   if (omp != other->omp) same = 0;
   if (intel != other->intel) same = 0;
+  if (ssa != other->ssa) same = 0;
 
   if (copy != other->copy_original) same = 0;
   if (same_skip(other) == 0) same = 0;
@@ -163,6 +166,7 @@ int NeighRequest::same_kind(NeighRequest *other)
   if (cudable != other->cudable) same = 0;
   if (omp != other->omp) same = 0;
   if (intel != other->intel) same = 0;
+  if (ssa != other->ssa) same = 0;
 
   return same;
 }
@@ -214,4 +218,5 @@ void NeighRequest::copy_request(NeighRequest *other)
   cudable = other->cudable;
   omp = other->omp;
   intel = other->intel;
+  ssa = other->ssa;
 }
diff --git a/src/neigh_request.h b/src/neigh_request.h
index e835230638..047f3d0bd6 100644
--- a/src/neigh_request.h
+++ b/src/neigh_request.h
@@ -90,6 +90,10 @@ class NeighRequest : protected Pointers {
   int kokkos_host;
   int kokkos_device;
 
+  // 1 if using Shardlow Splitting Algorithm (SSA) neighbor list build
+  
+  int ssa;
+  
   // set by neighbor and pair_hybrid after all requests are made
   // these settings do not change kind value
 
diff --git a/src/set.cpp b/src/set.cpp
index 182fe528a3..a0bd6c2b65 100644
--- a/src/set.cpp
+++ b/src/set.cpp
@@ -44,7 +44,8 @@ enum{ATOM_SELECT,MOL_SELECT,TYPE_SELECT,GROUP_SELECT,REGION_SELECT};
 enum{TYPE,TYPE_FRACTION,MOLECULE,X,Y,Z,CHARGE,MASS,SHAPE,LENGTH,TRI,
      DIPOLE,DIPOLE_RANDOM,QUAT,QUAT_RANDOM,THETA,THETA_RANDOM,ANGMOM,OMEGA,
      DIAMETER,DENSITY,VOLUME,IMAGE,BOND,ANGLE,DIHEDRAL,IMPROPER,
-     MESO_E,MESO_CV,MESO_RHO,SMD_MASS_DENSITY,SMD_CONTACT_RADIUS,INAME,DNAME};
+     MESO_E,MESO_CV,MESO_RHO,SMD_MASS_DENSITY,SMD_CONTACT_RADIUS,DPDTHETA,
+     INAME,DNAME};
 
 #define BIG INT_MAX
 
@@ -381,51 +382,61 @@ void Set::command(int narg, char **arg)
       topology(IMPROPER);
       iarg += 2;
 
-    } else if (strcmp(arg[iarg],"meso_e") == 0) {
+    } else if (strcmp(arg[iarg],"meso/e") == 0) {
       if (iarg+2 > narg) error->all(FLERR,"Illegal set command");
       if (strstr(arg[iarg+1],"v_") == arg[iarg+1]) varparse(arg[iarg+1],1);
       else dvalue = force->numeric(FLERR,arg[iarg+1]);
       if (!atom->e_flag)
-        error->all(FLERR,"Cannot set this attribute for this atom style");
+        error->all(FLERR,"Cannot set meso/e for this atom style");
       set(MESO_E);
       iarg += 2;
 
-    } else if (strcmp(arg[iarg],"meso_cv") == 0) {
+    } else if (strcmp(arg[iarg],"meso/cv") == 0) {
       if (iarg+2 > narg) error->all(FLERR,"Illegal set command");
       if (strstr(arg[iarg+1],"v_") == arg[iarg+1]) varparse(arg[iarg+1],1);
       else dvalue = force->numeric(FLERR,arg[iarg+1]);
       if (!atom->cv_flag)
-            error->all(FLERR,"Cannot set this attribute for this atom style");
+            error->all(FLERR,"Cannot set meso/cv for this atom style");
       set(MESO_CV);
       iarg += 2;
 
-    } else if (strcmp(arg[iarg],"meso_rho") == 0) {
+    } else if (strcmp(arg[iarg],"meso/rho") == 0) {
       if (iarg+2 > narg) error->all(FLERR,"Illegal set command");
       if (strstr(arg[iarg+1],"v_") == arg[iarg+1]) varparse(arg[iarg+1],1);
       else dvalue = force->numeric(FLERR,arg[iarg+1]);
       if (!atom->rho_flag)
-        error->all(FLERR,"Cannot set meso_rho for this atom style");
+        error->all(FLERR,"Cannot set meso/rho for this atom style");
       set(MESO_RHO);
       iarg += 2;
 
-    } else if (strcmp(arg[iarg],"smd_mass_density") == 0) {
+    } else if (strcmp(arg[iarg],"smd/mass/density") == 0) {
           if (iarg+2 > narg) error->all(FLERR,"Illegal set command");
           if (strstr(arg[iarg+1],"v_") == arg[iarg+1]) varparse(arg[iarg+1],1);
           else dvalue = force->numeric(FLERR,arg[iarg+1]);
           if (!atom->smd_flag)
-            error->all(FLERR,"Cannot set smd_mass_density for this atom style");
+            error->all(FLERR,"Cannot set smd/mass/density for this atom style");
           set(SMD_MASS_DENSITY);
           iarg += 2;
 
-    } else if (strcmp(arg[iarg],"smd_contact_radius") == 0) {
+    } else if (strcmp(arg[iarg],"smd/contact/radius") == 0) {
           if (iarg+2 > narg) error->all(FLERR,"Illegal set command");
           if (strstr(arg[iarg+1],"v_") == arg[iarg+1]) varparse(arg[iarg+1],1);
           else dvalue = force->numeric(FLERR,arg[iarg+1]);
           if (!atom->smd_flag)
-        	  error->all(FLERR,"Cannot set smd_contact_radius for this atom style");
+            error->all(FLERR,"Cannot set smd/contact/radius "
+                       "for this atom style");
           set(SMD_CONTACT_RADIUS);
           iarg += 2;
 
+    } else if (strcmp(arg[iarg],"dpd/theta") == 0) {
+      if (iarg+2 > narg) error->all(FLERR,"Illegal set command");
+      if (strstr(arg[iarg+1],"v_") == arg[iarg+1]) varparse(arg[iarg+1],1);
+      else dvalue = force->numeric(FLERR,arg[iarg+1]);
+      if (!atom->dpd_flag)
+        error->all(FLERR,"Cannot set dpd/theta for this atom style");
+      set(DPDTHETA);
+      iarg += 2;
+
     } else if (strstr(arg[iarg],"i_") == arg[iarg]) {
       if (iarg+2 > narg) error->all(FLERR,"Illegal set command");
       if (strstr(arg[iarg+1],"v_") == arg[iarg+1]) varparse(arg[iarg+1],1);
@@ -621,6 +632,7 @@ void Set::set(int keyword)
       atom->rmass[i] = atom->vfrac[i] * dvalue;
     }
     else if (keyword == SMD_CONTACT_RADIUS) atom->contact_radius[i] = dvalue;
+    else if (keyword == DPDTHETA) atom->dpdTheta[i] = dvalue;
 
     // set shape of ellipsoidal particle