diff --git a/src/.gitignore b/src/.gitignore
index ef74318fc7..c7e022797f 100644
--- a/src/.gitignore
+++ b/src/.gitignore
@@ -219,6 +219,8 @@
 /fix_rheo_pressure.h
 /fix_rheo_stress.cpp
 /fix_rheo_stress.h
+/fix_rheo_tension.cpp
+/fix_rheo_tension.h
 /fix_rheo_thermal.cpp
 /fix_rheo_thermal.h
 /fix_rheo_viscosity.cpp
@@ -227,8 +229,6 @@
 /pair_rheo.h
 /pair_rheo_react.cpp
 /pair_rheo_react.h
-/pair_rheo_tension.cpp
-/pair_rheo_tension.h
 
 /compute_grid.cpp
 /compute_grid.h
diff --git a/src/RHEO/compute_rheo_grad.cpp b/src/RHEO/compute_rheo_grad.cpp
index 7fdff41403..46c1500556 100644
--- a/src/RHEO/compute_rheo_grad.cpp
+++ b/src/RHEO/compute_rheo_grad.cpp
@@ -355,8 +355,6 @@ int ComputeRHEOGrad::pack_forward_comm(int n, int *list, double *buf,
         }
       }
 
-
-
       if (rho_flag)
         buf[m++] = rho[j];
 
diff --git a/src/RHEO/compute_rheo_surface.cpp b/src/RHEO/compute_rheo_surface.cpp
index 230bca8219..88b33c09af 100644
--- a/src/RHEO/compute_rheo_surface.cpp
+++ b/src/RHEO/compute_rheo_surface.cpp
@@ -87,27 +87,26 @@ void ComputeRHEOSurface::init()
 
   cutsq = cut * cut;
 
-  // Create rsurface, divr, nsurface arrays if they don't already exist
-  // Create a custom atom property so it works with compute property/atom
-  // Do not create grow callback as there's no reason to copy/exchange data
-  // Manually grow if nmax_store exceeded
+  // Create rsurface, divr, nsurface arrays as custom atom properties,
+  //   can print with compute property/atom
+  //   no grow callback as there's no reason to copy/exchange data, manually grow
   // For B and gradC, create a local array since they are unlikely to be printed
 
+  int dim = domain->dimension;
   int tmp1, tmp2;
-  int index = atom->find_custom("rheo_divr", tmp1, tmp2);
-  if (index == -1)  index = atom->add_custom("rheo_divr", 1, 0);
-  divr = atom->dvector[index];
+  index_divr = atom->find_custom("rheo_divr", tmp1, tmp2);
+  if (index_divr == -1) index_divr = atom->add_custom("rheo_divr", 1, 0);
+  divr = atom->dvector[index_divr];
 
-  index = atom->find_custom("rheo_rsurface", tmp1, tmp2);
-  if (index == -1)  index = atom->add_custom("rheo_rsurface", 1, 0);
-  rsurface = atom->dvector[index];
+  index_rsurf = atom->find_custom("rheo_rsurface", tmp1, tmp2);
+  if (index_rsurf == -1) index_rsurf = atom->add_custom("rheo_rsurface", 1, 0);
+  rsurface = atom->dvector[index_rsurf];
 
-  index = atom->find_custom("rheo_nsurface", tmp1, tmp2);
-  if (index == -1)  index = atom->add_custom("rheo_nsurface", 1, 3);
-  nsurface = atom->darray[index];
+  index_nsurf = atom->find_custom("rheo_nsurface", tmp1, tmp2);
+  if (index_nsurf == -1) index_nsurf = atom->add_custom("rheo_nsurface", 1, dim);
+  nsurface = atom->darray[index_nsurf];
 
   nmax_store = atom->nmax;
-  int dim = domain->dimension;
   memory->create(B, nmax_store, dim * dim, "rheo/surface:B");
   memory->create(gradC, nmax_store, dim * dim, "rheo/surface:gradC");
 
@@ -148,32 +147,21 @@ void ComputeRHEOSurface::compute_peratom()
   numneigh = list->numneigh;
   firstneigh = list->firstneigh;
 
