diff --git a/doc/src/fix_deform.rst b/doc/src/fix_deform.rst
index 829e4947eb..93f4136744 100644
--- a/doc/src/fix_deform.rst
+++ b/doc/src/fix_deform.rst
@@ -298,7 +298,8 @@ Note that it should return the "change" in box length, not the
 absolute box length.  This means it should evaluate to 0.0 when
 invoked on the initial timestep of the run following the definition of
 fix deform.  It should evaluate to a value > 0.0 to dilate the box at
-future times, or a value < 0.0 to compress the box.
+future times, or a value < 0.0 to compress the box. The exception
+would be if the run command uses the *pre no* option.
 
 The variable *name2* must also be an :doc:`equal-style variable
 <variable>` and should calculate the rate of box length change, in
diff --git a/src/RHEO/bond_rheo_shell.cpp b/src/RHEO/bond_rheo_shell.cpp
index d5b2dab1b0..b74a969cf2 100644
--- a/src/RHEO/bond_rheo_shell.cpp
+++ b/src/RHEO/bond_rheo_shell.cpp
@@ -147,7 +147,8 @@ void BondRHEOShell::store_data()
 
       // map to find index n
       j = atom->map(atom->bond_atom[i][m]);
-      if (j == -1) error->one(FLERR, "Atom missing in BPM bond");
+      if (j == -1) error->one(FLERR, Error::NOLASTLINE,
+                              "Atom {} missing in BPM bond", atom->bond_atom[i][m]);
 
       fix_bond_history->update_atom_value(i, m, 0, 0.0);
       fix_bond_history->update_atom_value(i, m, 1, 0.0);
@@ -217,12 +218,6 @@ void BondRHEOShell::compute(int eflag, int vflag)
       i2 = itmp;
     }
 
-    delx = x[i1][0] - x[i2][0];
-    dely = x[i1][1] - x[i2][1];
-    delz = x[i1][2] - x[i2][2];
-    rsq = delx * delx + dely * dely + delz * delz;
-    r = sqrt(rsq);
-
     // If bond hasn't been set - zero data
     if (t < EPSILON || std::isnan(t)) t = store_bond(n, i1, i2);
 
@@ -358,19 +353,20 @@ void BondRHEOShell::init_style()
   BondBPM::init_style();
 
   if (comm->ghost_velocity == 0)
-    error->all(FLERR, "Bond rheo/shell requires ghost atoms store velocity");
+    error->all(FLERR, Error::NOLASTLINE, "Bond rheo/shell requires ghost atoms store velocity");
 
   auto fixes = modify->get_fix_by_style("^rheo$");
-  if (fixes.size() == 0) error->all(FLERR, "Need to define fix rheo to use bond rheo/shell");
+  if (fixes.size() == 0)
+    error->all(FLERR, Error::NOLASTLINE, "Need to define fix rheo to use bond rheo/shell");
   class FixRHEO *fix_rheo = dynamic_cast<FixRHEO *>(fixes[0]);
 
   if (!fix_rheo->surface_flag)
-    error->all(FLERR, "Bond rheo/shell requires surface calculation in fix rheo");
+    error->all(FLERR, Error::NOLASTLINE, "Bond rheo/shell requires surface calculation in fix rheo");
   compute_surface = fix_rheo->compute_surface;
 
   fixes = modify->get_fix_by_style("^rheo/oxidation$");
   if (fixes.size() == 0)
-    error->all(FLERR, "Need to define fix rheo/oxidation to use bond rheo/shell");
+    error->all(FLERR, Error::NOLASTLINE, "Need to define fix rheo/oxidation to use bond rheo/shell");
   class FixRHEOOxidation *fix_rheo_oxidation = dynamic_cast<FixRHEOOxidation *>(fixes[0]);
 
   rsurf = fix_rheo_oxidation->rsurf;
