diff --git a/src/RHEO/compute_rheo_grad.cpp b/src/RHEO/compute_rheo_grad.cpp
index 03e76b194a..31f60550c1 100644
--- a/src/RHEO/compute_rheo_grad.cpp
+++ b/src/RHEO/compute_rheo_grad.cpp
@@ -38,7 +38,7 @@ enum{COMMGRAD, COMMFIELD};
 /* ---------------------------------------------------------------------- */
 
 ComputeRHEOGrad::ComputeRHEOGrad(LAMMPS *lmp, int narg, char **arg) :
-  Compute(lmp, narg, arg), compute_interface(nullptr), compute_kernel(nullptr),
+  Compute(lmp, narg, arg), fix_rheo(nullptr), compute_interface(nullptr), compute_kernel(nullptr),
   gradv(nullptr), gradr(nullptr), gradt(nullptr), gradn(nullptr)
 {
   if (narg < 4) error->all(FLERR,"Illegal compute rheo/grad command");
@@ -61,40 +61,26 @@ ComputeRHEOGrad::ComputeRHEOGrad(LAMMPS *lmp, int narg, char **arg) :
     ncomm_grad += dim * dim;
     ncomm_field += dim;
     comm_reverse += dim * dim;
-    indexv = atom->add_custom("rheo_grad_v", 1, dim * dim);
-    gradv = atom->darray[indexv];
   }
 
   if (rho_flag) {
     ncomm_grad += dim;
     ncomm_field += 1;
     comm_reverse += dim;
-    indexr = atom->add_custom("rheo_grad_rho", 1, dim);
-    gradr = atom->darray[indexr];
   }
 
   if (temperature_flag) {
     ncomm_grad += dim;
     ncomm_field += 1;
     comm_reverse += dim;
-    indext= atom->add_custom("rheo_grad_temp", 1, dim);
-    gradt = atom->darray[indext];
   }
 
   if (eta_flag) {
     ncomm_grad += dim;
     comm_reverse += dim;
-    indexn = atom->add_custom("rheo_grad_eta", 1, dim);
-    gradn = atom->darray[indexn];
   }
 
-  // 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_old exceeded
-
   comm_forward = ncomm_grad;
-  nmax_old = 0;
-  grow_arrays(atom->nmax);
 }
 
 /* ---------------------------------------------------------------------- */
@@ -102,14 +88,16 @@ ComputeRHEOGrad::ComputeRHEOGrad(LAMMPS *lmp, int narg, char **arg) :
 ComputeRHEOGrad::~ComputeRHEOGrad()
 {
   int dim = domain->dimension;
-  if (velocity_flag)
-    atom->remove_custom(indexv, 1, dim * dim);
-  if (rho_flag)
-    atom->remove_custom(indexr, 1, dim);
-  if (temperature_flag)
-    atom->remove_custom(indext, 1, dim);
-  if (eta_flag)
-    atom->remove_custom(indexn, 1, dim);
+  int tmp1, tmp2, index;
+
+  index = atom->find_custom("rheo_grad_v", tmp1, tmp2);
+  if (index != 1) atom->remove_custom(index, 1, dim * dim);
+  index = atom->find_custom("rheo_grad_rho", tmp1, tmp2);
+  if (index != 1)  atom->remove_custom(index, 1, dim);
+  index = atom->find_custom("rheo_grad_t", tmp1, tmp2);
+  if (index != 1)  atom->remove_custom(index, 1, dim);
+  index = atom->find_custom("rheo_grad_eta", tmp1, tmp2);
+  if (index != 1)  atom->remove_custom(index, 1, dim);
 }
 
 /* ---------------------------------------------------------------------- */
@@ -123,6 +111,36 @@ void ComputeRHEOGrad::init()
   rho0 = fix_rheo->rho0;
   compute_kernel = fix_rheo->compute_kernel;
   compute_interface = fix_rheo->compute_interface;
+
+  // Create coordination array if it doesn'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_old exceeded
+
+  int index;
+  int dim = domain->dimension;
+  if (velocity_flag) {
+    index = atom->add_custom("rheo_grad_v", 1, dim * dim);
+    gradv = atom->darray[index];
+  }
+
+  if (rho_flag) {
+    index = atom->add_custom("rheo_grad_rho", 1, dim);
+    gradr = atom->darray[index];
+  }
+
+  if (temperature_flag) {
+    index= atom->add_custom("rheo_grad_temp", 1, dim);
+    gradt = atom->darray[index];
+  }
+
+  if (eta_flag) {
+    index = atom->add_custom("rheo_grad_eta", 1, dim);
+    gradn = atom->darray[index];
+  }
+
+  nmax_old = 0;
+  grow_arrays(atom->nmax);
 }
 
 /* ---------------------------------------------------------------------- */
@@ -151,13 +169,17 @@ void ComputeRHEOGrad::compute_peratom()
   double **v = atom->v;
   double *rho = atom->rho;
   double *temperature = atom->temperature;
-  double *eta = atom->viscosity;
   int *status = atom->status;
   int *type = atom->type;
   double *mass = atom->mass;
   int newton = force->newton;
   int dim = domain->dimension;
 
