diff --git a/doc/src/fix_rheo.rst b/doc/src/fix_rheo.rst
index 7aaab48f54..5eae617419 100644
--- a/doc/src/fix_rheo.rst
+++ b/doc/src/fix_rheo.rst
@@ -23,7 +23,7 @@ Syntax
        *surface/detection* values = *sdstyle* *limit*
          *sdstyle* = *coordination* or *divergence*
          *limit* = threshold for surface particles (unitless)
-       *interface/reconstruction* values = none, reconstructs interfaces with solid particles
+       *interface/reconstruct* values = none, reconstructs interfaces with solid particles
        *rho/sum* values = none, uses the kernel to compute the density of particles
        *density* values = *rho0* (density)
        *sound/squared* values = *csq* (velocity\^2)
diff --git a/src/RHEO/compute_rheo_grad.cpp b/src/RHEO/compute_rheo_grad.cpp
index d71e7a24ac..7d0b1d5b14 100644
--- a/src/RHEO/compute_rheo_grad.cpp
+++ b/src/RHEO/compute_rheo_grad.cpp
@@ -179,9 +179,6 @@ 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];
     fluidi = !(status[i] & PHASECHECK);
     jlist = firstneigh[i];
@@ -203,6 +200,10 @@ void ComputeRHEOGrad::compute_peratom()
         rhoi = rho[i];
         rhoj = rho[j];
 
+        vi[0] = v[i][0];
+        vi[1] = v[i][1];
+        vi[2] = v[i][2];
+
         vj[0] = v[j][0];
         vj[1] = v[j][1];
         vj[2] = v[j][2];
diff --git a/src/RHEO/compute_rheo_interface.cpp b/src/RHEO/compute_rheo_interface.cpp
index ea4916087f..8cd69b49e3 100644
--- a/src/RHEO/compute_rheo_interface.cpp
+++ b/src/RHEO/compute_rheo_interface.cpp
@@ -170,9 +170,10 @@ void ComputeRHEOInterface::compute_peratom()
           chi[i] += w;
         } else {
           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;
+            dot = (-fp_store[j][0] + fp_store[i][0]) * delx;
+            dot += (-fp_store[j][1] + fp_store[i][1]) * dely;
+            dot += (-fp_store[j][2] + fp_store[i][2]) * delz;
+
             rho[i] += w * (csq * (rho[j] - rho0) - rho[j] * dot);
             normwf[i] += w;
           }
@@ -184,9 +185,10 @@ void ComputeRHEOInterface::compute_peratom()
             chi[j] += w;
           } else {
             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;
+              dot = (-fp_store[i][0] + fp_store[j][0]) * delx;
+              dot += (-fp_store[i][1] + fp_store[j][1]) * dely;
+              dot += (-fp_store[i][2] + fp_store[j][2]) * delz;
+
               rho[j] += w * (csq * (rho[i] - rho0) + rho[i] * dot);
               normwf[j] += w;
             }
@@ -342,7 +344,7 @@ void ComputeRHEOInterface::store_forces()
   // If forces are overwritten by a fix, there are no pressure forces
   // so just normalize
   auto fixlist = modify->get_fix_by_style("setforce");
