diff --git a/src/RHEO/compute_rheo_interface.cpp b/src/RHEO/compute_rheo_interface.cpp
index ae1bc13831..ff76b75f77 100644
--- a/src/RHEO/compute_rheo_interface.cpp
+++ b/src/RHEO/compute_rheo_interface.cpp
@@ -61,8 +61,8 @@ ComputeRHEOInterface::~ComputeRHEOInterface()
   if (index != -1) atom->remove_custom(index, 1, 0);
 
   if (id_fix_pa && modify->nfix) modify->delete_fix(id_fix_pa);
+  delete[] id_fix_pa;
 
-  memory->destroy(fx_m_norm);
   memory->destroy(norm);
   memory->destroy(normwf);
 }
@@ -72,7 +72,10 @@ ComputeRHEOInterface::~ComputeRHEOInterface()
 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;
   cutsq = cut * cut;
   wall_max = sqrt(3.0) / 12.0 * cut;
 
@@ -87,21 +90,25 @@ void ComputeRHEOInterface::init()
   int index = atom->find_custom("rheo_chi", tmp1, tmp2);
   if (index == -1) {
     index = atom->add_custom("rheo_chi", 1, 0);
+    memory->destroy(norm);
+    memory->destroy(normwf);
+    memory->create(norm, nmax, "rheo/interface:norm");
+    memory->create(normwf, nmax, "rheo/interface:normwf");
     nmax_old = nmax;
   }
   chi = atom->dvector[index];
 
-  // For fpressure, go ahead and create an instance of fix property atom
+  // For fp_store, go ahead and create an instance of fix property atom
   // Need restarts + exchanging with neighbors since it needs to persist
   // between timesteps (fix property atom will handle callbacks)
 
-  index = atom->find_custom("rheo_pressure", tmp1, tmp2);
+  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_f_pressure 3", id_fix_pa)));
-    index = atom->find_custom("rheo_pressure", tmp1, tmp2);
+    modify->add_fix(fmt::format("{} all property/atom d2_fp_store 3", id_fix_pa)));
+    index = atom->find_custom("fp_store", tmp1, tmp2);
   }
-  f_pressure = atom->darray[index];
+  fp_store = atom->darray[index];
 
   // need an occasional half neighbor list
   neighbor->add_request(this, NeighConst::REQ_HALF);
@@ -116,51 +123,37 @@ void ComputeRHEOInterface::init_list(int /*id*/, NeighList *ptr)
 
 /* ---------------------------------------------------------------------- */
 
