diff --git a/src/RHEO/atom_vec_rheo.cpp b/src/RHEO/atom_vec_rheo.cpp
index a8496a18e9..de2e2f77ad 100644
--- a/src/RHEO/atom_vec_rheo.cpp
+++ b/src/RHEO/atom_vec_rheo.cpp
@@ -120,7 +120,7 @@ void AtomVecRHEO::pack_property_atom(int index, double *buf, int nvalues, int gr
         buf[n] = 0.0;
       n += nvalues;
     }
-  } if else (index == 1) {
+  } else if (index == 1) {
     for (int i = 0; i < nlocal; i++) {
       if (mask[i] & groupbit)
         buf[n] = rho[i];
diff --git a/src/RHEO/compute_rheo_grad.cpp b/src/RHEO/compute_rheo_grad.cpp
index add3ed712d..606cec7dfc 100644
--- a/src/RHEO/compute_rheo_grad.cpp
+++ b/src/RHEO/compute_rheo_grad.cpp
@@ -21,11 +21,12 @@
 #include "atom.h"
 #include "comm.h"
 #include "compute_rheo_kernel.h"
-#include "compute_rheo_solids.h"
+#include "compute_rheo_interface.h"
 #include "domain.h"
 #include "error.h"
 #include "fix_rheo.h"
 #include "force.h"
+#include "memory.h"
 #include "neighbor.h"
 #include "neigh_list.h"
 #include "update.h"
@@ -33,6 +34,8 @@
 #include <cmath>
 
 using namespace LAMMPS_NS;
+using namespace RHEO_NS;
+
 enum{COMMGRAD, COMMFIELD};
 
 /* ---------------------------------------------------------------------- */
@@ -112,10 +115,13 @@ void ComputeRHEOGrad::init()
   compute_kernel = fix_rheo->compute_kernel;
   compute_interface = fix_rheo->compute_interface;
 
+  int tmp1, tmp2;
+  index_visc = atom->find_custom("rheo_viscosity", tmp1, tmp2);
+
   // Create coordination array if it doesn't already exist
   // Create a custom atom property so it works with compute property/atom
   // Do not create grow callback as there's no reason to copy/exchange data
-  // Manually grow if nmax_old exceeded
+  // Manually grow if nmax_store exceeded
 
   int index;
   int dim = domain->dimension;
@@ -139,7 +145,7 @@ void ComputeRHEOGrad::init()
     gradn = atom->darray[index];
   }
 
-  nmax_old = 0;
+  nmax_store = 0;
   grow_arrays(atom->nmax);
 }
 
@@ -158,7 +164,7 @@ void ComputeRHEOGrad::compute_peratom()
   double xtmp, ytmp, ztmp, delx, dely, delz;
   double rsq, imass, jmass;
   double rhoi, rhoj, Voli, Volj, drho, dT, deta;
-  double vij[3];
+  double vi[3], vj[3], vij[3];
   double wp, *dWij, *dWji;
 
   int inum, *ilist, *numneigh, **firstneigh;
@@ -169,26 +175,22 @@ void ComputeRHEOGrad::compute_peratom()
   double **v = atom->v;
   double *rho = atom->rho;
   double *temperature = atom->temperature;
+  double *viscosity = atom->dvector[index_visc];
   int *status = atom->status;
   int *type = atom->type;
   double *mass = atom->mass;
   int newton = force->newton;
   int dim = domain->dimension;
 
-  int tmp1, tmp2;
-  int index_visc = atom->find_custom("rheo_viscosity", tmp1, tmp2);
-  if (index_visc == -1) error->all(FLERR, "Cannot find rheo viscosity");
-  double *viscosity = atom->dvector[index_visc];
-
   inum = list->inum;
   ilist = list->ilist;
   numneigh = list->numneigh;
   firstneigh = list->firstneigh;
 
   // initialize arrays
-  if (nmax > nmax_old) grow_arrays(nmax);
+  if (atom->nmax > nmax_store) grow_arrays(atom->nmax);
 
