Ensuring fix order

src/RHEO/fix_rheo.cpp

@@ -197,8 +197,22 @@ void FixRHEO::init()
 
 void FixRHEO::setup_pre_force(int /*vflag*/)
 {
-  //Need to rethink and replan
+  // Check to confirm no accessory fixes are yet defined
+  // FixRHEO must be the first fix
+  // Note: these fixes set this flag in setup_pre_force()
+  if (viscosity_fix_defined || pressure_fix_defined || thermal_fix_defined || surface_fix_defined)
+    error->all(FLERR, "Fix RHEO must be defined before all other RHEO fixes");
+
+  pre_force(0);
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixRHEO::setup()
+{
   // Check to confirm all accessory fixes are defined
+  // Does not ensure fixes correctly cover all atoms (could be a subset group)
+  // Note: these fixes set this flag in setup_pre_force()
   if (!viscosity_fix_defined)
     error->all(FLERR, "Missing fix rheo/viscosity");
 
@@ -216,8 +230,6 @@ void FixRHEO::setup_pre_force(int /*vflag*/)
   viscosity_fix_defined = 0;
   pressure_fix_defined = 0;
   surface_fix_defined = 0;
-
-  pre_force(0);
 }
 
 /* ---------------------------------------------------------------------- */

src/RHEO/fix_rheo_thermal.cpp

@@ -149,7 +149,6 @@ void FixRHEOThermal::init()
   if (fixes.size() == 0) error->all(FLERR, "Need to define fix rheo to use fix rheo/viscosity");
   fix_rheo = dynamic_cast<FixRHEO *>(fixes[0]);
 
-  fix_rheo->thermal_fix_defined = 1;
   if (!fix_rheo->thermal_flag)
     error->all(FLERR, "Need to define thermal setting in fix rheo");
   compute_grad = fix_rheo->compute_grad;
@@ -162,6 +161,8 @@ void FixRHEOThermal::init()
 
 void FixRHEOThermal::setup_pre_force(int /*vflag*/)
 {
+  fix_rheo->thermal_fix_defined = 1;
+
   // Identify whether this is the first/last instance of fix thermal
   // First will handle growing arrays
   // Last will handle communication

src/RHEO/fix_rheo_viscosity.cpp

@@ -102,7 +102,6 @@ void FixRHEOViscosity::init()
   if (fixes.size() == 0) error->all(FLERR, "Need to define fix rheo to use fix rheo/viscosity");
   fix_rheo = dynamic_cast<FixRHEO *>(fixes[0]);
 
-  fix_rheo->viscosity_fix_defined = 1;
   compute_grad = fix_rheo->compute_grad;
 }
 
@@ -110,6 +109,8 @@ void FixRHEOViscosity::init()
 
 void FixRHEOViscosity::setup_pre_force(int /*vflag*/)
 {
+  fix_rheo->viscosity_fix_defined = 1;
+
   // Identify whether this is the first/last instance of fix viscosity
   // First will handle growing arrays
   // Last will handle communication