-  int nmax = atom->nmax;
-  if (nmax_store <= nmax) {
-    memory->grow(divr, nmax, "atom:rheo_divr");
-    memory->grow(rsurface, nmax, "atom:rheo_rsurface");
-    memory->grow(nsurface, nmax, 3, "atom:rheo_nsurface");
+  // Grow and zero arrays
+  if (nmax_store <= atom->nmax)
+    grow_arrays(atom->nmax);
 
-    memory->grow(B, nmax, dim * dim, "rheo/surface:B");
-    memory->grow(gradC, nmax, dim * dim, "rheo/surface:gradC");
-
-    nmax_store = atom->nmax;
-  }
+  size_t nbytes = nmax_store * sizeof(double);
+  memset(&divr, 0, nbytes);
+  memset(&rsurface, 0, nbytes);
+  memset(&nsurface, 0, 3 * nbytes);
+  memset(&gradC, 0, 3 * 3 * nbytes);
+  memset(&B, 0, 3 * 3 * nbytes);
 
+  // Remove surface settings
   int nall = nlocal + atom->nghost;
-  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;
-      }
-      nsurface[i][a] = 0.0;
-    }
-    divr[i] = 0.0;
-
-    // Remove surface settings
+  for (i = 0; i < nall; i++)
     status[i] &= SURFACEMASK;
-  }
 
   // loop over neighbors to calculate the average orientation of neighbors
   for (ii = 0; ii < inum; ii++) {
@@ -424,3 +412,25 @@ void ComputeRHEOSurface::unpack_forward_comm(int n, int first, double *buf)
     }
   }
 }
+
+/* ---------------------------------------------------------------------- */
+
+void ComputeRHEOSurface::grow_arrays(int nmax)
+{
+  int dim = domain->dimension;
+
+  // Grow atom variables and reassign pointers
+  memory->grow(atom->dvector[index_divr], nmax, "atom:rheo_divr");
+  memory->grow(atom->dvector[index_rsurf], nmax, "atom:rheo_rsurface");
+  memory->grow(atom->darray[index_nsurf], nmax, dim, "atom:rheo_nsurface");
+
+  divr = atom->dvector[index_divr];
+  rsurface = atom->dvector[index_rsurf];
+  nsurface = atom->darray[index_nsurf];
+
+  // Grow local variables
+  memory->grow(B, nmax, dim * dim, "rheo/surface:B");
+  memory->grow(gradC, nmax, dim * dim, "rheo/surface:gradC");
+
+  nmax_store = atom->nmax;
+}
diff --git a/src/RHEO/compute_rheo_surface.h b/src/RHEO/compute_rheo_surface.h
index 1d1cdba3fa..f1b1af7742 100644
--- a/src/RHEO/compute_rheo_surface.h
+++ b/src/RHEO/compute_rheo_surface.h
@@ -40,14 +40,18 @@ class ComputeRHEOSurface : public Compute {
   class FixRHEO *fix_rheo;
 
  private:
-  double cut, cutsq, rho0, threshold_divr;
   int surface_style, nmax_store, threshold_z, threshold_splash, interface_flag;
-  double **B, **gradC;
   int threshold_style, comm_stage;
+  int index_divr, index_rsurf, index_nsurf;
+
+  double cut, cutsq, rho0, threshold_divr;
+  double **B, **gradC;
 
   class NeighList *list;
   class ComputeRHEOKernel *compute_kernel;
   class ComputeRHEOInterface *compute_interface;
+
+  void grow_arrays(int);
 };
 
 }    // namespace LAMMPS_NS
