diff --git a/src/RHEO/compute_rheo_grad.h b/src/RHEO/compute_rheo_grad.h
index 7804ee72dd..ee4b4a5bd6 100644
--- a/src/RHEO/compute_rheo_grad.h
+++ b/src/RHEO/compute_rheo_grad.h
@@ -13,7 +13,7 @@
 
 #ifdef COMPUTE_CLASS
 // clang-format off
-ComputeStyle(rheo/grad,ComputeRHEOGrad)
+ComputeStyle(RHEO/GRAD,ComputeRHEOGrad)
 // clang-format on
 #else
 
diff --git a/src/RHEO/compute_rheo_interface.cpp b/src/RHEO/compute_rheo_interface.cpp
index ff76b75f77..55cb81382b 100644
--- a/src/RHEO/compute_rheo_interface.cpp
+++ b/src/RHEO/compute_rheo_interface.cpp
@@ -398,7 +398,7 @@ void ComputeRHEOInterface::store_forces()
 
 double ComputeRHEOInterface::memory_usage()
 {
-  double bytes = nmax * sizeof(double);
+  double bytes = 3 * nmax_old * sizeof(double);
   return bytes;
 }
 
diff --git a/src/RHEO/compute_rheo_interface.h b/src/RHEO/compute_rheo_interface.h
index 8e190c5430..cdb2eb6c54 100644
--- a/src/RHEO/compute_rheo_interface.h
+++ b/src/RHEO/compute_rheo_interface.h
@@ -13,7 +13,7 @@
 
 #ifdef COMPUTE_CLASS
 // clang-format off
-ComputeStyle(rheo/interface,ComputeRHEOInterface)
+ComputeStyle(RHEO/INTERFACE,ComputeRHEOInterface)
 // clang-format on
 #else
 
diff --git a/src/RHEO/compute_rheo_kernel.h b/src/RHEO/compute_rheo_kernel.h
index df659c47de..ed2b6ff5f2 100644
--- a/src/RHEO/compute_rheo_kernel.h
+++ b/src/RHEO/compute_rheo_kernel.h
@@ -13,7 +13,7 @@
 
 #ifdef COMPUTE_CLASS
 // clang-format off
-ComputeStyle(rheo/kernel,ComputeRHEOKernel)
+ComputeStyle(RHEO/KERNEL,ComputeRHEOKernel)
 // clang-format on
 #else
 