@@ -389,14 +385,13 @@ void BondRHEOShell::settings(int narg, char **arg)
     if (strcmp(arg[iarg], "t/form") == 0) {
       if (iarg + 1 > narg) utils::missing_cmd_args(FLERR, "bond rheo/shell t/form", error);
       tform = utils::numeric(FLERR, arg[iarg + 1], false, lmp);
+      if (tform < 0.0)
+        error->all(FLERR, iarg + 1, "Illegal bond rheo/shell command, must specify positive formation time");
       i += 1;
     } else {
-      error->all(FLERR, "Illegal bond rheo/shell command, invalid argument {}", arg[iarg]);
+      error->all(FLERR, iarg, "Illegal bond rheo/shell command, invalid argument {}", arg[iarg]);
     }
   }
-
-  if (tform < 0.0)
-    error->all(FLERR, "Illegal bond rheo/shell command, must specify positive formation time");
 }
 
 /* ----------------------------------------------------------------------
@@ -472,7 +467,7 @@ int BondRHEOShell::pack_reverse_comm(int n, int first, double *buf)
   m = 0;
   last = first + n;
 
-  for (i = first; i < last; i++) { buf[m++] = dbond[i]; }
+  for (i = first; i < last; i++) { buf[m++] = ubuf(dbond[i]).d; }
   return m;
 }
 
@@ -485,7 +480,7 @@ void BondRHEOShell::unpack_reverse_comm(int n, int *list, double *buf)
   m = 0;
   for (i = 0; i < n; i++) {
     j = list[i];
-    dbond[j] += buf[m++];
+    dbond[j] += (int) ubuf(buf[m++]).i;
   }
 }
 
@@ -495,14 +490,20 @@ double BondRHEOShell::single(int type, double rsq, int i, int j, double &fforce)
 {
   if (type <= 0) return 0.0;
 
-  double r0, t;
+  double r0 = -1;
+  double t = -1;
+  const auto *tag = atom->tag;
   for (int n = 0; n < atom->num_bond[i]; n++) {
-    if (atom->bond_atom[i][n] == atom->tag[j]) {
+    if (atom->bond_atom[i][n] == tag[j]) {
       r0 = fix_bond_history->get_atom_value(i, n, 0);
       t = fix_bond_history->get_atom_value(i, n, 1);
     }
   }
 
+  if (r0 == -1)
+    error->one(FLERR, Error::NOLASTLINE,
+               "Could not find bond between atoms {} and {}", tag[i], tag[j]);
+
   svector[1] = t;
   if (t < tform) return 0.0;
 
diff --git a/src/RHEO/compute_rheo_kernel.cpp b/src/RHEO/compute_rheo_kernel.cpp
index e44e356133..16ecca00e9 100644
--- a/src/RHEO/compute_rheo_kernel.cpp
+++ b/src/RHEO/compute_rheo_kernel.cpp
@@ -251,13 +251,14 @@ double ComputeRHEOKernel::calc_w_quintic(double r)
   double w, tmp1, tmp2, tmp3, tmp1sq, tmp2sq, tmp3sq, s;
   s = r * 3.0 * cutinv;
 
-  if (s > 3.0) { w = 0.0; }
-
   if (s <= 3.0) {
     tmp3 = 3.0 - s;
     tmp3sq = tmp3 * tmp3;
     w = tmp3sq * tmp3sq * tmp3;
+  } else {
+    w = 0.0;
   }
+
   if (s <= 2.0) {
     tmp2 = 2.0 - s;
     tmp2sq = tmp2 * tmp2;
@@ -285,12 +286,14 @@ double ComputeRHEOKernel::calc_dw_scalar_quintic(double r)
 
   s = r * 3.0 * cutinv;
 
-  if (s > 3.0) { wp = 0.0; }
   if (s <= 3.0) {
     tmp3 = 3.0 - s;
     tmp3sq = tmp3 * tmp3;
     wp = -5.0 * tmp3sq * tmp3sq;
+  } else {
+    wp = 0.0;
   }
+
   if (s <= 2.0) {
     tmp2 = 2.0 - s;
     tmp2sq = tmp2 * tmp2;
@@ -855,7 +858,7 @@ int ComputeRHEOKernel::pack_forward_comm(int n, int *list, double *buf, int /*pb
   for (int i = 0; i < n; i++) {
     int j = list[i];
     if (comm_stage == 0) {
-      buf[m++] = coordination[j];
+      buf[m++] = ubuf(coordination[j]).d;
     } else {
       if (kernel_style == RK0) {
         buf[m++] = C0[j];
@@ -878,7 +881,7 @@ void ComputeRHEOKernel::unpack_forward_comm(int n, int first, double *buf)
 
   for (int i = first; i < last; i++) {
     if (comm_stage == 0) {
-      coordination[i] = buf[m++];
+      coordination[i] = (int) ubuf(buf[m++]).i;
     } else {
       if (kernel_style == RK0) {
         C0[i] = buf[m++];
diff --git a/src/RHEO/compute_rheo_vshift.cpp b/src/RHEO/compute_rheo_vshift.cpp
index 5fbc780e91..03dc001ad2 100644
--- a/src/RHEO/compute_rheo_vshift.cpp
+++ b/src/RHEO/compute_rheo_vshift.cpp
@@ -271,6 +271,8 @@ void ComputeRHEOVShift::correct_surfaces()
   int dim = domain->dimension;
 
   double nx, ny, nz, vx, vy, vz, dot;
+  nz = 0.0;
+  vz = 0.0;
   for (int i = 0; i < nlocal; i++) {
     if (mask[i] & groupbit) {
 
diff --git a/src/RHEO/fix_rheo.cpp b/src/RHEO/fix_rheo.cpp
index c2ff2d8b09..423847c715 100644
--- a/src/RHEO/fix_rheo.cpp
+++ b/src/RHEO/fix_rheo.cpp
@@ -110,7 +110,7 @@ FixRHEO::FixRHEO(LAMMPS *lmp, int narg, char **arg) :
     kernel_style = RK2;
   } else
     error->all(FLERR, "Unknown kernel style {} in fix rheo", arg[4]);
-  zmin_kernel = utils::numeric(FLERR, arg[5], false, lmp);
+  zmin_kernel = utils::inumeric(FLERR, arg[5], false, lmp);
 
   int iarg = 6;
   while (iarg < narg) {
@@ -342,8 +342,6 @@ void FixRHEO::initial_integrate(int /*vflag*/)
   double *rmass = atom->rmass;
   double **gradr = compute_grad->gradr;
   double **gradv = compute_grad->gradv;