+  int tmp1, tmp2;
+  int index_visc = atom->find_custom("rheo_viscosity", tmp1, tmp2);
+  if (index_visc == -1) error->all(FLERR, "Cannot find rheo viscosity");
+  double *viscosity = atom->dvector[index_visc];
+
   inum = list->inum;
   ilist = list->ilist;
   numneigh = list->numneigh;
@@ -229,7 +251,7 @@ void ComputeRHEOGrad::compute_peratom()
 
         if (rho_flag) drho = rhoi - rhoj;
         if (temperature_flag) dT = temperature[i] - temperature[j];
-        if (eta_flag) deta = eta[i] - eta[j];
+        if (eta_flag) deta = viscosity[i] - viscosity[j];
 
         wp = compute_kernel->calc_dw(i, j, delx, dely, delz, sqrt(rsq));
         dWij = compute_kernel->dWij;
diff --git a/src/RHEO/compute_rheo_grad.h b/src/RHEO/compute_rheo_grad.h
index d86efe47ee..7804ee72dd 100644
--- a/src/RHEO/compute_rheo_grad.h
+++ b/src/RHEO/compute_rheo_grad.h
@@ -42,11 +42,9 @@ class ComputeRHEOGrad : public Compute {
   double **gradr;
   double **gradt;
   double **gradn;
-  int stage;
 
  private:
   int comm_stage, ncomm_grad, ncomm_field, nmax_old;
-  int indexv, indexr, indext, indexn;
   double cut, cutsq, rho0;
   class NeighList *list;
 
diff --git a/src/RHEO/compute_rheo_interface.cpp b/src/RHEO/compute_rheo_interface.cpp
index e76e46d422..ae1bc13831 100644
--- a/src/RHEO/compute_rheo_interface.cpp
+++ b/src/RHEO/compute_rheo_interface.cpp
@@ -18,24 +18,20 @@
 
 #include "compute_rheo_interface.h"
 
-#include "fix_rheo.h"
-#include "compute_rheo_kernel.h"
-#include "fix_store.h"
-#include "fix.h"
-#include <cmath>
-#include <cstring>
 #include "atom.h"
-#include "update.h"
-#include "modify.h"
+#include "comm.h"
 #include "domain.h"
+#include "compute_rheo_kernel.h"
+#include "error.h"
+#include "force.h"
+#include "fix_rheo.h"
+#include "memory.h"
+#include "modify.h"
 #include "neighbor.h"
 #include "neigh_list.h"
 #include "neigh_request.h"
-#include "force.h"
-#include "pair.h"
-#include "comm.h"
-#include "memory.h"
-#include "error.h"
+
+#include <cmath>
 
 using namespace LAMMPS_NS;
 
@@ -44,45 +40,29 @@ using namespace LAMMPS_NS;
 /* ---------------------------------------------------------------------- */
 
 ComputeRHEOInterface::ComputeRHEOInterface(LAMMPS *lmp, int narg, char **arg) :
-  Compute(lmp, narg, arg),
-  id_fix_chi(nullptr)
+  Compute(lmp, narg, arg), fix_rheo(nullptr), compute_kernel(nullptr), fx_m_norm(nullptr),
+  norm(nullptr), normwf(nullptr), chi(nullptr), f_pressure(nullptr), id_fix_pa(nullptr)
 {
-  if (narg != 4) error->all(FLERR,"Illegal compute RHEO/chi command");
-
-  cut = utils::numeric(FLERR,arg[3],false,lmp);
-  cutsq = cut*cut;
-
-  wall_max = sqrt(3)/12.0*cut;
+  if (narg != 3) error->all(FLERR,"Illegal compute rheo/interface command");
 
   nmax = 0;
 
   comm_forward = 3;
   comm_reverse = 4;
-
-  fix_chi = nullptr;
-  norm = nullptr;
-  normwf = nullptr;
-
-  // new id = fix-ID + FIX_STORE_ATTRIBUTE
-  // new fix group = group for this fix
-
-  id_fix_chi = nullptr;
-  std::string fixcmd = id + std::string("_chi");
-  id_fix_chi = new char[fixcmd.size()+1];
-  strcpy(id_fix_chi,fixcmd.c_str());
-  fixcmd += fmt::format(" all STORE peratom 0 {}", 1);
-  modify->add_fix(fixcmd);
-  fix_chi = (FixStore *) modify->fix[modify->nfix-1];
-  chi = fix_chi->vstore;
 }
 
 /* ---------------------------------------------------------------------- */
 
 ComputeRHEOInterface::~ComputeRHEOInterface()
 {
-  if (id_fix_chi && modify->nfix) modify->delete_fix(id_fix_chi);
-  if (modify->nfix) modify->delete_fix("PROPERTY_ATOM_COMP_RHEO_SOLIDS");
+  // Remove custom property if it exists
+  int tmp1, tmp2, index;
+  index = atom->find_custom("rheo_chi", tmp1, tmp2);
+  if (index != -1) atom->remove_custom(index, 1, 0);
 
+  if (id_fix_pa && modify->nfix) modify->delete_fix(id_fix_pa);
+
+  memory->destroy(fx_m_norm);
   memory->destroy(norm);
   memory->destroy(normwf);
 }
@@ -91,42 +71,40 @@ ComputeRHEOInterface::~ComputeRHEOInterface()
 
 void ComputeRHEOInterface::init()
 {
-  // need an occasional full neighbor list
-  int irequest = neighbor->request(this,instance_me);
-  neighbor->requests[irequest]->pair = 0;
-  neighbor->requests[irequest]->compute = 1;
-  neighbor->requests[irequest]->half = 1;
-  neighbor->requests[irequest]->full = 0;
-  //neighbor->requests[irequest]->occasional = 1; //Anticipate needing regulalry
+  compute_kernel = fix_rheo->compute_kernel;
+  cut = fix_rheo->cut;
+  cutsq = cut * cut;
+  wall_max = sqrt(3.0) / 12.0 * cut;
 
-  int icompute = modify->find_compute("rheo_kernel");
-  if (icompute == -1) error->all(FLERR, "Using compute/RHEO/chi without compute/RHEO/kernel");
+  // Create chi array if it doesn'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_old exceeded
 
-  compute_kernel = ((ComputeRHEOKernel *) modify->compute[icompute]);
+  int create_flag = 0;
+  int tmp1, tmp2;
+  int nmax = atom->nmax;
+  int index = atom->find_custom("rheo_chi", tmp1, tmp2);
+  if (index == -1) {
+    index = atom->add_custom("rheo_chi", 1, 0);
+    nmax_old = nmax;
+  }
+  chi = atom->dvector[index];
 
-  //Store persistent per atom quantities - need to be exchanged
-  char **fixarg = new char*[8];
-  fixarg[0] = (char *) "PROPERTY_ATOM_COMP_RHEO_SOLIDS";
-  fixarg[1] = (char *) "all";
-  fixarg[2] = (char *) "property/atom";
-  fixarg[3] = (char *) "d_fx";
-  fixarg[4] = (char *) "d_fy";
-  fixarg[5] = (char *) "d_fz";
-  fixarg[6] = (char *) "ghost";
-  fixarg[7] = (char *) "yes";
-  modify->add_fix(8,fixarg,1);
-  delete [] fixarg;
+  // For fpressure, 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)
 
-  int temp_flag;
-  index_fx = atom->find_custom("fx", temp_flag);
-  if ((index_fx < 0) || (temp_flag != 1))
-      error->all(FLERR, "Compute rheo/solids can't find fix property/atom fx");
-  index_fy = atom->find_custom("fy", temp_flag);
-  if ((index_fy < 0) || (temp_flag != 1))
-      error->all(FLERR, "Compute rheo/solids can't find fix property/atom fy");
-  index_fz = atom->find_custom("fz", temp_flag);
-  if ((index_fz < 0) || (temp_flag != 1))
-      error->all(FLERR, "Compute rheo/solids can't find fix property/atom fz");
+  index = atom->find_custom("rheo_pressure", 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);
+  }
+  f_pressure = atom->darray[index];
+
+  // need an occasional half neighbor list
+  neighbor->add_request(this, NeighConst::REQ_HALF);
 }
 
 /* ---------------------------------------------------------------------- */