diff --git a/src/RHEO/fix_rheo_tension.cpp b/src/RHEO/fix_rheo_tension.cpp
new file mode 100644
index 0000000000..44bc3bf580
--- /dev/null
+++ b/src/RHEO/fix_rheo_tension.cpp
@@ -0,0 +1,452 @@
+/* ----------------------------------------------------------------------
+ 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 "fix_rheo_tension.h"
+
+#include "atom.h"
+#include "comm.h"
+#include "compute_rheo_kernel.h"
+#include "compute_rheo_interface.h"
+#include "domain.h"
+#include "error.h"
+#include "fix_rheo.h"
+#include "force.h"
+#include "math_extra.h"
+#include "memory.h"
+#include "modify.h"
+#include "neighbor.h"
+#include "neigh_list.h"
+#include "update.h"
+#include "utils.h"
+
+#include <cmath>
+
+using namespace LAMMPS_NS;
+using namespace RHEO_NS;
+using namespace MathExtra;
+using namespace FixConst;
+
+/* ---------------------------------------------------------------------- */
+
+FixRHEOTension::FixRHEOTension(LAMMPS *lmp, int narg, char **arg) :
+  Fix(lmp, narg, arg), compute_kernel(nullptr), compute_interface(nullptr), fix_rheo(nullptr)
+{
+  if (narg != 4) error->all(FLERR,"Illegal fix command");
+  alpha = utils::numeric(FLERR,arg[3],false,lmp);
+
+  comm_forward = 3;
+  comm_reverse = 3;
+
+  // Create cgrad, n, and divr arrays as custom atom properties,
+  //   can print with compute property/atom
+  //   no grow callback as there's no reason to copy/exchange data, manually grow
+  // For norm, create a local array since they are unlikely to be printed
+
+  int tmp1, tmp2;
+  index_cgradt = atom->find_custom("cgrad_rheo_tension", tmp1, tmp2);
+  if (index_cgradt == -1)  index_cgradt = atom->add_custom("cgrad_rheo_tension", 1, 3);
+  cgradt = atom->darray[index_cgradt];
+
+  index_nt = atom->find_custom("n_rheo_tension", tmp1, tmp2);
+  if (index_nt == -1)  index_nt = atom->add_custom("n_rheo_tension", 1, 3);
+  nt = atom->darray[index_nt];
+
+  index_divnt = atom->find_custom("divn_rheo_tension", tmp1, tmp2);
+  if (index_divnt == -1) index_divnt = atom->add_custom("divn_rheo_tension", 1, 0);
+  divnt = atom->dvector[index_divnt];
+
+  norm = nullptr;
+  nmax_store = 0;
+}
+
+/* ---------------------------------------------------------------------- */
+
+FixRHEOTension::~FixRHEOTension()
+{
+  // Remove custom property if it exists
+  int tmp1, tmp2, index;
+
+  index = atom->find_custom("cgrad_rheo_tension", tmp1, tmp2);
+  if (index != -1) atom->remove_custom(index, 1, 3);
+
+  index = atom->find_custom("n_rheo_tension", tmp1, tmp2);
+  if (index != -1) atom->remove_custom(index, 1, 3);
+
+  index = atom->find_custom("divn_rheo_tension", tmp1, tmp2);
+  if (index != -1) atom->remove_custom(index, 1, 0);
+
+  memory->destroy(norm);
+}
+
+/* ---------------------------------------------------------------------- */
+
+int FixRHEOTension::setmask()
+{
+  int mask = 0;
+  mask |= POST_FORCE;
+  return mask;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixRHEOTension::init()
+{
+  auto fixes = modify->get_fix_by_style("^rheo$");
+  if (fixes.size() == 0) error->all(FLERR, "Need to define fix rheo to use fix rheo/tension");
+  fix_rheo = dynamic_cast<FixRHEO *>(fixes[0]);
+
+  compute_kernel = fix_rheo->compute_kernel;
+  compute_interface = fix_rheo->compute_interface;
+  interface_flag = fix_rheo->interface_flag;
+  h = fix_rheo->h;
+  rho0 = fix_rheo->rho0;
+
+  hsq = h * h;
+
+  neighbor->add_request(this, NeighConst::REQ_DEFAULT);
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixRHEOTension::init_list(int /*id*/, NeighList *ptr)
+{
+  list = ptr;
+}
+
+
+/* ---------------------------------------------------------------------- */
+
+void FixRHEOTension::setup(int vflag)
+{
+  // Grow and populate arrays
+  post_force(vflag);
+}
+
+/* ----------------------------------------------------------------------
+  Calculate and apply tension forces
+------------------------------------------------------------------------- */
+
+void FixRHEOTension::post_force(int vflag)
+{
+  int i, j, a, ii, jj, inum, jnum, itype, jtype;
+  int fluidi, fluidj;
+  double xtmp, ytmp, ztmp, w, wp, c;
+  double rhoi, rhoj, Voli, Volj;
+  double *dWij, *dWji;
+  double dx[3], ft[3];
+
+  int *ilist, *jlist, *numneigh, **firstneigh;
+  double imass, jmass, rsq, r, rinv;
+
+  int nlocal = atom->nlocal;
+  int newton = force->newton;
+  int dim = domain->dimension;
+
+  v_init(vflag);
+
+  double **x = atom->x;
+  double **f = atom->f;
+  double *rho = atom->rho;
+  double *mass = atom->mass;
+  imageint *image = atom->image;
+  int *type = atom->type;
+  int *status = atom->status;
+
+  inum = list->inum;
+  ilist = list->ilist;
+  numneigh = list->numneigh;
+  firstneigh = list->firstneigh;
+
+  if (nmax_store <= atom->nmax)
+    grow_arrays(atom->nmax);
+
+  for (i = 0; i < nlocal+atom->nghost; i++) {
+    cgradt[i][0] = 0.0;
+    cgradt[i][1] = 0.0;
+    cgradt[i][2] = 0.0;
+    norm[i] = 0.0;
+    divnt[i] = 0.0;
+  }
+
+  // Calculate color gradient
+  for (ii = 0; ii < inum; ii++) {
+    i = ilist[ii];
+    xtmp = x[i][0];
+    ytmp = x[i][1];
+    ztmp = x[i][2];
+    itype = type[i];
+    fluidi = !(status[i] & PHASECHECK);
+    jlist = firstneigh[i];
+    jnum = numneigh[i];
+    imass = mass[itype];
+
+    for (jj = 0; jj < jnum; jj++) {
+      j = jlist[jj];
+      j &= NEIGHMASK;
+
+      dx[0] = xtmp - x[j][0];
+      dx[1] = ytmp - x[j][1];
+      dx[2] = ztmp - x[j][2];
+      rsq = lensq3(dx);
+
+      if (rsq > hsq) continue;
+
+      fluidj = !(status[j] & PHASECHECK);
+      jtype = type[j];
+      r = sqrt(rsq);
+      rinv = 1 / r;
+
+      rhoi = rho[i];
+      rhoj = rho[j];
+
+      // Add corrections for walls
+      if (interface_flag) {
+        if (fluidi && (!fluidj)) {
+          rhoj = compute_interface->correct_rho(j, i);
+        } else if ((!fluidi) && fluidj) {
+          rhoi = compute_interface->correct_rho(i, j);
+        } else if ((!fluidi) && (!fluidj)) {
+          rhoi = rho0;
+          rhoj = rho0;
+        }
+      }
+
+      Voli = mass[itype] / rhoi;
+      Volj = mass[jtype] / rhoj;
+
+      wp = compute_kernel->calc_dw(i, j, dx[0], dx[1], dx[2],r);
+      dWij = compute_kernel->dWij;
+      dWji = compute_kernel->dWji;
+
+      c = 0;
+      if (itype == jtype) c += rhoi;
+      c /= (rhoi + rhoj);
+
+      for (a = 0; a < 3; a++) {
+        cgradt[i][a] -= c * Volj * dWij[a];
+        if (newton || j < nlocal)
+          cgradt[j][a] -= c * Voli * dWji[a];
+      }
+    }
+  }
+
+  comm_stage = 0;
+  comm_reverse = 3;
+  if (newton) comm->reverse_comm(this);
+
+  // Calculate normal direction
+  double minv;
+  for (i = 0; i < nlocal; i++) {
+    minv = 1.0 / sqrt(cgradt[i][0] * cgradt[i][0] + cgradt[i][1] * cgradt[i][1] + cgradt[i][2] * cgradt[i][2]);
+
+    for (a = 0; a < 3; a++)
+      nt[i][a] = cgradt[i][a] * minv;
+  }
+
+  comm->forward_comm(this);
+
+  // Calculate divergence
+  for (ii = 0; ii < inum; ii++) {
+    i = ilist[ii];
+    xtmp = x[i][0];
+    ytmp = x[i][1];
+    ztmp = x[i][2];
+    itype = type[i];
+    fluidi = !(status[i] & PHASECHECK);
+    jlist = firstneigh[i];
+    jnum = numneigh[i];
+    imass = mass[itype];
+
+    for (jj = 0; jj < jnum; jj++) {
+      j = jlist[jj];
+      j &= NEIGHMASK;
+
+      dx[0] = xtmp - x[j][0];
+      dx[1] = ytmp - x[j][1];
+      dx[2] = ztmp - x[j][2];
+      rsq = lensq3(dx);
+
+      if (rsq > hsq) continue;
+
+      fluidj = !(status[j] & PHASECHECK);
+      jtype = type[j];
+      r = sqrt(rsq);
+      rinv = 1 / r;
+
+      rhoi = rho[i];
+      rhoj = rho[j];
+
+      // Add corrections for walls
+      if (interface_flag) {
+        if (fluidi && (!fluidj)) {
+          rhoj = compute_interface->correct_rho(j, i);
+        } else if ((!fluidi) && fluidj) {
+          rhoi = compute_interface->correct_rho(i, j);
+        } else if ((!fluidi) && (!fluidj)) {
+          rhoi = rho0;
+          rhoj = rho0;
+        }
+      }
+
+      Voli = mass[itype] / rhoi;
+      Volj = mass[jtype] / rhoj;
+
+      wp = compute_kernel->calc_dw(i, j, dx[0], dx[1], dx[2],r);
+      dWij = compute_kernel->dWij;
+      dWji = compute_kernel->dWji;
+
+      for (a = 0; a < 3; a++) {
+        divnt[i] -= nt[i][a] * Volj * dWij[a];
+        norm[i] -= dx[a] * Volj * dWij[a];
+        if (newton || j < nlocal) {
+          divnt[j] -= nt[j][a] * Voli * dWji[a];
+          norm[j] += dx[a] * Voli * dWji[a];
+        }
+      }
+    }
+  }
+
+  comm_stage = 1;
+  comm_reverse = 2;
+  if (newton) comm->reverse_comm(this);
+
+  // Skip forces if it's setup
+  if (update->setupflag) return;
+
+  // apply force
+  int prefactor;
+  double unwrap[3];
+  double v[6];
+  for (i = 0; i < nlocal; i++) {
+    itype = type[i];
+    divnt[i] /= norm[i];
+
+    prefactor *= -alpha * divnt[i] / mass[itype];
+
+    for (a = 0; a < 3; a++)
+      f[i][a] += prefactor * cgradt[i][a];
+
+    if (evflag) {
+      domain->unmap(x[i], image[i], unwrap);
+      v[0] = prefactor * cgradt[i][0] * unwrap[0];
+      v[1] = prefactor * cgradt[i][1] * unwrap[1];
+      v[2] = prefactor * cgradt[i][2] * unwrap[2];
+      v[3] = prefactor * cgradt[i][0] * unwrap[1];
+      v[4] = prefactor * cgradt[i][0] * unwrap[2];
+      v[5] = prefactor * cgradt[i][1] * unwrap[2];
+      v_tally(i, v);
+    }
+  }
+
+
+
+  if (evflag) {
+
+        }
+
+}
+
+
+/* ---------------------------------------------------------------------- */
+
+int FixRHEOTension::pack_forward_comm(int n, int *list, double *buf, int /*pbc_flag*/, int * /*pbc*/)
+{
+  int i, j, a, m;
+  m = 0;
+
+  for (i = 0; i < n; i++) {
+    j = list[i];
+    for (a = 0; a < 3; a++)
+      buf[m++] = nt[j][a];
+  }
+  return m;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixRHEOTension::unpack_forward_comm(int n, int first, double *buf)
+{
+  int i, a, m, last;
+
+  m = 0;
+  last = first + n;
+  for (i = first; i < last; i++)
+    for (a = 0; a < 3; a++)
+      nt[i][a] = buf[m++];
+}
+
+
+/* ---------------------------------------------------------------------- */
+
+int FixRHEOTension::pack_reverse_comm(int n, int first, double *buf)
+{
+  int i, a, m, last;
+
+  m = 0;
+  last = first + n;
+  if (comm_stage == 0)
+    for (i = first; i < last; i++)
+      for (a = 0; a < 3; a++)
+        buf[m++] = cgradt[i][a];
+  else
+    for (i = first; i < last; i++) {
+      buf[m++] = norm[i];
+      buf[m++] = divnt[i];
+    }
+  return m;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixRHEOTension::unpack_reverse_comm(int n, int *list, double *buf)
+{
+  int i, j, a, m;
+
+  m = 0;
+  if (comm_stage == 0)
+    for (i = 0; i < n; i++) {
+      j = list[i];
+      for (a = 0; a < 3; a++)
+        cgradt[j][a] += buf[m++];
+    }
+  else
+    for (i = 0; i < n; i++) {
+      j = list[i];
+      norm[j] += buf[m++];
+      divnt[j] += buf[m++];
+    }
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixRHEOTension::grow_arrays(int nmax)
+{
+  // Grow atom variables and reassign pointers
+  memory->grow(atom->darray[index_cgradt], nmax, 3, "atom:rheo_cgradt");
+  memory->grow(atom->darray[index_nt], nmax, 3, "atom:rheo_nt");
+  memory->grow(atom->dvector[index_divnt], nmax, "atom:rheo_divnt");
+
+  cgradt = atom->darray[index_cgradt];
+  nt = atom->darray[index_nt];
+  divnt = atom->dvector[index_divnt];
+
+  // Grow local variables
+  memory->grow(norm, nmax, "rheo/tension:norm");
+
+  nmax_store = atom->nmax;
+}
\ No newline at end of file
diff --git a/src/RHEO/pair_rheo_tension.h b/src/RHEO/fix_rheo_tension.h
similarity index 60%
rename from src/RHEO/pair_rheo_tension.h
rename to src/RHEO/fix_rheo_tension.h
index 2a046ff324..74b8b71436 100644
--- a/src/RHEO/pair_rheo_tension.h
+++ b/src/RHEO/fix_rheo_tension.h
@@ -11,44 +11,45 @@
    See the README file in the top-level LAMMPS directory.
 ------------------------------------------------------------------------- */
 
-#ifdef PAIR_CLASS
+#ifdef FIX_CLASS
 // clang-format off
-PairStyle(rheo/tension,PairRHEOTension)
+FixStyle(rheo/tension,FixRHEOTension)
 // clang-format on
 #else
 
-#ifndef LMP_PAIR_RHEO_TENSION_H
-#define LMP_PAIR_RHEO_TENSION_H
+#ifndef LMP_FIX_RHEO_TENSION_H
+#define LMP_FIX_RHEO_TENSION_H
 
-#include "pair.h"
+#include "fix.h"
 
 namespace LAMMPS_NS {
 
-class PairRHEOTension : public Pair {
+class FixRHEOTension : public Fix {
  public:
-  PairRHEOTension(class LAMMPS *);
-  ~PairRHEOTension() override;
-  void compute(int, int) override;
-  void settings(int, char **) override;
-  void coeff(int, char **) override;
-  void setup() override;
-  void init_style() override;
-  double init_one(int, int) override;
+  FixRHEOTension(class LAMMPS *, int, char **);
+  ~FixRHEOTension() override;
+  int setmask() override;
+  void init() override;
+  void init_list(int, class NeighList *) override;
+  void setup(int) override;
+  void post_force(int) override;
   int pack_forward_comm(int, int *, double *, int, int *) override;
   void unpack_forward_comm(int, int, double *) override;
   int pack_reverse_comm(int, int, double *) override;
   void unpack_reverse_comm(int, int *, double *) override;
+  void grow_arrays(int) override;
 
- protected:
-  int nmax_store;
-  double **nt, *ct;
-  double **alpha;
-  double h, hsq, hinv, hinv3;
+ private:
+  int nmax_store, comm_stage, interface_flag;
+  int index_nt, index_cgradt, index_divnt;
 
-  void allocate();
+  double **nt, **cgradt, *divnt, *norm;
+  double alpha, h, hsq, hinv, hinv3, rho0;
 
   class ComputeRHEOKernel *compute_kernel;
+  class ComputeRHEOInterface *compute_interface;
   class FixRHEO *fix_rheo;
+  class NeighList *list;
 };
 
 }    // namespace LAMMPS_NS
diff --git a/src/RHEO/pair_rheo_tension.cpp b/src/RHEO/pair_rheo_tension.cpp
deleted file mode 100644
index ef0d0b60b4..0000000000
--- a/src/RHEO/pair_rheo_tension.cpp
+++ /dev/null
@@ -1,383 +0,0 @@
-/* ----------------------------------------------------------------------
- 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)
------------------------------------------------------------------------ */
-
-#include "pair_rheo_tension.h"
-
-#include "atom.h"
-#include "comm.h"
-#include "compute_rheo_kernel.h"
-#include "domain.h"
-#include "error.h"
-#include "fix_rheo.h"
-#include "force.h"
-#include "math_extra.h"
-#include "memory.h"
-#include "modify.h"
-#include "neighbor.h"
-#include "neigh_list.h"
-#include "update.h"
-#include "utils.h"
-
-#include <cmath>
-
-using namespace LAMMPS_NS;
-using namespace RHEO_NS;
-using namespace MathExtra;
-
-static constexpr double EPSILON = 1e-2;
-
-/* ---------------------------------------------------------------------- */
-
-PairRHEOTension::PairRHEOTension(LAMMPS *lmp) :
-  Pair(lmp), compute_kernel(nullptr), fix_rheo(nullptr)
-{
-  restartinfo = 0;
-  single_enable = 0;
-
-  comm_forward = 3;
-  comm_reverse = 3;
-}
-
-/* ---------------------------------------------------------------------- */
-
-PairRHEOTension::~PairRHEOTension()
-{
-  // Remove custom property if it exists
-  int tmp1, tmp2, index;
-
-  index = atom->find_custom("rheo_c_tension", tmp1, tmp2);
-  if (index != -1) atom->remove_custom(index, 1, 0);
-
-  index = atom->find_custom("rheo_n_tension", tmp1, tmp2);
-  if (index != -1) atom->remove_custom(index, 1, 3);
-
-  if (allocated) {
-    memory->destroy(alpha);
-    memory->destroy(setflag);
-    memory->destroy(cutsq);
-  }
-}
-
-/* ---------------------------------------------------------------------- */
-
-void PairRHEOTension::compute(int eflag, int vflag)
-{
-  int i, j, a, b, ii, jj, inum, jnum, itype, jtype;
-  int fluidi, fluidj;
-  double xtmp, ytmp, ztmp, w, wp;
-  double rhoi, rhoj, voli, volj;
-  double *dWij, *dWji;
-  double dx[3], ft[3];
-
-  int *ilist, *jlist, *numneigh, **firstneigh;
-  double imass, jmass, rsq, r, rinv;
-
-  int nlocal = atom->nlocal;
-  int newton_pair = force->newton_pair;
-  int dim = domain->dimension;
-
-  ev_init(eflag, vflag);
-
-  double **x = atom->x;
-  double **f = atom->f;
-  double *rho = atom->rho;
-  double *mass = atom->mass;
-  double *special_lj = force->special_lj;
-  int *type = atom->type;
-  int *status = atom->status;
-  tagint *tag = atom->tag;
-
-  inum = list->inum;
-  ilist = list->ilist;
-  numneigh = list->numneigh;
-  firstneigh = list->firstneigh;
-/*
-  int nmax = atom->nmax;
-  if (nmax_store <= nmax) {
-    memory->grow(ct, nmax, "atom:rheo_c_tension");
-    memory->grow(nnt_tension, nmax, 3, "atom:rheo_n_tension");
-    nmax_store = atom->nmax;
-  }
-
-  // loop over neighbors of my atoms
-
-  for (ii = 0; ii < inum; ii++) {
-    i = ilist[ii];
-    xtmp = x[i][0];
-    ytmp = x[i][1];
-    ztmp = x[i][2];
-    itype = type[i];
-    jlist = firstneigh[i];
-    jnum = numneigh[i];
-    imass = mass[itype];
-    rhoi = rho[i];
-    voli = imass / rhoi;
-    fluidi = !(status[i] & PHASECHECK);
-
-    for (jj = 0; jj < jnum; jj++) {
-      j = jlist[jj];
-      j &= NEIGHMASK;
-
-      dx[0] = xtmp - x[j][0];
-      dx[1] = ytmp - x[j][1];
-      dx[2] = ztmp - x[j][2];
-      rsq = lensq3(dx);
-      jtype = type[j];
-
-      if (rsq > hsq) continue;
-
-      r = sqrt(rsq);
-      rinv = 1 / r;
-
-      jmass = mass[jtype];
-      rhoj = rho[j];
-      volj = jmass / rhoj;
-      fluidj = !(status[j] & PHASECHECK);
-
-      wp = compute_kernel->calc_dw(i, j, dx[0], dx[1], dx[2],r);
-      dWij = compute_kernel->dWij;
-      dWji = compute_kernel->dWji;
-
-      f[i][0] += ft[0];
-      f[i][1] += ft[1];
-      f[i][2] += ft[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]);
-
-      if (newton_pair || j < nlocal) {
-
-        f[j][0] -= ft[0];
-        f[j][1] -= ft[1];
-        f[j][2] -= ft[2];
-      }
-    }
-  }
-
-  if (vflag_fdotr) virial_fdotr_compute();
-
-  comm->reverse_comm(this);
-  comm->forward_comm(this);
-*/
-}
-
-/* ----------------------------------------------------------------------
- allocate all arrays
- ------------------------------------------------------------------------- */
-
-void PairRHEOTension::allocate()
-{
-  allocated = 1;
-  int n = atom->ntypes;
-
-  memory->create(setflag, n + 1, n + 1, "pair:setflag");
-  for (int i = 1; i <= n; i++)
-    for (int j = i; j <= n; j++)
-      setflag[i][j] = 0;
-
-  memory->create(alpha, n + 1, n + 1, "pair:alpha");
-  memory->create(cutsq, n + 1, n + 1, "pair:cutsq");
-}
-
-/* ----------------------------------------------------------------------
- global settings
- ------------------------------------------------------------------------- */
-
-void PairRHEOTension::settings(int narg, char **arg)
-{
-}
-
-/* ----------------------------------------------------------------------
- set coeffs for one or more type pairs
- ------------------------------------------------------------------------- */
-
-void PairRHEOTension::coeff(int narg, char **arg)
-{
-  if (narg != 3)
-    error->all(FLERR,"Incorrect number of args for pair_style rheo coefficients");
-  if (!allocated)
-    allocate();
-
-  int ilo, ihi, jlo, jhi;
-  utils::bounds(FLERR, arg[0], 1, atom->ntypes, ilo, ihi,error);
-  utils::bounds(FLERR, arg[1], 1, atom->ntypes, jlo, jhi,error);
-
-  double alpha_one = utils::numeric(FLERR, arg[2], false, lmp);
-
-  int count = 0;
-  for (int i = ilo; i <= ihi; i++) {
-    for (int j = 0; j <= atom->ntypes; j++) {
-      alpha[i][j] = alpha_one;
-      setflag[i][j] = 1;
-      count++;
-    }
-  }
-
-  if (count == 0)
-    error->all(FLERR,"Incorrect args for pair rheo/tension coefficients");
-}
-
-/* ----------------------------------------------------------------------
- setup specific to this pair style
- ------------------------------------------------------------------------- */
-
-void PairRHEOTension::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");
-  fix_rheo = dynamic_cast<FixRHEO *>(fixes[0]);
-  /*
-  compute_kernel = fix_rheo->compute_kernel;
-  compute_grad = fix_rheo->compute_grad;
-  compute_interface = fix_rheo->compute_interface;
-  h = fix_rheo->h;
-  csq = fix_rheo->csq;
-  rho0 = fix_rheo->rho0;
-
-  hsq = h * h;
-  hinv = 1.0 / h;
-  hinv3 = hinv * 3.0;
-  */
-}
-
-/* ----------------------------------------------------------------------
-   init specific to this pair style
-------------------------------------------------------------------------- */
-
-void PairRHEOTension::init_style()
-{
-  neighbor->add_request(this);
-
-
-  // Create c_tension arrays n_tension arrays if they don't already exist
-  // Create a custom atom property so it works with compute property/atom
-  // Do not create grow callback as there's no reason to copy/exchange data
-  // Manually grow if nmax_store exceeded
-  // For B and gradC, create a local array since they are unlikely to be printed
-
-  int tmp1, tmp2;
-  int index = atom->find_custom("rheo_c_tension", tmp1, tmp2);
-  if (index == -1)  index = atom->add_custom("rheo_c_tension", 1, 0);
-  ct = atom->dvector[index];
-
-  index = atom->find_custom("rheo_n_tension", tmp1, tmp2);
-  if (index == -1)  index = atom->add_custom("rheo_n_tension", 1, 3);
-  nt = atom->darray[index];
-
-  nmax_store = atom->nmax;
-}
-
-/* ----------------------------------------------------------------------
- init for one type pair i,j and corresponding j,i
- ------------------------------------------------------------------------- */
-
-double PairRHEOTension::init_one(int i, int j)
-{
-  if (setflag[i][j] == 0)
-    error->all(FLERR,"All pair rheo/tension coeffs are not set");
-
-  alpha[j][i] = alpha[i][j];
-
-  return h;
-}
-
-
-/* ---------------------------------------------------------------------- */
-
-int PairRHEOTension::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;
-
-  for (i = 0; i < n; i++) {
-    j = list[i];
-    if (comm_stage == 0) {
-      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];
-    }
-  }
-  return m;
-  */
-}
-
-/* ---------------------------------------------------------------------- */
-
-void PairRHEOTension::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++) {
-    if (comm_stage == 0) {
-      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++];
-    }
-  }
-  */
-}
-
-
-/* ---------------------------------------------------------------------- */
-
-int PairRHEOTension::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 PairRHEOTension::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++];
-  }
-  */
-}