diff --git a/examples/PACKAGES/rheo/poiseuille/in.rheo.poiseuille b/examples/PACKAGES/rheo/poiseuille/in.rheo.poiseuille
new file mode 100644
index 0000000000..af5728c1a3
--- /dev/null
+++ b/examples/PACKAGES/rheo/poiseuille/in.rheo.poiseuille
@@ -0,0 +1,83 @@
+dimension          2
+units              lj
+atom_style         rheo
+boundary           p p p
+comm_modify        vel yes
+
+
+# ------ Particle Lattice/Resolution Parameters ------ #
+
+variable           L equal 10
+variable           sf equal 0.2
+variable           n equal 1.0/(${sf}^2)
+variable           cut equal 3.5*${sf}
+
+
+# ------ Create simulation box ------ #
+
+region             box block 0 20 -10 10 -0.01 0.01 units box
+create_box         2 box
+lattice            sq ${n}
+
+region             topwall block INF INF  7 10 INF INF units box
+region             block   block INF INF -6.99 6.99 INF INF units box
+region             botwall block INF INF -10 -7 INF INF units box
+
+create_atoms       2 region topwall
+create_atoms       2 region botwall
+create_atoms       1 region block
+
+group              fluid type 1
+group              rig type 2
+
+variable           dr equal 0.1*${sf}
+displace_atoms     fluid random $(0.1*v_sf) ${dr} 0 135414 units box
+
+
+# ------ Model parameters ------ #
+
+variable           rho0 equal 1.0
+variable           cs equal 1.0
+variable           mp equal ${rho0}/${n}
+variable           zeta equal 1.0
+variable	       kappa equal 1.0*${rho0}/${mp}
+variable           fext equal 1e-3/${n}
+variable           eta equal 0.1
+variable           dt_max equal  0.1*${cut}/${cs}/3
+variable           Dr equal 0.05*${cut}*${cs}
+
+mass               1 ${mp}
+mass               2 ${mp}
+set                group all rheo/rho ${rho0}
+set                group all rheo/status 0
+set                group rig rheo/status 2
+timestep           ${dt_max}
+
+pair_style         rheo ${cut} artificial/visc ${zeta} rho/damp ${Dr}
+pair_coeff         * *
+
+
+# ------ Fixes & computes ------ #
+
+fix                1 all rheo ${cut} Quintic 0 shift
+fix                2 all rheo/viscosity constant ${eta}
+fix                3 all rheo/pressure linear
+fix                4 rig setforce 0.0 0.0 0.0
+fix                5 fluid addforce ${fext} 0.0 0.0
+
+compute            1 all rheo/property/atom rho phase
+
+
+# ------ Output & Run ------ #
+
+thermo             200
+thermo_style       custom step time temp press
+
+variable           skin equal 0.2*${cut}
+neighbor           ${skin} bin
+neigh_modify       one 5000
+
+dump               1 all custom 200 atomDump id type x y vx vy fx fy c_1[*]
+
+run                20000
+
diff --git a/src/RHEO/compute_rheo_kernel.cpp b/src/RHEO/compute_rheo_kernel.cpp
index 9cfa86df7e..205ae6fb72 100644
--- a/src/RHEO/compute_rheo_kernel.cpp
+++ b/src/RHEO/compute_rheo_kernel.cpp
@@ -103,10 +103,6 @@ void ComputeRHEOKernel::init()
 {
   neighbor->add_request(this, NeighConst::REQ_FULL);
 
-  auto fixes = modify->get_fix_by_style("rheo");
-  if (fixes.size() == 0) error->all(FLERR, "Need to define fix rheo to use compute rheo/kernel");
-  fix_rheo = dynamic_cast<FixRHEO *>(fixes[0]);
-
   interface_flag = fix_rheo->interface_flag;
   compute_interface = fix_rheo->compute_interface;
 
@@ -147,17 +143,16 @@ void ComputeRHEOKernel::init_list(int /*id*/, NeighList *ptr)
 
 int ComputeRHEOKernel::check_corrections(int i)
 {
-  int corrections = 1;
-
-  if (gsl_error_flag) {
-    // If there were errors, check to see if it occured for this atom
+  // Skip if there were gsl errors for this atom
+  if (gsl_error_flag)
     if (gsl_error_tags.find(atom->tag[i]) != gsl_error_tags.end())
-      corrections = 0;
-  }
+      return 0;
 
-  if (coordination[i] < zmin) corrections = 0;
+  // Skip if undercoordinated
+  if (coordination[i] < zmin)
+    return 0;
 
-  return corrections;
+  return 1;
 }
 
 /* ---------------------------------------------------------------------- */
@@ -165,12 +160,17 @@ int ComputeRHEOKernel::check_corrections(int i)
 double ComputeRHEOKernel::calc_w(int i, int j, double delx, double dely, double delz, double r)
 {
   double w;
+  int corrections_i, corrections_j, corrections;
 
-  int corrections_i = check_corrections(i);
-  int corrections_j = check_corrections(j);
-  int corrections = corrections_i & corrections_j;
+  if (kernel_style != QUINTIC) {
+    corrections_i = check_corrections(i);
+    corrections_j = check_corrections(j);
+    corrections = corrections_i & corrections_j;
+  } else {
+    corrections = 0;
+  }
 
-  if (kernel_style == QUINTIC || !corrections) w = calc_w_quintic(i,j,delx,dely,delz,r);
+  if (!corrections) w = calc_w_quintic(i,j,delx,dely,delz,r);
   else if (kernel_style == CRK0) w = calc_w_crk0(i,j,delx,dely,delz,r);
   else if (kernel_style == CRK1) w = calc_w_crk1(i,j,delx,dely,delz,r);
   else if (kernel_style == CRK2) w = calc_w_crk2(i,j,delx,dely,delz,r);
@@ -183,17 +183,21 @@ double ComputeRHEOKernel::calc_w(int i, int j, double delx, double dely, double
 double ComputeRHEOKernel::calc_dw(int i, int j, double delx, double dely, double delz, double r)
 {
   double wp;
+  int corrections_i, corrections_j;
 
-  int corrections_i = check_corrections(i);
-  int corrections_j = check_corrections(j);
+  if (kernel_style != QUINTIC) {
+    corrections_i = check_corrections(i);
+    corrections_j = check_corrections(j);
+  }
 
   // Calc wp and default dW's, a bit inefficient but can redo later
   wp = calc_dw_quintic(i,j,delx,dely,delz,r,dWij,dWji);
-  if(kernel_style == CRK1) {
-    //check if kernel correction calculated successfully. If not, revert to quintic
+
+  // Overwrite if there are corrections
+  if (kernel_style == CRK1) {
     if (corrections_i) calc_dw_crk1(i,j,delx,dely,delz,r,dWij);
     if (corrections_j) calc_dw_crk1(j,i,-delx,-dely,-delz,r,dWji);
-  } else if(kernel_style == CRK2) {
+  } else if (kernel_style == CRK2) {
     if (corrections_i) calc_dw_crk2(i,j,delx,dely,delz,r,dWij);
     if (corrections_j) calc_dw_crk2(j,i,-delx,-dely,-delz,r,dWji);
   }
diff --git a/src/RHEO/compute_rheo_property_atom.cpp b/src/RHEO/compute_rheo_property_atom.cpp
index 7682552abe..3cd5d468b2 100644
--- a/src/RHEO/compute_rheo_property_atom.cpp
+++ b/src/RHEO/compute_rheo_property_atom.cpp
@@ -52,7 +52,7 @@ ComputeRHEOPropertyAtom::ComputeRHEOPropertyAtom(LAMMPS *lmp, int narg, char **a
   if (nvalues == 1) size_peratom_cols = 0;
   else size_peratom_cols = nvalues;
 
-  thermal_flag, interface_flag, surface_flag, shift_flag = 0;
+  thermal_flag = interface_flag = surface_flag = shift_flag = 0;
 
   // parse input values
   // customize a new keyword by adding to if statement
@@ -132,7 +132,7 @@ ComputeRHEOPropertyAtom::~ComputeRHEOPropertyAtom()
 
 void ComputeRHEOPropertyAtom::init()
 {
-  auto fixes = modify->get_fix_by_style("rheo");
+  auto fixes = modify->get_fix_by_style("^rheo$");
   if (fixes.size() == 0) error->all(FLERR, "Need to define fix rheo to use fix rheo/viscosity");
   fix_rheo = dynamic_cast<FixRHEO *>(fixes[0]);
 
diff --git a/src/RHEO/fix_rheo.cpp b/src/RHEO/fix_rheo.cpp
index ac870affd5..48eece239a 100644
--- a/src/RHEO/fix_rheo.cpp
+++ b/src/RHEO/fix_rheo.cpp
@@ -195,7 +195,7 @@ void FixRHEO::init()
   dtv = update->dt;
   dtf = 0.5 * update->dt * force->ftm2v;
 
-  if (modify->get_fix_by_style("rheo").size() > 1)
+  if (modify->get_fix_by_style("^rheo$").size() > 1)
     error->all(FLERR,"Can only specify one instance of fix rheo");
 }
 
@@ -258,6 +258,8 @@ void FixRHEO::setup(int /*vflag*/)
     error->one(FLERR, "Fix rheo/viscosity does not fully cover all atoms");
   if (!t_coverage_flag)
     error->one(FLERR, "Fix rheo/thermal does not fully cover all atoms");
+
+  pre_force(0);
 }
 
 /* ---------------------------------------------------------------------- */
diff --git a/src/RHEO/fix_rheo_pressure.cpp b/src/RHEO/fix_rheo_pressure.cpp
index 45557794cd..63a6995646 100644
--- a/src/RHEO/fix_rheo_pressure.cpp
+++ b/src/RHEO/fix_rheo_pressure.cpp
@@ -92,7 +92,7 @@ int FixRHEOPressure::setmask()
 
 void FixRHEOPressure::init()
 {
-  auto fixes = modify->get_fix_by_style("rheo");
+  auto fixes = modify->get_fix_by_style("^rheo$");
   if (fixes.size() == 0) error->all(FLERR, "Need to define fix rheo to use fix rheo/pressure");
   fix_rheo = dynamic_cast<FixRHEO *>(fixes[0]);
 
@@ -161,7 +161,7 @@ void FixRHEOPressure::pre_force(int /*vflag*/)
     }
   }
 
-  if (last_flag && comm_forward) comm->forward_comm(this);
+  if (comm_forward) comm->forward_comm(this);
 }
 
 /* ---------------------------------------------------------------------- */
diff --git a/src/RHEO/fix_rheo_pressure.h b/src/RHEO/fix_rheo_pressure.h
index c257f1dbfb..e8f7f3cb88 100644
--- a/src/RHEO/fix_rheo_pressure.h
+++ b/src/RHEO/fix_rheo_pressure.h
@@ -41,7 +41,6 @@ class FixRHEOPressure : public Fix {
   double c_cubic, csq, rho0, rho0inv;
   double *pressure;
   int pressure_style;
-  int first_flag, last_flag;
   int nmax_store;
 
   class FixRHEO *fix_rheo;
diff --git a/src/RHEO/fix_rheo_thermal.cpp b/src/RHEO/fix_rheo_thermal.cpp
index 0d8909b908..bd7a22ce1e 100644
--- a/src/RHEO/fix_rheo_thermal.cpp
+++ b/src/RHEO/fix_rheo_thermal.cpp
@@ -155,7 +155,7 @@ int FixRHEOThermal::setmask()
 
 void FixRHEOThermal::init()
 {
-  auto fixes = modify->get_fix_by_style("rheo");
+  auto fixes = modify->get_fix_by_style("^rheo$");
   if (fixes.size() == 0) error->all(FLERR, "Need to define fix rheo to use fix rheo/viscosity");
   fix_rheo = dynamic_cast<FixRHEO *>(fixes[0]);
 
diff --git a/src/RHEO/fix_rheo_viscosity.cpp b/src/RHEO/fix_rheo_viscosity.cpp
index 4d70ffaca3..7d915c9b93 100644
--- a/src/RHEO/fix_rheo_viscosity.cpp
+++ b/src/RHEO/fix_rheo_viscosity.cpp
@@ -102,7 +102,7 @@ int FixRHEOViscosity::setmask()
 
 void FixRHEOViscosity::init()
 {
-  auto fixes = modify->get_fix_by_style("rheo");
+  auto fixes = modify->get_fix_by_style("^rheo$");
   if (fixes.size() == 0) error->all(FLERR, "Need to define fix rheo to use fix rheo/viscosity");
   fix_rheo = dynamic_cast<FixRHEO *>(fixes[0]);
 
diff --git a/src/RHEO/pair_rheo.cpp b/src/RHEO/pair_rheo.cpp
index 03e859316f..8e76a6d413 100644
--- a/src/RHEO/pair_rheo.cpp
+++ b/src/RHEO/pair_rheo.cpp
@@ -356,7 +356,9 @@ void PairRHEO::settings(int narg, char **arg)
 {
   if (narg < 1) error->all(FLERR,"Illegal pair_style command");
 
-  int iarg = 0;
+  h = utils::numeric(FLERR,arg[0],false,lmp);
+printf("settings\n");
+  int iarg = 1;
   while (iarg < narg) {
     if (strcmp(arg[iarg], "rho/damp") == 0) {
       if (iarg + 1 >= narg) error->all(FLERR,"Illegal pair_style command");
@@ -419,10 +421,13 @@ void PairRHEO::setup()
   compute_grad = fix_rheo->compute_grad;
   compute_interface = fix_rheo->compute_interface;
   thermal_flag = fix_rheo->thermal_flag;
-  h = fix_rheo->h;
   csq = fix_rheo->csq;
   rho0 = fix_rheo->rho0;
+printf("setup\n");
+  if (h != fix_rheo->h)
+    error->all(FLERR, "Pair rheo cutoff {} does not agree with fix rheo cutoff {}", h, fix_rheo->h);
 
+  hsq = h * h;
   hinv = 1.0 / h;
   hinv3 = hinv * 3.0;
   cs = sqrt(csq);
@@ -450,6 +455,6 @@ double PairRHEO::init_one(int i, int j)
   if (setflag[i][j] == 0) {
       error->all(FLERR,"All pair rheo coeffs are not set");
   }
-
+printf("init one\n");
   return h;
 }
diff --git a/src/compute_rheo_property_atom.cpp b/src/compute_rheo_property_atom.cpp
index 7682552abe..3cd5d468b2 100644
--- a/src/compute_rheo_property_atom.cpp
+++ b/src/compute_rheo_property_atom.cpp
@@ -52,7 +52,7 @@ ComputeRHEOPropertyAtom::ComputeRHEOPropertyAtom(LAMMPS *lmp, int narg, char **a
   if (nvalues == 1) size_peratom_cols = 0;
   else size_peratom_cols = nvalues;
 
-  thermal_flag, interface_flag, surface_flag, shift_flag = 0;
+  thermal_flag = interface_flag = surface_flag = shift_flag = 0;
 
   // parse input values
   // customize a new keyword by adding to if statement
@@ -132,7 +132,7 @@ ComputeRHEOPropertyAtom::~ComputeRHEOPropertyAtom()
 
 void ComputeRHEOPropertyAtom::init()
 {
-  auto fixes = modify->get_fix_by_style("rheo");
+  auto fixes = modify->get_fix_by_style("^rheo$");
   if (fixes.size() == 0) error->all(FLERR, "Need to define fix rheo to use fix rheo/viscosity");
   fix_rheo = dynamic_cast<FixRHEO *>(fixes[0]);