-  if (fixlist.size() == 0) {
+  if (fixlist.size() != 0) {
     for (const auto &fix : fixlist) {
       for (int i = 0; i < atom->nlocal; i++) {
         minv = 1.0 / mass[type[i]];
diff --git a/src/RHEO/fix_rheo.cpp b/src/RHEO/fix_rheo.cpp
index fd436cab6e..cecd3183fd 100644
--- a/src/RHEO/fix_rheo.cpp
+++ b/src/RHEO/fix_rheo.cpp
@@ -105,7 +105,7 @@ FixRHEO::FixRHEO(LAMMPS *lmp, int narg, char **arg) :
       }
 
       iarg += 2;
-    } else if (strcmp(arg[iarg],"interface/reconstruction") == 0) {
+    } else if (strcmp(arg[iarg],"interface/reconstruct") == 0) {
       interface_flag = 1;
     } else if (strcmp(arg[iarg],"rho/sum") == 0) {
       rhosum_flag = 1;
@@ -360,6 +360,7 @@ void FixRHEO::initial_integrate(int /*vflag*/)
         }
 
         if (!rhosum_flag) {
+          if (status[i] & PHASECHECK) continue;
           for (a = 0; a < dim; a++) {
             rho[i] += dtv * vshift[i][a] * gradr[i][a];
           }
@@ -380,7 +381,10 @@ void FixRHEO::pre_force(int /*vflag*/)
 
   compute_grad->forward_fields(); // also forwards v and rho for chi
   compute_kernel->compute_peratom();
-  if (interface_flag) compute_interface->compute_peratom();
+  if (interface_flag) {
+    // Note on first setup, have no forces for pressure to reference
+    compute_interface->compute_peratom();
+  }
 
   compute_grad->compute_peratom();
   compute_grad->forward_gradients();
diff --git a/src/RHEO/fix_rheo_pressure.cpp b/src/RHEO/fix_rheo_pressure.cpp
index ff206937f4..84d21ee872 100644
--- a/src/RHEO/fix_rheo_pressure.cpp
+++ b/src/RHEO/fix_rheo_pressure.cpp
@@ -122,20 +122,9 @@ void FixRHEOPressure::pre_force(int /*vflag*/)
 
   int nlocal = atom->nlocal;
 
-  for (i = 0; i < nlocal; i++) {
-    if (mask[i] & groupbit) {
-      if (pressure_style == LINEAR) {
-        pressure[i] = csq * (rho[i] - rho0);
-      } else if (pressure_style == CUBIC) {
-        dr = rho[i] - rho0;
-        pressure[i] = csq * (dr + c_cubic * dr * dr * dr);
-      } else if (pressure_style == TAITWATER) {
-        rho_ratio = rho[i] / rho0inv;
-        rr3 = rho_ratio * rho_ratio * rho_ratio;
-        pressure[i] = csq * rho0 * SEVENTH * (rr3 * rr3 * rho_ratio - 1.0);
-      }
-    }
-  }
+  for (i = 0; i < nlocal; i++)
+    if (mask[i] & groupbit)
+      pressure[i] = calc_pressure(rho[i]);
 
   if (comm_forward) comm->forward_comm(this);
 }
@@ -186,5 +175,5 @@ double FixRHEOPressure::calc_pressure(double rho)
     rr3 = rho_ratio * rho_ratio * rho_ratio;
     p = csq * rho0 * SEVENTH * (rr3 * rr3 * rho_ratio - 1.0);
   }
-  return rho;
+  return p;
 }
diff --git a/src/RHEO/pair_rheo.cpp b/src/RHEO/pair_rheo.cpp
index 8b0e2c2df6..0feb6af445 100644
--- a/src/RHEO/pair_rheo.cpp
+++ b/src/RHEO/pair_rheo.cpp
@@ -56,6 +56,8 @@ PairRHEO::PairRHEO(LAMMPS *lmp) :
   artificial_visc_flag = 0;
   rho_damp_flag = 0;
   thermal_flag = 0;
+
+  comm_reverse = 3;
 }
 
 /* ---------------------------------------------------------------------- */
@@ -79,7 +81,7 @@ void PairRHEO::compute(int eflag, int vflag)
   double rhoi, rhoj, voli, volj, Pi, Pj, etai, etaj, kappai, kappaj;
   double mu, q, fp_prefactor, drho_damp, fmag, psi_ij, Fij;
   double *dWij, *dWji, *dW1ij, *dW1ji;
-  double dx[3], du[3], fv[3], dfp[3], fsolid[3], ft[3], vi[3], vj[3];
+  double dx[3], du[3], dv[3], fv[3], dfp[3], fsolid[3], ft[3], vi[3], vj[3];
 
   int *ilist, *jlist, *numneigh, **firstneigh;
   double imass, jmass, rsq, r, rinv;
