diff --git a/src/RHEO/bond_rheo_shell.cpp b/src/RHEO/bond_rheo_shell.cpp
index 258d047086..dbe65e2b2f 100644
--- a/src/RHEO/bond_rheo_shell.cpp
+++ b/src/RHEO/bond_rheo_shell.cpp
@@ -45,7 +45,8 @@ using namespace RHEO_NS;
 /* ---------------------------------------------------------------------- */
 
 BondRHEOShell::BondRHEOShell(LAMMPS *_lmp) :
-    BondBPM(_lmp), compute_surface(nullptr), k(nullptr), ecrit(nullptr), gamma(nullptr)
+  BondBPM(_lmp), k(nullptr), ecrit(nullptr), gamma(nullptr), dbond(nullptr), nbond(nullptr),
+  id_fix(nullptr), compute_surface(nullptr)
 {
   partial_flag = 1;
   comm_reverse = 1;
diff --git a/src/RHEO/compute_rheo_grad.cpp b/src/RHEO/compute_rheo_grad.cpp
index cdf90e1dc5..fb48d90d56 100644
--- a/src/RHEO/compute_rheo_grad.cpp
+++ b/src/RHEO/compute_rheo_grad.cpp
@@ -42,8 +42,9 @@ enum{COMMGRAD, COMMFIELD};
 /* ---------------------------------------------------------------------- */
 
 ComputeRHEOGrad::ComputeRHEOGrad(LAMMPS *lmp, int narg, char **arg) :
-  Compute(lmp, narg, arg), fix_rheo(nullptr), list(nullptr), rho0(nullptr), compute_interface(nullptr), compute_kernel(nullptr),
-  gradv(nullptr), gradr(nullptr), grade(nullptr), gradn(nullptr)
+  Compute(lmp, narg, arg), gradv(nullptr), gradr(nullptr), grade(nullptr), gradn(nullptr),
+  fix_rheo(nullptr), rho0(nullptr), compute_kernel(nullptr), compute_interface(nullptr),
+  list(nullptr)
 {
   if (narg < 4) error->all(FLERR,"Illegal compute rheo/grad command");
 
@@ -130,10 +131,9 @@ void ComputeRHEOGrad::compute_peratom()
 {
   int i, j, k, ii, jj, jnum, itype, jtype, a, b, fluidi, fluidj;
   double xtmp, ytmp, ztmp, delx, dely, delz;
-  double rsq, imass, jmass;
-  double rhoi, rhoj, Voli, Volj, drho, de, deta;
+  double rsq, rhoi, rhoj, Voli, Volj, drho, de, deta;
   double vi[3], vj[3], vij[3];
-  double wp, *dWij, *dWji;
+  double *dWij, *dWji;
 
   int inum, *ilist, *numneigh, **firstneigh;
   int *jlist;
@@ -236,7 +236,7 @@ void ComputeRHEOGrad::compute_peratom()
         if (energy_flag) de = energy[i] - energy[j];
         if (eta_flag) deta = viscosity[i] - viscosity[j];
 
-        wp = compute_kernel->calc_dw(i, j, delx, dely, delz, sqrt(rsq));
+        compute_kernel->calc_dw(i, j, delx, dely, delz, sqrt(rsq));
         dWij = compute_kernel->dWij;
         dWji = compute_kernel->dWji;
 
diff --git a/src/RHEO/compute_rheo_interface.cpp b/src/RHEO/compute_rheo_interface.cpp
index 8ccd4e6a3b..05821ebe22 100644
--- a/src/RHEO/compute_rheo_interface.cpp
+++ b/src/RHEO/compute_rheo_interface.cpp
@@ -45,8 +45,10 @@ static constexpr double EPSILON = 1e-1;
 /* ---------------------------------------------------------------------- */
 
 ComputeRHEOInterface::ComputeRHEOInterface(LAMMPS *lmp, int narg, char **arg) :
-  Compute(lmp, narg, arg), fix_rheo(nullptr), compute_kernel(nullptr), fp_store(nullptr),
-  rho0(nullptr), norm(nullptr), normwf(nullptr), chi(nullptr), id_fix_pa(nullptr)
+  Compute(lmp, narg, arg), chi(nullptr), fp_store(nullptr), fix_rheo(nullptr),
+  rho0(nullptr), norm(nullptr), normwf(nullptr), id_fix_pa(nullptr), list(nullptr),
+  compute_kernel(nullptr), fix_pressure(nullptr)
+
 {
   if (narg != 3) error->all(FLERR,"Illegal compute rheo/interface command");
 
@@ -228,12 +230,11 @@ void ComputeRHEOInterface::compute_peratom()
 
 int ComputeRHEOInterface::pack_forward_comm(int n, int *list, double *buf, int /*pbc_flag*/, int * /*pbc*/)
 {
-  int i, j, k, m;
-  m = 0;
+  int m = 0;
   double *rho = atom->rho;
 
-  for (i = 0; i < n; i++) {
-    j = list[i];
+  for (int i = 0; i < n; i++) {
+    int j = list[i];
     if (comm_stage == 0) {
       buf[m++] = fp_store[j][0];
       buf[m++] = fp_store[j][1];
@@ -250,11 +251,10 @@ int ComputeRHEOInterface::pack_forward_comm(int n, int *list, double *buf, int /
 
 void ComputeRHEOInterface::unpack_forward_comm(int n, int first, double *buf)
 {
-  int i, k, m, last;
   double *rho = atom->rho;
-  m = 0;
-  last = first + n;
-  for (i = first; i < last; i++) {
+  int m = 0;
+  int last = first + n;
+  for (int i = first; i < last; i++) {
     if (comm_stage == 0) {
       fp_store[i][0] = buf[m++];
       fp_store[i][1] = buf[m++];
@@ -270,12 +270,10 @@ void ComputeRHEOInterface::unpack_forward_comm(int n, int first, double *buf)
 
 int ComputeRHEOInterface::pack_reverse_comm(int n, int first, double *buf)
 {
-  int i, k, m, last;
   double *rho = atom->rho;
-
-  m = 0;
-  last = first + n;
-  for (i = first; i < last; i++) {
+  int m = 0;
+  int last = first + n;
+  for (int i = first; i < last; i++) {
     buf[m++] = norm[i];
     buf[m++] = chi[i];
     buf[m++] = normwf[i];
@@ -288,12 +286,11 @@ int ComputeRHEOInterface::pack_reverse_comm(int n, int first, double *buf)
 
 void ComputeRHEOInterface::unpack_reverse_comm(int n, int *list, double *buf)
 {
-  int i, k, j, m;
   double *rho = atom->rho;
   int *status = atom->rheo_status;
-  m = 0;
-  for (i = 0; i < n; i++) {
-    j = list[i];
+  int m = 0;
+  for (int i = 0; i < n; i++) {
+    int j = list[i];
     norm[j] += buf[m++];
     chi[j] += buf[m++];
     if (status[j] & PHASECHECK){
diff --git a/src/RHEO/compute_rheo_kernel.cpp b/src/RHEO/compute_rheo_kernel.cpp
index 4558ddccc8..1913bb1e53 100644
--- a/src/RHEO/compute_rheo_kernel.cpp
+++ b/src/RHEO/compute_rheo_kernel.cpp
@@ -53,8 +53,8 @@ static constexpr int DELTA = 2000;
 /* ---------------------------------------------------------------------- */
 
 ComputeRHEOKernel::ComputeRHEOKernel(LAMMPS *lmp, int narg, char **arg) :
-  Compute(lmp, narg, arg),
-  list(nullptr), C(nullptr), C0(nullptr), coordination(nullptr), compute_interface(nullptr)
+  Compute(lmp, narg, arg), coordination(nullptr), fix_rheo(nullptr), C(nullptr), C0(nullptr),
+  list(nullptr), compute_interface(nullptr)
 {
   if (narg != 4) error->all(FLERR,"Illegal compute rheo/kernel command");
 
@@ -190,7 +190,7 @@ double ComputeRHEOKernel::calc_w_self(int i, int j)
 
 double ComputeRHEOKernel::calc_w(int i, int j, double delx, double dely, double delz, double r)
 {
-  double w;
+  double w = 0.0;
   int corrections_i, corrections_j, corrections;
 
   if (kernel_style == WENDLANDC4)
@@ -282,8 +282,6 @@ double ComputeRHEOKernel::calc_w_quintic(int i, int j, double delx, double dely,
 double ComputeRHEOKernel::calc_dw_quintic(int i, int j, double delx, double dely, double delz, double r, double *dW1, double *dW2)
 {
   double wp, tmp1, tmp2, tmp3, tmp1sq, tmp2sq, tmp3sq, s, wprinv;
-  double *mass = atom->mass;
-  int *type = atom->type;
 
   s = r * 3.0 * cutinv;
 
@@ -326,9 +324,9 @@ double ComputeRHEOKernel::calc_w_wendlandc4(int i, int j, double delx, double de
   double w, tmp6, s;
   s = r * cutinv;
 
-        if (s > 1.0) {
-          w = 0.0;
-        } else {
+  if (s > 1.0) {
+    w = 0.0;
+  } else {
     tmp6 = (1.0 - s) * (1.0 - s);
     tmp6 *= tmp6 * tmp6;
     w = tmp6 * (1.0 + 6.0 * s + 35.0 * THIRD * s * s);
@@ -347,14 +345,12 @@ double ComputeRHEOKernel::calc_w_wendlandc4(int i, int j, double delx, double de
 double ComputeRHEOKernel::calc_dw_wendlandc4(int i, int j, double delx, double dely, double delz, double r, double *dW1, double *dW2)
 {
   double wp, tmp1, tmp5, tmp6, s, wprinv;
-  double *mass = atom->mass;
-  int *type = atom->type;
 
   s = r * cutinv;
 
-        if (s > 1.0) {
-          wp = 0.0;
-        } else {
+  if (s > 1.0) {
+    wp = 0.0;
+  } else {
     tmp1 = 1.0 - s;
     tmp5 = tmp1 * tmp1;
     tmp5 = tmp5 * tmp5 * tmp1;
@@ -395,7 +391,7 @@ double ComputeRHEOKernel::calc_w_rk0(int i, int j, double delx, double dely, dou
 double ComputeRHEOKernel::calc_w_rk1(int i, int j, double delx, double dely, double delz, double r)
 {
   int b;
-  double w, wR, dx[3], H[Mdim];
+  double w, dx[3], H[Mdim];
 
   dx[0] = delx;
   dx[1] = dely;
@@ -437,7 +433,7 @@ double ComputeRHEOKernel::calc_w_rk1(int i, int j, double delx, double dely, dou
 double ComputeRHEOKernel::calc_w_rk2(int i, int j, double delx, double dely, double delz, double r)
 {
   int b;
-  double w, wR, dx[3], H[Mdim];
+  double w, dx[3], H[Mdim];
   dx[0] = delx;
   dx[1] = dely;
   dx[2] = delz;
@@ -574,7 +570,7 @@ void ComputeRHEOKernel::compute_peratom()
   if (kernel_style == QUINTIC) return;
   corrections_calculated = 1;
 
-  int i, j, ii, jj, inum, jnum, itype, g, a, b, gsl_error;
+  int i, j, ii, jj, inum, jnum, a, b, gsl_error;
   double xtmp, ytmp, ztmp, r, rsq, w, vj, rhoj;
   double dx[3];
   gsl_matrix_view gM;
@@ -854,19 +850,17 @@ void ComputeRHEOKernel::grow_arrays(int nmax)
 int ComputeRHEOKernel::pack_forward_comm(int n, int *list, double *buf,
                                         int /*pbc_flag*/, int * /*pbc*/)
 {
-  int i,j,k,m,a,b;
-  m = 0;
-
-  for (i = 0; i < n; i++) {
-    j = list[i];
+  int m = 0;
+  for (int i = 0; i < n; i++) {
+    int j = list[i];
     if (comm_stage == 0) {
       buf[m++] = coordination[j];
     } else {
       if (kernel_style == RK0) {
         buf[m++] = C0[j];
       } else {
-        for (a = 0; a < ncor; a++)
-          for (b = 0; b < Mdim; b++)
+        for (int a = 0; a < ncor; a++)
+          for (int b = 0; b < Mdim; b++)
             buf[m++] = C[j][a][b];
       }
     }
@@ -878,19 +872,17 @@ int ComputeRHEOKernel::pack_forward_comm(int n, int *list, double *buf,
 
 void ComputeRHEOKernel::unpack_forward_comm(int n, int first, double *buf)
 {
-  int i, k, m, last,a,b;
-  m = 0;
-  last = first + n;
-
-  for (i = first; i < last; i++) {
+  int m = 0;
+  int last = first + n;
+  for (int i = first; i < last; i++) {
     if (comm_stage == 0) {
       coordination[i] = buf[m++];
     } else {
       if (kernel_style == RK0) {
         C0[i] = buf[m++];
       } else {
-        for (a = 0; a < ncor; a++)
-          for (b = 0; b < Mdim; b++)
+        for (int a = 0; a < ncor; a++)
+          for (int b = 0; b < Mdim; b++)
             C[i][a][b] = buf[m++];
       }
     }
diff --git a/src/RHEO/compute_rheo_property_atom.cpp b/src/RHEO/compute_rheo_property_atom.cpp
index 7a450e7708..cd0ff6692d 100644
--- a/src/RHEO/compute_rheo_property_atom.cpp
+++ b/src/RHEO/compute_rheo_property_atom.cpp
@@ -46,9 +46,10 @@ using namespace RHEO_NS;
 /* ---------------------------------------------------------------------- */
 
 ComputeRHEOPropertyAtom::ComputeRHEOPropertyAtom(LAMMPS *lmp, int narg, char **arg) :
-  Compute(lmp, narg, arg), fix_rheo(nullptr), fix_pressure(nullptr), fix_thermal(nullptr),  compute_interface(nullptr),
-  compute_kernel(nullptr), compute_surface(nullptr), compute_vshift(nullptr), compute_grad(nullptr),
-  avec_index(nullptr), pack_choice(nullptr), col_index(nullptr)
+  Compute(lmp, narg, arg), avec_index(nullptr), col_index(nullptr), col_t_index(nullptr),
+  buf(nullptr), pack_choice(nullptr), fix_rheo(nullptr), fix_pressure(nullptr),
+  fix_thermal(nullptr),  compute_interface(nullptr), compute_kernel(nullptr),
+  compute_surface(nullptr), compute_vshift(nullptr), compute_grad(nullptr)
 {
   if (narg < 4)  utils::missing_cmd_args(FLERR, "compute property/atom", error);
 
@@ -88,7 +89,6 @@ ComputeRHEOPropertyAtom::ComputeRHEOPropertyAtom(LAMMPS *lmp, int narg, char **a
   col_t_index = new int[nvalues];
 
   int i = 0;
-  int index, a, b;
   for (int iarg = 3; iarg < narg; iarg++) {
     if (strcmp(arg[iarg], "phase") == 0) {
       pack_choice[i] = &ComputeRHEOPropertyAtom::pack_phase;
@@ -536,7 +536,7 @@ int ComputeRHEOPropertyAtom::add_tensor_component(char* option, int i, FnPtrPack
     }
     shift = dim * dim;
   } else {
-    int index;
+    int index = -1;
     int dim_error = 0;
 
     if (utils::strmatch(option, "xx$")) {
@@ -592,7 +592,7 @@ int ComputeRHEOPropertyAtom::add_vector_component(char* option, int i, FnPtrPack
     }
     shift = dim;
   } else {
-    int index;
+    int index = -1;
     if (utils::strmatch(option, "x$")) {
       index = 0;
     } else if (utils::strmatch(option, "y$")) {
diff --git a/src/RHEO/compute_rheo_rho_sum.cpp b/src/RHEO/compute_rheo_rho_sum.cpp
index 6e25b35374..8067bae21b 100644
--- a/src/RHEO/compute_rheo_rho_sum.cpp
+++ b/src/RHEO/compute_rheo_rho_sum.cpp
@@ -72,7 +72,7 @@ void ComputeRHEORhoSum::init_list(int /*id*/, NeighList *ptr)
 
 void ComputeRHEORhoSum::compute_peratom()
 {
-  int i, j, ii, jj, inum, jnum, itype, jtype;
+  int i, j, ii, jj, inum, jnum;
   double xtmp, ytmp, ztmp, delx, dely, delz;
   int *ilist, *jlist, *numneigh, **firstneigh;
   double rsq, w;
@@ -81,12 +81,11 @@ void ComputeRHEORhoSum::compute_peratom()
 
   double **x = atom->x;
   double *rho = atom->rho;
-  int *type = atom->type;
   double *mass = atom->mass;
+  double *rmass = atom->rmass;
+  int *type = atom->type;
   int newton = force->newton;
 
-  double jmass;
-
   inum = list->inum;
   ilist = list->ilist;
   numneigh = list->numneigh;
@@ -96,7 +95,8 @@ void ComputeRHEORhoSum::compute_peratom()
   // initialize arrays, local with quintic self-contribution, ghosts are zeroed
   for (i = 0; i < nlocal; i++) {
     w = compute_kernel->calc_w_self(i, i);
-    rho[i] = w * mass[type[i]];
+    if (rmass) rho[i] = w * rmass[i];
+    else rho[i] = w * mass[type[i]];
   }
 
   for (i = nlocal; i < nall; i++) rho[i] = 0.0;
@@ -106,7 +106,6 @@ void ComputeRHEORhoSum::compute_peratom()
     xtmp = x[i][0];
     ytmp = x[i][1];
     ztmp = x[i][2];
-    itype = type[i];
     jlist = firstneigh[i];
     jnum = numneigh[i];
 
@@ -121,14 +120,26 @@ void ComputeRHEORhoSum::compute_peratom()
       if (rsq < cutsq) {
         w = compute_kernel->calc_w(i, j, delx, dely, delz, sqrt(rsq));
 
-        if (self_mass_flag) {
-          rho[i] += w * mass[type[i]];
-          if (newton || j < nlocal)
-            rho[j] += w * mass[type[j]];
+        if (rmass) {
+          if (self_mass_flag) {
+            rho[i] += w * rmass[i];
+            if (newton || j < nlocal)
+              rho[j] += w * rmass[j];
+          } else {
+            rho[i] += w * rmass[j];
+            if (newton || j < nlocal)
+              rho[j] += w * rmass[i];
+          }
         } else {
-          rho[i] += w * mass[type[j]];
-          if (newton || j < nlocal)
-            rho[j] += w * mass[type[i]];
+          if (self_mass_flag) {
+            rho[i] += w * mass[type[i]];
+            if (newton || j < nlocal)
+              rho[j] += w * mass[type[j]];
+          } else {
+            rho[i] += w * mass[type[j]];
+            if (newton || j < nlocal)
+              rho[j] += w * mass[type[i]];
+          }
         }
       }
     }
@@ -143,7 +154,7 @@ 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;
+  int i, j, m;
   double *rho = atom->rho;
   m = 0;
 
@@ -157,7 +168,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;
+  int i, m, last;
   double *rho = atom->rho;
 
   m = 0;
@@ -171,7 +182,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, m, last;
   double *rho = atom->rho;
 
   m = 0;
@@ -186,7 +197,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, 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 c3a3774cdc..83ef659e84 100644
--- a/src/RHEO/compute_rheo_surface.cpp
+++ b/src/RHEO/compute_rheo_surface.cpp
@@ -43,8 +43,9 @@ static constexpr double EPSILON = 1e-10;
 /* ---------------------------------------------------------------------- */
 
 ComputeRHEOSurface::ComputeRHEOSurface(LAMMPS *lmp, int narg, char **arg) :
-  Compute(lmp, narg, arg), fix_rheo(nullptr), list(nullptr), rho0(nullptr), compute_kernel(nullptr), compute_interface(nullptr),
-  B(nullptr), gradC(nullptr), nsurface(nullptr), divr(nullptr), rsurface(nullptr)
+  Compute(lmp, narg, arg), nsurface(nullptr), rsurface(nullptr), divr(nullptr), fix_rheo(nullptr),
+  rho0(nullptr), B(nullptr), gradC(nullptr), list(nullptr), compute_kernel(nullptr),
+  compute_interface(nullptr)
 {
   if (narg != 3) error->all(FLERR,"Illegal compute RHEO/SURFACE command");
 
@@ -98,8 +99,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, wp;
+  int i, j, ii, jj, inum, jnum, a, itype, jtype, fluidi, fluidj;
+  double xtmp, ytmp, ztmp, rsq, Voli, Volj, rhoi, rhoj;
   double dWij[3], dWji[3], dx[3];
   int *ilist, *jlist, *numneigh, **firstneigh;
 
@@ -175,7 +176,7 @@ void ComputeRHEOSurface::compute_peratom()
         Voli = mass[itype] / rhoi;
         Volj = mass[jtype] / rhoj;
 
-        wp = compute_kernel->calc_dw_quintic(i, j, dx[0], dx[1], dx[2], sqrt(rsq), dWij, dWji);
+        compute_kernel->calc_dw_quintic(i, j, dx[0], dx[1], dx[2], sqrt(rsq), dWij, dWji);
 
         for (a = 0; a < dim; a++) {
           divr[i] -= dWij[a] * dx[a] * Volj;
@@ -310,17 +311,15 @@ void ComputeRHEOSurface::compute_peratom()
 
 int ComputeRHEOSurface::pack_reverse_comm(int n, int first, double *buf)
 {
-  int i,a,b,k,m,last;
   int dim = domain->dimension;
   int *status = atom->rheo_status;
-
-  m = 0;
-  last = first + n;
-  for (i = first; i < last; i++) {
+  int m = 0;
+  int last = first + n;
+  for (int i = first; i < last; i++) {
     if (comm_stage == 0) {
       buf[m++] = divr[i];
-      for (a = 0; a < dim; a ++ )
-        for (b = 0; b < dim; b ++)
+      for (int a = 0; a < dim; a ++ )
+        for (int b = 0; b < dim; b ++)
           buf[m++] = gradC[i][a * dim + b];
     } else if (comm_stage == 1) {
       buf[m++] = (double) status[i];
@@ -334,27 +333,23 @@ int ComputeRHEOSurface::pack_reverse_comm(int n, int first, double *buf)
 
 void ComputeRHEOSurface::unpack_reverse_comm(int n, int *list, double *buf)
 {
-  int i,a,b,k,j,m;
   int dim = domain->dimension;
   int *status = atom->rheo_status;
-  int tmp1;
-  double tmp2;
-
-  m = 0;
-  for (i = 0; i < n; i++) {
-    j = list[i];
+  int m = 0;
+  for (int i = 0; i < n; i++) {
+    int j = list[i];
     if (comm_stage == 0) {
       divr[j] += buf[m++];
-      for (a = 0; a < dim; a ++ )
-        for (b = 0; b < dim; b ++)
+      for (int a = 0; a < dim; a ++ )
+        for (int b = 0; b < dim; b ++)
           gradC[j][a * dim + b] += buf[m++];
     } else if (comm_stage == 1) {
-      tmp1 = (int) buf[m++];
+      auto tmp1 = (int) buf[m++];
       if ((status[j] & STATUS_BULK) && (tmp1 & STATUS_LAYER)) {
         status[j] &= SURFACEMASK;
         status[j] |= STATUS_LAYER;
       }
-      tmp2 = buf[m++];
+      auto tmp2 = buf[m++];
       rsurface[j] = MIN(rsurface[j], tmp2);
     }
   }
@@ -366,12 +361,10 @@ void ComputeRHEOSurface::unpack_reverse_comm(int n, int *list, double *buf)
 int ComputeRHEOSurface::pack_forward_comm(int n, int *list, double *buf,
                                         int /*pbc_flag*/, int * /*pbc*/)
 {
-  int i,j,a,b,k,m;
   int *status = atom->rheo_status;
-  m = 0;
-
-  for (i = 0; i < n; i++) {
-    j = list[i];
+  int m = 0;
+  for (int i = 0; i < n; i++) {
+    int j = list[i];
     if (comm_stage == 0) {
       buf[m++] = divr[j];
     } else if (comm_stage == 1) {
@@ -386,12 +379,10 @@ int ComputeRHEOSurface::pack_forward_comm(int n, int *list, double *buf,
 
 void ComputeRHEOSurface::unpack_forward_comm(int n, int first, double *buf)
 {
-  int i, k, a, b, m, last;
   int *status = atom->rheo_status;
-
-  m = 0;
-  last = first + n;
-  for (i = first; i < last; i++) {
+  int m = 0;
+  int last = first + n;
+  for (int i = first; i < last; i++) {
     if (comm_stage == 0) {
       divr[i] = buf[m++];
     } else if (comm_stage == 1) {
diff --git a/src/RHEO/compute_rheo_vshift.cpp b/src/RHEO/compute_rheo_vshift.cpp
index c06ef533ac..d030c32924 100644
--- a/src/RHEO/compute_rheo_vshift.cpp
+++ b/src/RHEO/compute_rheo_vshift.cpp
@@ -41,8 +41,8 @@ using namespace RHEO_NS;
 /* ---------------------------------------------------------------------- */
 
 ComputeRHEOVShift::ComputeRHEOVShift(LAMMPS *lmp, int narg, char **arg) :
-  Compute(lmp, narg, arg), list(nullptr), vshift(nullptr), fix_rheo(nullptr),
-  compute_kernel(nullptr), compute_interface(nullptr), compute_surface(nullptr)
+  Compute(lmp, narg, arg), vshift(nullptr), fix_rheo(nullptr), rho0(nullptr), list(nullptr),
+  compute_interface(nullptr), compute_kernel(nullptr), compute_surface(nullptr)
 {
   if (narg != 3) error->all(FLERR,"Illegal compute RHEO/VShift command");
 
@@ -90,7 +90,7 @@ void ComputeRHEOVShift::init_list(int /*id*/, NeighList *ptr)
 
 void ComputeRHEOVShift::compute_peratom()
 {
-  int i, j, a, b, ii, jj, jnum, itype, jtype;
+  int i, j, a, ii, jj, jnum, itype, jtype;
   int fluidi, fluidj;
   double xtmp, ytmp, ztmp, rsq, r, rinv;
   double w, wp, dr, w0, w4, vmag, prefactor;
@@ -229,8 +229,6 @@ void ComputeRHEOVShift::correct_surfaces()
 {
   if (!surface_flag) return;
 
-  int i, a, b;
-
   int *status = atom->rheo_status;
   int *mask = atom->mask;
   double **nsurface = compute_surface->nsurface;
@@ -239,7 +237,7 @@ void ComputeRHEOVShift::correct_surfaces()
   int dim = domain->dimension;
 
   double nx, ny, nz, vx, vy, vz, dot;
-  for (i = 0; i < nlocal; i++) {
+  for (int i = 0; i < nlocal; i++) {
     if (mask[i] & groupbit) {
 
       if (status[i] & PHASECHECK) continue;
@@ -283,11 +281,11 @@ void ComputeRHEOVShift::correct_surfaces()
 
 int ComputeRHEOVShift::pack_reverse_comm(int n, int first, double *buf)
 {
-  int i,m,last;
+  int m, last;
 
   m = 0;
   last = first + n;
-  for (i = first; i < last; i++) {
+  for (int i = first; i < last; i++) {
     buf[m++] = vshift[i][0];
     buf[m++] = vshift[i][1];
     buf[m++] = vshift[i][2];
diff --git a/src/RHEO/fix_rheo.cpp b/src/RHEO/fix_rheo.cpp
index f70b9e121f..2156a82649 100644
--- a/src/RHEO/fix_rheo.cpp
+++ b/src/RHEO/fix_rheo.cpp
@@ -50,8 +50,9 @@ static const char cite_rheo[] =
 /* ---------------------------------------------------------------------- */
 
 FixRHEO::FixRHEO(LAMMPS *lmp, int narg, char **arg) :
-  Fix(lmp, narg, arg), compute_grad(nullptr), compute_kernel(nullptr), compute_surface(nullptr),
-  compute_interface(nullptr), compute_rhosum(nullptr), compute_vshift(nullptr), rho0(nullptr), csq(nullptr)
+  Fix(lmp, narg, arg), rho0(nullptr), csq(nullptr), compute_grad(nullptr), compute_kernel(nullptr),
+  compute_interface(nullptr), compute_surface(nullptr), compute_rhosum(nullptr),
+  compute_vshift(nullptr)
 {
   time_integrate = 1;
 
@@ -457,7 +458,6 @@ void FixRHEO::final_integrate()
   double dtfm, divu;
   int i, a;
 
-  double **x = atom->x;
   double **v = atom->v;
   double **f = atom->f;
   double **gradv = compute_grad->gradv;
diff --git a/src/RHEO/fix_rheo_oxidation.cpp b/src/RHEO/fix_rheo_oxidation.cpp
index a51f2feb95..699f3428c8 100644
--- a/src/RHEO/fix_rheo_oxidation.cpp
+++ b/src/RHEO/fix_rheo_oxidation.cpp
@@ -155,7 +155,6 @@ void FixRHEOOxidation::post_integrate()
   double delx, dely, delz, rsq;
   tagint tagi, tagj;
 
-  int nlocal = atom->nlocal;
   int newton_bond = force->newton_bond;
 
   tagint *tag = atom->tag;
@@ -267,7 +266,7 @@ void FixRHEOOxidation::post_force(int /*vflag*/)
 int FixRHEOOxidation::pack_forward_comm(int n, int *list, double *buf,
                                         int /*pbc_flag*/, int * /*pbc*/)
 {
-  int i, j, k, m;
+  int i, j, m;
   double **x = atom->x;
   m = 0;
 
@@ -284,7 +283,7 @@ int FixRHEOOxidation::pack_forward_comm(int n, int *list, double *buf,
 
 void FixRHEOOxidation::unpack_forward_comm(int n, int first, double *buf)
 {
-  int i, k, m, last;
+  int i, m, last;
   double **x = atom->x;
   m = 0;
   last = first + n;
diff --git a/src/RHEO/fix_rheo_pressure.cpp b/src/RHEO/fix_rheo_pressure.cpp
index 82adf52012..84b0825570 100644
--- a/src/RHEO/fix_rheo_pressure.cpp
+++ b/src/RHEO/fix_rheo_pressure.cpp
@@ -39,7 +39,9 @@ static constexpr double SEVENTH = 1.0 / 7.0;
 /* ---------------------------------------------------------------------- */
 
 FixRHEOPressure::FixRHEOPressure(LAMMPS *lmp, int narg, char **arg) :
-  Fix(lmp, narg, arg), fix_rheo(nullptr), rho0(nullptr), csq(nullptr), rho0inv(nullptr), csqinv(nullptr), c_cubic(nullptr), tpower(nullptr), pbackground(nullptr), pressure_style(nullptr)
+  Fix(lmp, narg, arg), c_cubic(nullptr), csq(nullptr), csqinv(nullptr), rho0(nullptr),
+  rho0inv(nullptr), tpower(nullptr), pbackground(nullptr), pressure_style(nullptr),
+  fix_rheo(nullptr)
 {
   if (narg < 4) error->all(FLERR,"Illegal fix command");
 
@@ -160,9 +162,6 @@ void FixRHEOPressure::setup_pre_force(int /*vflag*/)
 
 void FixRHEOPressure::pre_force(int /*vflag*/)
 {
-  int i;
-  double dr, rr3, rho_ratio;
-
   int *mask = atom->mask;
   int *type = atom->type;
   double *rho = atom->rho;
@@ -170,7 +169,7 @@ void FixRHEOPressure::pre_force(int /*vflag*/)
 
   int nlocal = atom->nlocal;
 
-  for (i = 0; i < nlocal; i++)
+  for (int i = 0; i < nlocal; i++)
     if (mask[i] & groupbit)
       pressure[i] = calc_pressure(rho[i], type[i]);
 
@@ -182,12 +181,10 @@ void FixRHEOPressure::pre_force(int /*vflag*/)
 int FixRHEOPressure::pack_forward_comm(int n, int *list, double *buf,
                                         int /*pbc_flag*/, int * /*pbc*/)
 {
-  int i,j,k,m;
   double *pressure = atom->pressure;
-  m = 0;
-
-  for (i = 0; i < n; i++) {
-    j = list[i];
+  int m = 0;
+  for (int i = 0; i < n; i++) {
+    int j = list[i];
     buf[m++] = pressure[j];
   }
   return m;
@@ -197,12 +194,10 @@ int FixRHEOPressure::pack_forward_comm(int n, int *list, double *buf,
 
 void FixRHEOPressure::unpack_forward_comm(int n, int first, double *buf)
 {
-  int i, k, m, last;
   double *pressure = atom->pressure;
-
-  m = 0;
-  last = first + n;
-  for (i = first; i < last; i++) {
+  int m = 0;
+  int last = first + n;
+  for (int i = first; i < last; i++) {
     pressure[i] = buf[m++];
   }
 }
@@ -211,7 +206,8 @@ void FixRHEOPressure::unpack_forward_comm(int n, int first, double *buf)
 
 double FixRHEOPressure::calc_pressure(double rho, int type)
 {
-  double p, dr, rr3, rho_ratio;
+  double p = 0.0;
+  double dr, rr3, rho_ratio;
 
   if (pressure_style[type] == LINEAR) {
     p = csq[type] * (rho - rho0[type]);
@@ -234,7 +230,7 @@ double FixRHEOPressure::calc_pressure(double rho, int type)
 
 double FixRHEOPressure::calc_rho(double p, int type)
 {
-  double rho, dr, rr3, rho_ratio;
+  double rho = 0.0;
 
   if (pressure_style[type] == LINEAR) {
     rho = csqinv[type] * p + rho0[type];
diff --git a/src/RHEO/fix_rheo_thermal.cpp b/src/RHEO/fix_rheo_thermal.cpp
index c8c1877e32..a7d9fc8df9 100644
--- a/src/RHEO/fix_rheo_thermal.cpp
+++ b/src/RHEO/fix_rheo_thermal.cpp
@@ -60,10 +60,9 @@ static const char cite_rheo_oxide[] =
 /* ---------------------------------------------------------------------- */
 
 FixRHEOThermal::FixRHEOThermal(LAMMPS *lmp, int narg, char **arg) :
-  Fix(lmp, narg, arg), fix_rheo(nullptr), compute_grad(nullptr), compute_vshift(nullptr),
-  Tc(nullptr), kappa(nullptr), cv(nullptr), L(nullptr),
-  Tc_style(nullptr), kappa_style(nullptr), cv_style(nullptr), L_style(nullptr),
-  fix_update_special_bonds(nullptr)
+  Fix(lmp, narg, arg), cv(nullptr), Tc(nullptr), kappa(nullptr), L(nullptr), cv_style(nullptr),
+  Tc_style(nullptr), kappa_style(nullptr),  L_style(nullptr), list(nullptr), fix_rheo(nullptr),
+  compute_grad(nullptr), compute_vshift(nullptr), fix_update_special_bonds(nullptr)
 {
   if (narg < 4) error->all(FLERR,"Illegal fix command");
 
@@ -189,13 +188,14 @@ FixRHEOThermal::FixRHEOThermal(LAMMPS *lmp, int narg, char **arg) :
       cut_bond = utils::numeric(FLERR, arg[iarg + 1], false, lmp);
       btype = utils::numeric(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) error->all(FLERR, "Illegal value for bond type");
+      if (cut_bond <= 0.0) error->all(FLERR, "Illegal max bond length must be greater than zero");
+      if ((btype < 1) || (btype > atom->nbondtypes))
+        error->all(FLERR, "Illegal value {} for bond type", btype);
 
       cutsq_bond = cut_bond * cut_bond;
       iarg += 3;
     } else {
-      error->all(FLERR,"Illegal fix command, {}", arg[iarg]);
+      error->all(FLERR,"Unknown fix rheo/thermal keyword: {}", arg[iarg]);
     }
   }
 
@@ -460,20 +460,16 @@ void FixRHEOThermal::post_neighbor()
 
 void FixRHEOThermal::pre_force(int /*vflag*/)
 {
-  int i, itype;
   double cvi, Tci, Ti, Li;
 
   double *energy = atom->esph;
   double *temperature = atom->temperature;
   int *type = atom->type;
-  int *status = atom->rheo_status;
-
-  int nlocal = atom->nlocal;
-  int nall = nlocal + atom->nghost;
+  int nall = atom->nlocal + atom->nghost;
 
   // Calculate temperature
-  for (i = 0; i < nall; i++) {
-    itype = type[i];
+  for (int i = 0; i < nall; i++) {
+    int itype = type[i];
     cvi = calc_cv(i, itype);
     temperature[i] = energy[i] / cvi;
 
@@ -529,7 +525,6 @@ void FixRHEOThermal::break_bonds()
   int nbondlist = neighbor->nbondlist;
 
   int nlocal = atom->nlocal;
-  int nall = nlocal + atom->nghost;
 
   // Delete all bonds for local atoms that melt of a given type
   for (int i = 0; i < nlocal; i++) {
@@ -753,13 +748,11 @@ double FixRHEOThermal::calc_L(int i, int itype)
 int FixRHEOThermal::pack_forward_comm(int n, int *list, double *buf,
                                         int /*pbc_flag*/, int * /*pbc*/)
 {
-  int i, j, k, m;
   int *status = atom->rheo_status;
   double **x = atom->x;
-  m = 0;
-
-  for (i = 0; i < n; i++) {
-    j = list[i];
+  int m = 0;
+  for (int i = 0; i < n; i++) {
+    int j = list[i];
     buf[m++] = ubuf(status[j]).d;
     buf[m++] = x[j][0];
     buf[m++] = x[j][1];
@@ -772,12 +765,11 @@ int FixRHEOThermal::pack_forward_comm(int n, int *list, double *buf,
 
 void FixRHEOThermal::unpack_forward_comm(int n, int first, double *buf)
 {
-  int i, k, m, last;
   int *status = atom->rheo_status;
   double **x = atom->x;
-  m = 0;
-  last = first + n;
-  for (i = first; i < last; i++) {
+  int m = 0;
+  int last = first + n;
+  for (int i = first; i < last; i++) {
     status[i] = (int) ubuf(buf[m++]).i;
     x[i][0] = buf[m++];
     x[i][1] = buf[m++];
diff --git a/src/RHEO/fix_rheo_viscosity.cpp b/src/RHEO/fix_rheo_viscosity.cpp
index 3ca7fd8d13..5f7a1208c5 100644
--- a/src/RHEO/fix_rheo_viscosity.cpp
+++ b/src/RHEO/fix_rheo_viscosity.cpp
@@ -38,8 +38,8 @@ enum {NONE, CONSTANT, POWER};
 /* ---------------------------------------------------------------------- */
 
 FixRHEOViscosity::FixRHEOViscosity(LAMMPS *lmp, int narg, char **arg) :
-  Fix(lmp, narg, arg), fix_rheo(nullptr), compute_grad(nullptr), eta(nullptr),
-  npow(nullptr), K(nullptr), gd0(nullptr), tau0(nullptr), viscosity_style(nullptr)
+  Fix(lmp, narg, arg), eta(nullptr), npow(nullptr), K(nullptr), gd0(nullptr), tau0(nullptr),
+  viscosity_style(nullptr), fix_rheo(nullptr), compute_grad(nullptr)
 {
   if (narg < 4) error->all(FLERR, "Illegal fix command");
 
@@ -225,12 +225,10 @@ void FixRHEOViscosity::pre_force(int /*vflag*/)
 int FixRHEOViscosity::pack_forward_comm(int n, int *list, double *buf,
                                         int /*pbc_flag*/, int * /*pbc*/)
 {
-  int i, j, k, m;
   double *viscosity = atom->viscosity;
-  m = 0;
-
-  for (i = 0; i < n; i++) {
-    j = list[i];
+  int m = 0;
+  for (int i = 0; i < n; i++) {
+    int j = list[i];
     buf[m++] = viscosity[j];
   }
   return m;
@@ -240,12 +238,10 @@ int FixRHEOViscosity::pack_forward_comm(int n, int *list, double *buf,
 
 void FixRHEOViscosity::unpack_forward_comm(int n, int first, double *buf)
 {
-  int i, k, m, last;
   double *viscosity = atom->viscosity;
-
-  m = 0;
-  last = first + n;
-  for (i = first; i < last; i++) {
+  int m = 0;
+  int last = first + n;
+  for (int i = first; i < last; i++) {
     viscosity[i] = buf[m++];
   }
 }
diff --git a/src/RHEO/pair_rheo.cpp b/src/RHEO/pair_rheo.cpp
index 9ebf884b6e..c682474bb2 100644
--- a/src/RHEO/pair_rheo.cpp
+++ b/src/RHEO/pair_rheo.cpp
@@ -47,8 +47,9 @@ static constexpr double EPSILON = 1e-2;
 /* ---------------------------------------------------------------------- */
 
 PairRHEO::PairRHEO(LAMMPS *lmp) :
-  Pair(lmp), compute_kernel(nullptr), compute_grad(nullptr), compute_interface(nullptr), fix_rheo(nullptr),
-  fix_pressure(nullptr), rho0(nullptr), csq(nullptr), cs(nullptr)
+  Pair(lmp), csq(nullptr), rho0(nullptr), cs(nullptr), compute_kernel(nullptr),
+  compute_grad(nullptr), compute_interface(nullptr), fix_rheo(nullptr),
+  fix_pressure(nullptr)
 {
   restartinfo = 0;
   single_enable = 0;
@@ -80,10 +81,10 @@ 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, w, wp, Ti, Tj, dT, csq_ave, cs_ave;
+  double xtmp, ytmp, ztmp, wp, Ti, Tj, dT, csq_ave, cs_ave;
   double rhoi, rhoj, rho0i, rho0j, voli, volj, Pi, Pj, etai, etaj, kappai, kappaj, eta_ave, kappa_ave, dT_prefactor;
   double mu, q, fp_prefactor, drho_damp, fmag, psi_ij, Fij;
-  double *dWij, *dWji, *dW1ij, *dW1ji;
+  double *dWij, *dWji;
   double dx[3], du[3], dv[3], fv[3], dfp[3], fsolid[3], ft[3], vi[3], vj[3];
 
   int *ilist, *jlist, *numneigh, **firstneigh;
@@ -108,10 +109,8 @@ void PairRHEO::compute(int eflag, int vflag)
   double *conductivity = atom->conductivity;
   double *temperature = atom->temperature;
   double *heatflow = atom->heatflow;
-  double *special_lj = force->special_lj;
   int *type = atom->type;
   int *status = atom->rheo_status;
-  tagint *tag = atom->tag;
 
   double **fp_store, *chi;
   if (compute_interface) {
@@ -461,11 +460,13 @@ void PairRHEO::setup()
 {
   auto fixes = modify->get_fix_by_style("rheo");
   if (fixes.size() == 0) error->all(FLERR, "Need to define fix rheo to use pair rheo");
+  if (fixes.size() > 1) error->all(FLERR, "Must have only one fix rheo defined");
   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");
+  if (fixes.size() > 1) error->all(FLERR, "Must have only one fix rheo/pressure defined");
   fix_pressure = dynamic_cast<FixRHEOPressure *>(fixes[0]);
 
   compute_kernel = fix_rheo->compute_kernel;
@@ -518,12 +519,10 @@ double PairRHEO::init_one(int i, int j)
 
 int PairRHEO::pack_reverse_comm(int n, int first, double *buf)
 {
-  int i, k, m, last;
   double **fp_store = compute_interface->fp_store;
-
-  m = 0;
-  last = first + n;
-  for (i = first; i < last; i++) {
+  int m = 0;
+  int last = first + n;
+  for (int i = first; i < last; i++) {
     buf[m++] = fp_store[i][0];
     buf[m++] = fp_store[i][1];
     buf[m++] = fp_store[i][2];
@@ -536,12 +535,10 @@ int PairRHEO::pack_reverse_comm(int n, int first, double *buf)
 
 void PairRHEO::unpack_reverse_comm(int n, int *list, double *buf)
 {
-  int i, j, k, m;
   double **fp_store = compute_interface->fp_store;
-
-  m = 0;
-  for (i = 0; i < n; i++) {
-    j = list[i];
+  int m = 0;
+  for (int i = 0; i < n; i++) {
+    int j = list[i];
     fp_store[j][0] += buf[m++];
     fp_store[j][1] += buf[m++];
     fp_store[j][2] += buf[m++];