-  for (i = 0; i < nmax; i++) {
+  for (i = 0; i < nmax_store; i++) {
     if (velocity_flag) {
       for (k = 0; k < dim * dim; k++)
         gradv[i][k] = 0.0;
@@ -212,6 +214,9 @@ void ComputeRHEOGrad::compute_peratom()
     xtmp = x[i][0];
     ytmp = x[i][1];
     ztmp = x[i][2];
+    vi[0] = v[i][0];
+    vi[1] = v[i][1];
+    vi[2] = v[i][2];
     itype = type[i];
     jlist = firstneigh[i];
     jnum = numneigh[i];
@@ -230,14 +235,18 @@ void ComputeRHEOGrad::compute_peratom()
         rhoi = rho[i];
         rhoj = rho[j];
 
+        vj[0] = v[j][0];
+        vj[1] = v[j][1];
+        vj[2] = v[j][2];
+
         // Add corrections for walls
-        if ((status[i] & FixRHEO::STATUS_FLUID) && !(status[j] & FixRHEO::STATUS_FLUID)) {
-          compute_interface->correct_v(v[i], v[j], vi, i, j);
-          rhoj = compute_interface->correct_rho(j,i);
-        } else if (!(status[i] & FixRHEO::STATUS_FLUID) && (status[j] & FixRHEO::STATUS_FLUID)) {
-          compute_interface->correct_v(v[j], v[i], vj, j, i);
-          rhoi = compute_interface->correct_rho(i,j);
-        } else if (!(status[i] & FixRHEO::STATUS_FLUID) && !(status[j] & FixRHEO::STATUS_FLUID)) {
+        if ((status[i] & STATUS_FLUID) && !(status[j] & STATUS_FLUID)) {
+          compute_interface->correct_v(vi, vj, i, j);
+          rhoj = compute_interface->correct_rho(j, i);
+        } else if (!(status[i] & STATUS_FLUID) && (status[j] & STATUS_FLUID)) {
+          compute_interface->correct_v(vj, vi, j, i);
+          rhoi = compute_interface->correct_rho(i, j);
+        } else if (!(status[i] & STATUS_FLUID) && !(status[j] & STATUS_FLUID)) {
           rhoi = rho0;
           rhoj = rho0;
         }
@@ -324,7 +333,6 @@ int ComputeRHEOGrad::pack_forward_comm(int n, int *list, double *buf,
   int i,j,k,m;
   double *rho = atom->rho;
   double *temperature = atom->temperature;
-  double *eta = atom->viscosity;
   double **v = atom->v;
   int dim = domain->dimension;
 
@@ -371,9 +379,9 @@ int ComputeRHEOGrad::pack_forward_comm(int n, int *list, double *buf,
 void ComputeRHEOGrad::unpack_forward_comm(int n, int first, double *buf)
 {
   int i, k, m, last;
-  double * rho = atom->rho;
-  double * temperature = atom->temperature;
-  double ** v = atom->v;
+  double *rho = atom->rho;
+  double *temperature = atom->temperature;
+  double **v = atom->v;
   int dim = domain->dimension;
 
   m = 0;
@@ -483,25 +491,27 @@ void ComputeRHEOGrad::grow_arrays(int nmax)
 
   if (eta_flag)
     memory->grow(gradn, nmax, dim, "atom:rheo_grad_eta");
-  nmax_old = nmax;
+  nmax_store = nmax;
 }
 
 /* ---------------------------------------------------------------------- */
 
-double ComputeRHEOKernel::memory_usage()
+double ComputeRHEOGrad::memory_usage()
 {
   double bytes = 0.0;
+  int dim = domain->dimension;
+
   if (velocity_flag)
-    bytes = (size_t) nmax_old * dim * dim * sizeof(double);
+    bytes = (size_t) nmax_store * dim * dim * sizeof(double);
 
   if (rho_flag)
-    bytes = (size_t) nmax_old * dim * sizeof(double);
+    bytes = (size_t) nmax_store * dim * sizeof(double);
 
   if (temperature_flag)
-    bytes = (size_t) nmax_old * dim * sizeof(double);
+    bytes = (size_t) nmax_store * dim * sizeof(double);
 
   if (eta_flag)
-    bytes = (size_t) nmax_old * dim * sizeof(double);
+    bytes = (size_t) nmax_store * dim * sizeof(double);
 
   return bytes;
 }
diff --git a/src/RHEO/compute_rheo_grad.h b/src/RHEO/compute_rheo_grad.h
index 4c2461c73f..8c7962a978 100644
--- a/src/RHEO/compute_rheo_grad.h
+++ b/src/RHEO/compute_rheo_grad.h
@@ -45,7 +45,8 @@ class ComputeRHEOGrad : public Compute {
   class FixRHEO *fix_rheo;
 
  private:
-  int comm_stage, ncomm_grad, ncomm_field, nmax_old;
+  int comm_stage, ncomm_grad, ncomm_field, nmax_store;
+  int index_visc;
   double cut, cutsq, rho0;
   class NeighList *list;
 
@@ -53,6 +54,8 @@ class ComputeRHEOGrad : public Compute {
   class ComputeRHEOInterface *compute_interface;
 
   int velocity_flag, temperature_flag, rho_flag, eta_flag;
+
+  void grow_arrays(int);
 };
 
 }    // namespace LAMMPS_NS
diff --git a/src/RHEO/compute_rheo_interface.cpp b/src/RHEO/compute_rheo_interface.cpp
index 059f4057f4..72dcdc17d7 100644
--- a/src/RHEO/compute_rheo_interface.cpp
+++ b/src/RHEO/compute_rheo_interface.cpp
@@ -34,18 +34,19 @@
 #include <cmath>
 
 using namespace LAMMPS_NS;
+using namespace RHEO_NS;
 
-#define EPSILON 1e-1
+static constexpr double EPSILON = 1e-1;
 
 /* ---------------------------------------------------------------------- */
 
 ComputeRHEOInterface::ComputeRHEOInterface(LAMMPS *lmp, int narg, char **arg) :
-  Compute(lmp, narg, arg), fix_rheo(nullptr), compute_kernel(nullptr), fx_m_norm(nullptr),
-  norm(nullptr), normwf(nullptr), chi(nullptr), f_pressure(nullptr), id_fix_pa(nullptr)
+  Compute(lmp, narg, arg), fix_rheo(nullptr), compute_kernel(nullptr), fp_store(nullptr),
+  norm(nullptr), normwf(nullptr), chi(nullptr), id_fix_pa(nullptr)
 {
   if (narg != 3) error->all(FLERR,"Illegal compute rheo/interface command");
 
-  nmax = 0;
+  nmax_store = 0;
 
   comm_forward = 3;
   comm_reverse = 4;
@@ -74,15 +75,15 @@ void ComputeRHEOInterface::init()
   compute_kernel = fix_rheo->compute_kernel;
   rho0 = fix_rheo->rho0;
   cut = fix_rheo->cut;
-  cs = fix_rheo->cs;
-  cs_inv = 1.0 / cs;
+  csq = fix_rheo->csq;
+  csq_inv = 1.0 / csq;
   cutsq = cut * cut;
   wall_max = sqrt(3.0) / 12.0 * cut;
 
   // Create chi array if it doesn't already exist
   // Create a custom atom property so it works with compute property/atom
   // Do not create grow callback as there's no reason to copy/exchange data
-  // Manually grow if nmax_old exceeded
+  // Manually grow if nmax_store exceeded
 
   int tmp1, tmp2;
   int nmax = atom->nmax;
@@ -93,7 +94,7 @@ void ComputeRHEOInterface::init()
     memory->destroy(normwf);
     memory->create(norm, nmax, "rheo/interface:norm");
     memory->create(normwf, nmax, "rheo/interface:normwf");
-    nmax_old = nmax;
+    nmax_store = nmax;
   }
   chi = atom->dvector[index];
 
@@ -104,13 +105,13 @@ void ComputeRHEOInterface::init()
   index = atom->find_custom("fp_store", tmp1, tmp2);
   if (index == -1) {
     id_fix_pa = utils::strdup(id + std::string("_fix_property_atom"));
-    modify->add_fix(fmt::format("{} all property/atom d2_fp_store 3", id_fix_pa)));
+    modify->add_fix(fmt::format("{} all property/atom d2_fp_store 3", id_fix_pa));
     index = atom->find_custom("fp_store", tmp1, tmp2);
   }
   fp_store = atom->darray[index];
 
   // need an occasional half neighbor list
-  neighbor->add_request(this, NeighConst::REQ_HALF);
+  neighbor->add_request(this, NeighConst::REQ_DEFAULT);
 }
 
 /* ---------------------------------------------------------------------- */
@@ -142,17 +143,17 @@ void ComputeRHEOInterface::compute_peratom()
   numneigh = list->numneigh;
   firstneigh = list->firstneigh;
 
-  if (atom->nmax > nmax_old) {
-    nmax_old = atom->nmax;
+  if (atom->nmax > nmax_store) {
+    nmax_store = atom->nmax;
     memory->destroy(norm);
     memory->destroy(normwf);
-    memory->create(norm, nmax_old, "rheo/interface:norm");
-    memory->create(normwf, nmax_old, "rheo/interface:normwf");
-    memory->grow(chi, nmax_old, "rheo/interface:chi");
+    memory->create(norm, nmax_store, "rheo/interface:norm");
+    memory->create(normwf, nmax_store, "rheo/interface:normwf");
+    memory->grow(chi, nmax_store, "rheo/interface:chi");
   }
 
   for (i = 0; i < nall; i++) {
-    if (!(status[i] & FixRHEO::STATUS_FLUID)) rho[i] = 0.0;
+    if (!(status[i] & STATUS_FLUID)) rho[i] = 0.0;
     normwf[i] = 0.0;
     norm[i] = 0.0;
     chi[i] = 0.0;
@@ -164,7 +165,7 @@ void ComputeRHEOInterface::compute_peratom()
     ytmp = x[i][1];
     ztmp = x[i][2];
     itype = type[i];
-    fluidi = status[i] & FixRHEO::STATUS_FLUID;
+    fluidi = status[i] & STATUS_FLUID;
     jlist = firstneigh[i];
     jnum = numneigh[i];
 
@@ -179,12 +180,12 @@ void ComputeRHEOInterface::compute_peratom()
 
       if (rsq < cutsq) {
         jtype = type[j];
-        fluidj = status[j] & FixRHEO::STATUS_FLUID;
+        fluidj = status[j] & STATUS_FLUID;
         w = compute_kernel->calc_w_quintic(i, j, delx, dely, delz, sqrt(rsq));
 
         status_match = 0;
         norm[i] += w;
-        if ((fluidi && fluidj) || ((!fluid) && (!fluidj)))
+        if ((fluidi && fluidj) || ((!fluidi) && (!fluidj)))
           status_match = 1;
 
         if (status_match) {
@@ -194,7 +195,7 @@ void ComputeRHEOInterface::compute_peratom()
             dot = (-fp_store[0][j] + fp_store[0][i]) * delx;
             dot += (-fp_store[1][j] + fp_store[1][i]) * dely;
             dot += (-fp_store[2][j] + fp_store[2][i]) * delz;
-            rho[i] += w * (cs * (rho[j] - rho0) - rho[j] * dot);
+            rho[i] += w * (csq * (rho[j] - rho0) - rho[j] * dot);
             normwf[i] += w;
           }
         }
@@ -208,7 +209,7 @@ void ComputeRHEOInterface::compute_peratom()
               dot = (-fp_store[0][i] + fp_store[0][j]) * delx;
               dot += (-fp_store[1][i] + fp_store[1][j]) * dely;
               dot += (-fp_store[2][i] + fp_store[2][j]) * delz;
-              rho[j] += w * (cs * (rho[i] - rho0) + rho[i] * dot);
+              rho[j] += w * (csq * (rho[i] - rho0) + rho[i] * dot);
               normwf[j] += w;
             }
           }
@@ -223,10 +224,10 @@ void ComputeRHEOInterface::compute_peratom()
     if (norm[i] != 0.0) chi[i] /= norm[i];
 
     // Recalculate rho for non-fluid particles
-    if (!(status[i] & FixRHEO::STATUS_FLUID)) {
+    if (!(status[i] & STATUS_FLUID)) {
       if (normwf[i] != 0.0) {
         // Stores rho for solid particles 1+Pw in Adami Adams 2012
-        rho[i] = MAX(EPSILON, rho0 + (rho[i] / normwf[i]) * cs_inv);
+        rho[i] = MAX(EPSILON, rho0 + (rho[i] / normwf[i]) * csq_inv);
       } else {
         rho[i] = rho0;
       }
@@ -310,7 +311,7 @@ void ComputeRHEOInterface::unpack_reverse_comm(int n, int *list, double *buf)
     j = list[i];
     norm[j] += buf[m++];
     chi[j] += buf[m++];
-    if (!(status[j] & FixRHEO::STATUS_FLUID)){
+    if (!(status[j] & STATUS_FLUID)){
       normwf[j] += buf[m++];
       rho[j] += buf[m++];
     } else {
@@ -322,7 +323,7 @@ void ComputeRHEOInterface::unpack_reverse_comm(int n, int *list, double *buf)
 
 /* ---------------------------------------------------------------------- */
 
-void ComputeRHEOInterface::correct_v(double *vi, double *vj, double *vi_out, int i, int j)
+void ComputeRHEOInterface::correct_v(double *vi, double *vj, int i, int j)
 {
   double wall_prefactor, wall_denom, wall_numer;
 
@@ -333,9 +334,9 @@ void ComputeRHEOInterface::correct_v(double *vi, double *vj, double *vi_out, int
 
   wall_prefactor = wall_numer / wall_denom;
 
-  vi_out[0] = (vi[0] - vj[0]) * wall_prefactor + vi[0];
-  vi_out[1] = (vi[1] - vj[1]) * wall_prefactor + vi[1];
-  vi_out[2] = (vi[2] - vj[2]) * wall_prefactor + vi[2];
+  vi[0] = (vi[0] - vj[0]) * wall_prefactor + vi[0];
+  vi[1] = (vi[1] - vj[1]) * wall_prefactor + vi[1];
+  vi[2] = (vi[2] - vj[2]) * wall_prefactor + vi[2];
 }
 
 /* ---------------------------------------------------------------------- */
@@ -352,28 +353,30 @@ double ComputeRHEOInterface::correct_rho(int i, int j)
 void ComputeRHEOInterface::store_forces()
 {
   double minv;
-  double mass = atom->mass;
-  double type = atom->type;
-  double **f = atom->f;
+  int *type = atom->type;
   int *mask = atom->mask;
+  double *mass = atom->mass;
+  double **f = atom->f;
 
   // When this is called, fp_store stores the pressure force
   // After this method, fp_store instead stores non-pressure forces
   // and is also normalized by the particles mass
   // If forces are overwritten by a fix, there are no pressure forces
   // so just normalize
-  int ifix = modify->find_fix_by_style("setforce");
-  if (ifix != -1) {
-    for (int i = 0; i < atom->nlocal; i++) {
-      minv = 1.0 / mass[type[i]];
-      if (mask[i] & modify->fix[ifix]->groupbit) {
-        fp_store[i][0] = f[i][0] * minv;
-        fp_store[i][1] = f[i][1] * minv;
-        fp_store[i][2] = f[i][2] * minv;
-      } else {
-        fp_store[i][0] = (f[i][0] - fp_store[i][0]) * minv;
-        fp_store[i][1] = (f[i][1] - fp_store[i][1]) * minv;
-        fp_store[i][2] = (f[i][2] - fp_store[i][2]) * minv;
+  auto fixlist = modify->get_fix_by_style("setforce");
+  if (fixlist.size() == 0) {
+    for (const auto &fix : fixlist) {
+      for (int i = 0; i < atom->nlocal; i++) {
+        minv = 1.0 / mass[type[i]];
+        if (mask[i] & fix->groupbit) {
+          fp_store[i][0] = f[i][0] * minv;
+          fp_store[i][1] = f[i][1] * minv;
+          fp_store[i][2] = f[i][2] * minv;
+        } else {
+          fp_store[i][0] = (f[i][0] - fp_store[i][0]) * minv;
+          fp_store[i][1] = (f[i][1] - fp_store[i][1]) * minv;
+          fp_store[i][2] = (f[i][2] - fp_store[i][2]) * minv;
+        }
       }
     }
   } else {
@@ -397,7 +400,7 @@ void ComputeRHEOInterface::store_forces()
 
 double ComputeRHEOInterface::memory_usage()
 {
-  double bytes = 3 * nmax_old * sizeof(double);
+  double bytes = 3 * nmax_store * sizeof(double);
   return bytes;
 }
 
diff --git a/src/RHEO/compute_rheo_interface.h b/src/RHEO/compute_rheo_interface.h
index 04733ff334..50cec97790 100644
--- a/src/RHEO/compute_rheo_interface.h
+++ b/src/RHEO/compute_rheo_interface.h
@@ -36,17 +36,17 @@ class ComputeRHEOInterface : public Compute {
   int pack_reverse_comm(int, int, double *) override;
   void unpack_reverse_comm(int, int *, double *) override;
   double memory_usage() override;
-  void correct_v(double *, double *, double *, int, int);
+  void correct_v(double *, double *, int, int);
   double correct_rho(int, int);
   void store_forces();
 
-  double *chi, **f_pressure;
+  double *chi, **fp_store;
   class FixRHEO *fix_rheo;
 
  private:
-  int nmax_old, comm_stage;
-  double rho0, cut, cutsq, cs, cs_inv, wall_max;
-  double *norm, *normwf, **fom_store;
+  int nmax_store, comm_stage;
+  double rho0, cut, cutsq, csq, csq_inv, wall_max;
+  double *norm, *normwf;
 
   char *id_fix_pa;
 
diff --git a/src/RHEO/compute_rheo_kernel.cpp b/src/RHEO/compute_rheo_kernel.cpp
index c0b61daae3..d96d8d234e 100644
--- a/src/RHEO/compute_rheo_kernel.cpp
+++ b/src/RHEO/compute_rheo_kernel.cpp
@@ -42,10 +42,10 @@
 #include <gsl/gsl_cblas.h>
 
 using namespace LAMMPS_NS;
+using namespace RHEO_NS;
 using namespace MathExtra;
 
-enum {QUINTIC, CRK0, CRK1, CRK2};
-#define DELTA 2000
+static constexpr int DELTA = 2000;
 
 /* ---------------------------------------------------------------------- */
 
@@ -55,16 +55,16 @@ ComputeRHEOKernel::ComputeRHEOKernel(LAMMPS *lmp, int narg, char **arg) :
 {
   if (narg != 4) error->all(FLERR,"Illegal compute rheo/kernel command");
 
-  kernel_style = (SubModelType) utils::inumeric(FLERR,arg[3],false,lmp);
+  kernel_style = utils::inumeric(FLERR,arg[3],false,lmp);
 
 
-  if (kernel_style == FixRHEO::QUINTIC) {
+  if (kernel_style == QUINTIC) {
     correction_order = -1;
-  } else if (kernel_style == FixRHEO::CRK0) {
+  } else if (kernel_style == CRK0) {
     correction_order = 0;
-  } else if (kernel_style == FixRHEO::CRK1) {
+  } else if (kernel_style == CRK1) {
     correction_order = 1;
-  } else if (kernel_style == FixRHEO::CRK2) {
+  } else if (kernel_style == CRK2) {
     correction_order = 2;
   }
 
@@ -74,11 +74,11 @@ ComputeRHEOKernel::ComputeRHEOKernel(LAMMPS *lmp, int narg, char **arg) :
   comm_forward = 1;
   ncor = 0;
   Mdim = 0;
-  if (kernel_type == CRK1) {
+  if (kernel_style == CRK1) {
     Mdim = 1 + dim;
     ncor = 1 + dim;
     comm_forward = ncor * Mdim;
-  } else if (kernel_type == CRK2) {
+  } else if (kernel_style == CRK2) {
     //Polynomial basis size (up to quadratic order)
     Mdim = 1 + dim + dim * (dim + 1) / 2;
     //Number of sets of correction coefficients  (1 x y xx yy)  + z zz (3D)
@@ -126,35 +126,35 @@ void ComputeRHEOKernel::init()
 
 
   if (dim == 3) {
-    pre_w = 0.002652582384864922 * 27.0 * ihsq * ih;
-    pre_wp = pre_w * 3.0 * ih;
+    pre_w = 0.002652582384864922 * 27.0 * hsqinv * hinv;
+    pre_wp = pre_w * 3.0 * hinv;
   } else {
-    pre_w = 0.004661441847879780 * 9 * ihsq;
-    pre_wp = pre_w * 3.0 * ih;
+    pre_w = 0.004661441847879780 * 9 * hsqinv;
+    pre_wp = pre_w * 3.0 * hinv;
   }
 
   // Create coordination array if it doesn't already exist
   // Create a custom atom property so it works with compute property/atom
   // Do not create grow callback as there's no reason to copy/exchange data
-  // Manually grow if nmax_old exceeded
+  // Manually grow if nmax_store exceeded
 
   int tmp1, tmp2;
   int nmax = atom->nmax;
   int index = atom->find_custom("rheo_coordination", tmp1, tmp2);
   if (index == -1) {
     index = atom->add_custom("rheo_coordination", 0, 0);
-    nmax_old = nmax;
+    nmax_store = nmax;
   }
   coordination = atom->ivector[index];
 
   // Create local arrays for kernel arrays, I can't foresee a reason to print
 
-  if (kernel_type == CRK0) {
-    memory->create(C0, nmax_old, "rheo/kernel:C0");
-  } else if (kernel_type == CRK1) {
-    memory->create(C, nmax_old, ncor, Mdim, "rheo/kernel:C");
-  } else if (kernel_type == CRK2) {
-    memory->create(C, nmax_old, ncor, Mdim, "rheo/kernel:C");
+  if (kernel_style == CRK0) {
+    memory->create(C0, nmax_store, "rheo/kernel:C0");
+  } else if (kernel_style == CRK1) {
+    memory->create(C, nmax_store, ncor, Mdim, "rheo/kernel:C");
+  } else if (kernel_style == CRK2) {
+    memory->create(C, nmax_store, ncor, Mdim, "rheo/kernel:C");
   }
 }
 
@@ -192,10 +192,10 @@ double ComputeRHEOKernel::calc_w(int i, int j, double delx, double dely, double
   int corrections_j = check_corrections(j);
   int corrections = corrections_i & corrections_j;
 
-  if (kernel_type == QUINTIC || !corrections) w = calc_w_quintic(i,j,delx,dely,delz,r);
-  else if (kernel_type == CRK0) w = calc_w_crk0(i,j,delx,dely,delz,r);
-  else if (kernel_type == CRK1) w = calc_w_crk1(i,j,delx,dely,delz,r);
-  else if (kernel_type == CRK2) w = calc_w_crk2(i,j,delx,dely,delz,r);
+  if (kernel_style == QUINTIC || !corrections) w = calc_w_quintic(i,j,delx,dely,delz,r);
+  else if (kernel_style == CRK0) w = calc_w_crk0(i,j,delx,dely,delz,r);
+  else if (kernel_style == CRK1) w = calc_w_crk1(i,j,delx,dely,delz,r);
+  else if (kernel_style == CRK2) w = calc_w_crk2(i,j,delx,dely,delz,r);
 
   return w;
 }
@@ -211,11 +211,11 @@ double ComputeRHEOKernel::calc_dw(int i, int j, double delx, double dely, double
 
   // Calc wp and default dW's, a bit inefficient but can redo later
   wp = calc_dw_quintic(i,j,delx,dely,delz,r,dWij,dWji);
-  if(kernel_type == CRK1) {
+  if(kernel_style == CRK1) {
     //check if kernel correction calculated successfully. If not, revert to quintic
     if (corrections_i) calc_dw_crk1(i,j,delx,dely,delz,r,dWij);
     if (corrections_j) calc_dw_crk1(j,i,-delx,-dely,-delz,r,dWji);
-  } else if(kernel_type == CRK2) {
+  } else if(kernel_style == CRK2) {
     if (corrections_i) calc_dw_crk2(i,j,delx,dely,delz,r,dWij);
     if (corrections_j) calc_dw_crk2(j,i,-delx,-dely,-delz,r,dWji);
   }
@@ -228,7 +228,7 @@ double ComputeRHEOKernel::calc_dw(int i, int j, double delx, double dely, double
 double ComputeRHEOKernel::calc_w_quintic(int i, int j, double delx, double dely, double delz, double r)
 {
   double w, tmp1, tmp2, tmp3, tmp1sq, tmp2sq, tmp3sq, s;
-  s = r * 3.0 * ih;
+  s = r * 3.0 * hinv;
 
 	if (s > 3.0) {
 	  w = 0.0;
@@ -266,7 +266,7 @@ double ComputeRHEOKernel::calc_dw_quintic(int i, int j, double delx, double dely
   double *mass = atom->mass;
   int *type = atom->type;
 
-  s = r * 3.0 * ih;
+  s = r * 3.0 * hinv;
 
   if (s > 3.0) {
     wp = 0.0;
@@ -496,9 +496,9 @@ void ComputeRHEOKernel::compute_peratom()
   gsl_error_flag = 0;
   gsl_error_tags.clear();
 
-  if (kernel_type == FixRHEO::QUINTIC) return;
+  if (kernel_style == QUINTIC) return;
 
-  int i, j, ii, jj, inum, jnum, g, a, b, gsl_error;
+  int i, j, ii, jj, inum, jnum, itype, g, a, b, gsl_error;
   double xtmp, ytmp, ztmp, r, rsq, w, vj;
   double dx[3];
   gsl_matrix_view gM;
@@ -520,9 +520,9 @@ void ComputeRHEOKernel::compute_peratom()
   firstneigh = list->firstneigh;
 
   // Grow arrays if necessary
-  if (nmax_old < atom->nmax) grow_arrays(atom->nmax);
+  if (nmax_store < atom->nmax) grow_arrays(atom->nmax);
 
-  if (kernel_type == FixRHEO::CRK0) {
+  if (kernel_style == CRK0) {
 
     double M;
     for (ii = 0; ii < inum; ii++) {
@@ -544,12 +544,12 @@ void ComputeRHEOKernel::compute_peratom()
         dx[0] = xtmp - x[j][0];
         dx[1] = ytmp - x[j][1];
         dx[2] = ztmp - x[j][2];
-        rsq = lensq(dx);
+        rsq = lensq3(dx);
 
         if (rsq < hsq) {
           r = sqrt(rsq);
           w = calc_w_quintic(i,j,dx[0],dx[1],dx[2],r);
-          if (!(status[j] & FixRHEO::STATUS_FLUID) && solid_flag) {
+          if (!(status[j] & STATUS_FLUID) && solid_flag) {
             vj = mass[type[j]] / compute_interface->correct_rho(j,i);
           } else vj = mass[type[j]] / rho[j];
 
@@ -590,22 +590,24 @@ void ComputeRHEOKernel::compute_peratom()
         dx[1] = ytmp - x[j][1];
         dx[2] = ztmp - x[j][2];
 
-        rsq = lensq(dx);
+        rsq = lensq3(dx);
 
-        if (rsq < cutsq) {
+        if (rsq < hsq) {
           r = sqrt(rsq);
           w = calc_w_quintic(i,j,dx[0],dx[1],dx[2],r);
 
-          if (status[j] > FixRHEO::FLUID_MAX && solid_flag)
-            vj = mass[type[j]]/compute_interface->correct_rho(j,i);
-          else vj = mass[type[j]]/rho[j];
+          if (solid_flag)
+            if (!(status[j] & STATUS_FLUID))
+              vj = mass[type[j]]/compute_interface->correct_rho(j,i);
+          else
+            vj = mass[type[j]]/rho[j];
 
           //Populate the H-vector of polynomials (2D)
           if (dim == 2) {
             H[0] = 1.0;
             H[1] = dx[0] * hinv;
             H[2] = dx[1] * hinv;
-            if (kernel_type == FixRHEO::CRK2) {
+            if (kernel_style == CRK2) {
               H[3] = 0.5 * dx[0] * dx[0] * hsqinv;
               H[4] = 0.5 * dx[1] * dx[1] * hsqinv;
               H[5] = dx[0] * dx[1] * hsqinv;
@@ -615,7 +617,7 @@ void ComputeRHEOKernel::compute_peratom()
             H[1] = dx[0] * hinv;
             H[2] = dx[1] * hinv;
             H[3] = dx[2] * hinv;
-            if (kernel_type == FixRHEO::CRK2) {
+            if (kernel_style == CRK2) {
               H[4] = 0.5 * dx[0] * dx[0] * hsqinv;
               H[5] = 0.5 * dx[1] * dx[1] * hsqinv;
               H[6] = 0.5 * dx[2] * dx[2] * hsqinv;
@@ -642,7 +644,7 @@ void ComputeRHEOKernel::compute_peratom()
       }
 
       // Skip if undercoordinated
-      if (coordination[i] < zmin) continue
+      if (coordination[i] < zmin) continue;
 
       // Use gsl to get Minv, use Cholesky decomposition since the
       // polynomials are independent, M is symmetrix & positive-definite
@@ -656,7 +658,7 @@ void ComputeRHEOKernel::compute_peratom()
 
         //check if not positive-definite
         if (gsl_error != GSL_EDOM)
-          error->warn(FLERR, "Failed decomposition in rheo/kernel, gsl_error = {}", gsl_error);
+          error->warning(FLERR, "Failed decomposition in rheo/kernel, gsl_error = {}", gsl_error);
 
         continue;
       }
@@ -689,7 +691,7 @@ void ComputeRHEOKernel::compute_peratom()
           // columns 1-2 (2D)  or 1-3 (3D)
 
           //Second derivatives
-          if (kernel_type == FixRHEO::CRK2)
+          if (kernel_style == CRK2)
             C[i][1 + dim + b][a] = M[a * Mdim + b + 1 + dim] * hsqinv;
             // columns 3-4 (2D) or 4-6 (3D)
         }
@@ -721,7 +723,7 @@ void ComputeRHEOKernel::compute_coordination()
   firstneigh = list->firstneigh;
 
   // Grow arrays if necessary
-  if (nmax_old < atom->nmax) grow_arrays(atom->nmax);
+  if (nmax_store < atom->nmax) grow_arrays(atom->nmax);
 
   for (ii = 0; ii < inum; ii++) {
     i = ilist[ii];
@@ -740,7 +742,7 @@ void ComputeRHEOKernel::compute_coordination()
       dx[0] = xtmp - x[j][0];
       dx[1] = ytmp - x[j][1];
       dx[2] = ztmp - x[j][2];
-      rsq = lensq(dx);
+      rsq = lensq3(dx);
 
       if (rsq < hsq)
         coordination[i] += 1;
@@ -759,13 +761,13 @@ void ComputeRHEOKernel::grow_arrays(int nmax)
 {
   memory->grow(coordination, nmax, "atom:rheo_coordination");
 
-  if (kernel_type == FixRHEO::CRK0) {
+  if (kernel_style == CRK0) {
     memory->grow(C0, nmax, "rheo/kernel:C0");
   } else if (correction_order > 0) {
     memory->grow(C, nmax, ncor, Mdim, "rheo/kernel:C");
   }
 
-  nmax_old = nmax;
+  nmax_store = nmax;
 }
 
 /* ---------------------------------------------------------------------- */
@@ -781,7 +783,7 @@ int ComputeRHEOKernel::pack_forward_comm(int n, int *list, double *buf,
     if (comm_stage == 0) {
       buf[m++] = coordination[j];
     } else {
-      if (kernel_type == FixRHEO::CRK0) {
+      if (kernel_style == CRK0) {
         buf[m++] = C0[j];
       } else {
         for (a = 0; a < ncor; a++)
@@ -805,7 +807,7 @@ void ComputeRHEOKernel::unpack_forward_comm(int n, int first, double *buf)
     if (comm_stage == 0) {
       coordination[i] = buf[m++];
     } else {
-      if (kernel_type == FixRHEO::CRK0) {
+      if (kernel_style == CRK0) {
         C0[i] = buf[m++];
       } else {
         for (a = 0; a < ncor; a++)
@@ -821,12 +823,12 @@ void ComputeRHEOKernel::unpack_forward_comm(int n, int first, double *buf)
 double ComputeRHEOKernel::memory_usage()
 {
   double bytes = 0.0;
-  bytes = (size_t) nmax_old * sizeof(int);
+  bytes = (size_t) nmax_store * sizeof(int);
 
-  if (kernel_type == FixRHEO::CRK0) {
-    bytes += (size_t) nmax_old * sizeof(double);
+  if (kernel_style == CRK0) {
+    bytes += (size_t) nmax_store * sizeof(double);
   } else if (correction_order > 0) {
-    bytes += (size_t) nmax_old * ncor * Mdim * sizeof(double);
+    bytes += (size_t) nmax_store * ncor * Mdim * sizeof(double);
   }
   return bytes;
 }
diff --git a/src/RHEO/compute_rheo_kernel.h b/src/RHEO/compute_rheo_kernel.h
index 19062e483b..1842406977 100644
--- a/src/RHEO/compute_rheo_kernel.h
+++ b/src/RHEO/compute_rheo_kernel.h
@@ -54,7 +54,7 @@ class ComputeRHEOKernel : public Compute {
   std::unordered_set<tagint> gsl_error_tags;
 
   int kernel_style, zmin, dim, Mdim, ncor;
-  int nmax_old;
+  int nmax_store;
   double h, hsq, hinv, hsqinv, pre_w, pre_wp;
   double ***C;
   double *C0;
diff --git a/src/RHEO/compute_rheo_rho_sum.cpp b/src/RHEO/compute_rheo_rho_sum.cpp
index fafd948538..726d876ea1 100644
--- a/src/RHEO/compute_rheo_rho_sum.cpp
+++ b/src/RHEO/compute_rheo_rho_sum.cpp
@@ -49,7 +49,7 @@ void ComputeRHEORhoSum::init()
   cutsq = cut * cut;
 
   // need an occasional half neighbor list
-  neighbor->add_request(this, NeighConst::REQ_HALF);
+  neighbor->add_request(this, NeighConst::REQ_DEFAULT);
 }
 
 /* ---------------------------------------------------------------------- */
@@ -130,8 +130,9 @@ void ComputeRHEORhoSum::compute_peratom()
 int ComputeRHEORhoSum::pack_forward_comm(int n, int *list, double *buf,
                                         int /*pbc_flag*/, int * /*pbc*/)
 {
-  int i,j,k,m;
-  double * rho = atom->rho;
+  int i, j, k, m;
+  double *rho = atom->rho;
+  int *coordination = compute_kernel->coordination;
   m = 0;
 
   for (i = 0; i < n; i++) {
@@ -145,7 +146,7 @@ int ComputeRHEORhoSum::pack_forward_comm(int n, int *list, double *buf,
 void ComputeRHEORhoSum::unpack_forward_comm(int n, int first, double *buf)
 {
   int i, k, m, last;
-  double * rho = atom->rho;
+  double *rho = atom->rho;
 
   m = 0;
   last = first + n;
@@ -158,7 +159,7 @@ void ComputeRHEORhoSum::unpack_forward_comm(int n, int first, double *buf)
 
 int ComputeRHEORhoSum::pack_reverse_comm(int n, int first, double *buf)
 {
-  int i,k,m,last;
+  int i, k, m, last;
   double *rho = atom->rho;
 
   m = 0;
@@ -173,7 +174,7 @@ int ComputeRHEORhoSum::pack_reverse_comm(int n, int first, double *buf)
 
 void ComputeRHEORhoSum::unpack_reverse_comm(int n, int *list, double *buf)
 {
-  int i,k,j,m;
+  int i, k, j, m;
   double *rho = atom->rho;
 
   m = 0;
diff --git a/src/RHEO/compute_rheo_surface.cpp b/src/RHEO/compute_rheo_surface.cpp
index 9a94ea6a76..180c430dd1 100644
--- a/src/RHEO/compute_rheo_surface.cpp
+++ b/src/RHEO/compute_rheo_surface.cpp
@@ -20,8 +20,8 @@
 
 #include "atom.h"
 #include "comm.h"
+#include "compute_rheo_interface.h"
 #include "compute_rheo_kernel.h"
-#include "compute_rheo_solids.h"
 #include "domain.h"
 #include "error.h"
 #include "fix_rheo.h"
@@ -33,18 +33,19 @@
 #include "neigh_request.h"
 
 using namespace LAMMPS_NS;
+using namespace RHEO_NS;
 using namespace FixConst;
 using namespace MathExtra;
 
-#define EPSILON 1e-10;
+static constexpr double EPSILON = 1e-10;
 
 /* ---------------------------------------------------------------------- */
 
 ComputeRHEOSurface::ComputeRHEOSurface(LAMMPS *lmp, int narg, char **arg) :
-  Fix(lmp, narg, arg), fix_rheo(nullptr), list(nullptr), compute_kernel(nullptr), compute_solids(nullptr),
+  Compute(lmp, narg, arg), fix_rheo(nullptr), list(nullptr), compute_kernel(nullptr), compute_interface(nullptr),
   B(nullptr), gradC(nullptr), nsurface(nullptr), divr(nullptr), rsurface(nullptr)
 {
-  if (narg != 3) error->all(FLERR,"Illegal fix RHEO/SURFACE command");
+  if (narg != 3) error->all(FLERR,"Illegal compute RHEO/SURFACE command");
 
   int dim = domain->dimension;
   comm_forward = 2;
@@ -75,19 +76,19 @@ ComputeRHEOSurface::~ComputeRHEOSurface()
 void ComputeRHEOSurface::init()
 {
   compute_kernel = fix_rheo->compute_kernel;
-  compute_solids = fix_rheo->compute_solids;
+  compute_interface = fix_rheo->compute_interface;
   cut = fix_rheo->cut;
   rho0 = fix_rheo->rho0;
   threshold_style = fix_rheo->surface_style;
-  threshold_divr = fix_rheo->divrsurface;
-  threshold_z = fix_rheo->zminsurface;
+  threshold_divr = fix_rheo->divr_surface;
+  threshold_z = fix_rheo->zmin_surface;
 
   cutsq = cut * cut;
 
   // Create rsurface, divr, nsurface arrays if they don't already exist
   // Create a custom atom property so it works with compute property/atom
   // Do not create grow callback as there's no reason to copy/exchange data
-  // Manually grow if nmax_old exceeded
+  // Manually grow if nmax_store exceeded
   // For B and gradC, create a local array since they are unlikely to be printed
 
   int tmp1, tmp2;
@@ -103,12 +104,13 @@ void ComputeRHEOSurface::init()
   if (index == -1)  index = atom->add_custom("rheo_nsurface", 1, 3);
   nsurface = atom->darray[index];
 
-  nmax_old = atom->nmax;
-  memory->create(B, nmax_old, dim * dim, "rheo/surface:B");
-  memory->create(gradC, nmax_old, dim * dim, "rheo/surface:gradC");
+  nmax_store = atom->nmax;
+  int dim = domain->dimension;
+  memory->create(B, nmax_store, dim * dim, "rheo/surface:B");
+  memory->create(gradC, nmax_store, dim * dim, "rheo/surface:gradC");
 
   // need an occasional half neighbor list
-  neighbor->add_request(this, NeighConst::REQ_HALF);
+  neighbor->add_request(this, NeighConst::REQ_DEFAULT);
 }
 
 /* ---------------------------------------------------------------------- */
@@ -123,9 +125,8 @@ void ComputeRHEOSurface::init_list(int /*id*/, NeighList *ptr)
 void ComputeRHEOSurface::compute_peratom()
 {
   int i, j, ii, jj, inum, jnum, a, b, itype, jtype, fluidi, fluidj;
-  double xtmp, ytmp, ztmp, rsq, Voli, Volj, rhoi, rhoj;
-  double *dWij, *dWji;
-  double dx[3];
+  double xtmp, ytmp, ztmp, rsq, Voli, Volj, rhoi, rhoj, wp;
+  double *dWij, *dWji, dx[3];
   int *ilist, *jlist, *numneigh, **firstneigh;
 
   int nlocal = atom->nlocal;
@@ -146,7 +147,7 @@ void ComputeRHEOSurface::compute_peratom()
   firstneigh = list->firstneigh;
 
   int nmax = atom->nmax;
-  if (nmax_old <= nmax) {
+  if (nmax_store <= nmax) {
     memory->grow(divr, nmax, "atom:rheo_divr");
     memory->grow(rsurface, nmax, "atom:rheo_rsurface");
     memory->grow(nsurface, nmax, 3, "atom:rheo_nsurface");
@@ -154,7 +155,7 @@ void ComputeRHEOSurface::compute_peratom()
     memory->grow(B, nmax, dim * dim, "rheo/surface:B");
     memory->grow(gradC, nmax, dim * dim, "rheo/surface:gradC");
 
-    nmax_old = atom->nmax;
+    nmax_store = atom->nmax;
   }
 
   int nall = nlocal + atom->nghost;
@@ -169,7 +170,7 @@ void ComputeRHEOSurface::compute_peratom()
     divr[i] = 0.0;
 
     // Remove surface settings
-    status[i] &= FixRHEO::surfacemask;
+    status[i] &= SURFACEMASK;
   }
 
   // loop over neighbors to calculate the average orientation of neighbors
@@ -182,7 +183,7 @@ void ComputeRHEOSurface::compute_peratom()
     jlist = firstneigh[i];
     jnum = numneigh[i];
     itype = type[i];
-    fluidi = status[i] & FixRHEO::STATUS_FLUID;
+    fluidi = status[i] & STATUS_FLUID;
 
     for (jj = 0; jj < jnum; jj++) {
       j = jlist[jj];
@@ -192,19 +193,19 @@ void ComputeRHEOSurface::compute_peratom()
       dx[1] = ytmp - x[j][1];
       dx[2] = ztmp - x[j][2];
 
-      rsq = lensq(dx);
+      rsq = lensq3(dx);
       if (rsq < cutsq) {
         jtype = type[j];
-        fluidj = status[j] & FixRHEO::STATUS_FLUID;
+        fluidj = status[j] & STATUS_FLUID;
 
         rhoi = rho[i];
         rhoj = rho[j];
 
         // Add corrections for walls
         if (fluidi && (!fluidj)) {
-          rhoj = compute_solids->correct_rho(j, i);
+          rhoj = compute_interface->correct_rho(j, i);
         } else if ((!fluidi) && fluidj) {
-          rhoi = compute_solids->correct_rho(i, j);
+          rhoi = compute_interface->correct_rho(i, j);
         } else if ((!fluidi) && (!fluidj)) {
           rhoi = rho0;
           rhoj = rho0;
@@ -253,29 +254,29 @@ void ComputeRHEOSurface::compute_peratom()
   }
 
   // Find the free-surface
-  if (threshold_style == FixRHEO::DIVR) {
+  if (threshold_style == DIVR) {
     for (i = 0; i < nall; i++) {
       if (mask[i] & groupbit) {
-        status[i] |= FixRHEO::STATUS_BULK;
+        status[i] |= STATUS_BULK;
         rsurface[i] = cut;
         if (divr[i] < threshold_divr) {
-          status[i] |= FixRHEO::STATUS_SURFACE;
+          status[i] |= STATUS_SURFACE;
           rsurface[i] = 0.0;
           if (coordination[i] < threshold_z)
-            status[i] |= FixRHEO::STATUS_SPLASH;
+            status[i] |= STATUS_SPLASH;
         }
       }
     }
   } else {
     for (i = 0; i < nall; i++) {
       if (mask[i] & groupbit) {
-        status[i] |= FixRHEO::STATUS_BULK;
+        status[i] |= STATUS_BULK;
         rsurface[i] = cut;
-        if (coordination[i] < divR_limit) {
-          status[i] |= FixRHEO::STATUS_SURFACE;
+        if (coordination[i] < threshold_divr) {
+          status[i] |= STATUS_SURFACE;
           rsurface[i] = 0.0;
           if (coordination[i] < threshold_z)
-            status[i] |= FixRHEO::STATUS_SPLASH;
+            status[i] |= STATUS_SPLASH;
         }
       }
     }
@@ -297,23 +298,23 @@ void ComputeRHEOSurface::compute_peratom()
       dx[0] = xtmp - x[j][0];
       dx[1] = ytmp - x[j][1];
       dx[2] = ztmp - x[j][2];
-      rsq = lensq(dx);
+      rsq = lensq3(dx);
       if (rsq < cutsq) {
-        if ((status[i] & FixRHEO::STATUS_BULK) && (status[j] & FixRHEO::STATUS_SURFACE)) {
-          status[i] &= FixRHEO::surfacemask;
-          status[i] |= FixRHEO::STATUS_LAYER;
+        if ((status[i] & STATUS_BULK) && (status[j] & STATUS_SURFACE)) {
+          status[i] &= SURFACEMASK;
+          status[i] |= STATUS_LAYER;
         }
 
-        if (status[j] & FixRHEO::STATUS_SURFACE) rsurface[i] = MIN(rsurface[i], sqrt(rsq));
+        if (status[j] & STATUS_SURFACE) rsurface[i] = MIN(rsurface[i], sqrt(rsq));
 
 
         if (j < nlocal || newton) {
-          if ((status[j] & FixRHEO::STATUS_BULK) && (status[i] & FixRHEO::STATUS_SURFACE)) {
-            status[j] &= FixRHEO::surfacemask;
-            status[j] |= FixRHEO::STATUS_LAYER;
+          if ((status[j] & STATUS_BULK) && (status[i] & STATUS_SURFACE)) {
+            status[j] &= SURFACEMASK;
+            status[j] |= STATUS_LAYER;
           }
 
-          if (status[i] & FixRHEO::STATUS_SURFACE) rsurface[j] = MIN(rsurface[j], sqrt(rsq));
+          if (status[i] & STATUS_SURFACE) rsurface[j] = MIN(rsurface[j], sqrt(rsq));
         }
       }
     }
@@ -371,7 +372,7 @@ void ComputeRHEOSurface::unpack_reverse_comm(int n, int *list, double *buf)
     } else if (comm_stage == 1) {
 
       temp = (int) buf[m++];
-      if ((status[j] & FixRHEO::STATUS_BULK) && (temp & FixRHEO::STATUS_LAYER))
+      if ((status[j] & STATUS_BULK) && (temp & STATUS_LAYER))
         status[j] = temp;
 
       rsurface[j] = MIN(rsurface[j], buf[m++]);
diff --git a/src/RHEO/compute_rheo_surface.h b/src/RHEO/compute_rheo_surface.h
index 224b2594a1..58a3e3b9c4 100644
--- a/src/RHEO/compute_rheo_surface.h
+++ b/src/RHEO/compute_rheo_surface.h
@@ -17,8 +17,8 @@ ComputeStyle(RHEO/SURFACE,ComputeRHEOSurface)
 // clang-format on
 #else
 
-#ifndef LMP_COMPUTE_RHEO_INTERFACE_H
-#define LMP_COMPUTE_RHEO_INTERFACE_H
+#ifndef LMP_COMPUTE_RHEO_SURFACE_H
+#define LMP_COMPUTE_RHEO_SURFACE_H
 
 #include "compute.h"
 
@@ -36,18 +36,18 @@ class ComputeRHEOSurface : public Compute {
   int pack_forward_comm(int, int *, double *, int, int *) override;
   void unpack_forward_comm(int, int, double *) override;
 
-  double **nsurface, **rsurface;
+  double **nsurface, *rsurface;
   class FixRHEO *fix_rheo;
 
  private:
   double cut, cutsq, rho0, threshold_divr;
-  int surface_style, nmax_old, threshold_z;
+  int surface_style, nmax_store, threshold_z;
   double **B, **gradC, *divr;
   int threshold_style, comm_stage;
 
   class NeighList *list;
   class ComputeRHEOKernel *compute_kernel;
-  class ComputeRHEOSolids *compute_solids;
+  class ComputeRHEOInterface *compute_interface;
 };
 
 }    // namespace LAMMPS_NS
diff --git a/src/RHEO/compute_rheo_vshift.cpp b/src/RHEO/compute_rheo_vshift.cpp
index 1b3edcffd8..440bfc7fc7 100644
--- a/src/RHEO/compute_rheo_vshift.cpp
+++ b/src/RHEO/compute_rheo_vshift.cpp
@@ -32,6 +32,7 @@
 #include "neigh_request.h"
 
 using namespace LAMMPS_NS;
+using namespace RHEO_NS;
 
 /* ---------------------------------------------------------------------- */
 
@@ -47,15 +48,15 @@ ComputeRHEOVShift::ComputeRHEOVShift(LAMMPS *lmp, int narg, char **arg) :
   // Create vshift array if it doesn't already exist
   // Create a custom atom property so it works with compute property/atom
   // Do not create grow callback as there's no reason to copy/exchange data
-  // Manually grow if nmax_old exceeded
+  // Manually grow if nmax_store exceeded
 
   int tmp1, tmp2;
   int index = atom->find_custom("rheo_vshift", tmp1, tmp2);
   if (index == -1) {
     index = atom->add_custom("rheo_vshift", 1, 3);
-    nmax_old = atom->nmax;
+    nmax_store = atom->nmax;
   }
-  vshift = atom->dvector[index];
+  vshift = atom->darray[index];
 }
 
 /* ---------------------------------------------------------------------- */
@@ -108,16 +109,17 @@ void ComputeRHEOVShift::compute_peratom()
 
   int *jlist;
   int inum, *ilist, *numneigh, **firstneigh;
-  int nlocal = atom->nlocal;
-  int nall = nlocal + atom->nghost;
 
-  double **x = atom->x;
-  double **v = atom->v;
   int *type = atom->type;
   int *status = atom->status;
-  int *surface = atom->surface;
+  int *mask = atom->mask;
+  double **x = atom->x;
+  double **v = atom->v;
   double *rho = atom->rho;
   double *mass = atom->mass;
+
+  int nlocal = atom->nlocal;
+  int nall = nlocal + atom->nghost;
   int newton_pair = force->newton_pair;
 
   inum = list->inum;
@@ -125,9 +127,9 @@ void ComputeRHEOVShift::compute_peratom()
   numneigh = list->numneigh;
   firstneigh = list->firstneigh;
 
-  if (nmax_old < atom->nmax) {
+  if (nmax_store < atom->nmax) {
     memory->grow(vshift, atom->nmax, 3, "atom:rheo_vshift");
-    nmax_old = atom->nmax;
+    nmax_store = atom->nmax;
   }
 
   for (i = 0; i < nall; i++)
@@ -143,15 +145,15 @@ void ComputeRHEOVShift::compute_peratom()
     jlist = firstneigh[i];
     jnum = numneigh[i];
     imass = mass[itype];
-    fluidi = status[i] & FixRHEO::STATUS_FLUID;
+    fluidi = status[i] & STATUS_FLUID;
 
     for (jj = 0; jj < jnum; jj++) {
       j = jlist[jj];
       j &= NEIGHMASK;
 
-      fluidj = status[j] & FixRHEO::STATUS_FLUID;
+      fluidj = status[j] & STATUS_FLUID;
       if ((!fluidi) && (!fluidj)) continue;
-      if (!(status[i] & FixRHEO::STATUS_SHIFT) && !(status[j] & FixRHEO::STATUS_SHIFT)) continue;
+      if (!(status[i] & STATUS_SHIFT) && !(status[j] & STATUS_SHIFT)) continue;
 
       dx[0] = xtmp - x[j][0];
       dx[1] = ytmp - x[j][1];
@@ -175,10 +177,10 @@ void ComputeRHEOVShift::compute_peratom()
 
         // Add corrections for walls
         if (fluidi && (!fluidj)) {
-          compute_interface->correct_v(v[i], v[j], vi, i, j);
+          compute_interface->correct_v(vi, vj, i, j);
           rhoj = compute_interface->correct_rho(j,i);
         } else if ((!fluidi) && fluidj) {
-          compute_interface->correct_v(v[j], v[i], vj, j, i);
+          compute_interface->correct_v(vj, vi, j, i);
           rhoi = compute_interface->correct_rho(i,j);
         } else if ((!fluidi) && (!fluidj)) {
           rhoi = 1.0;
@@ -215,7 +217,7 @@ void ComputeRHEOVShift::compute_peratom()
     }
   }
 
-  if (newton_pair) comm->reverse_comm_compute(this);
+  if (newton_pair) comm->reverse_comm(this);
 }
 
 
@@ -239,7 +241,7 @@ void ComputeRHEOVShift::correct_surfaces()
   double nx,ny,nz,vx,vy,vz;
   for (i = 0; i < nlocal; i++) {
     if (mask[i] & groupbit) {
-      if ((status[i] & FixRHEO::STATUS_SURFACE) || (status[i] & FixRHEO::STATUS_LAYER)) {
+      if ((status[i] & STATUS_SURFACE) || (status[i] & STATUS_LAYER)) {
         nx = nsurf[i][0];
         ny = nsurf[i][1];
         vx = vshift[i][0];
@@ -297,6 +299,6 @@ void ComputeRHEOVShift::unpack_reverse_comm(int n, int *list, double *buf)
 
 double ComputeRHEOVShift::memory_usage()
 {
-  double bytes = 3 * nmax_old * sizeof(double);
+  double bytes = 3 * nmax_store * sizeof(double);
   return bytes;
 }
diff --git a/src/RHEO/compute_rheo_vshift.h b/src/RHEO/compute_rheo_vshift.h
index 88a9cdcd1d..8611e177d1 100644
--- a/src/RHEO/compute_rheo_vshift.h
+++ b/src/RHEO/compute_rheo_vshift.h
@@ -40,7 +40,7 @@ class ComputeRHEOVShift : public Compute {
   class FixRHEO *fix_rheo;
 
  private:
-  int nmax_old;
+  int nmax_store;
   double dtv, cut, cutsq, cutthird;
   int surface_flag;
 
diff --git a/src/RHEO/fix_rheo.cpp b/src/RHEO/fix_rheo.cpp
index 2b55320c4e..ff804fe007 100644
--- a/src/RHEO/fix_rheo.cpp
+++ b/src/RHEO/fix_rheo.cpp
@@ -33,6 +33,7 @@
 #include "utils.h"
 
 using namespace LAMMPS_NS;
+using namespace RHEO_NS;
 using namespace FixConst;
 
 /* ---------------------------------------------------------------------- */
@@ -68,6 +69,7 @@ FixRHEO::FixRHEO(LAMMPS *lmp, int narg, char **arg) :
     error->all(FLERR,"Insufficient arguments for fix rheo command");
 
   h = utils::numeric(FLERR,arg[3],false,lmp);
+  cut = h;
   if (strcmp(arg[4],"Quintic") == 0) {
       kernel_style = QUINTIC;
   } else if (strcmp(arg[4],"CRK0") == 0) {
@@ -101,7 +103,7 @@ FixRHEO::FixRHEO(LAMMPS *lmp, int narg, char **arg) :
       iarg += 2;
     } else if (strcmp(arg[iarg],"interface/reconstruction") == 0) {
       interface_flag = 1;
-    } else if (strcmp(arg[iarg],"rhosum") == 0) {
+    } else if (strcmp(arg[iarg],"rho/sum") == 0) {
       rhosum_flag = 1;
     } else if (strcmp(arg[iarg],"rho0") == 0) {
       if(iarg + 1 >= narg) error->all(FLERR,"Illegal rho0 option in fix rheo");
@@ -137,7 +139,8 @@ FixRHEO::~FixRHEO()
 
 void FixRHEO::post_constructor()
 {
-  compute_kernel = dynamic_cast<ComputeRHEOKernel *>(modify->add_compute("rheo_kernel all RHEO/KERNEL {}", kernel_style));
+  compute_kernel = dynamic_cast<ComputeRHEOKernel *>(modify->add_compute(
+    fmt::format("rheo_kernel all RHEO/KERNEL {}", kernel_style)));
   compute_kernel->fix_rheo = this;
 
   std::string cmd = "rheo_grad all RHEO/GRAD velocity rho viscosity";
@@ -146,22 +149,26 @@ void FixRHEO::post_constructor()
   compute_grad->fix_rheo = this;
 
   if (rhosum_flag) {
-    compute_rhosum = dynamic_cast<ComputeRHEORhoSum *>(modify->add_compute("rheo_rho_sum all RHEO/RHO/SUM"));
+    compute_rhosum = dynamic_cast<ComputeRHEORhoSum *>(modify->add_compute(
+      "rheo_rho_sum all RHEO/RHO/SUM"));
     compute_rhosum->fix_rheo = this;
   }
 
   if (shift_flag) {
-    compute_vshift = dynamic_cast<ComputeRHEOVShift *>(modify->add_compute("rheo_vshift all RHEO/VSHIFT"));
+    compute_vshift = dynamic_cast<ComputeRHEOVShift *>(modify->add_compute(
+      "rheo_vshift all RHEO/VSHIFT"));
     compute_vshift->fix_rheo = this;
   }
 
   if (interface_flag) {
-    compute_interface = dynamic_cast<ComputeRHEOInterface *>(modify->add_compute(fmt::format("rheo_interface all RHEO/INTERFACE")));
+    compute_interface = dynamic_cast<ComputeRHEOInterface *>(modify->add_compute(
+      "rheo_interface all RHEO/INTERFACE"));
     compute_interface->fix_rheo = this;
   }
 
   if (surface_flag) {
-    compute_surface = dynamic_cast<ComputeRHEOSurface *>(modify->add_compute(fmt::format("rheo_surface all RHEO/SURFACE")));
+    compute_surface = dynamic_cast<ComputeRHEOSurface *>(modify->add_compute(
+      "rheo_surface all RHEO/SURFACE"));
     compute_surface->fix_rheo = this;
   }
 }
@@ -193,7 +200,7 @@ void FixRHEO::init()
 void FixRHEO::setup_pre_force(int /*vflag*/)
 {
   // Check to confirm accessory fixes do not preceed FixRHEO
-  // Note: these fixes set this flag in setup_pre_force()
+  // Note: fixes set this flag in setup_pre_force()
   if (viscosity_fix_defined || pressure_fix_defined || thermal_fix_defined)
     error->all(FLERR, "Fix RHEO must be defined before all other RHEO fixes");
 
@@ -206,23 +213,23 @@ void FixRHEO::setup_pre_force(int /*vflag*/)
 void FixRHEO::setup(int /*vflag*/)
 {
   // Confirm all accessory fixes are defined
-  // Note: these fixes set this flag in setup_pre_force()
+  // Note: fixes set this flag in setup_pre_force()
   if (!viscosity_fix_defined)
     error->all(FLERR, "Missing fix rheo/viscosity");
 
   if (!pressure_fix_defined)
     error->all(FLERR, "Missing fix rheo/pressure");
 
-  if(!thermal_fix_defined && thermal_flag)
+  if((!thermal_fix_defined) && thermal_flag)
     error->all(FLERR, "Missing fix rheo/thermal");
 
-  // Reset to zero for next run
+  // Reset to zero for future runs
   thermal_fix_defined = 0;
   viscosity_fix_defined = 0;
   pressure_fix_defined = 0;
 
-  // Check fixes cover all atoms (doesnt ensure user covers atoms created midrun)
-  // (pressure is currently required to be group all)
+  // Check fixes cover all atoms (may still fail if atoms are created)
+  // FixRHEOPressure currently requires group all
   auto visc_fixes = modify->get_fix_by_style("rheo/viscosity");
   auto therm_fixes = modify->get_fix_by_style("rheo/thermal");
 
@@ -232,12 +239,12 @@ void FixRHEO::setup(int /*vflag*/)
   int covered;
   for (int i = 0; i < atom->nlocal; i++) {
     covered = 0;
-    for (auto fix in visc_fixes)
+    for (auto fix : visc_fixes)
       if (mask[i] & fix->groupbit) covered = 1;
     if (!covered) v_coverage_flag = 0;
     if (thermal_flag) {
       covered = 0;
-      for (auto fix in therm_fixes)
+      for (auto fix : therm_fixes)
         if (mask[i] & fix->groupbit) covered = 1;
       if (!covered) v_coverage_flag = 0;
     }
@@ -253,11 +260,12 @@ void FixRHEO::setup(int /*vflag*/)
 
 void FixRHEO::initial_integrate(int /*vflag*/)
 {
-  // update v and x and rho of atoms in group
+  // update v, x and rho of atoms in group
   int i, a, b;
   double dtfm, divu;
-  int dim = domain->dimension;
 
+  int *type = atom->type;
+  int *mask = atom->mask;
   int *status = atom->status;
   double **x = atom->x;
   double **v = atom->v;
@@ -266,16 +274,14 @@ void FixRHEO::initial_integrate(int /*vflag*/)
   double *drho = atom->drho;
   double *mass = atom->mass;
   double *rmass = atom->rmass;
-  int rmass_flag = atom->rmass_flag;
-
   double **gradr = compute_grad->gradr;
   double **gradv = compute_grad->gradv;
   double **vshift;
   if (shift_flag) compute_vshift->vshift;
 
-  int *type = atom->type;
-  int *mask = atom->mask;
   int nlocal = atom->nlocal;
+  int rmass_flag = atom->rmass_flag;
+  int dim = domain->dimension;
 
   if (igroup == atom->firstgroup)
     nlocal = atom->nfirst;
@@ -333,7 +339,7 @@ void FixRHEO::initial_integrate(int /*vflag*/)
 
   // Shifting atoms
   if (shift_flag) {
-    compute_vshift->correct_surfaces(); // COuld this be moved to preforce after the surface fix runs?
+    compute_vshift->correct_surfaces(); // Could this be moved to preforce after the surface fix runs?
     for (i = 0; i < nlocal; i++) {
 
       if (!(status[i] & STATUS_SHIFT)) continue;
@@ -376,6 +382,7 @@ void FixRHEO::pre_force(int /*vflag*/)
     compute_vshift->compute_peratom();
 
   // Remove extra shifting/no force options
+  int *mask = atom->mask;
   int *status = atom->status;
   int nall = atom->nlocal + atom->nghost;
   for (int i = 0; i < nall; i++) {
@@ -393,26 +400,28 @@ void FixRHEO::pre_force(int /*vflag*/)
 
 /* ---------------------------------------------------------------------- */
 
-void FixRHEO::final_integrate() {
-  int *status = atom->status;
-  double **gradv = compute_grad->gradv;
-  double **x = atom->x;
-  double **v = atom->v;
-  double **f = atom->f;
-
-  double *rho = atom->rho;
-  double *drho = atom->drho;
-  int *type = atom->type;
-  int *mask = atom->mask;
-  double *mass = atom->mass;
+void FixRHEO::final_integrate()
+{
   int nlocal = atom->nlocal;
   if (igroup == atom->firstgroup)
     nlocal = atom->nfirst;
+
   double dtfm, divu;
-  double *rmass = atom->rmass;
-  int rmass_flag = atom->rmass_flag;
   int i, a;
 
+  double **x = atom->x;
+  double **v = atom->v;
+  double **f = atom->f;
+  double **gradv = compute_grad->gradv;
+  double *rho = atom->rho;
+  double *drho = atom->drho;
+  double *mass = atom->mass;
+  double *rmass = atom->rmass;
+  int *type = atom->type;
+  int *mask = atom->mask;
+  int *status = atom->status;
+
+  int rmass_flag = atom->rmass_flag;
   int dim = domain->dimension;
 
   // Update velocity
diff --git a/src/RHEO/fix_rheo.h b/src/RHEO/fix_rheo.h
index 0064f4c90b..1d8ae06159 100644
--- a/src/RHEO/fix_rheo.h
+++ b/src/RHEO/fix_rheo.h
@@ -39,33 +39,10 @@ class FixRHEO : public Fix {
   void reset_dt() override;
 
   // Model parameters
-  double h, rho0, csq;
+  double h, cut, rho0, csq;
   int zmin_kernel, zmin_surface;
   int kernel_style, surface_style;
   double divr_surface;
-  enum {QUINTIC, CRK0, CRK1, CRK2};
-  enum {COORDINATION, DIVR};
-
-  // Status variables
-  enum {
-    // Phase status
-    STATUS_FLUID = 1 << 0,
-    STATUS_REACTIVE = 1 << 1,
-    STATUS_SOLID = 1 << 2,
-    STATUS_FREEZING = 1 << 3,
-
-    // Surface status
-    STATUS_BULK = 1 << 4,
-    STATUS_LAYER = 1 << 5,
-    STATUS_SURFACE = 1 << 6,
-    STATUS_SPLASH = 1 << 7,
-
-    // Temporary status options - reset in preforce
-    STATUS_SHIFT = 1 << 8,
-    STATUS_NO_FORCE = 1 << 9
-  };
-  int phasemask = 0xFFFFFFF0;
-  int surfacemask = 0xFFFFFF0F;
 
   // Accessory fixes/computes
   int thermal_flag;
@@ -89,6 +66,32 @@ class FixRHEO : public Fix {
   double dtv, dtf;
 };
 
+namespace RHEO_NS {
+
+  enum {QUINTIC, CRK0, CRK1, CRK2};
+  enum {COORDINATION, DIVR};
+
+  // Status variables
+  enum Status{
+    // Phase status
+    STATUS_FLUID = 1 << 0,
+    STATUS_REACTIVE = 1 << 1,
+    STATUS_SOLID = 1 << 2,
+    STATUS_FREEZING = 1 << 3,
+    // Surface status
+    STATUS_BULK = 1 << 4,
+    STATUS_LAYER = 1 << 5,
+    STATUS_SURFACE = 1 << 6,
+    STATUS_SPLASH = 1 << 7,
+    // Temporary status options - reset in preforce
+    STATUS_SHIFT = 1 << 8,
+    STATUS_NO_FORCE = 1 << 9
+  };
+
+  int PHASEMASK = 0xFFFFFFF0;
+  int SURFACEMASK = 0xFFFFFF0F;
+
+}    // namespace RHEO_NS
 }    // namespace LAMMPS_NS
 
 #endif
diff --git a/src/RHEO/fix_rheo_pressure.cpp b/src/RHEO/fix_rheo_pressure.cpp
index 726c44c32b..75edf42572 100644
--- a/src/RHEO/fix_rheo_pressure.cpp
+++ b/src/RHEO/fix_rheo_pressure.cpp
@@ -33,6 +33,8 @@ using namespace LAMMPS_NS;
 using namespace FixConst;
 enum {NONE, LINEAR, CUBIC, TAITWATER};
 
+static constexpr double SEVENTH = 1.0 / 7.0;
+
 /* ---------------------------------------------------------------------- */
 
 FixRHEOPressure::FixRHEOPressure(LAMMPS *lmp, int narg, char **arg) :
@@ -43,7 +45,7 @@ FixRHEOPressure::FixRHEOPressure(LAMMPS *lmp, int narg, char **arg) :
   pressure_style = NONE;
 
   comm_forward = 1;
-  nmax_old = 0;
+  nmax_store = 0;
 
   // Currently can only have one instance of fix rheo/pressure
   if (igroup != 0)
@@ -112,13 +114,13 @@ void FixRHEOPressure::setup_pre_force(int /*vflag*/)
   // Create pressure array if it doesn't already exist
   // Create a custom atom property so it works with compute property/atom
   // Do not create grow callback as there's no reason to copy/exchange data
-  // Manually grow if nmax_old exceeded
+  // Manually grow if nmax_store exceeded
 
   int tmp1, tmp2;
   int index = atom->find_custom("rheo_pressure", tmp1, tmp2);
   if (index == -1) {
     index = atom->add_custom("rheo_pressure", 1, 0);
-    nmax_old = atom->nmax;
+    nmax_store = atom->nmax;
   }
   pressure = atom->dvector[index];
 
@@ -139,13 +141,11 @@ void FixRHEOPressure::pre_force(int /*vflag*/)
 
   int nlocal = atom->nlocal;
 
-  if (nmax_old < atom->nmax) {
+  if (nmax_store < atom->nmax) {
     memory->grow(pressure, atom->nmax, "atom:rheo_pressure");
-    nmax_old = atom->nmax;
+    nmax_store = atom->nmax;
   }
 
-  if (pressure_style == TAITWATER) inv7 = 1.0 / 7.0;
-
   for (i = 0; i < nlocal; i++) {
     if (mask[i] & groupbit) {
       if (pressure_style == LINEAR) {
@@ -156,7 +156,7 @@ void FixRHEOPressure::pre_force(int /*vflag*/)
       } else if (pressure_style == TAITWATER) {
         rho_ratio = rho[i] / rho0inv;
         rr3 = rho_ratio * rho_ratio * rho_ratio;
-        pressure[i] = csq * rho0 * inv7 * (rr3 * rr3 * rho_ratio - 1.0);
+        pressure[i] = csq * rho0 * SEVENTH * (rr3 * rr3 * rho_ratio - 1.0);
       }
     }
   }
@@ -194,9 +194,10 @@ void FixRHEOPressure::unpack_forward_comm(int n, int first, double *buf)
 
 /* ---------------------------------------------------------------------- */
 
-double FixRHEOPressure::calculate_p(double rho)
+double FixRHEOPressure::calc_pressure(double rho)
 {
-  double rho;
+  double p, dr, rr3, rho_ratio;
+
   if (pressure_style == LINEAR) {
     p = csq * (rho - rho0);
   } else if (pressure_style == CUBIC) {
@@ -205,7 +206,7 @@ double FixRHEOPressure::calculate_p(double rho)
   } else if (pressure_style == TAITWATER) {
     rho_ratio = rho / rho0inv;
     rr3 = rho_ratio * rho_ratio * rho_ratio;
-    p = csq * rho0 * inv7 * (rr3 * rr3 * rho_ratio - 1.0);
+    p = csq * rho0 * SEVENTH * (rr3 * rr3 * rho_ratio - 1.0);
   }
   return rho;
 }
@@ -215,6 +216,6 @@ double FixRHEOPressure::calculate_p(double rho)
 double FixRHEOPressure::memory_usage()
 {
   double bytes = 0.0;
-  bytes += (size_t) nmax_old * sizeof(double);
+  bytes += (size_t) nmax_store * sizeof(double);
   return bytes;
 }
diff --git a/src/RHEO/fix_rheo_pressure.h b/src/RHEO/fix_rheo_pressure.h
index 197cab6e5c..c257f1dbfb 100644
--- a/src/RHEO/fix_rheo_pressure.h
+++ b/src/RHEO/fix_rheo_pressure.h
@@ -35,13 +35,14 @@ class FixRHEOPressure : public Fix {
   int pack_forward_comm(int, int *, double *, int, int *) override;
   void unpack_forward_comm(int, int, double *) override;
   double memory_usage() override;
-  double calculate_p(double);
+  double calc_pressure(double);
+
  private:
   double c_cubic, csq, rho0, rho0inv;
   double *pressure;
   int pressure_style;
   int first_flag, last_flag;
-  int nmax_old;
+  int nmax_store;
 
   class FixRHEO *fix_rheo;
 };
diff --git a/src/RHEO/fix_rheo_thermal.cpp b/src/RHEO/fix_rheo_thermal.cpp
index fd08f39fd7..5df8c1c506 100644
--- a/src/RHEO/fix_rheo_thermal.cpp
+++ b/src/RHEO/fix_rheo_thermal.cpp
@@ -22,6 +22,7 @@
 #include "comm.h"
 #include "compute_rheo_grad.h"
 #include "compute_rheo_vshift.h"
+#include "domain.h"
 #include "error.h"
 #include "fix_rheo.h"
 #include "force.h"
@@ -31,14 +32,15 @@
 #include "update.h"
 
 using namespace LAMMPS_NS;
+using namespace RHEO_NS;
 using namespace FixConst;
 enum {NONE, CONSTANT, TYPE};
 
 /* ---------------------------------------------------------------------- */
 
 FixRHEOThermal::FixRHEOThermal(LAMMPS *lmp, int narg, char **arg) :
-  Fix(lmp, narg, arg), Tc_type(nullptr), kappa_type(nullptr), cv_type(nullptr),
-  conductivity(nullptr)
+  Fix(lmp, narg, arg), fix_rheo(nullptr), compute_grad(nullptr), compute_vshift(nullptr),
+  Tc_type(nullptr), kappa_type(nullptr), cv_type(nullptr), conductivity(nullptr)
 {
   if (narg < 4) error->all(FLERR,"Illegal fix command");
 
@@ -47,7 +49,7 @@ FixRHEOThermal::FixRHEOThermal(LAMMPS *lmp, int narg, char **arg) :
   conductivity_style = NONE;
 
   comm_forward = 1;
-  nmax_old = 0;
+  nmax_store = 0;
 
   int ntypes = atom->ntypes;
   int iarg = 3;
@@ -181,13 +183,13 @@ void FixRHEOThermal::setup_pre_force(int /*vflag*/)
 
   // Identify whether this is the first/last instance of fix thermal
   // First will grow arrays, last will communicate
-  first_flag = 0
+  first_flag = 0;
   last_flag = 0;
 
   int i = 0;
   auto fixlist = modify->get_fix_by_style("rheo/thermal");
-  for (const auto &ifix : fixlist) {
-    if (strcmp(ifix->id, id) == 0) break;
+  for (const auto &fix : fixlist) {
+    if (strcmp(fix->id, id) == 0) break;
     i++;
   }
 
@@ -197,13 +199,13 @@ void FixRHEOThermal::setup_pre_force(int /*vflag*/)
   // Create conductivity array if it doesn't already exist
   // Create a custom atom property so it works with compute property/atom
   // Do not create grow callback as there's no reason to copy/exchange data
-  // Manually grow if nmax_old exceeded
+  // Manually grow if nmax_store exceeded
 
   int tmp1, tmp2;
-  index = atom->find_custom("rheo_conductivity", tmp1, tmp2);
+  int index = atom->find_custom("rheo_conductivity", tmp1, tmp2);
   if (index== -1) {
     index = atom->add_custom("rheo_conductivity", 1, 0);
-    nmax_old = atom->nmax;
+    nmax_store = atom->nmax;
   }
   conductivity = atom->dvector[index];
 
@@ -217,13 +219,16 @@ void FixRHEOThermal::initial_integrate(int /*vflag*/)
 {
   // update temperature from shifting
   if (!fix_rheo->shift_flag) return;
-  int i;
+  int i, a;
+
   int *status = atom->status;
+  int *mask = atom->mask;
+  double *temperature = atom->temperature;
   double **gradt = compute_grad->gradt;
   double **vshift = compute_vshift->array_atom;
 
-  int *mask = atom->mask;
   int nlocal = atom->nlocal;
+  int dim = domain->dimension;
 
   if (igroup == atom->firstgroup)
     nlocal = atom->nfirst;
@@ -248,14 +253,14 @@ void FixRHEOThermal::post_integrate()
   double *heatflow = atom->heatflow;
   double *rho = atom->rho;
   int *mask = atom->mask;
-  int *type = aotm->type;
+  int *type = atom->type;
 
   double cvi, Tci, Ti;
 
   //Integrate temperature and check status
   for (int i = 0; i < atom->nlocal; i++) {
     if (mask[i] & groupbit) {
-      if (status[i] == FixRHEO::FLUID_NO_FORCE) continue;
+      if (status[i] & STATUS_NO_FORCE) continue;
 
       cvi = calc_cv(i);
       temperature[i] += dtf * heatflow[i] / cvi;
@@ -265,15 +270,15 @@ void FixRHEOThermal::post_integrate()
         if (Tc_style == CONSTANT) {
           Tci = Tc;
         } else if (Tc_style == TYPE) {
-          Tci = Tc_type[type[i]]);
+          Tci = Tc_type[type[i]];
         }
 
         if (Ti > Tci) {
-          status[i] &= FixRHEO::phasemask;
-          status[i] |= FixRHEO::STATUS_FLUID;
-        } else if (!(status[i] & FixRHEO::STATUS_SOLID))
-          status[i] &= FixRHEO::phasemask;
-          status[i] |= FixRHEO::STATUS_FREEZING;
+          status[i] &= PHASEMASK;
+          status[i] |= STATUS_FLUID;
+        } else if (!(status[i] & STATUS_SOLID)) {
+          status[i] &= PHASEMASK;
+          status[i] |= STATUS_FREEZING;
         }
       }
     }
@@ -288,14 +293,13 @@ void FixRHEOThermal::post_neighbor()
 {
   int i;
   int *type = atom->type;
-  double *conductivity = atom->dvector[index_cond];
   int *mask = atom->mask;
   int nlocal = atom->nlocal;
   int nall = nlocal + atom->nghost;
 
-  if (first_flag && (nmax_old < atom->nmax)) {
+  if (first_flag && (nmax_store < atom->nmax)) {
     memory->grow(conductivity, atom->nmax, "atom:rheo_conductivity");
-    nmax_old = atom->nmax;
+    nmax_store = atom->nmax;
   }
 
   if (conductivity_style == CONSTANT) {
@@ -304,7 +308,6 @@ void FixRHEOThermal::post_neighbor()
   } else if (conductivity_style == TYPE) {
     for (i = 0; i < nall; i++)
       if (mask[i] & groupbit) conductivity[i] = kappa_type[type[i]];
-    }
   }
 }
 
@@ -329,9 +332,9 @@ void FixRHEOThermal::pre_force(int /*vflag*/)
   //int *mask = atom->mask;
   //int nlocal = atom->nlocal;
 
-  //if (first_flag && (nmax_old < atom->nmax)) {
+  //if (first_flag && (nmax_store < atom->nmax)) {
   //  memory->grow(conductivity, atom->nmax, "atom:rheo_conductivity");
-  //  nmax_old = atom->nmax;
+  //  nmax_store = atom->nmax;
   //}
 
   //if (conductivity_style == TBD) {
@@ -358,7 +361,7 @@ void FixRHEOThermal::final_integrate()
   //Integrate temperature and check status
   for (int i = 0; i < atom->nlocal; i++) {
     if (mask[i] & groupbit) {
-      if (status[i] & FixRHEO::STATUS_NO_FORCE) continue;
+      if (status[i] & STATUS_NO_FORCE) continue;
 
       cvi = calc_cv(i);
       temperature[i] += dtf * heatflow[i] / cvi;
@@ -447,6 +450,6 @@ void FixRHEOThermal::unpack_reverse_comm(int n, int *list, double *buf)
 double FixRHEOThermal::memory_usage()
 {
   double bytes = 0.0;
-  bytes += (size_t) nmax_old * sizeof(double);
+  bytes += (size_t) nmax_store * sizeof(double);
   return bytes;
 }
diff --git a/src/RHEO/fix_rheo_thermal.h b/src/RHEO/fix_rheo_thermal.h
index 4f0e89f17c..cf64c0b8d1 100644
--- a/src/RHEO/fix_rheo_thermal.h
+++ b/src/RHEO/fix_rheo_thermal.h
@@ -31,7 +31,7 @@ class FixRHEOThermal : public Fix {
   int setmask() override;
   void init() override;
   void setup_pre_force(int) override;
-  void initial_integrate() override;
+  void initial_integrate(int) override;
   void post_integrate() override;
   void post_neighbor() override;
   void pre_force(int) override;
@@ -53,9 +53,11 @@ class FixRHEOThermal : public Fix {
   int cv_style;
   int conductivity_style;
   int first_flag, last_flag;
-  int nmax_old;
+  int nmax_store;
 
   class FixRHEO *fix_rheo;
+  class ComputeRHEOGrad *compute_grad;
+  class ComputeRHEOVShift *compute_vshift;
 
   double calc_cv(int);
 };
diff --git a/src/RHEO/fix_rheo_viscosity.cpp b/src/RHEO/fix_rheo_viscosity.cpp
index 5ae1b95529..b15f488370 100644
--- a/src/RHEO/fix_rheo_viscosity.cpp
+++ b/src/RHEO/fix_rheo_viscosity.cpp
@@ -44,7 +44,7 @@ FixRHEOViscosity::FixRHEOViscosity(LAMMPS *lmp, int narg, char **arg) :
   viscosity_style = NONE;
 
   comm_forward = 0;
-  nmax_old = 0;
+  nmax_store = 0;
 
   int ntypes = atom->ntypes;
   int iarg = 3;
@@ -89,7 +89,7 @@ FixRHEOViscosity::~FixRHEOViscosity()
   // Remove custom property if it exists
   int tmp1, tmp2, index;
   index = atom->find_custom("rheo_viscosity", tmp1, tmp2);
-  if (index != -1) atom->remove_custom(index_visc, 1, 0);
+  if (index != -1) atom->remove_custom(index, 1, 0);
 
   memory->destroy(eta_type);
 }
@@ -123,13 +123,13 @@ void FixRHEOViscosity::setup_pre_force(int /*vflag*/)
 
   // Identify whether this is the first/last instance of fix viscosity
   // First will grow arrays, last will communicate
-  first_flag = 0
+  first_flag = 0;
   last_flag = 0;
 
   int i = 0;
   auto fixlist = modify->get_fix_by_style("rheo/viscosity");
-  for (const auto &ifix : fixlist) {
-    if (strcmp(ifix->id, id) == 0) break;
+  for (const auto &fix : fixlist) {
+    if (strcmp(fix->id, id) == 0) break;
     i++;
   }
 
@@ -139,13 +139,13 @@ void FixRHEOViscosity::setup_pre_force(int /*vflag*/)
   // Create viscosity array if it doesn't already exist
   // Create a custom atom property so it works with compute property/atom
   // Do not create grow callback as there's no reason to copy/exchange data
-  // Manually grow if nmax_old exceeded
+  // Manually grow if nmax_store exceeded
 
   int tmp1, tmp2;
   int index = atom->find_custom("rheo_viscosity", tmp1, tmp2);
-  if (index_visc == -1) {
+  if (index == -1) {
     index = atom->add_custom("rheo_viscosity", 1, 0);
-    nmax_old = atom->nmax;
+    nmax_store = atom->nmax;
   }
   viscosity = atom->dvector[index];
 
@@ -167,9 +167,9 @@ void FixRHEOViscosity::post_neighbor()
   int nlocal = atom->nlocal;
   int nall = nlocal + atom->nghost;
 
-  if (first_flag && (nmax_old < atom->nmax)) {
+  if (first_flag && (nmax_store < atom->nmax)) {
     memory->grow(viscosity, atom->nmax, "atom:rheo_viscosity");
-    nmax_old = atom->nmax;
+    nmax_store = atom->nmax;
   }
 
   if (viscosity_style == CONSTANT) {
@@ -178,7 +178,6 @@ void FixRHEOViscosity::post_neighbor()
   } else if (viscosity_style == TYPE) {
     for (i = 0; i < nall; i++)
       if (mask[i] & groupbit) viscosity[i] = eta_type[type[i]];
-    }
   }
 }
 
@@ -197,9 +196,9 @@ void FixRHEOViscosity::pre_force(int /*vflag*/)
   int nlocal = atom->nlocal;
   int dim = domain->dimension;
 
-  if (first_flag && (nmax_old < atom->nmax)) {
+  if (first_flag && (nmax_store < atom->nmax)) {
     memory->grow(viscosity, atom->nmax, "atom:rheo_viscosity");
-    nmax_old = atom->nmax;
+    nmax_store = atom->nmax;
   }
 
   if (viscosity_style == POWER) {
@@ -260,6 +259,6 @@ void FixRHEOViscosity::unpack_forward_comm(int n, int first, double *buf)
 double FixRHEOViscosity::memory_usage()
 {
   double bytes = 0.0;
-  bytes += (size_t) nmax_old * sizeof(double);
+  bytes += (size_t) nmax_store * sizeof(double);
   return bytes;
 }
diff --git a/src/RHEO/fix_rheo_viscosity.h b/src/RHEO/fix_rheo_viscosity.h
index 14f8b70de9..66df51601e 100644
--- a/src/RHEO/fix_rheo_viscosity.h
+++ b/src/RHEO/fix_rheo_viscosity.h
@@ -43,7 +43,7 @@ class FixRHEOViscosity : public Fix {
   double *viscosity;
   int viscosity_style;
   int first_flag, last_flag;
-  int nmax_old;
+  int nmax_store;
 
   class FixRHEO *fix_rheo;
   class ComputeRHEOGrad *compute_grad;
diff --git a/src/RHEO/pair_rheo.cpp b/src/RHEO/pair_rheo.cpp
index 5974b9b756..c61d613d82 100644
--- a/src/RHEO/pair_rheo.cpp
+++ b/src/RHEO/pair_rheo.cpp
@@ -39,9 +39,10 @@
 #include <cmath>
 
 using namespace LAMMPS_NS;
+using namespace RHEO_NS;
 using namespace MathExtra;
 
-#define EPSILON 1e-2
+static constexpr double EPSILON = 1e-2;
 
 /* ---------------------------------------------------------------------- */
 
@@ -83,6 +84,10 @@ void PairRHEO::compute(int eflag, int vflag)
   int *ilist, *jlist, *numneigh, **firstneigh;
   double imass, jmass, rsq, r, rinv;
 
+  int nlocal = atom->nlocal;
+  int newton_pair = force->newton_pair;
+  int dim = domain->dimension;
+
   ev_init(eflag, vflag);
 
   double **gradv = compute_grad->gradv;
@@ -123,11 +128,6 @@ void PairRHEO::compute(int eflag, int vflag)
     conductivity = atom->dvector[index];
   }
 
-  int *ilist, *jlist, *numneigh, **firstneigh;
-  int nlocal = atom->nlocal;
-  int newton_pair = force->newton_pair;
-  int dim = domain->dimension;
-
   inum = list->inum;
   ilist = list->ilist;
   numneigh = list->numneigh;
@@ -145,7 +145,7 @@ void PairRHEO::compute(int eflag, int vflag)
     jnum = numneigh[i];
     imass = mass[itype];
     etai = viscosity[i];
-    fluidi = status[i] & FixRHEO::STATUS_FLUID;
+    fluidi = status[i] & STATUS_FLUID;
     if (thermal_flag) {
       kappai = conductivity[i];
       Ti = temperature[i];
@@ -167,7 +167,7 @@ void PairRHEO::compute(int eflag, int vflag)
 
         jmass = mass[jtype];
         etaj = viscosity[j];
-        fluidj = status[j] & FixRHEO::STATUS_FLUID;
+        fluidj = status[j] & STATUS_FLUID;
         if (thermal_flag) {
           Tj = temperature[j];
           kappaj = conductivity[j];
@@ -202,7 +202,7 @@ void PairRHEO::compute(int eflag, int vflag)
         if (fluidi && (!fluidj)) {
           compute_interface->correct_v(vi, vj, i, j);
           rhoj = compute_interface->correct_rho(j, i);
-          Pj = fix_pressure->calculate_p(rhoj);
+          Pj = fix_pressure->calc_pressure(rhoj);
 
           if ((chi[j] > 0.9) && (r < (h * 0.5)))
             fmag = (chi[j] - 0.9) * (h * 0.5 - r) * rho0 * csq * h * rinv;
@@ -210,9 +210,9 @@ void PairRHEO::compute(int eflag, int vflag)
         } else if ((!fluidi) && fluidj) {
           compute_interface->correct_v(vj, vi, j, i);
           rhoi = compute_interface->correct_rho(i, j);
-          Pi = calc_pressure(rhoi, itype);
+          Pi = fix_pressure->calc_pressure(rhoi);
 
-          if (chi[i] > 0.9 && r < (h * 0.5)) {
+          if (chi[i] > 0.9 && r < (h * 0.5))
             fmag = (chi[i] - 0.9) * (h * 0.5 - r) * rho0 * csq * h * rinv;
 
         } else if ((!fluidi) && (!fluidj)) {
@@ -244,7 +244,7 @@ void PairRHEO::compute(int eflag, int vflag)
 
           //Hydrostatic pressure forces
           fp_prefactor = voli * volj * (Pj + Pi);
-          sub3(v1, vj, du);
+          sub3(vi, vj, du);
 
           //Add artificial viscous pressure if required
           if (artificial_visc_flag && pair_avisc_flag){
@@ -423,15 +423,11 @@ void PairRHEO::setup()
   if (fixes.size() == 0) error->all(FLERR, "Need to define fix rheo to use pair rheo");
   fix_rheo = dynamic_cast<FixRHEO *>(fixes[0]);
 
+  // Currently only allow one instance of fix rheo/pressure
   fixes = modify->get_fix_by_style("rheo/pressure");
   if (fixes.size() == 0) error->all(FLERR, "Need to define fix rheo/pressure to use pair rheo");
   fix_pressure = dynamic_cast<FixRHEOPressure *>(fixes[0]);
 
-  int tmp1, tmp2;
-  index_pressure = atom->find_custom("rheo_pressure", tmp1, tmp2);
-  if (index_pressure == -1) index_pressure = atom->add_custom("rheo_pressure", 1, 0);
-  else error->all(FLERR, "Cannot find pressure value in pair rheo");
-
   compute_kernel = fix_rheo->compute_kernel;
   compute_grad = fix_rheo->compute_grad;
   compute_interface = fix_rheo->compute_interface;
@@ -449,9 +445,9 @@ void PairRHEO::setup()
     error->all(FLERR,"Pair RHEO requires ghost atoms store velocity");
 
   if (laplacian_order == -1) {
-    if (fix_rheo->kernel_type == FixRHEO::CRK2)
+    if (fix_rheo->kernel_style == CRK2)
       laplacian_order = 2;
-    else if (fix_rheo->kernel_type == FixRHEO::CRK1)
+    else if (fix_rheo->kernel_style == CRK1)
       laplacian_order = 1;
     else
       laplacian_order = 0;
@@ -468,8 +464,5 @@ double PairRHEO::init_one(int i, int j)
       error->all(FLERR,"All pair rheo coeffs are not set");
   }
 
-  cut[i][j] = h;
-  cut[j][i] = cut[i][j];
-
-  return cut[i][j];
+  return h;
 }