-  double **vshift;
-  if (shift_flag) vshift = compute_vshift->vshift;
 
   int nlocal = atom->nlocal;
   int rmass_flag = atom->rmass_flag;
@@ -400,6 +398,7 @@ void FixRHEO::initial_integrate(int /*vflag*/)
 
   // Shifting atoms
   if (shift_flag) {
+    double **vshift = compute_vshift->vshift;
     for (i = 0; i < nlocal; i++) {
       if (status[i] & STATUS_NO_SHIFT) continue;
       if (status[i] & PHASECHECK) continue;
diff --git a/src/RHEO/fix_rheo_thermal.cpp b/src/RHEO/fix_rheo_thermal.cpp
index e9842e19e2..1d3eaf6f19 100644
--- a/src/RHEO/fix_rheo_thermal.cpp
+++ b/src/RHEO/fix_rheo_thermal.cpp
@@ -189,7 +189,7 @@ FixRHEOThermal::FixRHEOThermal(LAMMPS *lmp, int narg, char **arg) :
     } else if (strcmp(arg[iarg], "react") == 0) {
       if (iarg + 2 >= narg) utils::missing_cmd_args(FLERR, "fix rheo/thermal react", error);
       cut_bond = utils::numeric(FLERR, arg[iarg + 1], false, lmp);
-      btype = utils::numeric(FLERR, arg[iarg + 2], false, lmp);
+      btype = utils::inumeric(FLERR, arg[iarg + 2], false, lmp);
       comm_forward = 4;
       if (cut_bond <= 0.0) error->all(FLERR, "Illegal max bond length must be greater than zero");
       if ((btype < 1) || (btype > atom->nbondtypes))
diff --git a/src/RHEO/pair_rheo.cpp b/src/RHEO/pair_rheo.cpp
index db83b01539..af9802e07c 100644
--- a/src/RHEO/pair_rheo.cpp
+++ b/src/RHEO/pair_rheo.cpp
@@ -80,7 +80,7 @@ void PairRHEO::compute(int eflag, int vflag)
   int i, j, a, b, ii, jj, inum, jnum, itype, jtype;
   int pair_force_flag, pair_rho_flag, pair_avisc_flag;
   int fluidi, fluidj;
-  double xtmp, ytmp, ztmp, wp, Ti, Tj, dT, csq_ave, cs_ave;
+  double xtmp, ytmp, ztmp, wp, Ti, Tj, dT, cs_ave;
   double rhoi, rhoj, rho0i, rho0j, voli, volj, Pi, Pj, etai, etaj, kappai, kappaj, csqi, csqj;
   double eta_ave, kappa_ave, dT_prefactor;
   double mu, q, fp_prefactor, drho_damp, fmag, psi_ij, Fij;
@@ -189,11 +189,10 @@ void PairRHEO::compute(int eflag, int vflag)
         if (!variable_csq) {
           cs_ave = 0.5 * (cs[itype] + cs[jtype]);
         } else {
-          csqi = fix_pressure->calc_csq(rhoi, i);
-          csqj = fix_pressure->calc_csq(rhoj, j);
+          csqi = fix_pressure->calc_csq(rho[i], i);
+          csqj = fix_pressure->calc_csq(rho[j], j);
           cs_ave = 0.5 * (sqrt(csqi) + sqrt(csqj));
         }
-        csq_ave = cs_ave * cs_ave;
 
         pair_rho_flag = 0;
         pair_force_flag = 0;
@@ -231,7 +230,7 @@ void PairRHEO::compute(int eflag, int vflag)
             Pj = fix_pressure->calc_pressure(rhoj, j);
 
             if ((chi[j] > 0.9) && (r < (cutk * 0.5)))
-              fmag = (chi[j] - 0.9) * (cutk * 0.5 - r) * rho0j * csq_ave * cutk * rinv;
+              fmag = (chi[j] - 0.9) * (cutk * 0.5 - r) * rho0j * cs_ave * cs_ave * cutk * rinv;
 
           } else if ((!fluidi) && fluidj) {
             compute_interface->correct_v(vi, vj, i, j);
@@ -239,7 +238,7 @@ void PairRHEO::compute(int eflag, int vflag)
             Pi = fix_pressure->calc_pressure(rhoi, i);
 
             if (chi[i] > 0.9 && r < (cutk * 0.5))
-              fmag = (chi[i] - 0.9) * (cutk * 0.5 - r) * rho0i * csq_ave * cutk * rinv;
+              fmag = (chi[i] - 0.9) * (cutk * 0.5 - r) * rho0i * cs_ave * cs_ave * cutk * rinv;
 
           } else if ((!fluidi) && (!fluidj)) {
             rhoi = rho0i;
@@ -251,7 +250,6 @@ void PairRHEO::compute(int eflag, int vflag)
             csqi = fix_pressure->calc_csq(rhoi, i);
             csqj = fix_pressure->calc_csq(rhoj, j);
             cs_ave = 0.5 * (sqrt(csqi) + sqrt(csqj));
-            csq_ave = cs_ave * cs_ave;
           }
         }
 