diff --git a/src/RHEO/compute_rheo_surface.cpp b/src/RHEO/compute_rheo_surface.cpp
new file mode 100644
index 0000000000..dff5db4cfa
--- /dev/null
+++ b/src/RHEO/compute_rheo_surface.cpp
@@ -0,0 +1,595 @@
+/* ----------------------------------------------------------------------
+ LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+ https://www.lammps.org/, Sandia National Laboratories
+ LAMMPS development team: developers@lammps.org
+
+ Copyright (2003) Sandia Corporation.  Under the terms of Contract
+ DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+ certain rights in this software.  This software is distributed under
+ the GNU General Public License.
+
+ See the README file in the top-level LAMMPS directory.
+ ------------------------------------------------------------------------- */
+
+/* ----------------------------------------------------------------------
+   Contributing authors:
+   Joel Clemmer (SNL), Thomas O'Connor (CMU), Eric Palermo (CMU)
+----------------------------------------------------------------------- */
+
+#include "compute_rheo_surface.h"
+
+#include "fix_rheo.h"
+#include "compute_rheo_kernel.h"
+#include "compute_rheo_solids.h"
+#include "atom.h"
+#include "memory.h"
+#include "atom.h"
+#include "comm.h"
+#include "modify.h"
+#include "neighbor.h"
+#include "neigh_list.h"
+#include "neigh_request.h"
+#include "error.h"
+#include "force.h"
+#include "domain.h"
+
+using namespace LAMMPS_NS;
+using namespace FixConst;
+
+/* ---------------------------------------------------------------------- */
+
+ComputeRHEOSurface::ComputeRHEOSurface(LAMMPS *lmp, int narg, char **arg) :
+  Fix(lmp, narg, arg)
+{
+  if (narg < 6) error->all(FLERR,"Illegal fix rheo/surface command");
+
+  cut = utils::numeric(FLERR,arg[3],false,lmp);
+  divR_limit = utils::numeric(FLERR,arg[4],false,lmp);
+  coord_limit = utils::inumeric(FLERR,arg[5],false,lmp);
+
+  divr_flag = 1;
+  if (narg == 7) {
+    divr_flag = 0;
+  }
+
+  int dim = domain->dimension;
+
+  peratom_flag = 1;
+  size_peratom_cols = dim;
+  peratom_freq = 1;
+
+
+  comm_forward = 2;
+  comm_reverse = dim*dim + 1;
+
+  cutsq = cut*cut;
+
+  B = nullptr;
+  gradC = nullptr;
+  n_surface = nullptr;
+
+  int nall = atom->nlocal + atom->nghost;
+  nmax = nall;
+  memory->create(B,nmax,dim*dim,"fix/rheo/surface:B");
+  memory->create(gradC,nmax,dim*dim,"fix/rheo/surface:gradC");
+  memory->create(n_surface,nmax,dim,"fix/rheo/surface:B");
+  array_atom = n_surface;
+
+  compute_kernel = nullptr;
+  compute_solids = NULL;
+  fix_rheo = nullptr;
+}
+
+/* ---------------------------------------------------------------------- */
+
+ComputeRHEOSurface::~ComputeRHEOSurface()
+{
+  if (modify->nfix) modify->delete_fix("PROPERTY_ATOM_RHEO_SURFACE");
+
+  memory->destroy(B);
+  memory->destroy(gradC);
+  memory->destroy(n_surface);
+}
+
+void ComputeRHEOSurface::post_constructor()
+{
+  //Store persistent per atom quantities
+  char **fixarg = new char*[5];
+  fixarg[0] = (char *) "PROPERTY_ATOM_RHEO_SURFACE";
+  fixarg[1] = (char *) "all";
+  fixarg[2] = (char *) "property/atom";
+  fixarg[3] = (char *) "d_divr";
+  fixarg[4] = (char *) "d_rsurf";
+  modify->add_fix(5,fixarg,1);
+
+  int temp_flag;
+  index_divr = atom->find_custom("divr", temp_flag);
+  if ((index_divr < 0) || (temp_flag != 1))
+      error->all(FLERR, "Pair rheo/surface can't find fix property/atom divr");
+
+  index_rsurf = atom->find_custom("rsurf", temp_flag);
+  if ((index_rsurf < 0) || (temp_flag != 1))
+      error->all(FLERR, "Pair rheo/surface can't find fix property/atom rsurf");
+
+  delete [] fixarg;
+  divr = atom->dvector[index_divr];
+}
+/* ---------------------------------------------------------------------- */
+
+int ComputeRHEOSurface::setmask()
+{
+  int mask = 0;
+  mask |= PRE_FORCE;
+  return mask;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void ComputeRHEOSurface::init()
+{
+  // need an occasional full neighbor list
+  int irequest = neighbor->request(this,instance_me);
+  neighbor->requests[irequest]->pair = 0;
+  neighbor->requests[irequest]->fix = 1;
+  neighbor->requests[irequest]->half = 1;
+  neighbor->requests[irequest]->full = 0;
+
+  int flag;
+  int ifix = modify->find_fix_by_style("rheo");
+  if (ifix == -1) error->all(FLERR, "Need to define fix rheo to use fix rheo/surface");
+  fix_rheo = ((FixRHEO *) modify->fix[ifix]);
+  compute_kernel = fix_rheo->compute_kernel;
+  compute_solids = fix_rheo->compute_solids;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void ComputeRHEOSurface::setup_pre_force(int /*vflag*/)
+{
+  pre_force(0);
+}
+
+/* ---------------------------------------------------------------------- */
+
+void ComputeRHEOSurface::init_list(int /*id*/, NeighList *ptr)
+{
+  list = ptr;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void ComputeRHEOSurface::pre_force(int /*vflag*/)
+{
+  int i, j, ii, jj, jnum, a, b, itype, jtype;
+  double xtmp, ytmp, ztmp, delx, dely, delz, rsq, r, wp, Voli, Volj, rhoi, rhoj;
+  int *jlist;
+  double *dWij, *dWji;
+  double dx[3];
+
+  divr = atom->dvector[index_divr];
+
+  // neighbor list variables
+  int inum, *ilist, *numneigh, **firstneigh;
+  int nlocal = atom->nlocal;
+  int nall = nlocal + atom->nghost;
+
+  double **x = atom->x;
+  int *surface = atom->surface;
+  int *phase = atom->phase;
+  double *rsurf = atom->dvector[index_rsurf];
+  int newton = force->newton;
+  int dim = domain->dimension;
+  int *mask = atom->mask;
+  int *type = atom->type;
+  double *mass = atom->mass;
+  double *rho = atom->rho;
+  double *temp = atom->temp;
+
+  inum = list->inum;
+  ilist = list->ilist;
+  numneigh = list->numneigh;
+  firstneigh = list->firstneigh;
+
+  if (nmax <= nall) {
+    nmax = nall;
+    memory->destroy(B);
+    memory->destroy(gradC);
+    memory->destroy(n_surface);
+
+    memory->create(B,nmax,dim*dim,"fix/rheo/surface:B");
+    memory->create(gradC,nmax,dim*dim,"fix/rheo/surface:gradC");
+    memory->create(n_surface,nmax,dim,"fix/rheo/surface:n_surface");
+    array_atom = n_surface;
+  }
+
+  for (i = 0; i < nall; i++) {
+    for (a = 0; a < dim; a++) {
+      for (b = 0; b < dim; b++) {
+        B[i][a*dim + b] = 0.0;
+        gradC[i][a*dim + b] = 0.0;
+      }
+      n_surface[i][a] = 0.0;
+    }
+    divr[i] = 0.0;
+    surface[i] = 0;
+  }
+
+  // loop over neighbors to calculate the average orientation of neighbors
+  for (ii = 0; ii < inum; ii++) {
+    i = ilist[ii];
+    xtmp = x[i][0];
+    ytmp = x[i][1];
+    ztmp = x[i][2];
+
+    jlist = firstneigh[i];
+    jnum = numneigh[i];
+    itype = type[i];
+
+    for (jj = 0; jj < jnum; jj++) {
+      j = jlist[jj];
+      j &= NEIGHMASK;
+
+      delx = xtmp - x[j][0];
+      dely = ytmp - x[j][1];
+      delz = ztmp - x[j][2];
+      dx[0] = delx;
+      dx[1] = dely;
+      dx[2] = delz;
+      rsq = delx * delx + dely * dely + delz * delz;
+      if (rsq < cutsq) {
+        jtype = type[j];
+
+        rhoi = rho[i];
+        rhoj = rho[j];
+
+        // Add corrections for walls
+        if (phase[i] <= FixRHEO::FLUID_MAX && phase[j] > FixRHEO::FLUID_MAX) {
+          rhoj = compute_solids->correct_rho(j,i);
+        } else if (phase[i] > FixRHEO::FLUID_MAX && phase[j] <= FixRHEO::FLUID_MAX) {
+          rhoi = compute_solids->correct_rho(i,j);
+        } else if (phase[i] > FixRHEO::FLUID_MAX && phase[j] > FixRHEO::FLUID_MAX) {
+          rhoi = 1.0;
+          rhoj = 1.0;
+        }
+
+        Voli = mass[itype]/rhoi;
+        Volj = mass[jtype]/rhoj;
+
+        //compute kernel gradient
+        wp = compute_kernel->calc_dw_quintic(i, j, delx, dely, delz, sqrt(rsq),compute_kernel->dWij,compute_kernel->dWji);
+        //wp = compute_kernel->calc_dw(i, j, delx, dely, delz, sqrt(rsq));//,compute_kernel->dWij,compute_kernel->dWji);
+
+        dWij = compute_kernel->dWij;
+        dWji = compute_kernel->dWji;
+
+        for (a=0; a<dim; a++){
+          divr[i] -= dWij[a]*dx[a]*Volj; // dx = xi-xj = xji = -xij
+          gradC[i][a] +=  dWij[a]*Volj;
+        }
+
+        if (j < nlocal || newton) {
+          for (a=0; a<dim; a++){
+            divr[j] += dWji[a]*dx[a]*Voli;
+            gradC[j][a] +=  dWji[a]*Voli;
+          }
+        }
+      }
+    }
+  }
+
+  comm_stage = 0;
+  comm_reverse = dim*dim + 1;   // gradC and divr
+  comm_forward = 1;             // divr
+  if (newton) comm->reverse_comm_fix(this);
+  comm->forward_comm_fix(this);
+
+  int *coordination = compute_kernel->coordination;
+  // Find the free-surface
+  //0-bulk 1-surf vicinity  2-surface 3-splash
+  if (divr_flag) {
+    for (i = 0; i < nall; i++) {
+      if (mask[i] & groupbit) {
+        surface[i] = 0;
+        rsurf[i] = cut; //Maximum range that can be seen
+        if (divr[i] < divR_limit) {
+          surface[i] = 2;
+          rsurf[i] = 0.0;
+          if (coordination[i] < coord_limit) surface[i] = 3;
+        }
+      }
+    }
+  } else {
+    for (i = 0; i < nall; i++) {
+      if (mask[i] & groupbit) {
+        surface[i] = 0;
+        rsurf[i] = cut; //Maximum range that can be seen
+        if (coordination[i] < divR_limit) {
+          surface[i] = 2;
+          rsurf[i] = 0.0;
+          if (coordination[i] < coord_limit) surface[i] = 3;
+        }
+      }
+    }
+  }
+
+  //comm_stage = 1;
+  //comm_forward = 1;
+  //comm->forward_comm_fix(this);  // communicate free surface particles
+
+  for (ii = 0; ii < inum; ii++) {
+    i = ilist[ii];
+    xtmp = x[i][0];
+    ytmp = x[i][1];
+    ztmp = x[i][2];
+
+    jlist = firstneigh[i];
+    jnum = numneigh[i];
+
+    for (jj = 0; jj < jnum; jj++) {
+      j = jlist[jj];
+      j &= NEIGHMASK;
+
+      delx = xtmp - x[j][0];
+      dely = ytmp - x[j][1];
+      delz = ztmp - x[j][2];
+      rsq = delx * delx + dely * dely + delz * delz;
+      if (rsq < cutsq) {
+        r = sqrt(rsq);
+        if (surface[i] == 0 && surface[j] == 2) surface[i] = 1;
+        if (surface[j] == 0 && surface[i] == 2) surface[j] = 1;
+        if (surface[j] == 2) rsurf[i] = MIN(rsurf[i], r);
+        if (surface[i] == 2) rsurf[j] = MIN(rsurf[j], r);
+      }
+    }
+  }
+
+  comm_stage = 1;
+  comm_reverse = 2;
+  comm_forward = 2;
+  if (newton) comm->reverse_comm_fix(this);
+  comm->forward_comm_fix(this);
+
+  //Now loop again and for each surface particle (2)
+  // find its neighbors that are bulk (0) and convert to surface vicinity (1)
+  // if the surface particle has no (0) or (1) neighbors then it is a spash (3)
+
+  //for (ii = 0; ii < inum; ii++) {  // is this the right i and j loop for this?
+  //  i = ilist[ii];
+  //
+  //  if (surface[i]!=2) continue; //Only consider surface particles
+  //
+  //  bool nobulkneigh = true; // whether we have no bulk neighbors
+  //  xtmp = x[i][0];
+  //  ytmp = x[i][1];
+  //  ztmp = x[i][2];
+  //  jlist = firstneigh[i];
+  //  jnum = numneigh[i];
+  //
+  //  for (jj = 0; jj < jnum; jj++) {
+  //    j = jlist[jj];
+  //    j &= NEIGHMASK;
+  //
+  //    //other surface or splash neighbors do not need labeling
+  //    if (surface[j]>=2){
+  //      continue;
+  //    }
+  //
+  //    //check distance criterion rij < h = cutsq/9 for quintic kernel
+  //    delx = xtmp - x[j][0];
+  //    dely = ytmp - x[j][1];
+  //    delz = ztmp - x[j][2];
+  //    dx[0] = 3.0*delx;   // multiplied by three here to make criterion r<h instead of r<3*h
+  //    dx[1] = 3.0*dely;
+  //    dx[2] = 3.0*delz;
+  //    rsq = delx * delx + dely * dely + delz * delz;
+  //    if (rsq < cutsq) {
+  //      //We have identified 1 bulk fluid neighbor
+  //      nobulkneigh = false;
+  //      //that bulk fluid neighbor is in the vicinity of hte surface
+  //      surface[j] = 1;
+  //    }
+  //  }
+  //  if (nobulkneigh){
+  //    surface[i] = 3;
+  //  }
+  //}
+//
+//  //Reverse comm surface?
+//
+//  // loop over neighbors to calculate the average orientation
+//  // skip for bulk or splash
+//  for (ii = 0; ii < inum; ii++) {
+//    i = ilist[ii];
+//    if ((surface[i]==0)||(surface[i]==3)){
+//      continue;
+//    }
+//
+//    itype = type[i];
+//    rhoi = rho[i];
+//    Voli = mass[itype]/rhoi;
+//
+//    xtmp = x[i][0];
+//    ytmp = x[i][1];
+//    ztmp = x[i][2];
+//
+//    jlist = firstneigh[i];
+//    jnum = numneigh[i];
+//    for (jj = 0; jj < jnum; jj++) {
+//      j = jlist[jj];
+//      j &= NEIGHMASK;
+//
+//      delx = xtmp - x[j][0];
+//      dely = ytmp - x[j][1];
+//      delz = ztmp - x[j][2];
+//      dx[0] = delx;   // multiplied by three here to make criterion r<h instead of r<3*h
+//      dx[1] = dely;
+//      dx[2] = delz;
+//      rsq = delx * delx + dely * dely + delz * delz;
+//      if (rsq < cutsq) {
+//
+//        jtype = type[j];
+//        rhoj = rho[j];
+//        Volj = mass[jtype]/rhoj;
+//
+//        for (a=0; a<dim; a++){
+//          for (b=0; b<dim; b++){
+//            B[i][a*dim+b] -= dx[a]*dWij[b]*Volj;
+//          }
+//        }
+//
+//        if (j < nlocal || newton) {
+//          for (a=0; a<dim; a++){
+//            for (b=0; b<dim; b++){
+//              B[j][a*dim+b] += dx[a]*dWji[b]*Voli;
+//            }
+//          }
+//        }
+//      }
+//    }
+//  }
+//  //reverse comm to populate B[j] if Newton is on
+//  comm_stage = 2;
+//  comm_reverse = dim*dim;        // B
+//  if (newton) comm->reverse_comm_fix(this);
+//
+//
+//  // Now need to invert each B
+//  int status, s;
+//  //LU requires a permuation matrix
+//  gsl_permutation * p = gsl_permutation_alloc(dim);
+//  for (ii = 0; ii < inum; ii++) {
+//    i = ilist[ii];
+//    if ((surface[i]==0)||(surface[i]==3)){
+//      continue;
+//    }
+//
+//    //Use gsl to get Binv
+//    //B is not symmteric so we will use a LU decomp
+//    gsl_matrix_view gB = gsl_matrix_view_array(B[i],dim,dim);
+//    status = 0;
+//    status = gsl_linalg_LU_decomp(&gB.matrix,p,&s);  //B[i] is now the LU decomp
+//    // check if decomposition failure
+//    if (status) {
+//      fprintf(stderr, "failed, gsl_errno=%d.n", status);
+//      continue;
+//    } else {
+//      gsl_linalg_LU_invx(&gB.matrix,p); //B[i] is now inv(B[i])
+//    }
+//  }
+//  gsl_permutation_free(p);
+  double maggC = 0.0;
+  for (i = 0; i < nlocal; i++) {
+    if (mask[i] & groupbit) {
+      maggC=0;
+      for (a=0;a<dim;a++){
+        maggC += gradC[i][a]*gradC[i][a];
+      }
+      maggC = sqrt(maggC) + 1e-10;
+      for (a=0;a<dim;a++){
+        n_surface[i][a] = -gradC[i][a]/maggC;
+      }//dr can then be calculated by fix vshift
+    }
+  }
+}
+
+/* ---------------------------------------------------------------------- */
+
+int ComputeRHEOSurface::pack_reverse_comm(int n, int first, double *buf)
+{
+  int i,a,b,k,m,last;
+  int dim = domain->dimension;
+  int *surface = atom->surface;
+  double *rsurf = atom->dvector[index_rsurf];
+
+  m = 0;
+  last = first + n;
+  for (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 ++)
+          buf[m++] = gradC[i][a*dim + b];
+    } else if (comm_stage == 1) {
+      buf[m++] = (double) surface[i];
+      buf[m++] = rsurf[i];
+    } else if (comm_stage == 2) {
+      for (a = 0; a < dim; a ++ )
+        for (b = 0; b < dim; b ++)
+          buf[m++] = B[i][a*dim + b];
+    }
+  }
+  return m;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void ComputeRHEOSurface::unpack_reverse_comm(int n, int *list, double *buf)
+{
+  int i,a,b,k,j,m;
+  int dim = domain->dimension;
+  int *surface = atom->surface;
+  double *rsurf = atom->dvector[index_rsurf];
+
+  m = 0;
+  for (i = 0; i < n; i++) {
+    j = list[i];
+    if (comm_stage == 0) {
+      divr[j] += buf[m++];
+      for (a = 0; a < dim; a ++ )
+        for (b = 0; b < dim; b ++)
+          gradC[j][a*dim + b] += buf[m++];
+    } else if (comm_stage == 1) {
+      int temp = (int) buf[m++];
+      surface[j] = MAX(surface[j], temp);
+      double temp2 = buf[m++];
+      rsurf[j] = MIN(rsurf[j], temp2);
+    } else if (comm_stage == 2) {
+      for (a = 0; a < dim; a ++ )
+        for (b = 0; b < dim; b ++)
+          B[j][a*dim + b] += buf[m++];
+    }
+  }
+}
+
+
+/* ---------------------------------------------------------------------- */
+
+int ComputeRHEOSurface::pack_forward_comm(int n, int *list, double *buf,
+                                        int /*pbc_flag*/, int * /*pbc*/)
+{
+  int i,j,a,b,k,m;
+  int *surface = atom->surface;
+  double *rsurf = atom->dvector[index_rsurf];
+  m = 0;
+
+  for (i = 0; i < n; i++) {
+    j = list[i];
+    if (comm_stage == 0) {
+      buf[m++] = divr[j];
+    } else if (comm_stage == 1) {
+      buf[m++] = (double) surface[j];
+      buf[m++] = rsurf[j];
+    }
+  }
+  return m;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void ComputeRHEOSurface::unpack_forward_comm(int n, int first, double *buf)
+{
+  int i, k, a, b, m, last;
+  int *surface = atom->surface;
+  double *rsurf = atom->dvector[index_rsurf];
+
+  m = 0;
+  last = first + n;
+  for (i = first; i < last; i++) {
+    if (comm_stage == 0) {
+      divr[i] = buf[m++];
+    } else if (comm_stage == 1) {
+      surface[i] = (int) buf[m++];
+      rsurf[i] = buf[m++];
+    }
+  }
+}
diff --git a/src/RHEO/compute_rheo_surface.h b/src/RHEO/compute_rheo_surface.h
new file mode 100644
index 0000000000..480a8cd8ec
--- /dev/null
+++ b/src/RHEO/compute_rheo_surface.h
@@ -0,0 +1,72 @@
+/* -*- c++ -*- ----------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   https://www.lammps.org/, Sandia National Laboratories
+   LAMMPS development team: developers@lammps.org
+
+   Copyright (2003) Sandia Corporation.  Under the terms of Contract
+   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+   certain rights in this software.  This software is distributed under
+   the GNU General Public License.
+
+   See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+#ifdef COMPUTE_CLASS
+// clang-format off
+ComputeStyle(RHEO/SURFACE,ComputeRHEOSurface)
+// clang-format on
+#else
+
+#ifndef LMP_COMPUTE_RHEO_INTERFACE_H
+#define LMP_COMPUTE_RHEO_INTERFACE_H
+
+#include "compute.h"
+
+namespace LAMMPS_NS {
+
+class ComputeRHEOSurface : public Compute {
+ public:
+  ComputeRHEOSurface(class LAMMPS *, int, char **);
+  ~ComputeRHEOSurface();
+  void init() override;
+  void init_list(int, class NeighList *) override;
+  void compute_peratom() override;
+  int pack_reverse_comm(int, int, double *) override;
+  void unpack_reverse_comm(int, int *, double *) override;
+  int pack_forward_comm(int, int *, double *, int, int *) override;
+  void unpack_forward_comm(int, int, double *) override;
+
+  double **gradC, **n_surface;
+
+ private:
+  double cut, cutsq, threshold;
+  int surface_style, nmax_old;
+  double **B, *divr;
+  int comm_stage;
+
+  int index_divr;
+  int index_rsurf;
+
+  double divR_limit;
+  int coord_limit;
+
+  class NeighList *list;
+  class FixRHEO *fix_rheo;
+  class ComputeRHEOKernel *compute_kernel;
+  class ComputeRHEOSolids *compute_solids;
+};
+
+}
+
+#endif
+#endif
+
+/* ERROR/WARNING messages:
+
+E: Illegal ... command
+
+Self-explanatory.  Check the input script syntax and compare to the
+documentation for the command.  You can use -echo screen as a
+command-line option when running LAMMPS to see the offending line.
+
+*/
diff --git a/src/RHEO/compute_rheo_vshift.h b/src/RHEO/compute_rheo_vshift.h
index af5675fd8d..e76476e7fd 100644
--- a/src/RHEO/compute_rheo_vshift.h
+++ b/src/RHEO/compute_rheo_vshift.h
@@ -13,7 +13,7 @@
 
 #ifdef COMPUTE_CLASS
 // clang-format off
-ComputeStyle(rheo/vshift,ComputeRHEOVShift)
+ComputeStyle(RHEO/VSHIFT,ComputeRHEOVShift)
 // clang-format on
 #else
 
diff --git a/src/RHEO/fix_rheo.cpp b/src/RHEO/fix_rheo.cpp
index e83a95751a..8fc48c0b3b 100644
--- a/src/RHEO/fix_rheo.cpp
+++ b/src/RHEO/fix_rheo.cpp
@@ -88,12 +88,24 @@ FixRHEO::FixRHEO(LAMMPS *lmp, int narg, char **arg) :
       thermal_flag = 1;
     } else if (strcmp(arg[iarg],"surface/detection") == 0) {
       surface_flag = 1;
+      if(iarg + 2 >= narg) error->all(FLERR,"Illegal surface/detection option in fix rheo");
+      if (strcmp(arg[iarg + 1], "coordination")) {
+        surface_style = COORDINATION;
+        zmin_surface = utils::inumeric(FLERR,arg[iarg + 2],false,lmp);
+      } else if (strcmp(arg[iarg + 1], "divergence")) {
+        surface_style = DIVR;
+        divr_surface = utils::numeric(FLERR,arg[iarg + 2],false,lmp);
+      } else {
+        error->all(FLERR,"Illegal surface/detection option in fix rheo, {}", arg[iarg + 1]);
+      }
+
+      iarg += 2;
     } else if (strcmp(arg[iarg],"interface/reconstruction") == 0) {
       interface_flag = 1;
     } else if (strcmp(arg[iarg],"rhosum") == 0) {
       rhosum_flag = 1;
       if(iarg + 1 >= narg) error->all(FLERR,"Illegal rhosum option in fix rheo");
-      zmin_rhosum = utils::inumeric(FLERR,arg[iarg + 1],false,lmp);
+      rhosum_zmin = utils::inumeric(FLERR,arg[iarg + 1],false,lmp);
       iarg += 1;
     } else if (strcmp(arg[iarg],"rho0") == 0) {
       if(iarg + 1 >= narg) error->all(FLERR,"Illegal rho0 option in fix rheo");
@@ -117,6 +129,7 @@ FixRHEO::~FixRHEO()
   if (compute_kernel) modify->delete_compute("rheo_kernel");
   if (compute_grad) modify->delete_compute("rheo_grad");
   if (compute_interface) modify->delete_compute("rheo_interface");
+  if (compute_surface) modify->delete_compute("compute_surface");
   if (compute_rhosum) modify->delete_compute("rheo_rhosum");
   if (compute_vshift) modify->delete_compute("rheo_vshift");
 }
@@ -128,28 +141,33 @@ FixRHEO::~FixRHEO()
 
 void FixRHEO::post_constructor()
 {
-  compute_kernel = dynamic_cast<ComputeRHEOKernel *>(modify->add_compute("rheo_kernel all rheo/kernel"));
+  compute_kernel = dynamic_cast<ComputeRHEOKernel *>(modify->add_compute("rheo_kernel all RHEO/KERNEL"));
   compute_kernel->fix_rheo = this;
 
-  std::string cmd = "rheo_grad all rheo/grad velocity rho viscosity";
+  std::string cmd = "rheo_grad all RHEO/GRAD velocity rho viscosity";
   if (thermal_flag) cmd += "temperature";
   compute_grad = dynamic_cast<ComputeRHEOGrad *>(modify->add_compute(cmd));
   compute_grad->fix_rheo = this;
 
   if (rhosum_flag) {
-    compute_rhosum = dynamic_cast<ComputeRHEORhoSum *>(modify->add_compute("rheo_rhosum all rheo/rho/sum"));
+    compute_rhosum = dynamic_cast<ComputeRHEORhoSum *>(modify->add_compute("rheo_rhosum 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(fmt::format("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->fix_rheo = this;
+  }
 }
 
 /* ---------------------------------------------------------------------- */
@@ -180,10 +198,11 @@ 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()
-  if (viscosity_fix_defined || pressure_fix_defined || thermal_fix_defined || surface_fix_defined)
+  if (viscosity_fix_defined || pressure_fix_defined || thermal_fix_defined)
     error->all(FLERR, "Fix RHEO must be defined before all other RHEO fixes");
 
-  pre_force(0);
+  // Calculate surfaces
+  if (surface_flag) compute_surface->compute_peratom();
 }
 
 /* ---------------------------------------------------------------------- */
@@ -201,14 +220,10 @@ void FixRHEO::setup()
   if(!thermal_fix_defined && thermal_flag)
     error->all(FLERR, "Missing fix rheo/thermal");
 
-  if(!surface_fix_defined && surface_flag)
-    error->all(FLERR, "Missing fix rheo/surface");
-
   // Reset to zero for next run
   thermal_fix_defined = 0;
   viscosity_fix_defined = 0;
   pressure_fix_defined = 0;
-  surface_fix_defined = 0;
 
   // Check fixes cover all atoms (doesnt ensure user covers atoms created midrun)
   // (pressure is currently required to be group all)
@@ -259,7 +274,8 @@ void FixRHEO::initial_integrate(int /*vflag*/)
 
   double **gradr = compute_grad->gradr;
   double **gradv = compute_grad->gradv;
-  double **vshift = compute_vshift->vshift;
+  double **vshift;
+  if (shift_flag) compute_vshift->vshift;
 
   int *type = atom->type;
   int *mask = atom->mask;
@@ -287,8 +303,11 @@ void FixRHEO::initial_integrate(int /*vflag*/)
 
   // Update gradients and interpolate solid properties
   compute_grad->forward_fields(); // also forwards v and rho for chi
-  compute_interface->store_forces(); // Need to save, wiped in exchange
-  compute_interface->compute_peratom();
+  if (interface_flag) {
+    // Need to save, wiped in exchange
+    compute_interface->store_forces();
+    compute_interface->compute_peratom();
+  }
   compute_grad->compute_peratom();
 
   // Position half-step
@@ -350,7 +369,7 @@ void FixRHEO::pre_force(int /*vflag*/)
 
   compute_grad->forward_fields(); // also forwards v and rho for chi
   compute_kernel->compute_peratom();
-  compute_interface->compute_peratom();
+  if (interface_flag) compute_interface->compute_peratom();
 
   compute_grad->compute_peratom();
   compute_grad->forward_gradients();
@@ -369,6 +388,9 @@ void FixRHEO::pre_force(int /*vflag*/)
         status[i] &= ~STATUS_SHIFT;
     }
   }
+
+  // Calculate surfaces, update status
+  if (surface_flag) compute_surface->compute_peratom();
 }
 
 /* ---------------------------------------------------------------------- */
diff --git a/src/RHEO/fix_rheo.h b/src/RHEO/fix_rheo.h
index 89304fe22f..d2097ade71 100644
--- a/src/RHEO/fix_rheo.h
+++ b/src/RHEO/fix_rheo.h
@@ -40,9 +40,11 @@ class FixRHEO : public Fix {
 
   // Model parameters
   double h, rho0, csq;
-  int zmin_kernel, rhosum_zmin;
-  int kernel_style;
+  int zmin_kernel, zmin_rhosum, zmin_surface;
+  int kernel_style, surface_style;
+  double divr_surface;
   enum {QUINTIC, CRK0, CRK1, CRK2};
+  enum {COORDINATION, DIVR};
 
   // Status variables
   enum {
@@ -75,8 +77,6 @@ class FixRHEO : public Fix {
   int viscosity_fix_defined;
   int pressure_fix_defined;
   int thermal_fix_defined;
-  int interface_fix_defined;
-  int surface_fix_defined;
 
   class ComputeRHEOGrad *compute_grad;
   class ComputeRHEOKernel *compute_kernel;