-// Left off here
 void ComputeRHEOInterface::compute_peratom()
 {
-  int i, j, ii, jj, jnum, itype, jtype, phase_match;
-  double xtmp, ytmp, ztmp, delx, dely, delz, rsq;
-  int *jlist;
-  double w;
+  int i, j, ii, jj, jnum, itype, jtype, fluidi, fluidj, status_match;
+  double xtmp, ytmp, ztmp, delx, dely, delz, rsq, w, dot;
 
-  // neighbor list ariables
-  int inum, *ilist, *numneigh, **firstneigh;
+  int inum, *ilist, *jlist, *numneigh, **firstneigh;
   int nlocal = atom->nlocal;
   int nall = nlocal + atom->nghost;
 
   double **x = atom->x;
   int *type = atom->type;
   int newton = force->newton;
-  int *phase = atom->phase;
+  int *status = atom->status;
   double *rho = atom->rho;
-  double *fx = atom->dvector[index_fx];
-  double *fy = atom->dvector[index_fy];
-  double *fz = atom->dvector[index_fz];
-  double mi, mj;
-
-  //Declare mass pointer to calculate acceleration from force
-  double *mass = atom->mass;
-
-  cs2 = 1.0;
 
   inum = list->inum;
   ilist = list->ilist;
   numneigh = list->numneigh;
   firstneigh = list->firstneigh;
 
-  if (atom->nmax > nmax) {
-    nmax = atom->nmax;
-    fix_chi->grow_arrays(nmax);
-    chi = fix_chi->vstore;
+  if (atom->nmax > nmax_old) {
+    nmax_old = atom->nmax;
     memory->destroy(norm);
     memory->destroy(normwf);
-    memory->create(norm, nmax, "RHEO/chi:norm");
-    memory->create(normwf, nmax, "RHEO/chi: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");
   }
 
   for (i = 0; i < nall; i++) {
-    if (phase[i] > FixRHEO::FLUID_MAX) rho[i] = 0.0;
+    if (!(status[i] & FixRHEO::STATUS_FLUID)) rho[i] = 0.0;
     normwf[i] = 0.0;
     norm[i] = 0.0;
     chi[i] = 0.0;
@@ -172,9 +165,9 @@ void ComputeRHEOInterface::compute_peratom()
     ytmp = x[i][1];
     ztmp = x[i][2];
     itype = type[i];
+    fluidi = status[i] & FixRHEO::STATUS_FLUID;
     jlist = firstneigh[i];
     jnum = numneigh[i];
-    mi = mass[type[i]];
 
     for (jj = 0; jj < jnum; jj++) {
       j = jlist[jj];
@@ -187,37 +180,36 @@ void ComputeRHEOInterface::compute_peratom()
 
       if (rsq < cutsq) {
         jtype = type[j];
-        mj = mass[type[j]];
+        fluidj = status[j] & FixRHEO::STATUS_FLUID;
         w = compute_kernel->calc_w_quintic(i, j, delx, dely, delz, sqrt(rsq));
 
-        phase_match = 0;
+        status_match = 0;
         norm[i] += w;
-        if ((phase[i] <= FixRHEO::FLUID_MAX && phase[j] <= FixRHEO::FLUID_MAX)
-            || (phase[i] > FixRHEO::FLUID_MAX && phase[j] > FixRHEO::FLUID_MAX)) {
-          phase_match = 1;
-        }
+        if ((fluidi && fluidj) || ((!fluid) && (!fluidj)))
+          status_match = 1;
 
-        if (phase_match) {
+        if (status_match) {
           chi[i] += w;
         } else {
-          if (phase[i] > FixRHEO::FLUID_MAX) {
-            //speed of sound and rho0 assumed to = 1 (units in density not pressure)
-            //In general, rho is calculated using the force vector on the wall particle
-            //fx stores f-fp
-            rho[i] += w*(cs2*(rho[j] - 1.0) - rho[j]*((-fx[j]/mj+fx[i]/mi)*delx + (-fy[j]/mj+fy[i]/mi)*dely + (-fz[j]/mj+fz[i]/mi)*delz));
-            //For the specific taste case whre force on wall particles = 0
-            //rho[i] += w*(1.0*(rho[j] - 1.0) + rho[j]*(1e-3*delx));
+          if (!fluidi) {
+            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);
             normwf[i] += w;
           }
         }
 
         if (newton || j < nlocal) {
           norm[j] += w;
-          if (phase_match) {
+          if (status_match) {
             chi[j] += w;
           } else {
-            if (phase[j] > FixRHEO::FLUID_MAX) {
-              rho[j] += w*(cs2*(rho[i] - 1.0) + rho[i]*((-fx[i]/mi+fx[j]/mj)*delx + (-fy[i]/mi+fy[j]/mj)*dely + (-fz[i]/mi+fz[j]/mj)*delz));
+            if (!fluidj) {
+              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);
               normwf[j] += w;
             }
           }
@@ -226,41 +218,41 @@ void ComputeRHEOInterface::compute_peratom()
     }
   }
 
-  if (newton) comm->reverse_comm_compute(this);
+  if (newton) comm->reverse_comm(this);
 
   for (i = 0; i < nlocal; i++) {
     if (norm[i] != 0.0) chi[i] /= norm[i];
-    if (normwf[i] != 0.0) { // Only if it's a wall particle
-      rho[i] = 1.0 + (rho[i] / normwf[i])/cs2; // Stores rho for solid particles 1+Pw in Adami Adams 2012
-      if (rho[i] < EPSILON) rho[i] = EPSILON;
-    }
 
-    if (normwf[i] == 0.0 && phase[i] > FixRHEO::FLUID_MAX) rho[i] = 1.0;
+    // Recalculate rho for non-fluid particles
+    if (!(status[i] & FixRHEO::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);
+      } else {
+        rho[i] = rho0;
+      }
+    }
   }
 
   comm_stage = 1;
   comm_forward = 2;
-  comm->forward_comm_compute(this);
+  comm->forward_comm(this);
 }
 
 /* ---------------------------------------------------------------------- */
 
-int ComputeRHEOInterface::pack_forward_comm(int n, int *list, double *buf,
-                                        int /*pbc_flag*/, int * /*pbc*/)
+int ComputeRHEOInterface::pack_forward_comm(int n, int *list, double *buf, int /*pbc_flag*/, int * /*pbc*/)
 {
   int i,j,k,m;
   m = 0;
   double *rho = atom->rho;
-  double *fx = atom->dvector[index_fx];
-  double *fy = atom->dvector[index_fy];
-  double *fz = atom->dvector[index_fz];
 
   for (i = 0; i < n; i++) {
     j = list[i];
     if (comm_stage == 0) {
-      buf[m++] = fx[j];
-      buf[m++] = fy[j];
-      buf[m++] = fz[j];
+      buf[m++] = fp_store[j][0];
+      buf[m++] = fp_store[j][1];
+      buf[m++] = fp_store[j][2];
     } else {
       buf[m++] = chi[j];
       buf[m++] = rho[j];
@@ -275,20 +267,16 @@ void ComputeRHEOInterface::unpack_forward_comm(int n, int first, double *buf)
 {
   int i, k, m, last;
   double *rho = atom->rho;
-  double *fx = atom->dvector[index_fx];
-  double *fy = atom->dvector[index_fy];
-  double *fz = atom->dvector[index_fz];
-
   m = 0;
   last = first + n;
   for (i = first; i < last; i++) {
     if (comm_stage == 0) {
-      fx[i] = buf[m++];
-      fy[i] = buf[m++];
-      fz[i] = buf[m++];
+      fp_store[i][0] = buf[m++];
+      fp_store[i][1] = buf[m++];
+      fp_store[i][2] = buf[m++];
     } else {
       chi[i] = buf[m++];
-      rho[i] = buf[m++]; // Won't do anything for fluids
+      rho[i] = buf[m++];
     }
   }
 }
@@ -317,13 +305,13 @@ void ComputeRHEOInterface::unpack_reverse_comm(int n, int *list, double *buf)
 {
   int i,k,j,m;
   double *rho = atom->rho;
-  int *phase = atom->phase;
+  int *status = atom->status;
   m = 0;
   for (i = 0; i < n; i++) {
     j = list[i];
     norm[j] += buf[m++];
     chi[j] += buf[m++];
-    if (phase[j] > FixRHEO::FLUID_MAX){
+    if (!(status[j] & FixRHEO::STATUS_FLUID)){
       normwf[j] += buf[m++];
       rho[j] += buf[m++];
     } else {
@@ -339,22 +327,23 @@ void ComputeRHEOInterface::correct_v(double *vi, double *vj, double *vi_out, int
 {
   double wall_prefactor, wall_denom, wall_numer;
 
-  wall_numer = 2.0*cut*(chi[i]-0.5);
+  wall_numer = 2.0 * cut * (chi[i] - 0.5);
   if (wall_numer < 0) wall_numer = 0;
-  wall_denom = 2.0*cut*(chi[j]-0.5);
+  wall_denom = 2.0 * cut * (chi[j] - 0.5);
   if (wall_denom < wall_max) wall_denom = wall_max;
 
   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_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];
 }
 
 /* ---------------------------------------------------------------------- */
 
-double ComputeRHEOInterface::correct_rho(int i, int j) // i is wall, j is fluid
+double ComputeRHEOInterface::correct_rho(int i, int j)
 {
+  // i is wall, j is fluid
   //In future may depend on atom type j's pressure equation
   return atom->rho[i];
 }
@@ -363,42 +352,46 @@ double ComputeRHEOInterface::correct_rho(int i, int j) // i is wall, j is fluid
 
 void ComputeRHEOInterface::store_forces()
 {
-  double *fx = atom->dvector[index_fx];
-  double *fy = atom->dvector[index_fy];
-  double *fz = atom->dvector[index_fz];
+  double minv;
+  double mass = atom->mass;
+  double type = atom->type;
   double **f = atom->f;
-  double **fp = atom->fp;
   int *mask = atom->mask;
 
-  int flag;
+  // 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) {
-        fx[i] = f[i][0];
-        fy[i] = f[i][1];
-        fz[i] = f[i][2];
+        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 {
-        fx[i] = f[i][0] - fp[i][0];
-        fy[i] = f[i][1] - fp[i][1];
-        fz[i] = f[i][2] - fp[i][2];
+        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 {
     for (int i = 0; i < atom->nlocal; i++) {
-      fx[i] = f[i][0] - fp[i][0];
-      fy[i] = f[i][1] - fp[i][1];
-      fz[i] = f[i][2] - fp[i][2];
+      minv = 1.0 / mass[type[i]];
+      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;
     }
   }
 
-  //Forward comm forces -note only needed here b/c property atom will forward otherwise
+  // Forward comm forces
   comm_forward = 3;
   comm_stage = 0;
-  comm->forward_comm_compute(this);
+  comm->forward_comm(this);
 }
 
-
 /* ----------------------------------------------------------------------
    memory usage of local atom-based array
 ------------------------------------------------------------------------- */
diff --git a/src/RHEO/compute_rheo_interface.h b/src/RHEO/compute_rheo_interface.h
index f314c16aa9..8e190c5430 100644
--- a/src/RHEO/compute_rheo_interface.h
+++ b/src/RHEO/compute_rheo_interface.h
@@ -44,8 +44,8 @@ class ComputeRHEOInterface : public Compute {
 
  private:
   int nmax_old, comm_stage;
-  double cut, cutsq, cs, wall_max;
-  double **fx_m_norm, *norm, *normwf;
+  double rho0, cut, cutsq, cs, cs_inv, wall_max;
+  double *norm, *normwf, **fom_store;
 
   char *id_fix_pa;
 
diff --git a/src/RHEO/compute_rheo_vshift.cpp b/src/RHEO/compute_rheo_vshift.cpp
index dccd3867cd..6cbd6e96da 100644
--- a/src/RHEO/compute_rheo_vshift.cpp
+++ b/src/RHEO/compute_rheo_vshift.cpp
@@ -125,9 +125,10 @@ void ComputeRHEOVShift::compute_peratom()
   numneigh = list->numneigh;
   firstneigh = list->firstneigh;
 
-  if (nmax_old < atom->nmax)
+  if (nmax_old < atom->nmax) {
     memory->grow(vshift, atom->nmax, 3, "atom:rheo_vshift");
-  nmax_old = atom->nmax;
+    nmax_old = atom->nmax;
+  }
 
   for (i = 0; i < nall; i++)
     for (a = 0; a < dim; a++)
@@ -231,7 +232,6 @@ void ComputeRHEOVShift::correct_surfaces()
   int dim = domain->dimension;
 
   int tmp1, tmp2;
-  define after surf
   int index_nsurf = atom->find_custom("rheo_nsurf", tmp1, tmp2);
   if (index_nsurf == -1) error->all(FLERR, "Cannot find rheo nsurf");
   double **nsurf = atom->darray[index_nsurf];
@@ -239,7 +239,7 @@ void ComputeRHEOVShift::correct_surfaces()
   double nx,ny,nz,vx,vy,vz;
   for (i = 0; i < nlocal; i++) {
     if (mask[i] & groupbit) {
-      if ((surface[i] & FixRHEO::STATUS_SURFACE) || (surface[i] & FixRHEO::STATUS_LAYER)) {
+      if ((status[i] & FixRHEO::STATUS_SURFACE) || (status[i] & FixRHEO::STATUS_LAYER)) {
         nx = nsurf[i][0];
         ny = nsurf[i][1];
         vx = vshift[i][0];
diff --git a/src/RHEO/fix_rheo_thermal.cpp b/src/RHEO/fix_rheo_thermal.cpp
index 1a5d894ced..1912dc9f8c 100644
--- a/src/RHEO/fix_rheo_thermal.cpp
+++ b/src/RHEO/fix_rheo_thermal.cpp
@@ -293,9 +293,10 @@ void FixRHEOThermal::post_neighbor()
   int nlocal = atom->nlocal;
   int nall = nlocal + atom->nghost;
 
-  if (first_flag && (nmax_old < atom->nmax))
+  if (first_flag && (nmax_old < atom->nmax)) {
     memory->grow(conductivity, atom->nmax, "atom:rheo_conductivity");
-  nmax_old = atom->nmax;
+    nmax_old = atom->nmax;
+  }
 
   if (conductivity_style == CONSTANT) {
     for (i = 0; i < nall; i++)
@@ -328,9 +329,10 @@ 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_old < atom->nmax)) {
   //  memory->grow(conductivity, atom->nmax, "atom:rheo_conductivity");
-  //nmax_old = atom->nmax;
+  //  nmax_old = atom->nmax;
+  //}
 
   //if (conductivity_style == TBD) {
   //  for (i = 0; i < nlocal; i++) {
diff --git a/src/RHEO/fix_rheo_viscosity.cpp b/src/RHEO/fix_rheo_viscosity.cpp
index bb7094c43e..5ae1b95529 100644
--- a/src/RHEO/fix_rheo_viscosity.cpp
+++ b/src/RHEO/fix_rheo_viscosity.cpp
@@ -167,9 +167,10 @@ void FixRHEOViscosity::post_neighbor()
   int nlocal = atom->nlocal;
   int nall = nlocal + atom->nghost;
 
-  if (first_flag && (nmax_old < atom->nmax))
+  if (first_flag && (nmax_old < atom->nmax)) {
     memory->grow(viscosity, atom->nmax, "atom:rheo_viscosity");
-  nmax_old = atom->nmax;
+    nmax_old = atom->nmax;
+  }
 
   if (viscosity_style == CONSTANT) {
     for (i = 0; i < nall; i++)
@@ -196,9 +197,10 @@ 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_old < atom->nmax)) {
     memory->grow(viscosity, atom->nmax, "atom:rheo_viscosity");
-  nmax_old = atom->nmax;
+    nmax_old = atom->nmax;
+  }
 
   if (viscosity_style == POWER) {
     for (i = 0; i < nlocal; i++) {
diff --git a/src/RHEO/pair_rheo.cpp b/src/RHEO/pair_rheo.cpp
index 7b9425cf20..1eb2a43e1a 100644
--- a/src/RHEO/pair_rheo.cpp
+++ b/src/RHEO/pair_rheo.cpp
@@ -91,17 +91,21 @@ void PairRHEO::compute(int eflag, int vflag)
   double **v = atom->v;
   double **x = atom->x;
   double **f = atom->f;
-  double **f_pressure = compute_interface->f_pressure;
   double *rho = atom->rho;
   double *mass = atom->mass;
   double *drho = atom->drho;
   double *temperature = atom->temperature;
   double *heatflow = atom->heatflow;
   double *special_lj = force->special_lj;
-  tagint *tag = atom->tag;
-  int *chi = compute_interface->chi;
   int *type = atom->type;
   int *status = atom->status;
+  tagint *tag = atom->tag;
+
+  double **fp_store, *chi;
+  if (compute_interface) {
+    fp_store = compute_interface->fp_store;
+    chi = compute_interface->chi;
+  }
 
   int tmp1, tmp2;
   int index = atom->find_custom("rheo_viscosity", tmp1, tmp2);
@@ -273,9 +277,6 @@ void PairRHEO::compute(int eflag, int vflag)
           f[i][0] += ft[0];
           f[i][1] += ft[1];
           f[i][2] += ft[2];
-          fp[i][0] += dfp[0];
-          fp[i][1] += dfp[1];
-          fp[i][2] += dfp[2];
 
           if (evflag) // Does not account for unbalanced forces
             ev_tally_xyz(i, j, nlocal, newton_pair, 0.0, 0.0, ft[0], ft[1], ft[2], dx[0], dx[1], dx[2]);
@@ -289,7 +290,7 @@ void PairRHEO::compute(int eflag, int vflag)
               // flip sign here b/c -= at accummulator
             }
 
-            scale3(fp_prefact,r dWji, dfp);
+            scale3(fp_prefactor, dWji, dfp);
 
             add3(fv, dfp, ft);
             add3(fsolid, ft, ft);
@@ -297,10 +298,18 @@ void PairRHEO::compute(int eflag, int vflag)
             f[j][0] -= ft[0];
             f[j][1] -= ft[1];
             f[j][2] -= ft[2];
+          }
 
-            fp[j][0] -= dfp[0];
-            fp[j][1] -= dfp[1];
-            fp[j][2] -= dfp[2];
+          if (compute_interface) {
+            fp_store[i][0] += dfp[0];
+            fp_store[i][1] += dfp[1];
+            fp_store[i][2] += dfp[2];
+
+            if (newton_pair || j < nlocal) {
+              fp_store[j][0] -= dfp[0];
+              fp_store[j][1] -= dfp[1];
+              fp_store[j][2] -= dfp[2];
+            }
           }
         }