diff --git a/src/set.cpp b/src/set.cpp
index 8e7ef9203f..73bc7ccbfe 100644
--- a/src/set.cpp
+++ b/src/set.cpp
@@ -48,7 +48,7 @@ enum{ATOM_SELECT,MOL_SELECT,TYPE_SELECT,GROUP_SELECT,REGION_SELECT};
 
 enum{ANGLE,ANGMOM,APIP_LAMBDA,BOND,CC,CHARGE,DENSITY,DIAMETER,DIHEDRAL,DIPOLE,
   DIPOLE_RANDOM,DPD_THETA,EDPD_CV,EDPD_TEMP,EPSILON,IMAGE,IMPROPER,LENGTH,
-  MASS,MOLECULE,OMEGA,QUAT,QUAT_RANDOM,RADIUS_ELECTRON,SHAPE,
+  MASS,MOLECULE,OMEGA,QUAT,QUAT_RANDOM,RADIUS_ELECTRON,RHEO_STATUS,SHAPE,
   SMD_CONTACT_RADIUS,SMD_MASS_DENSITY,SPH_CV,SPH_E,SPH_RHO,
   SPIN_ATOM,SPIN_ATOM_RANDOM,SPIN_ELECTRON,TEMPERATURE,THETA,THETA_RANDOM,
   TRI,TYPE,TYPE_FRACTION,TYPE_RATIO,TYPE_SUBSET,VOLUME,VX,VY,VZ,X,Y,Z,
@@ -296,6 +296,14 @@ void Set::process_args(int caller_flag, int narg, char **arg)
       action->keyword = RADIUS_ELECTRON;
       process_radius_election(iarg,narg,arg,action);
       invoke_choice[naction++] = &Set::invoke_radius_election;
+    } else if (strcmp(arg[iarg],"rheo/rho") == 0) {
+      action->keyword = SPH_RHO;
+      process_sph_rho(iarg,narg,arg,action);
+      invoke_choice[naction++] = &Set::invoke_sph_rho;
+    } else if (strcmp(arg[iarg],"rheo/status") == 0) {
+      action->keyword = RHEO_STATUS;
+      process_rheo_status(iarg,narg,arg,action);
+      invoke_choice[naction++] = &Set::invoke_rheo_status;
     } else if (strcmp(arg[iarg],"shape") == 0) {
       action->keyword = SHAPE;
       process_shape(iarg,narg,arg,action);
@@ -2029,6 +2037,44 @@ void Set::invoke_radius_election(Action *action)
 
 /* ---------------------------------------------------------------------- */
 
+void Set::process_rheo_status(int &iarg, int narg, char **arg, Action *action)
+{
+  if (!atom->rheo_status_flag)
+    error->all(FLERR,"Cannot set attribute {} for atom style {}", arg[iarg], atom->get_style());
+  if (iarg+2 > narg) utils::missing_cmd_args(FLERR, "set rheo/status", error);
+
+  if (utils::strmatch(arg[iarg+1],"^v_")) varparse(arg[iarg+1],1,action);
+  else {
+    action->ivalue1 = utils::inumeric(FLERR,arg[iarg+1],false,lmp);
+    if (action->ivalue1 != 0 && action->ivalue1 != 1)
+      error->one(FLERR,"Invalid rheo/status {} in set command", action->ivalue1);
+  }
+
+  iarg += 2;
+}
+
+void Set::invoke_rheo_status(Action *action)
+{
+  int nlocal = atom->nlocal;
+  int *status = atom->rheo_status;
+
+  int varflag = action->varflag;
+  int rheo_status;
+  if (!action->varflag1) rheo_status = action->ivalue1;
+
+  for (int i = 0; i < nlocal; i++) {
+    if (!select[i]) continue;
+    if (varflag) {
+      rheo_status = static_cast<int> (vec1[i]);
+      if (rheo_status != 0 && rheo_status != 1)
+        error->one(FLERR,"Invalid rheo/status in set command");
+    }
+    status[i] = rheo_status;
+  }
+}
+
+/* ---------------------------------------------------------------------- */
+
 void Set::process_shape(int &iarg, int narg, char **arg, Action *action)
 {
   if (!atom->ellipsoid_flag)
diff --git a/src/set.h b/src/set.h
index 5e410bdd8c..856ffd12ac 100644
--- a/src/set.h
+++ b/src/set.h
@@ -111,6 +111,7 @@ class Set : public Command {
   void process_quat(int &, int, char **, Action *);
   void process_quat_random(int &, int, char **, Action *);
   void process_radius_election(int &, int, char **, Action *);
+  void process_rheo_status(int &, int, char **, Action *);
   void process_shape(int &, int, char **, Action *);
   void process_smd_contact_radius(int &, int, char **, Action *);
   void process_smd_mass_density(int &, int, char **, Action *);
@@ -164,6 +165,7 @@ class Set : public Command {
   void invoke_quat(Action *);
   void invoke_quat_random(Action *);
   void invoke_radius_election(Action *);
+  void invoke_rheo_status(Action *);
   void invoke_shape(Action *);
   void invoke_smd_contact_radius(Action *);
   void invoke_smd_mass_density(Action *);