diff --git a/src/RHEO/compute_rheo_grad.cpp b/src/RHEO/compute_rheo_grad.cpp
index fb48d90d56..fdddea18b5 100644
--- a/src/RHEO/compute_rheo_grad.cpp
+++ b/src/RHEO/compute_rheo_grad.cpp
@@ -147,6 +147,7 @@ void ComputeRHEOGrad::compute_peratom()
   int *status = atom->rheo_status;
   int *type = atom->type;
   double *mass = atom->mass;
+  double *rmass = atom->rmass;
   int newton = force->newton;
   int dim = domain->dimension;
 
@@ -225,8 +226,13 @@ void ComputeRHEOGrad::compute_peratom()
           }
         }
 
-        Voli = mass[itype] / rhoi;
-        Volj = mass[jtype] / rhoj;
+        if (rmass) {
+          Voli = rmass[i] / rhoi;
+          Volj = rmass[j] / rhoj;
+        } else {
+          Voli = mass[itype] / rhoi;
+          Volj = mass[jtype] / rhoj;
+        }
 
         vij[0] = vi[0] - vj[0];
         vij[1] = vi[1] - vj[1];
diff --git a/src/RHEO/compute_rheo_interface.cpp b/src/RHEO/compute_rheo_interface.cpp
index 05821ebe22..2d2a368926 100644
--- a/src/RHEO/compute_rheo_interface.cpp
+++ b/src/RHEO/compute_rheo_interface.cpp
@@ -335,6 +335,7 @@ void ComputeRHEOInterface::store_forces()
   int *type = atom->type;
   int *mask = atom->mask;
   double *mass = atom->mass;
+  double *rmass = atom->rmass;
   double **f = atom->f;
 
   // When this is called, fp_store stores the pressure force
@@ -346,7 +347,8 @@ void ComputeRHEOInterface::store_forces()
   if (fixlist.size() != 0) {
     for (const auto &fix : fixlist) {
       for (int i = 0; i < atom->nlocal; i++) {
-        minv = 1.0 / mass[type[i]];
+        if (rmass) minv = 1.0 / rmass[i];
+        else minv = 1.0 / mass[type[i]];
         if (mask[i] & fix->groupbit)
           for (int a = 0; a < 3; a++)
             fp_store[i][a] = f[i][a] * minv;
@@ -356,10 +358,18 @@ void ComputeRHEOInterface::store_forces()
       }
     }
   } else {
-    for (int i = 0; i < atom->nlocal; i++) {
-      minv = 1.0 / mass[type[i]];
-      for (int a = 0; a < 3; a++)
-        fp_store[i][a] = (f[i][a] - fp_store[i][a]) * minv;
+    if (rmass) {
+      for (int i = 0; i < atom->nlocal; i++) {
+        minv = 1.0 / rmass[i];
+        for (int a = 0; a < 3; a++)
+          fp_store[i][a] = (f[i][a] - fp_store[i][a]) * minv;
+      }
+    } else {
+      for (int i = 0; i < atom->nlocal; i++) {
+        minv = 1.0 / mass[type[i]];
+        for (int a = 0; a < 3; a++)
+          fp_store[i][a] = (f[i][a] - fp_store[i][a]) * minv;
+      }
     }
   }
 
diff --git a/src/RHEO/compute_rheo_kernel.cpp b/src/RHEO/compute_rheo_kernel.cpp
index 1913bb1e53..108ec58b09 100644
--- a/src/RHEO/compute_rheo_kernel.cpp
+++ b/src/RHEO/compute_rheo_kernel.cpp
@@ -48,8 +48,6 @@ using namespace RHEO_NS;
 using namespace MathConst;
 using namespace MathExtra;
 
-static constexpr int DELTA = 2000;
-
 /* ---------------------------------------------------------------------- */
 
 ComputeRHEOKernel::ComputeRHEOKernel(LAMMPS *lmp, int narg, char **arg) :
@@ -581,6 +579,7 @@ void ComputeRHEOKernel::compute_peratom()
   double **x = atom->x;
   int *type = atom->type;
   double *mass = atom->mass;
+  double *rmass = atom->rmass;
   double *rho = atom->rho;
   int *status = atom->rheo_status;
   tagint *tag = atom->tag;
@@ -626,7 +625,8 @@ void ComputeRHEOKernel::compute_peratom()
             if (status[j] & PHASECHECK)
               rhoj = compute_interface->correct_rho(j, i);
 
-          vj = mass[type[j]] / rhoj;
+          if (rmass) vj = rmass[j] / rhoj;
+          else vj = mass[type[j]] / rhoj;
           M += w * vj;
         }
       }
@@ -647,7 +647,6 @@ void ComputeRHEOKernel::compute_peratom()
 
       jlist = firstneigh[i];
       jnum = numneigh[i];