@@ -138,6 +116,7 @@ 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;
diff --git a/src/RHEO/compute_rheo_interface.h b/src/RHEO/compute_rheo_interface.h
index 635f190aed..f314c16aa9 100644
--- a/src/RHEO/compute_rheo_interface.h
+++ b/src/RHEO/compute_rheo_interface.h
@@ -27,35 +27,34 @@ namespace LAMMPS_NS {
 class ComputeRHEOInterface : public Compute {
  public:
   ComputeRHEOInterface(class LAMMPS *, int, char **);
-  ~ComputeRHEOInterface();
-  void init();
-  void init_list(int, class NeighList *);
-  void compute_peratom();
-  int pack_forward_comm(int, int *, double *, int, int *);
-  void unpack_forward_comm(int, int, double *);
-  int pack_reverse_comm(int, int, double *);
-  void unpack_reverse_comm(int, int *, double *);
+  ~ComputeRHEOInterface() override;
+  void init() override;
+  void init_list(int, class NeighList *) override;
+  void compute_peratom() 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;
+  double memory_usage() override;
   void correct_v(double *, double *, double *, int, int);
   double correct_rho(int, int);
-  double memory_usage();
   void store_forces();
-  double *chi;
+
+  double *chi, **f_pressure;
 
  private:
-  int nmax;
-  double cut, cutsq, cs2;
+  int nmax_old, comm_stage;
+  double cut, cutsq, cs, wall_max;
+  double **fx_m_norm, *norm, *normwf;
+
+  char *id_fix_pa;
+
   class NeighList *list;
-  double *norm, *normwf, wall_max;
-
+  class FixRHEO *fix_rheo;
   class ComputeRHEOKernel *compute_kernel;
-  char *id_fix_chi;
-  class FixStore *fix_chi;
-
-  int index_fx, index_fy, index_fz;
-  int comm_stage;
 };
 
-}
+}    // namespace LAMMPS_NS
 
 #endif
 #endif
diff --git a/src/RHEO/compute_rheo_kernel.cpp b/src/RHEO/compute_rheo_kernel.cpp
index 15b5f9568b..ac9ea75fc4 100644
--- a/src/RHEO/compute_rheo_kernel.cpp
+++ b/src/RHEO/compute_rheo_kernel.cpp
@@ -55,7 +55,7 @@ Move away from h notation, use cut?
 
 ComputeRHEOKernel::ComputeRHEOKernel(LAMMPS *lmp, int narg, char **arg) :
   Compute(lmp, narg, arg),