@@ -113,6 +115,12 @@ void PairRHEO::compute(int eflag, int vflag)
   if (compute_interface) {
     fp_store = compute_interface->fp_store;
     chi = compute_interface->chi;
+
+    for (i = 0; i < atom->nmax; i++) {
+      fp_store[i][0] = 0.0;
+      fp_store[i][1] = 0.0;
+      fp_store[i][2] = 0.0;
+    }
   }
 
   inum = list->inum;
@@ -243,11 +251,12 @@ void PairRHEO::compute(int eflag, int vflag)
 
           //Hydrostatic pressure forces
           fp_prefactor = voli * volj * (Pj + Pi);
-          sub3(vi, vj, du);
+          sub3(vi, vj, dv);
 
           //Add artificial viscous pressure if required
           if (artificial_visc_flag && pair_avisc_flag) {
             //Interpolate velocities to midpoint and use this difference for artificial viscosity
+            copy3(dv, du);
             for (a = 0; a < dim; a++)
               for (b = 0; b < dim; b++)
                 du[a] -= 0.5 * (gradv[i][a * dim + b] + gradv[j][a * dim + b]) * dx[b];
@@ -257,10 +266,6 @@ void PairRHEO::compute(int eflag, int vflag)
             mu = MIN(0.0, mu);
             q = av * (-2.0 * cs * mu + mu * mu);
             fp_prefactor += voli * volj * q * (rhoj + rhoi);
-
-  if (fabs(fp_prefactor*dWij[0]) > 1e-9)
-            if (atom->tag[i] == 643 or atom->tag[j] == 643 or atom->tag[i] == 963 or atom->tag[j] == 963)
-              printf("%d-%d (%d %d) fp_prefactor %g %g %g\n", atom->tag[i], atom->tag[j], i, j, fp_prefactor, dWij[0], -fp_prefactor*dWij[0]);
           }
 
           // -Grad[P + Q]
@@ -270,7 +275,7 @@ void PairRHEO::compute(int eflag, int vflag)
           for (a = 0; a < dim; a++) {
             fv[a] = 0.0;
             for (b = 0; b < dim; b++)
-              fv[a] += du[a] * dx[b] * dWij[b];
+              fv[a] += dv[a] * dx[b] * dWij[b];
             fv[a] *= (etai + etaj) * voli * volj * rinv * rinv;
           }
 
@@ -346,6 +351,11 @@ void PairRHEO::compute(int eflag, int vflag)
   }
 
   if (vflag_fdotr) virial_fdotr_compute();
+
+  if (compute_interface) {
+    comm->reverse_comm(this);
+    comm->forward_comm(this);
+  }
 }
 
 /* ----------------------------------------------------------------------
@@ -476,3 +486,37 @@ double PairRHEO::init_one(int i, int j)
 
   return h;
 }
+
+/* ---------------------------------------------------------------------- */
+
+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++) {
+    buf[m++] = fp_store[i][0];
+    buf[m++] = fp_store[i][1];
+    buf[m++] = fp_store[i][2];
+  }
+
+  return m;
+}
+
+/* ---------------------------------------------------------------------- */
+
+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];
+    fp_store[j][0] += buf[m++];
+    fp_store[j][1] += buf[m++];
+    fp_store[j][2] += buf[m++];
+  }
+}
diff --git a/src/RHEO/pair_rheo.h b/src/RHEO/pair_rheo.h
index b30b0c3c04..cb2227c8d6 100644
--- a/src/RHEO/pair_rheo.h
+++ b/src/RHEO/pair_rheo.h
@@ -33,6 +33,8 @@ class PairRHEO : public Pair {
   void coeff(int, char **) override;
   void setup() override;
   double init_one(int, int) override;
+  int pack_reverse_comm(int, int, double *) override;
+  void unpack_reverse_comm(int, int *, double *) override;
 
  protected:
   double h, csq, rho0;        // From fix RHEO