diff --git a/src/RHEO/compute_rheo_kernel.cpp b/src/RHEO/compute_rheo_kernel.cpp
index 453f2b9028..658abd3189 100644
--- a/src/RHEO/compute_rheo_kernel.cpp
+++ b/src/RHEO/compute_rheo_kernel.cpp
@@ -137,6 +137,9 @@ void ComputeRHEOKernel::init()
     }
   }
 
+  if (correction_order != -1)
+    fix_rheo->coordination_flag = 1;
+
   nmax_store = atom->nmax;
   memory->create(coordination, nmax_store, "rheo:coordination");
   if (kernel_style == RK0) {
diff --git a/src/RHEO/compute_rheo_property_atom.cpp b/src/RHEO/compute_rheo_property_atom.cpp
index fb5848ac75..fd05c0373e 100644
--- a/src/RHEO/compute_rheo_property_atom.cpp
+++ b/src/RHEO/compute_rheo_property_atom.cpp
@@ -79,7 +79,7 @@ ComputeRHEOPropertyAtom::ComputeRHEOPropertyAtom(LAMMPS *lmp, int narg, char **a
     size_peratom_cols = nvalues;
 
   pressure_flag = thermal_flag = interface_flag = 0;
-  surface_flag = shift_flag = shell_flag = 0;
+  surface_flag = shift_flag = shell_flag = coordination_flag = 0;
 
   // parse input values
   // customize a new keyword by adding to if statement
@@ -109,6 +109,7 @@ ComputeRHEOPropertyAtom::ComputeRHEOPropertyAtom(LAMMPS *lmp, int narg, char **a
       surface_flag = 1;
       pack_choice[i] = &ComputeRHEOPropertyAtom::pack_surface_divr;
     } else if (strcmp(arg[iarg], "coordination") == 0) {
+      coordination_flag = 1;
       pack_choice[i] = &ComputeRHEOPropertyAtom::pack_coordination;
     } else if (strcmp(arg[iarg], "pressure") == 0) {
       pressure_flag = 1;
@@ -223,6 +224,11 @@ void ComputeRHEOPropertyAtom::compute_peratom()
 {
   invoked_peratom = update->ntimestep;
 
+  // calculate optional values, if needed
+
+  if (coordination_flag && !(fix_rheo->coordination_flag))
+    compute_kernel->compute_coordination();
+
   // grow vector or array if necessary
 
   if (atom->nmax > nmax) {
diff --git a/src/RHEO/compute_rheo_property_atom.h b/src/RHEO/compute_rheo_property_atom.h
index 14f28075a3..745423b4ae 100644
--- a/src/RHEO/compute_rheo_property_atom.h
+++ b/src/RHEO/compute_rheo_property_atom.h
@@ -36,7 +36,7 @@ class ComputeRHEOPropertyAtom : public Compute {
  private:
   int nvalues, nmax;
   int pressure_flag, thermal_flag, interface_flag;
-  int surface_flag, shift_flag, shell_flag;
+  int surface_flag, shift_flag, shell_flag, coordination_flag;
   int *avec_index;
   int *col_index, *col_t_index;
   double *buf;
diff --git a/src/RHEO/compute_rheo_surface.cpp b/src/RHEO/compute_rheo_surface.cpp
index 4ff9587f6c..b4f7fe7c1c 100644
--- a/src/RHEO/compute_rheo_surface.cpp
+++ b/src/RHEO/compute_rheo_surface.cpp
@@ -81,6 +81,8 @@ void ComputeRHEOSurface::init()
   threshold_splash = fix_rheo->zmin_splash;
   interface_flag = fix_rheo->interface_flag;
 
+  fix_rheo->coordination_flag = 1;
+
   cutsq = cut * cut;
 
   // need an occasional half neighbor list
diff --git a/src/RHEO/fix_rheo.cpp b/src/RHEO/fix_rheo.cpp
index 4a25a46501..a978333aba 100644
--- a/src/RHEO/fix_rheo.cpp
+++ b/src/RHEO/fix_rheo.cpp
@@ -70,6 +70,7 @@ FixRHEO::FixRHEO(LAMMPS *lmp, int narg, char **arg) :
   surface_flag = 0;
   oxidation_flag = 0;
   self_mass_flag = 0;
+  coordination_flag = 0;
 
   int i;
   int n = atom->ntypes;
@@ -398,34 +399,39 @@ void FixRHEO::initial_integrate(int /*vflag*/)
 
 void FixRHEO::pre_force(int /*vflag*/)
 {
-  compute_kernel->compute_coordination();    // Needed for rho sum
+  if (coordination_flag)
+    compute_kernel->compute_coordination();
 
-  if (rhosum_flag) compute_rhosum->compute_peratom();
+  if (rhosum_flag)
+    compute_rhosum->compute_peratom();
 
   compute_kernel->compute_peratom();
 
-  if (interface_flag) {
-    // Note on first setup, have no forces for pressure to reference
+  // Note on first setup, have no forces for pressure to reference
+  if (interface_flag)
     compute_interface->compute_peratom();
-  }
 
   // No need to forward v, rho, or T for compute_grad since already done
   compute_grad->compute_peratom();
   compute_grad->forward_gradients();
 
-  if (shift_flag) compute_vshift->compute_peratom();
+  // Depends on NO_SHIFT status
+  if (shift_flag)
+    compute_vshift->compute_peratom();
 
   // Remove temporary options
   int *mask = atom->mask;
   int *status = atom->rheo_status;
   int nall = atom->nlocal + atom->nghost;
   for (int i = 0; i < nall; i++)
-    if (mask[i] & groupbit) status[i] &= OPTIONSMASK;
+    if (mask[i] & groupbit)
+      status[i] &= OPTIONSMASK;
 
   // Calculate surfaces, update status
   if (surface_flag) {
     compute_surface->compute_peratom();
-    if (shift_flag) compute_vshift->correct_surfaces();
+    if (shift_flag)
+      compute_vshift->correct_surfaces();
   }
 }
 
diff --git a/src/RHEO/fix_rheo.h b/src/RHEO/fix_rheo.h
index 24e86b21d6..4b75299f27 100644
--- a/src/RHEO/fix_rheo.h
+++ b/src/RHEO/fix_rheo.h
@@ -53,6 +53,7 @@ class FixRHEO : public Fix {
   int interface_flag;
   int surface_flag;
   int oxidation_flag;
+  int coordination_flag;
 
   int viscosity_fix_defined;
   int pressure_fix_defined;