-  C(nullptr), C0(nullptr), compute_interface(nullptr);
+  C(nullptr), C0(nullptr), coordination(nullptr), compute_interface(nullptr)
 {
   if (narg != 3) error->all(FLERR,"Illegal compute rheo/kernel command");
 
@@ -71,7 +71,7 @@ ComputeRHEOKernel::~ComputeRHEOKernel()
   // Remove custom property if it exists
   int tmp1, tmp2, index;
   index = atom->find_custom("rheo_coordination", tmp1, tmp2);
-  if (index != -1) atom->remove_custom(index_coord, 1, 0);
+  if (index != -1) atom->remove_custom(index, 0, 0);
 
   memory->destroy(C);
   memory->destroy(C0);
@@ -125,12 +125,14 @@ void ComputeRHEOKernel::init()
 
   int tmp1, tmp2;
   int nmax = atom->nmax;
-  index_coord = atom->find_custom("rheo_coordination", tmp1, tmp2);
-  if (index_coord == -1) {
-    index_coord = atom->add_custom("rheo_coordination", 0, 0);
+  int index = atom->find_custom("rheo_coordination", tmp1, tmp2);
+  if (index == -1) {
+    index = atom->add_custom("rheo_coordination", 0, 0);
     nmax_old = nmax;
   }
+  coordination = atom->ivector[index];
 
+  // Create local arrays for kernel arrays, I can't foresee a reason to print
   comm_forward = 1;
   ncor = 0;
   Mdim = 0;
@@ -170,7 +172,6 @@ int ComputeRHEOKernel::check_corrections(int i)
       corrections = 0;
   }
 
-  int *coordination = atom->ivector[index_coord];
   if (coordination[i] < zmin) corrections = 0;
 
   return corrections;
@@ -503,7 +504,6 @@ void ComputeRHEOKernel::compute_peratom()
   double *mass = atom->mass;
   double *rho = atom->rho;
   int *status = atom->status;
-  int *coordination = atom->ivector[index_coord];
   tagint *tag = atom->tag;
 
   int inum, *ilist, *jlist, *numneigh, **firstneigh;
@@ -742,7 +742,6 @@ int ComputeRHEOKernel::pack_forward_comm(int n, int *list, double *buf,
                                         int /*pbc_flag*/, int * /*pbc*/)
 {
   int i,j,k,m,a,b;
-  coordination = atom->ivector[index_coord];
   m = 0;
   if (correction_order > 0) {
     for (i = 0; i < n; i++) {
@@ -774,7 +773,6 @@ int ComputeRHEOKernel::pack_forward_comm(int n, int *list, double *buf,
 void ComputeRHEOKernel::unpack_forward_comm(int n, int first, double *buf)
 {
   int i, k, m, last,a,b;
-  coordination = atom->ivector[index_coord];
   m = 0;
   last = first + n;
   if (correction_order > 0) {
diff --git a/src/RHEO/compute_rheo_kernel.h b/src/RHEO/compute_rheo_kernel.h
index 346df5a20c..df659c47de 100644
--- a/src/RHEO/compute_rheo_kernel.h
+++ b/src/RHEO/compute_rheo_kernel.h
@@ -42,6 +42,7 @@ class ComputeRHEOKernel : public Compute {
 
   double dWij[3], dWji[3], Wij, Wji;
   int correction_order;
+  int *coordination;
 
  private:
   int solid_flag;
@@ -49,7 +50,7 @@ class ComputeRHEOKernel : public Compute {
   std::unordered_set<tagint> gsl_error_tags;
 
   int kernel_style, zmin, dim, Mdim, ncor;
-  int nmax_old, index_coord;
+  int nmax_old;
   double h, hsq, hinv, hsqinv, pre_w, pre_wp;
   double ***C;
   double *C0;
diff --git a/src/RHEO/compute_rheo_vshift.cpp b/src/RHEO/compute_rheo_vshift.cpp
index 2b7c9394d3..dccd3867cd 100644
--- a/src/RHEO/compute_rheo_vshift.cpp
+++ b/src/RHEO/compute_rheo_vshift.cpp
@@ -21,7 +21,6 @@
 #include "atom.h"
 #include "comm.h"
 #include "compute_rheo_interface.h"
-#include "compute_rheo_grad.h"
 #include "compute_rheo_kernel.h"
 #include "domain.h"
 #include "error.h"
@@ -37,8 +36,8 @@ using namespace LAMMPS_NS;
 /* ---------------------------------------------------------------------- */
 
 ComputeRHEOVShift::ComputeRHEOVShift(LAMMPS *lmp, int narg, char **arg) :
-  Compute(lmp, narg, arg), vshift(nullptr), fix_rheo(nullptr), compute_kernel(nullptr),
-  compute_grad(nullptr), compute_surface(nullptr), compute_interface(nullptr)
+  Compute(lmp, narg, arg), vshift(nullptr), fix_rheo(nullptr), fix_rheo(nullptr),
+  compute_kernel(nullptr), compute_interface(nullptr)
 {
   if (narg != 3) error->all(FLERR,"Illegal compute RHEO/VShift command");
 
@@ -51,11 +50,12 @@ ComputeRHEOVShift::ComputeRHEOVShift(LAMMPS *lmp, int narg, char **arg) :
   // Manually grow if nmax_old exceeded
 
   int tmp1, tmp2;
-  index_vshift = atom->find_custom("rheo_vshift", tmp1, tmp2);
-  if (index_vshift == -1) {
-    index_vshift = atom->add_custom("rheo_vshift", 1, 3);
+  int index = atom->find_custom("rheo_vshift", tmp1, tmp2);
+  if (index == -1) {
+    index = atom->add_custom("rheo_vshift", 1, 3);
     nmax_old = atom->nmax;
   }
+  vshift = atom->dvector[index];
 }
 
 /* ---------------------------------------------------------------------- */
@@ -65,7 +65,7 @@ ComputeRHEOVShift::~ComputeRHEOVShift()
   // Remove custom property if it exists
   int tmp1, tmp2, index;
   index = atom->find_custom("rheo_vshift", tmp1, tmp2);
-  if (index != -1) atom->remove_custom(index_vshift, 1, 3);
+  if (index != -1) atom->remove_custom(index, 1, 3);
 
 }
 
@@ -80,7 +80,6 @@ void ComputeRHEOVShift::init()
     surface_flag = 1;
 
   compute_kernel = fix_rheo->compute_kernel;
-  compute_grad = fix_rheo->compute_grad;
   compute_interface = fix_rheo->compute_interface;
 
   cut = fix_rheo->cut;
@@ -119,7 +118,6 @@ void ComputeRHEOVShift::compute_peratom()
   int *surface = atom->surface;
   double *rho = atom->rho;
   double *mass = atom->mass;
-  double **vshift = atom->darray[index_vshift];
   int newton_pair = force->newton_pair;
 
   inum = list->inum;
@@ -127,7 +125,6 @@ void ComputeRHEOVShift::compute_peratom()
   numneigh = list->numneigh;
   firstneigh = list->firstneigh;
 
-
   if (nmax_old < atom->nmax)
     memory->grow(vshift, atom->nmax, 3, "atom:rheo_vshift");
   nmax_old = atom->nmax;
@@ -232,9 +229,9 @@ void ComputeRHEOVShift::correct_surfaces()
   int nlocal = atom->nlocal;
   int i, a, b;
   int dim = domain->dimension;
-  double **vshift = atom->darray[index_vshift];
 
   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];
@@ -268,7 +265,6 @@ void ComputeRHEOVShift::correct_surfaces()
 int ComputeRHEOVShift::pack_reverse_comm(int n, int first, double *buf)
 {
   int i,m,last;
-  double **vshift = atom->darray[index_vshift];
 
   m = 0;
   last = first + n;
@@ -285,7 +281,6 @@ int ComputeRHEOVShift::pack_reverse_comm(int n, int first, double *buf)
 void ComputeRHEOVShift::unpack_reverse_comm(int n, int *list, double *buf)
 {
   int i,j,m;
-  double **vshift = atom->darray[index_vshift];
 
   m = 0;
   for (i = 0; i < n; i++) {
diff --git a/src/RHEO/compute_rheo_vshift.h b/src/RHEO/compute_rheo_vshift.h
index f93cced31c..af5675fd8d 100644
--- a/src/RHEO/compute_rheo_vshift.h
+++ b/src/RHEO/compute_rheo_vshift.h
@@ -35,17 +35,17 @@ class ComputeRHEOVShift : public Compute {
   void unpack_reverse_comm(int, int *, double *) override;
   double memory_usage() override;
   void correct_surfaces();
+  double **vshift;
 
  private:
   int nmax_old;
   double dtv, cut, cutsq, cutthird;
-  int surface_flag, index_vshift;
+  int surface_flag;
 
   class NeighList *list;
   class FixRHEO *fix_rheo;
   class ComputeRHEOInterface *compute_interface ;
   class ComputeRHEOKernel *compute_kernel;
-  class ComputeRHEOGrad *compute_grad;
 };
 
 }    // namespace LAMMPS_NS
diff --git a/src/RHEO/fix_rheo.cpp b/src/RHEO/fix_rheo.cpp
index 25a97bd52a..e83a95751a 100644
--- a/src/RHEO/fix_rheo.cpp
+++ b/src/RHEO/fix_rheo.cpp
@@ -21,6 +21,7 @@
 #include "atom.h"
 #include "compute_rheo_grad.h"
 #include "compute_rheo_interface.h"
+#include "compute_rheo_surface.h"
 #include "compute_rheo_kernel.h"
 #include "compute_rheo_rhosum.h"
 #include "compute_rheo_vshift.h"
@@ -37,7 +38,7 @@ using namespace FixConst;
 /* ---------------------------------------------------------------------- */
 
 FixRHEO::FixRHEO(LAMMPS *lmp, int narg, char **arg) :
-  Fix(lmp, narg, arg), compute_grad(nullptr), compute_kernel(nullptr),
+  Fix(lmp, narg, arg), compute_grad(nullptr), compute_kernel(nullptr), compute_surface(nullptr),
   compute_interface(nullptr), compute_rhosum(nullptr), compute_vshift(nullptr)
 {
   time_integrate = 1;
@@ -258,7 +259,7 @@ void FixRHEO::initial_integrate(int /*vflag*/)
 
   double **gradr = compute_grad->gradr;
   double **gradv = compute_grad->gradv;
-  double **vshift = compute_vshift->array_atom;
+  double **vshift = compute_vshift->vshift;
 
   int *type = atom->type;
   int *mask = atom->mask;
diff --git a/src/RHEO/fix_rheo.h b/src/RHEO/fix_rheo.h
index 0132f32bcc..89304fe22f 100644
--- a/src/RHEO/fix_rheo.h
+++ b/src/RHEO/fix_rheo.h
@@ -81,6 +81,7 @@ class FixRHEO : public Fix {
   class ComputeRHEOGrad *compute_grad;
   class ComputeRHEOKernel *compute_kernel;
   class ComputeRHEOInterface *compute_interface;
+  class ComputeRHEOSurface *compute_surface;
   class ComputeRHEORhoSum *compute_rhosum;
   class ComputeRHEOVShift *compute_vshift;
 
diff --git a/src/RHEO/fix_rheo_pressure.cpp b/src/RHEO/fix_rheo_pressure.cpp
index 4fa4ece0f9..726c44c32b 100644
--- a/src/RHEO/fix_rheo_pressure.cpp
+++ b/src/RHEO/fix_rheo_pressure.cpp
@@ -36,7 +36,7 @@ enum {NONE, LINEAR, CUBIC, TAITWATER};
 /* ---------------------------------------------------------------------- */
 
 FixRHEOPressure::FixRHEOPressure(LAMMPS *lmp, int narg, char **arg) :
-  Fix(lmp, narg, arg)
+  Fix(lmp, narg, arg), fix_rheo(nullptr), pressure(nullptr)
 {
   if (narg < 4) error->all(FLERR,"Illegal fix command");
 
@@ -74,7 +74,7 @@ FixRHEOPressure::~FixRHEOPressure()
   // Remove custom property if it exists
   int tmp1, tmp2, index;
   index = atom->find_custom("rheo_pressure", tmp1, tmp2);
-  if (index != -1) atom->remove_custom(index_pres, 1, 0);
+  if (index != -1) atom->remove_custom(index, 1, 0);
 }
 
 /* ---------------------------------------------------------------------- */
@@ -115,11 +115,12 @@ void FixRHEOPressure::setup_pre_force(int /*vflag*/)
   // Manually grow if nmax_old exceeded
 
   int tmp1, tmp2;
-  index_pres = atom->find_custom("rheo_pressure", tmp1, tmp2);
-  if (index_pres == -1) {
-    index_pres = atom->add_custom("rheo_pressure", 1, 0);
+  int index = atom->find_custom("rheo_pressure", tmp1, tmp2);
+  if (index == -1) {
+    index = atom->add_custom("rheo_pressure", 1, 0);
     nmax_old = atom->nmax;
   }
+  pressure = atom->dvector[index];
 
   pre_force(0);
 }
@@ -133,7 +134,6 @@ void FixRHEOPressure::pre_force(int /*vflag*/)
   int i;
   double dr, rr3, rho_ratio;
 
-  double *pressure = atom->dvector[index_pres];
   int *mask = atom->mask;
   double *rho = atom->rho;
 
@@ -170,7 +170,6 @@ int FixRHEOPressure::pack_forward_comm(int n, int *list, double *buf,
                                         int /*pbc_flag*/, int * /*pbc*/)
 {
   int i,j,k,m;
-  double *pressure = atom->dvector[index_pres];
   m = 0;
 
   for (i = 0; i < n; i++) {
@@ -185,7 +184,6 @@ int FixRHEOPressure::pack_forward_comm(int n, int *list, double *buf,
 void FixRHEOPressure::unpack_forward_comm(int n, int first, double *buf)
 {
   int i, k, m, last;
-  double *pressure = atom->dvector[index_pres];
 
   m = 0;
   last = first + n;
@@ -194,6 +192,8 @@ void FixRHEOPressure::unpack_forward_comm(int n, int first, double *buf)
   }
 }
 
+/* ---------------------------------------------------------------------- */
+
 double FixRHEOPressure::calculate_p(double rho)
 {
   double rho;
diff --git a/src/RHEO/fix_rheo_pressure.h b/src/RHEO/fix_rheo_pressure.h
index 8272cc38f5..197cab6e5c 100644
--- a/src/RHEO/fix_rheo_pressure.h
+++ b/src/RHEO/fix_rheo_pressure.h
@@ -38,9 +38,10 @@ class FixRHEOPressure : public Fix {
   double calculate_p(double);
  private:
   double c_cubic, csq, rho0, rho0inv;
+  double *pressure;
   int pressure_style;
   int first_flag, last_flag;
-  int nmax_old, index_pressure;
+  int nmax_old;
 
   class FixRHEO *fix_rheo;
 };
diff --git a/src/RHEO/fix_rheo_thermal.cpp b/src/RHEO/fix_rheo_thermal.cpp
index 3ef9fab853..1a5d894ced 100644
--- a/src/RHEO/fix_rheo_thermal.cpp
+++ b/src/RHEO/fix_rheo_thermal.cpp
@@ -37,7 +37,8 @@ enum {NONE, CONSTANT, TYPE};
 /* ---------------------------------------------------------------------- */
 
 FixRHEOThermal::FixRHEOThermal(LAMMPS *lmp, int narg, char **arg) :
-  Fix(lmp, narg, arg), Tc_type(nullptr), kappa_type(nullptr), cv_type(nullptr)
+  Fix(lmp, narg, arg), Tc_type(nullptr), kappa_type(nullptr), cv_type(nullptr),
+  conductivity(nullptr)
 {
   if (narg < 4) error->all(FLERR,"Illegal fix command");
 
@@ -131,7 +132,7 @@ FixRHEOThermal::~FixRHEOThermal()
   // Remove custom property if it exists
   int tmp1, tmp2, index;
   index = atom->find_custom("rheo_conductivity", tmp1, tmp2);
-  if (index != -1) atom->remove_custom(index_cond, 1, 0);
+  if (index != -1) atom->remove_custom(index, 1, 0);
 
   memory->destroy(cv_type);
   memory->destroy(Tc_type);
@@ -199,11 +200,12 @@ void FixRHEOThermal::setup_pre_force(int /*vflag*/)
   // Manually grow if nmax_old exceeded
 
   int tmp1, tmp2;
-  index_cond = atom->find_custom("rheo_conductivity", tmp1, tmp2);
-  if (index_cond == -1) {
-    index_cond = atom->add_custom("rheo_conductivity", 1, 0);
+  index = atom->find_custom("rheo_conductivity", tmp1, tmp2);
+  if (index== -1) {
+    index = atom->add_custom("rheo_conductivity", 1, 0);
     nmax_old = atom->nmax;
   }
+  conductivity = atom->dvector[index];
 
   post_neighbor();
   pre_force(0);
diff --git a/src/RHEO/fix_rheo_thermal.h b/src/RHEO/fix_rheo_thermal.h
index 3b62e86fa8..4f0e89f17c 100644
--- a/src/RHEO/fix_rheo_thermal.h
+++ b/src/RHEO/fix_rheo_thermal.h
@@ -48,11 +48,12 @@ class FixRHEOThermal : public Fix {
   double *Tc_type, Tc;
   double *kappa_type, kappa;
   double dtf, dtv;
+  double *conductivity;
   int Tc_style;
   int cv_style;
   int conductivity_style;
   int first_flag, last_flag;
-  int nmax_old, index_cond;
+  int nmax_old;
 
   class FixRHEO *fix_rheo;
 
diff --git a/src/RHEO/fix_rheo_viscosity.cpp b/src/RHEO/fix_rheo_viscosity.cpp
index 3f779aa3c6..bb7094c43e 100644
--- a/src/RHEO/fix_rheo_viscosity.cpp
+++ b/src/RHEO/fix_rheo_viscosity.cpp
@@ -37,7 +37,7 @@ enum {NONE, CONSTANT, TYPE, POWER};
 /* ---------------------------------------------------------------------- */
 
 FixRHEOViscosity::FixRHEOViscosity(LAMMPS *lmp, int narg, char **arg) :
-  Fix(lmp, narg, arg), eta_type(nullptr)
+  Fix(lmp, narg, arg), fix_rheo(nullptr), eta_type(nullptr), viscosity(nullptr), compute_grad(nullptr)
 {
   if (narg < 4) error->all(FLERR,"Illegal fix command");
 
@@ -142,11 +142,12 @@ void FixRHEOViscosity::setup_pre_force(int /*vflag*/)
   // Manually grow if nmax_old exceeded
 
   int tmp1, tmp2;
-  index_visc = atom->find_custom("rheo_viscosity", tmp1, tmp2);
+  int index = atom->find_custom("rheo_viscosity", tmp1, tmp2);
   if (index_visc == -1) {
-    index_visc = atom->add_custom("rheo_viscosity", 1, 0);
+    index = atom->add_custom("rheo_viscosity", 1, 0);
     nmax_old = atom->nmax;
   }
+  viscosity = atom->dvector[index];
 
   post_neighbor();
   pre_force(0);
@@ -161,7 +162,6 @@ void FixRHEOViscosity::post_neighbor()
   int i;
 
   int *type = atom->type;
-  double *viscosity = atom->dvector[index_visc];
   int *mask = atom->mask;
 
   int nlocal = atom->nlocal;
@@ -190,7 +190,6 @@ void FixRHEOViscosity::pre_force(int /*vflag*/)
   int i, a, b;
   double tmp, gdot;
 
-  double *viscosity = atom->dvector[index_visc];
   int *mask = atom->mask;
   double **gradv = compute_grad->gradv;
 
@@ -232,7 +231,6 @@ int FixRHEOViscosity::pack_forward_comm(int n, int *list, double *buf,
                                         int /*pbc_flag*/, int * /*pbc*/)
 {
   int i,j,k,m;
-  double *viscosity = atom->dvector[index_visc];
   m = 0;
 
   for (i = 0; i < n; i++) {
@@ -247,7 +245,6 @@ int FixRHEOViscosity::pack_forward_comm(int n, int *list, double *buf,
 void FixRHEOViscosity::unpack_forward_comm(int n, int first, double *buf)
 {
   int i, k, m, last;
-  double *viscosity = atom->dvector[index_visc];
 
   m = 0;
   last = first + n;
diff --git a/src/RHEO/fix_rheo_viscosity.h b/src/RHEO/fix_rheo_viscosity.h
index 407170f668..14f8b70de9 100644
--- a/src/RHEO/fix_rheo_viscosity.h
+++ b/src/RHEO/fix_rheo_viscosity.h
@@ -40,9 +40,10 @@ class FixRHEOViscosity : public Fix {
  private:
   double *eta_type, eta;
   double npow, K, gd0, tau0;
+  double *viscosity;
   int viscosity_style;
   int first_flag, last_flag;
-  int nmax_old, index_visc;
+  int nmax_old;
 
   class FixRHEO *fix_rheo;
   class ComputeRHEOGrad *compute_grad;
diff --git a/src/RHEO/pair_rheo.cpp b/src/RHEO/pair_rheo.cpp
index 7d5d0f425f..7b9425cf20 100644
--- a/src/RHEO/pair_rheo.cpp
+++ b/src/RHEO/pair_rheo.cpp
@@ -47,7 +47,7 @@ using namespace MathExtra;
 
 PairRHEO::PairRHEO(LAMMPS *lmp) :
   Pair(lmp), compute_kernel(nullptr), compute_grad(nullptr),
-  compute_interface(nullptr), fix_rheo(nullptr), fix_rheo_pressure(nullptr)
+  compute_interface(nullptr), fix_rheo(nullptr), fix_pressure(nullptr)
 {
   restartinfo = 0;
   single_enable = 0;
@@ -91,8 +91,7 @@ void PairRHEO::compute(int eflag, int vflag)
   double **v = atom->v;
   double **x = atom->x;
   double **f = atom->f;
-  double **f_pressure = fix_rheo->f_pressure; // rewrite later
-  double *pressure = atom->pressure; // rewrite later
+  double **f_pressure = compute_interface->f_pressure;
   double *rho = atom->rho;
   double *mass = atom->mass;
   double *drho = atom->drho;
@@ -105,15 +104,19 @@ void PairRHEO::compute(int eflag, int vflag)
   int *status = atom->status;
 
   int tmp1, tmp2;
-  int index_visc = atom->find_custom("rheo_viscosity", tmp1, tmp2);
-  if (index_visc == -1) error->all(FLERR, "Cannot find rheo viscosity");
-  double *viscosity = atom->dvector[index_visc];
+  int index = atom->find_custom("rheo_viscosity", tmp1, tmp2);
+  if (index == -1) error->all(FLERR, "Cannot find rheo viscosity");
+  double *viscosity = atom->dvector[index];
+
+  index = atom->find_custom("rheo_pressure", tmp1, tmp2);
+  if (index == -1) error->all(FLERR, "Cannot find rheo pressure");
+  double *pressure = atom->dvector[index];
 
   double *conductivity;
   if (thermal_flag) {
-    int index_cond = atom->find_custom("rheo_conductivity", tmp1, tmp2);
-    if (index_cond == -1) error->all(FLERR, "Cannot find rheo conductivity");
-    conductivity = atom->dvector[index_cond];
+    index = atom->find_custom("rheo_conductivity", tmp1, tmp2);
+    if (index == -1) error->all(FLERR, "Cannot find rheo conductivity");
+    conductivity = atom->dvector[index];
   }
 
   int *ilist, *jlist, *numneigh, **firstneigh;
@@ -195,7 +198,7 @@ void PairRHEO::compute(int eflag, int vflag)
         if (fluidi && (!fluidj)) {
           compute_interface->correct_v(vi, vj, i, j);
           rhoj = compute_interface->correct_rho(j, i);
-          Pj = fix_rheo_pressure->calculate_p(rhoj);
+          Pj = fix_pressure->calculate_p(rhoj);
 
           if ((chi[j] > 0.9) && (r < (h * 0.5)))
             fmag = (chi[j] - 0.9) * (h * 0.5 - r) * rho0 * csq * h * rinv;
@@ -413,7 +416,12 @@ void PairRHEO::setup()
 
   fixes = modify->get_fix_by_style("rheo/pressure");
   if (fixes.size() == 0) error->all(FLERR, "Need to define fix rheo/pressure to use pair rheo");
-  fix_rheo_pressure = dynamic_cast<FixRHEOPressure *>(fixes[0]);
+  fix_pressure = dynamic_cast<FixRHEOPressure *>(fixes[0]);
+
+  int tmp1, tmp2;
+  index_pressure = atom->find_custom("rheo_pressure", tmp1, tmp2);
+  if (index_pressure == -1) index_pressure = atom->add_custom("rheo_pressure", 1, 0);
+  else error->all(FLERR, "Cannot find pressure value in pair rheo");
 
   compute_kernel = fix_rheo->compute_kernel;
   compute_grad = fix_rheo->compute_grad;
diff --git a/src/RHEO/pair_rheo.h b/src/RHEO/pair_rheo.h
index aa98e63ddc..49aa1ad025 100644
--- a/src/RHEO/pair_rheo.h
+++ b/src/RHEO/pair_rheo.h
@@ -49,7 +49,7 @@ class PairRHEO : public Pair {
   class ComputeRHEOGrad *compute_grad;
   class ComputeRHEOInterface *compute_interface;
   class FixRHEO *fix_rheo;
-  class FixRHEOPressure *fix_rheo_pressure;
+  class FixRHEOPressure *fix_pressure;
 };
 
 }    // namespace LAMMPS_NS