-      itype = type[i];
 
       // Zero upper-triangle M and cut (will be symmetric):
       for (a = 0; a < Mdim; a++) {
@@ -675,7 +674,8 @@ void ComputeRHEOKernel::compute_peratom()
             if (status[j] & PHASECHECK)
               rhoj = compute_interface->correct_rho(j, i);
 
-          vj = mass[type[j]] / rhoj;
+          if (rmass) vj = rmass[j] / rhoj;
+          else vj = mass[type[j]] / rhoj;
 
           //Populate the H-vector of polynomials (2D)
           if (dim == 2) {
diff --git a/src/RHEO/compute_rheo_surface.cpp b/src/RHEO/compute_rheo_surface.cpp
index 83ef659e84..314d2d2ef5 100644
--- a/src/RHEO/compute_rheo_surface.cpp
+++ b/src/RHEO/compute_rheo_surface.cpp
@@ -113,6 +113,7 @@ void ComputeRHEOSurface::compute_peratom()
   int *mask = atom->mask;
   int *type = atom->type;
   double *mass = atom->mass;
+  double *rmass = atom->rmass;
   double *rho = atom->rho;
   int *coordination = compute_kernel->coordination;
 
@@ -173,9 +174,13 @@ void ComputeRHEOSurface::compute_peratom()
           }
         }
 
-        Voli = mass[itype] / rhoi;
-        Volj = mass[jtype] / rhoj;
-
+        if (rmass) {
+          Voli = rmass[i] / rhoi;
+          Volj = rmass[j] / rhoj;
+        } else {
+          Voli = mass[itype] / rhoi;
+          Volj = mass[jtype] / rhoj;
+        }
         compute_kernel->calc_dw_quintic(i, j, dx[0], dx[1], dx[2], sqrt(rsq), dWij, dWji);
 
         for (a = 0; a < dim; a++) {
diff --git a/src/RHEO/compute_rheo_vshift.cpp b/src/RHEO/compute_rheo_vshift.cpp
index d030c32924..8d2f04d7f6 100644
--- a/src/RHEO/compute_rheo_vshift.cpp
+++ b/src/RHEO/compute_rheo_vshift.cpp
@@ -108,6 +108,7 @@ void ComputeRHEOVShift::compute_peratom()
   double **v = atom->v;
   double *rho = atom->rho;
   double *mass = atom->mass;
+  double *rmass = atom->rmass;
 
   int nlocal = atom->nlocal;
   int nall = nlocal + atom->nghost;
@@ -140,7 +141,8 @@ void ComputeRHEOVShift::compute_peratom()
     itype = type[i];
     jlist = firstneigh[i];
     jnum = numneigh[i];
-    imass = mass[itype];
+    if (rmass) imass = rmass[i];
+    else imass = mass[itype];
     fluidi = !(status[i] & PHASECHECK);
 
     for (jj = 0; jj < jnum; jj++) {
@@ -160,7 +162,8 @@ void ComputeRHEOVShift::compute_peratom()
 
       if (rsq < cutsq) {
         jtype = type[j];
-        jmass = mass[jtype];
+        if (rmass) jmass = rmass[j];
+        else jmass = mass[jtype];
 
         r = sqrt(rsq);
         rinv = 1 / r;
diff --git a/src/RHEO/fix_rheo.cpp b/src/RHEO/fix_rheo.cpp
index 2156a82649..b093eff681 100644
--- a/src/RHEO/fix_rheo.cpp
+++ b/src/RHEO/fix_rheo.cpp
@@ -469,7 +469,6 @@ void FixRHEO::final_integrate()
   int *mask = atom->mask;
   int *status = atom->rheo_status;
 
-  int rmass_flag = atom->rmass_flag;
   int dim = domain->dimension;
 
   // Update velocity
@@ -477,7 +476,7 @@ void FixRHEO::final_integrate()
     if (mask[i] & groupbit) {
       if (status[i] & STATUS_NO_INTEGRATION) continue;
 
-      if (rmass_flag) {
+      if (rmass) {
         dtfm = dtf / rmass[i];
       } else {
         dtfm = dtf / mass[type[i]];
diff --git a/src/RHEO/pair_rheo.cpp b/src/RHEO/pair_rheo.cpp
index c682474bb2..520a1c4470 100644
--- a/src/RHEO/pair_rheo.cpp
+++ b/src/RHEO/pair_rheo.cpp
@@ -103,6 +103,7 @@ void PairRHEO::compute(int eflag, int vflag)
   double **f = atom->f;
   double *rho = atom->rho;
   double *mass = atom->mass;
+  double *rmass = atom->rmass;
   double *drho = atom->drho;
   double *pressure = atom->pressure;
   double *viscosity = atom->viscosity;
@@ -149,7 +150,8 @@ void PairRHEO::compute(int eflag, int vflag)
     itype = type[i];
     jlist = firstneigh[i];
     jnum = numneigh[i];
-    imass = mass[itype];
+    if (rmass) imass = rmass[i];
+    else imass = mass[itype];
     etai = viscosity[i];
     fluidi = !(status[i] & PHASECHECK);
     if (thermal_flag) {
@@ -171,7 +173,8 @@ void PairRHEO::compute(int eflag, int vflag)
         r = sqrt(rsq);
         rinv = 1 / r;
 
-        jmass = mass[jtype];
+        if (rmass) jmass = rmass[i];
+        else jmass = mass[jtype];
         etaj = viscosity[j];
         fluidj = !(status[j] & PHASECHECK);
         if (thermal_flag) {