diff --git a/lib/Grid.h b/lib/Grid.h
index c14be568..6dbd5f04 100644
--- a/lib/Grid.h
+++ b/lib/Grid.h
@@ -47,9 +47,11 @@
 #include <parallelIO/BinaryIO.h>
 #include <qcd/QCD.h>
 #include <parallelIO/NerscIO.h>
-#include <qcd/hmc/NerscCheckpointer.h>
-
 #include <Init.h>
 
+#include <qcd/hmc/NerscCheckpointer.h>
+#include <qcd/hmc/HmcRunner.h>
+
+
 
 #endif
diff --git a/lib/qcd/hmc/HmcRunner.h b/lib/qcd/hmc/HmcRunner.h
new file mode 100644
index 00000000..4f3be698
--- /dev/null
+++ b/lib/qcd/hmc/HmcRunner.h
@@ -0,0 +1,113 @@
+#ifndef HMC_RUNNER
+#define HMC_RUNNER
+
+namespace Grid{
+  namespace QCD{
+
+
+class NerscHmcRunner {
+public:
+
+  enum StartType_t { ColdStart, HotStart, TepidStart, CheckpointStart };
+
+  ActionSet<LatticeGaugeField> TheAction;
+
+  GridCartesian         * UGrid   ;
+  GridCartesian         * FGrid   ;
+  GridRedBlackCartesian * UrbGrid ;
+  GridRedBlackCartesian * FrbGrid ;
+
+  virtual void BuildTheAction (int argc, char **argv) = 0;
+
+  
+  void Run (int argc, char  **argv){
+
+    StartType_t StartType = HotStart;
+
+    std::string arg;
+
+    if( GridCmdOptionExists(argv,argv+argc,"--StartType") ){
+      arg = GridCmdOptionPayload(argv,argv+argc,"--StartType");
+      if ( arg == "HotStart" ) { StartType = HotStart; }
+      else if ( arg == "ColdStart" ) { StartType = ColdStart; }
+      else if ( arg == "TepidStart" ) { StartType = TepidStart; }
+      else if ( arg == "CheckpointStart" ) { StartType = CheckpointStart; }
+      else assert(0);
+    }
+
+    int StartTraj = 0;
+    if( GridCmdOptionExists(argv,argv+argc,"--StartTrajectory") ){
+      arg= GridCmdOptionPayload(argv,argv+argc,"--StartTrajectory");
+      std::vector<int> ivec(0);
+      GridCmdOptionIntVector(arg,ivec);
+      StartTraj = ivec[0];
+    }    
+
+    int NumTraj = 1;
+    if( GridCmdOptionExists(argv,argv+argc,"--Trajectories") ){
+      arg= GridCmdOptionPayload(argv,argv+argc,"--Trajectories");
+      std::vector<int> ivec(0);
+      GridCmdOptionIntVector(arg,ivec);
+      NumTraj = ivec[0];
+    }
+
+    // Create integrator
+    typedef MinimumNorm2<LatticeGaugeField>  IntegratorType;// change here to change the algorithm
+    IntegratorParameters MDpar(20);
+    IntegratorType MDynamics(UGrid,MDpar, TheAction);
+
+    // Checkpoint strategy
+    NerscHmcCheckpointer<LatticeGaugeField> Checkpoint(std::string("ckpoint_lat"),std::string("ckpoint_rng"),1);
+    PlaquetteLogger<LatticeGaugeField> PlaqLog(std::string("plaq"));
+
+    HMCparameters HMCpar;
+    HMCpar.StartTrajectory = StartTraj;
+    HMCpar.Trajectories    = NumTraj;
+    
+    GridSerialRNG    sRNG;
+    GridParallelRNG  pRNG(UGrid);
+    LatticeGaugeField  U(UGrid);
+
+    std::vector<int> SerSeed({1,2,3,4,5});
+    std::vector<int> ParSeed({6,7,8,9,10});
+
+    if ( StartType == HotStart ) {
+      // Hot start
+      HMCpar.NoMetropolisUntil =0;
+      HMCpar.MetropolisTest = true;
+      sRNG.SeedFixedIntegers(SerSeed);
+      pRNG.SeedFixedIntegers(ParSeed);
+      SU3::HotConfiguration(pRNG, U);
+    } else if ( StartType == ColdStart ) { 
+      // Cold start
+      HMCpar.NoMetropolisUntil =0;
+      HMCpar.MetropolisTest = true;
+      sRNG.SeedFixedIntegers(SerSeed);
+      pRNG.SeedFixedIntegers(ParSeed);
+      SU3::ColdConfiguration(pRNG, U);
+    } else if ( StartType == TepidStart ) {       
+      // Tepid start
+      HMCpar.NoMetropolisUntil =0;
+      HMCpar.MetropolisTest = true;
+      sRNG.SeedFixedIntegers(SerSeed);
+      pRNG.SeedFixedIntegers(ParSeed);
+      SU3::TepidConfiguration(pRNG, U);
+    } else if ( StartType == CheckpointStart ) { 
+      HMCpar.NoMetropolisUntil =0;
+      HMCpar.MetropolisTest = true;
+      // CheckpointRestart
+      Checkpoint.CheckpointRestore(StartTraj, U, sRNG, pRNG);
+    }
+
+    HybridMonteCarlo<LatticeGaugeField,IntegratorType>  HMC(HMCpar, MDynamics,sRNG,pRNG,U);
+    HMC.AddObservable(&Checkpoint);
+    HMC.AddObservable(&PlaqLog);
+    
+    // Run it
+    HMC.evolve();
+    
+  }
+  
+};
+}}
+#endif
diff --git a/tests/Test_hmc_EODWFRatio.cc b/tests/Test_hmc_EODWFRatio.cc
index 435757a6..a7eecb69 100644
--- a/tests/Test_hmc_EODWFRatio.cc
+++ b/tests/Test_hmc_EODWFRatio.cc
@@ -1,4 +1,5 @@
 #include "Grid.h"
+
 using namespace std;
 using namespace Grid;
 using namespace Grid::QCD;
@@ -6,20 +7,17 @@ using namespace Grid::QCD;
 namespace Grid { 
   namespace QCD { 
 
-class NerscHmcRunner {
+
+class HmcRunner : public NerscHmcRunner {
 public:
 
-  enum StartType_t { ColdStart, HotStart, TepidStart, CheckpointStart };
-
-  ActionSet<LatticeGaugeField> TheAction;
-
-  GridCartesian         * UGrid   ;
-  GridCartesian         * FGrid   ;
-  GridRedBlackCartesian * UrbGrid ;
-  GridRedBlackCartesian * FrbGrid ;
-
   void BuildTheAction (int argc, char **argv)
+
   {
+    typedef WilsonImplR ImplPolicy;
+    typedef DomainWallFermionR FermionAction;
+    typedef typename FermionAction::FermionField FermionField;
+
     const int Ls = 8;
 
     UGrid   = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd,vComplex::Nsimd()),GridDefaultMpi());
@@ -37,11 +35,11 @@ public:
     Real mass=0.04;
     Real pv  =1.0;
     RealD M5=1.5;
-    DomainWallFermionR DenOp(U,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5);
-    DomainWallFermionR NumOp(U,*FGrid,*FrbGrid,*UGrid,*UrbGrid,pv,M5);
+    FermionAction DenOp(U,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5);
+    FermionAction NumOp(U,*FGrid,*FrbGrid,*UGrid,*UrbGrid,pv,M5);
   
-    ConjugateGradient<LatticeFermion>  CG(1.0e-8,10000);
-    TwoFlavourEvenOddRatioPseudoFermionAction<WilsonImplR> Nf2(NumOp, DenOp,CG,CG);
+    ConjugateGradient<FermionField>  CG(1.0e-8,10000);
+    TwoFlavourEvenOddRatioPseudoFermionAction<ImplPolicy> Nf2(NumOp, DenOp,CG,CG);
   
     //Collect actions
     ActionLevel<LatticeGaugeField> Level1;
@@ -51,95 +49,7 @@ public:
 
     Run(argc,argv);
   };
-  
-  void Run (int argc, char  **argv){
 
-    StartType_t StartType = HotStart;
-
-    std::string arg;
-
-    if( GridCmdOptionExists(argv,argv+argc,"--StartType") ){
-      arg = GridCmdOptionPayload(argv,argv+argc,"--StartType");
-      if ( arg == "HotStart" ) { StartType = HotStart; }
-      else if ( arg == "ColdStart" ) { StartType = ColdStart; }
-      else if ( arg == "TepidStart" ) { StartType = TepidStart; }
-      else if ( arg == "CheckpointStart" ) { StartType = CheckpointStart; }
-      else assert(0);
-    }
-
-    int StartTraj = 0;
-    if( GridCmdOptionExists(argv,argv+argc,"--StartTrajectory") ){
-      arg= GridCmdOptionPayload(argv,argv+argc,"--StartTrajectory");
-      std::vector<int> ivec(0);
-      GridCmdOptionIntVector(arg,ivec);
-      StartTraj = ivec[0];
-    }    
-
-    int NumTraj = 1;
-    if( GridCmdOptionExists(argv,argv+argc,"--Trajectories") ){
-      arg= GridCmdOptionPayload(argv,argv+argc,"--Trajectories");
-      std::vector<int> ivec(0);
-      GridCmdOptionIntVector(arg,ivec);
-      NumTraj = ivec[0];
-    }
-
-    // Create integrator
-    typedef MinimumNorm2<LatticeGaugeField>  IntegratorType;// change here to change the algorithm
-    IntegratorParameters MDpar(20);
-    IntegratorType MDynamics(UGrid,MDpar, TheAction);
-
-    // Checkpoint strategy
-    NerscHmcCheckpointer<LatticeGaugeField> Checkpoint(std::string("ckpoint_lat"),std::string("ckpoint_rng"),1);
-    PlaquetteLogger<LatticeGaugeField> PlaqLog(std::string("plaq"));
-
-    HMCparameters HMCpar;
-    HMCpar.StartTrajectory = StartTraj;
-    HMCpar.Trajectories    = NumTraj;
-    
-    GridSerialRNG    sRNG;
-    GridParallelRNG  pRNG(UGrid);
-    LatticeGaugeField  U(UGrid);
-
-    std::vector<int> SerSeed({1,2,3,4,5});
-    std::vector<int> ParSeed({6,7,8,9,10});
-
-    if ( StartType == HotStart ) {
-      // Hot start
-      HMCpar.NoMetropolisUntil =0;
-      HMCpar.MetropolisTest = true;
-      sRNG.SeedFixedIntegers(SerSeed);
-      pRNG.SeedFixedIntegers(ParSeed);
-      SU3::HotConfiguration(pRNG, U);
-    } else if ( StartType == ColdStart ) { 
-      // Cold start
-      HMCpar.NoMetropolisUntil =0;
-      HMCpar.MetropolisTest = true;
-      sRNG.SeedFixedIntegers(SerSeed);
-      pRNG.SeedFixedIntegers(ParSeed);
-      SU3::ColdConfiguration(pRNG, U);
-    } else if ( StartType == TepidStart ) {       
-      // Tepid start
-      HMCpar.NoMetropolisUntil =0;
-      HMCpar.MetropolisTest = true;
-      sRNG.SeedFixedIntegers(SerSeed);
-      pRNG.SeedFixedIntegers(ParSeed);
-      SU3::TepidConfiguration(pRNG, U);
-    } else if ( StartType == CheckpointStart ) { 
-      HMCpar.NoMetropolisUntil =0;
-      HMCpar.MetropolisTest = true;
-      // CheckpointRestart
-      Checkpoint.CheckpointRestore(StartTraj, U, sRNG, pRNG);
-    }
-
-    HybridMonteCarlo<LatticeGaugeField,IntegratorType>  HMC(HMCpar, MDynamics,sRNG,pRNG,U);
-    HMC.AddObservable(&Checkpoint);
-    HMC.AddObservable(&PlaqLog);
-    
-    // Run it
-    HMC.evolve();
-    
-  }
-  
 };
 
 }}
@@ -151,7 +61,7 @@ int main (int argc, char ** argv)
   int threads = GridThread::GetThreads();
   std::cout<<GridLogMessage << "Grid is setup to use "<<threads<<" threads"<<std::endl;
 
-  NerscHmcRunner TheHMC;
+  HmcRunner TheHMC;
   
   TheHMC.BuildTheAction(argc,argv);
 
diff --git a/tests/Test_hmc_EODWFRatio_Gparity.cc b/tests/Test_hmc_EODWFRatio_Gparity.cc
index b5166362..0ea2ca3f 100644
--- a/tests/Test_hmc_EODWFRatio_Gparity.cc
+++ b/tests/Test_hmc_EODWFRatio_Gparity.cc
@@ -1,73 +1,69 @@
 #include "Grid.h"
 
-
 using namespace std;
 using namespace Grid;
 using namespace Grid::QCD;
 
+namespace Grid { 
+  namespace QCD { 
+
+
+class HmcRunner : public NerscHmcRunner {
+public:
+
+  void BuildTheAction (int argc, char **argv)
+
+  {
+    typedef GparityWilsonImplR ImplPolicy;
+    typedef GparityDomainWallFermionR FermionAction;
+
+    typedef typename FermionAction::FermionField FermionField;
+
+    const int Ls = 8;
+
+    UGrid   = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd,vComplex::Nsimd()),GridDefaultMpi());
+    UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid);
+  
+    FGrid   = SpaceTimeGrid::makeFiveDimGrid(Ls,UGrid);
+    FrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,UGrid);
+
+    // temporarily need a gauge field
+    LatticeGaugeField  U(UGrid);
+
+    // Gauge action
+    WilsonGaugeActionR Waction(5.6);
+
+    Real mass=0.04;
+    Real pv  =1.0;
+    RealD M5=1.5;
+    FermionAction DenOp(U,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5);
+    FermionAction NumOp(U,*FGrid,*FrbGrid,*UGrid,*UrbGrid,pv,M5);
+  
+    ConjugateGradient<FermionField>  CG(1.0e-8,10000);
+    TwoFlavourEvenOddRatioPseudoFermionAction<ImplPolicy> Nf2(NumOp, DenOp,CG,CG);
+  
+    //Collect actions
+    ActionLevel<LatticeGaugeField> Level1;
+    Level1.push_back(&Nf2);
+    Level1.push_back(&Waction);
+    TheAction.push_back(Level1);
+
+    Run(argc,argv);
+  };
+
+};
+
+}}
+
 int main (int argc, char ** argv)
 {
   Grid_init(&argc,&argv);
 
-  const int Ls = 8;
+  int threads = GridThread::GetThreads();
+  std::cout<<GridLogMessage << "Grid is setup to use "<<threads<<" threads"<<std::endl;
 
-  GridCartesian         * UGrid   = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd,vComplex::Nsimd()),GridDefaultMpi());
-  GridRedBlackCartesian * UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid);
-
-  GridCartesian         * FGrid   = SpaceTimeGrid::makeFiveDimGrid(Ls,UGrid);
-  GridRedBlackCartesian * FrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,UGrid);
-
-  GridSerialRNG    sRNG;
-  GridParallelRNG  pRNG(UGrid);
-  sRNG.SeedRandomDevice();
-  pRNG.SeedRandomDevice();
-
-  LatticeLorentzColourMatrix     U(UGrid);
-
-  SU3::HotConfiguration(pRNG, U);
-
-  // simplify template declaration? Strip the lorentz from the second template
-  WilsonGaugeActionR Waction(5.6);
-
-  Real mass=0.04;
-  Real pv  =1.0;
-  RealD M5=1.5;
-
-  typedef typename GparityDomainWallFermionR::FermionField FermionField;
-
-  const int nu = 3;
-  std::vector<int> twists(Nd,0);
-  twists[nu] = 1;
-  GparityDomainWallFermionR::ImplParams params;
-  params.twists = twists;
-
-  GparityDomainWallFermionR DenOp(U,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,params);
-  GparityDomainWallFermionR NumOp(U,*FGrid,*FrbGrid,*UGrid,*UrbGrid,pv,M5,params);
+  HmcRunner TheHMC;
   
-  ConjugateGradient<FermionField>  CG(1.0e-8,10000);
-
-  TwoFlavourEvenOddRatioPseudoFermionAction<GparityWilsonImplR> Nf2(NumOp, DenOp,CG,CG);
-  
-  //Collect actions
-  ActionLevel<LatticeGaugeField> Level1;
-  Level1.push_back(&Nf2);
-  Level1.push_back(&Waction);
-  ActionSet<LatticeGaugeField> FullSet;
-  FullSet.push_back(Level1);
-
-  // Create integrator
-  typedef MinimumNorm2<LatticeGaugeField>  IntegratorType;// change here to change the algorithm
-  IntegratorParameters MDpar(20);
-  IntegratorType MDynamics(UGrid,MDpar, FullSet);
-
-  // Create HMC
-  NerscHmcCheckpointer<LatticeGaugeField> Checkpoint(std::string("ckpoint_lat"),std::string("ckpoint_rng"),1);
-  HMCparameters HMCpar;
-  HybridMonteCarlo<LatticeGaugeField,IntegratorType>  HMC(HMCpar, MDynamics,sRNG,pRNG,U);
-  HMC.AddObservable(&Checkpoint);
-
-  // Run it
-  HMC.evolve();
-
+  TheHMC.BuildTheAction(argc,argv);
 
 }
diff --git a/tests/Test_hmc_EOWilsonFermionGauge.cc b/tests/Test_hmc_EOWilsonFermionGauge.cc
index 9f2ba82c..0a6f7aaa 100644
--- a/tests/Test_hmc_EOWilsonFermionGauge.cc
+++ b/tests/Test_hmc_EOWilsonFermionGauge.cc
@@ -1,69 +1,72 @@
 #include "Grid.h"
 
-
 using namespace std;
 using namespace Grid;
 using namespace Grid::QCD;
 
+namespace Grid { 
+  namespace QCD { 
+
+
+class HmcRunner : public NerscHmcRunner {
+public:
+
+  void BuildTheAction (int argc, char **argv)
+
+  {
+    typedef WilsonImplR ImplPolicy;
+    typedef WilsonFermionR FermionAction;
+    typedef typename FermionAction::FermionField FermionField;
+
+    UGrid   = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd,vComplex::Nsimd()),GridDefaultMpi());
+    UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid);
+  
+    FGrid   = UGrid;
+    FrbGrid = UrbGrid;
+
+    // temporarily need a gauge field
+    LatticeGaugeField  U(UGrid);
+
+    // Gauge action
+    WilsonGaugeActionR Waction(5.6);
+
+    Real mass=-0.77;
+    FermionAction FermOp(U,*FGrid,*FrbGrid,mass);
+  
+    ConjugateGradient<FermionField>  CG(1.0e-8,10000);
+
+    TwoFlavourEvenOddPseudoFermionAction<ImplPolicy> Nf2(FermOp,CG,CG);
+  
+    //Collect actions
+    ActionLevel<LatticeGaugeField> Level1(1);
+    Level1.push_back(&Nf2);
+
+    ActionLevel<LatticeGaugeField> Level2(4);
+    Level2.push_back(&Waction);
+
+    TheAction.push_back(Level1);
+    TheAction.push_back(Level2);
+
+    Run(argc,argv);
+  };
+
+};
+
+}}
+
 int main (int argc, char ** argv)
 {
   Grid_init(&argc,&argv);
 
-  std::vector<int> latt_size   = GridDefaultLatt();
-  std::vector<int> simd_layout = GridDefaultSimd(4,vComplex::Nsimd());
-  std::vector<int> mpi_layout  = GridDefaultMpi();
+  int threads = GridThread::GetThreads();
+  std::cout<<GridLogMessage << "Grid is setup to use "<<threads<<" threads"<<std::endl;
+
+  HmcRunner TheHMC;
   
-  GridCartesian            Fine(latt_size,simd_layout,mpi_layout);
-  GridRedBlackCartesian  RBFine(latt_size,simd_layout,mpi_layout);
-
-  std::vector<int> seeds({6,7,8,80});
-  GridParallelRNG  pRNG(&Fine);
-  pRNG.SeedFixedIntegers(seeds);
-
-  std::vector<int> seedsS({1,2,3,4});
-  GridSerialRNG    sRNG;
-  sRNG.SeedFixedIntegers(seedsS);
-
-  LatticeLorentzColourMatrix     U(&Fine);
-
-  SU3::HotConfiguration(pRNG, U);
-
-  // simplify template declaration? Strip the lorentz from the second template
-  WilsonGaugeActionR Waction(5.6);
-
-  Real mass=-0.77;
-  WilsonFermionR FermOp(U,Fine,RBFine,mass);
-  
-  ConjugateGradient<LatticeFermion>  CG(1.0e-8,10000);
-  
-  TwoFlavourEvenOddPseudoFermionAction<WilsonImplR> WilsonNf2(FermOp,CG,CG);
-  
-  //Collect actions
-  ActionLevel<LatticeGaugeField> Level1(1);
-  ActionLevel<LatticeGaugeField> Level2(4);
-  Level1.push_back(&WilsonNf2);
-  Level2.push_back(&Waction);
-
-  ActionSet<LatticeGaugeField> FullSet;
-  FullSet.push_back(Level1);
-  FullSet.push_back(Level2);
-
-  // Create integrator
-  //typedef LeapFrog<LatticeGaugeField>  IntegratorAlgorithm;// change here to change the algorithm
-  //typedef MinimumNorm2<LatticeGaugeField>  IntegratorAlgorithm;// change here to change the algorithm
-  typedef ForceGradient<LatticeGaugeField>  IntegratorAlgorithm;// change here to change the algorithm
-
-  IntegratorParameters MDpar(16);
-  IntegratorAlgorithm MDynamics(&Fine,MDpar, FullSet);
-
-
-  // Create HMC
-  NerscHmcCheckpointer<LatticeGaugeField> Checkpoint(std::string("ckpoint_lat"),std::string("ckpoint_rng"),1);
-  HMCparameters HMCpar;
-  HybridMonteCarlo<LatticeGaugeField,IntegratorAlgorithm>  HMC(HMCpar, MDynamics,sRNG,pRNG,U);
-  HMC.AddObservable(&Checkpoint);
-
-  // Run it
-  HMC.evolve();
+  TheHMC.BuildTheAction(argc,argv);
 
 }
+
+
+
+
diff --git a/tests/Test_hmc_EOWilsonRatio.cc b/tests/Test_hmc_EOWilsonRatio.cc
index c6efbab3..f6fb0f72 100644
--- a/tests/Test_hmc_EOWilsonRatio.cc
+++ b/tests/Test_hmc_EOWilsonRatio.cc
@@ -1,67 +1,66 @@
 #include "Grid.h"
 
-
 using namespace std;
 using namespace Grid;
 using namespace Grid::QCD;
 
+namespace Grid { 
+  namespace QCD { 
+
+
+class HmcRunner : public NerscHmcRunner {
+public:
+
+  void BuildTheAction (int argc, char **argv)
+
+  {
+    typedef WilsonImplR ImplPolicy;
+    typedef WilsonFermionR FermionAction;
+    typedef typename FermionAction::FermionField FermionField;
+
+    UGrid   = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd,vComplex::Nsimd()),GridDefaultMpi());
+    UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid);
+  
+    FGrid   = UGrid;
+    FrbGrid = UrbGrid;
+
+    // temporarily need a gauge field
+    LatticeGaugeField  U(UGrid);
+
+    // Gauge action
+    WilsonGaugeActionR Waction(5.6);
+
+    RealD mass=-0.77;
+    RealD pv  =0.0;
+    FermionAction DenOp(U,*FGrid,*FrbGrid,mass);
+    FermionAction NumOp(U,*FGrid,*FrbGrid,pv);
+  
+    ConjugateGradient<FermionField>  CG(1.0e-8,10000);
+    TwoFlavourEvenOddRatioPseudoFermionAction<ImplPolicy> Nf2(NumOp, DenOp,CG,CG);
+  
+    //Collect actions
+    ActionLevel<LatticeGaugeField> Level1;
+    Level1.push_back(&Nf2);
+    Level1.push_back(&Waction);
+    TheAction.push_back(Level1);
+
+    Run(argc,argv);
+  };
+
+};
+
+}}
+
 int main (int argc, char ** argv)
 {
   Grid_init(&argc,&argv);
 
-  std::vector<int> latt_size   = GridDefaultLatt();
-  std::vector<int> simd_layout = GridDefaultSimd(4,vComplex::Nsimd());
-  std::vector<int> mpi_layout  = GridDefaultMpi();
+  int threads = GridThread::GetThreads();
+  std::cout<<GridLogMessage << "Grid is setup to use "<<threads<<" threads"<<std::endl;
+
+  HmcRunner TheHMC;
   
-  GridCartesian            Fine(latt_size,simd_layout,mpi_layout);
-  GridRedBlackCartesian  RBFine(latt_size,simd_layout,mpi_layout);
-
-  std::vector<int> seeds({6,7,8,80});
-  GridParallelRNG  pRNG(&Fine);
-  pRNG.SeedFixedIntegers(seeds);
-
-  std::vector<int> seedsS({1,2,3,4});
-  GridSerialRNG    sRNG;
-  sRNG.SeedFixedIntegers(seedsS);
-
-  LatticeLorentzColourMatrix     U(&Fine);
-
-  SU3::HotConfiguration(pRNG, U);
-
-  // simplify template declaration? Strip the lorentz from the second template
-  WilsonGaugeActionR Waction(5.6);
-
-  Real mass=-0.77;
-  Real pv  =0.0;
-  WilsonFermionR  DenOp(U,Fine,RBFine,mass);
-  WilsonFermionR  NumOp(U,Fine,RBFine,pv);
-  
-  ConjugateGradient<LatticeFermion>  CG(1.0e-8,10000);
-  
-  TwoFlavourEvenOddRatioPseudoFermionAction<WilsonImplR> WilsonNf2(NumOp, DenOp,CG,CG);
-  
-  //Collect actions
-  ActionLevel<LatticeGaugeField> Level1;
-  Level1.push_back(&WilsonNf2);
-  Level1.push_back(&Waction);
-  ActionSet<LatticeGaugeField> FullSet;
-  FullSet.push_back(Level1);
-
-
-  // Create integrator
-  typedef MinimumNorm2<LatticeGaugeField>  IntegratorAlgorithm;// change here to change the algorithm
-  //  typedef LeapFrog  IntegratorAlgorithm;// change here to change the algorithm
-  IntegratorParameters MDpar(20);
-  IntegratorAlgorithm MDynamics(&Fine,MDpar, FullSet);
-
-
-  // Create HMC
-  NerscHmcCheckpointer<LatticeGaugeField> Checkpoint(std::string("ckpoint_lat"),std::string("ckpoint_rng"),1);
-  HMCparameters HMCpar;
-  HybridMonteCarlo<LatticeGaugeField,IntegratorAlgorithm>  HMC(HMCpar, MDynamics,sRNG,pRNG,U);
-  HMC.AddObservable(&Checkpoint);
-
-  // Run it
-  HMC.evolve();
+  TheHMC.BuildTheAction(argc,argv);
 
 }
+
diff --git a/tests/Test_hmc_WilsonFermionGauge.cc b/tests/Test_hmc_WilsonFermionGauge.cc
index bb383ad4..4fc088ba 100644
--- a/tests/Test_hmc_WilsonFermionGauge.cc
+++ b/tests/Test_hmc_WilsonFermionGauge.cc
@@ -1,68 +1,70 @@
 #include "Grid.h"
 
-
 using namespace std;
 using namespace Grid;
 using namespace Grid::QCD;
 
+namespace Grid { 
+  namespace QCD { 
+
+
+class HmcRunner : public NerscHmcRunner {
+public:
+
+  void BuildTheAction (int argc, char **argv)
+
+  {
+    typedef WilsonImplR ImplPolicy;
+    typedef WilsonFermionR FermionAction;
+    typedef typename FermionAction::FermionField FermionField;
+
+    UGrid   = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd,vComplex::Nsimd()),GridDefaultMpi());
+    UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid);
+  
+    FGrid   = UGrid;
+    FrbGrid = UrbGrid;
+
+    // temporarily need a gauge field
+    LatticeGaugeField  U(UGrid);
+
+    // Gauge action
+    WilsonGaugeActionR Waction(5.6);
+
+    Real mass=-0.77;
+    FermionAction FermOp(U,*FGrid,*FrbGrid,mass);
+  
+    ConjugateGradient<FermionField>  CG(1.0e-8,10000);
+
+    TwoFlavourPseudoFermionAction<ImplPolicy> Nf2(FermOp,CG,CG);
+  
+    //Collect actions
+    ActionLevel<LatticeGaugeField> Level1(1);
+    Level1.push_back(&Nf2);
+
+    ActionLevel<LatticeGaugeField> Level2(4);
+    Level2.push_back(&Waction);
+
+    TheAction.push_back(Level1);
+    TheAction.push_back(Level2);
+
+    Run(argc,argv);
+  };
+
+};
+
+}}
+
 int main (int argc, char ** argv)
 {
   Grid_init(&argc,&argv);
 
-  std::vector<int> latt_size   = GridDefaultLatt();
-  std::vector<int> simd_layout = GridDefaultSimd(4,vComplex::Nsimd());
-  std::vector<int> mpi_layout  = GridDefaultMpi();
+  int threads = GridThread::GetThreads();
+  std::cout<<GridLogMessage << "Grid is setup to use "<<threads<<" threads"<<std::endl;
+
+  HmcRunner TheHMC;
   
-  GridCartesian            Fine(latt_size,simd_layout,mpi_layout);
-  GridRedBlackCartesian  RBFine(latt_size,simd_layout,mpi_layout);
-
-  std::vector<int> seeds({6,7,8,80});
-  GridParallelRNG  pRNG(&Fine);
-  pRNG.SeedFixedIntegers(seeds);
-
-  std::vector<int> seedsS({1,2,3,4});
-  GridSerialRNG    sRNG;
-  sRNG.SeedFixedIntegers(seedsS);
-
-  LatticeLorentzColourMatrix     U(&Fine);
-
-  SU3::HotConfiguration(pRNG, U);
-
-  // simplify template declaration? Strip the lorentz from the second template
-  WilsonGaugeActionR Waction(5.6);
-
-  Real mass=-0.77;
-  WilsonFermionR FermOp(U,Fine,RBFine,mass);
-  
-  ConjugateGradient<LatticeFermion>  CG(1.0e-8,10000);
-  
-  TwoFlavourPseudoFermionAction<WilsonImplR> WilsonNf2(FermOp,CG,CG);
-  
-  //Collect actions
-  ActionLevel<LatticeGaugeField> Level1(1);
-  ActionLevel<LatticeGaugeField> Level2(4);
-  Level1.push_back(&WilsonNf2);
-  Level1.push_back(&Waction);
-  //  Level1.push_back(&Waction);
-
-  ActionSet<LatticeGaugeField> FullSet;
-  FullSet.push_back(Level1);
-  //  FullSet.push_back(Level2);
-
-  // Create integrator
-  typedef MinimumNorm2<LatticeGaugeField>  IntegratorAlgorithm;// change here to change the algorithm
-  //  typedef LeapFrog  IntegratorAlgorithm;// change here to change the algorithm
-  IntegratorParameters MDpar(20);
-  IntegratorAlgorithm MDynamics(&Fine,MDpar, FullSet);
-
-
-  // Create HMC
-  NerscHmcCheckpointer<LatticeGaugeField> Checkpoint(std::string("ckpoint_lat"),std::string("ckpoint_rng"),1);
-  HMCparameters HMCpar;
-  HybridMonteCarlo<LatticeGaugeField,IntegratorAlgorithm>  HMC(HMCpar, MDynamics,sRNG,pRNG,U);
-  HMC.AddObservable(&Checkpoint);
-
-  // Run it
-  HMC.evolve();
+  TheHMC.BuildTheAction(argc,argv);
 
 }
+
+
diff --git a/tests/Test_hmc_WilsonGauge.cc b/tests/Test_hmc_WilsonGauge.cc
index a89ff3da..8f112cdc 100644
--- a/tests/Test_hmc_WilsonGauge.cc
+++ b/tests/Test_hmc_WilsonGauge.cc
@@ -1,62 +1,57 @@
 #include "Grid.h"
 
-
 using namespace std;
 using namespace Grid;
 using namespace Grid::QCD;
 
+namespace Grid { 
+  namespace QCD { 
+
+
+class HmcRunner : public NerscHmcRunner {
+public:
+
+  void BuildTheAction (int argc, char **argv)
+
+  {
+    typedef WilsonImplR ImplPolicy;
+    typedef WilsonFermionR FermionAction;
+    typedef typename FermionAction::FermionField FermionField;
+
+    UGrid   = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd,vComplex::Nsimd()),GridDefaultMpi());
+    UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid);
+  
+    FGrid   = UGrid;
+    FrbGrid = UrbGrid;
+
+    // temporarily need a gauge field
+    LatticeGaugeField  U(UGrid);
+
+    // Gauge action
+    WilsonGaugeActionR Waction(5.6);
+
+    //Collect actions
+    ActionLevel<LatticeGaugeField> Level1(1);
+    Level1.push_back(&Waction);
+    TheAction.push_back(Level1);
+
+    Run(argc,argv);
+  };
+
+};
+
+}}
+
 int main (int argc, char ** argv)
 {
   Grid_init(&argc,&argv);
 
-  std::vector<int> latt_size   = GridDefaultLatt();
-  std::vector<int> simd_layout = GridDefaultSimd(4,vComplex::Nsimd());
-  std::vector<int> mpi_layout  = GridDefaultMpi();
+  int threads = GridThread::GetThreads();
+  std::cout<<GridLogMessage << "Grid is setup to use "<<threads<<" threads"<<std::endl;
 
-  latt_size.resize(4);
-
-  latt_size[0] = 8;
-  latt_size[1] = 8;
-  latt_size[2] = 8;
-  latt_size[3] = 8;
-  double volume = latt_size[0]*latt_size[1]*latt_size[2]*latt_size[3];
+  HmcRunner TheHMC;
   
-  GridCartesian           Fine(latt_size,simd_layout,mpi_layout);
-
-
-  std::vector<int> seeds({6,7,8,80});
-  GridParallelRNG  pRNG(&Fine);
-  pRNG.SeedFixedIntegers(seeds);
-
-  std::vector<int> seedsS({1,2,3,4});
-  GridSerialRNG    sRNG;
-  sRNG.SeedFixedIntegers(seedsS);
-
-  LatticeGaugeField U(&Fine);
-
-  SU3::HotConfiguration(pRNG, U);
-
-  // simplify template declaration? Strip the lorentz from the second template
-  WilsonGaugeActionR Waction(6.0);
-
-  //Collect actions
-  ActionLevel<LatticeGaugeField> Level1;
-  Level1.push_back(&Waction);
-  ActionSet<LatticeGaugeField> FullSet;
-  FullSet.push_back(Level1);
-
-  // Create integrator
-  typedef MinimumNorm2<LatticeGaugeField>  IntegratorAlgorithm;// change here to modify the algorithm
-  IntegratorParameters MDpar(20);
-  IntegratorAlgorithm  MDynamics(&Fine,MDpar, FullSet);
-
-  // Create HMC
-  NerscHmcCheckpointer<LatticeGaugeField> Checkpoint(std::string("ckpoint_lat"),std::string("ckpoint_rng"),1);
-  HMCparameters HMCpar;
-  HybridMonteCarlo<LatticeGaugeField,IntegratorAlgorithm>  HMC(HMCpar, MDynamics,sRNG,pRNG,U);
-  HMC.AddObservable(&Checkpoint);
-
-  // Run it
-  HMC.evolve();
+  TheHMC.BuildTheAction(argc,argv);
 
 }
+
diff --git a/tests/Test_hmc_WilsonRatio.cc b/tests/Test_hmc_WilsonRatio.cc
index d3792493..29128716 100644
--- a/tests/Test_hmc_WilsonRatio.cc
+++ b/tests/Test_hmc_WilsonRatio.cc
@@ -1,66 +1,67 @@
 #include "Grid.h"
 
-
 using namespace std;
 using namespace Grid;
 using namespace Grid::QCD;
 
+namespace Grid { 
+  namespace QCD { 
+
+
+class HmcRunner : public NerscHmcRunner {
+public:
+
+  void BuildTheAction (int argc, char **argv)
+
+  {
+    typedef WilsonImplR ImplPolicy;
+    typedef WilsonFermionR FermionAction;
+    typedef typename FermionAction::FermionField FermionField;
+
+    UGrid   = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd,vComplex::Nsimd()),GridDefaultMpi());
+    UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid);
+  
+    FGrid   = UGrid;
+    FrbGrid = UrbGrid;
+
+    // temporarily need a gauge field
+    LatticeGaugeField  U(UGrid);
+
+    // Gauge action
+    WilsonGaugeActionR Waction(5.6);
+
+    RealD mass=-0.77;
+    RealD pv  =0.0;
+    FermionAction DenOp(U,*FGrid,*FrbGrid,mass);
+    FermionAction NumOp(U,*FGrid,*FrbGrid,pv);
+  
+    ConjugateGradient<FermionField>  CG(1.0e-8,10000);
+    TwoFlavourRatioPseudoFermionAction<ImplPolicy> Nf2(NumOp, DenOp,CG,CG);
+  
+    //Collect actions
+    ActionLevel<LatticeGaugeField> Level1;
+    Level1.push_back(&Nf2);
+    Level1.push_back(&Waction);
+    TheAction.push_back(Level1);
+
+    Run(argc,argv);
+  };
+
+};
+
+}}
+
 int main (int argc, char ** argv)
 {
   Grid_init(&argc,&argv);
 
-  std::vector<int> latt_size   = GridDefaultLatt();
-  std::vector<int> simd_layout = GridDefaultSimd(4,vComplex::Nsimd());
-  std::vector<int> mpi_layout  = GridDefaultMpi();
+  int threads = GridThread::GetThreads();
+  std::cout<<GridLogMessage << "Grid is setup to use "<<threads<<" threads"<<std::endl;
+
+  HmcRunner TheHMC;
   
-  GridCartesian            Fine(latt_size,simd_layout,mpi_layout);
-  GridRedBlackCartesian  RBFine(latt_size,simd_layout,mpi_layout);
-
-  std::vector<int> seeds({6,7,8,80});
-  GridParallelRNG  pRNG(&Fine);
-  pRNG.SeedFixedIntegers(seeds);
-
-  std::vector<int> seedsS({1,2,3,4});
-  GridSerialRNG    sRNG;
-  sRNG.SeedFixedIntegers(seedsS);
-
-  LatticeLorentzColourMatrix     U(&Fine);
-
-  SU3::HotConfiguration(pRNG, U);
-
-  // simplify template declaration? Strip the lorentz from the second template
-  WilsonGaugeActionR Waction(5.6);
-
-  Real mass=-0.77;
-  Real pv  =0.0;
-  WilsonFermionR FermOp(U,Fine,RBFine,mass);
-  WilsonFermionR NumOp(U,Fine,RBFine,pv);
-  
-  ConjugateGradient<LatticeFermion>  CG(1.0e-8,10000);
-  
-  TwoFlavourRatioPseudoFermionAction<WilsonImplR> WilsonNf2(NumOp,FermOp,CG,CG);
-  
-  //Collect actions
-  ActionLevel<LatticeGaugeField> Level1;
-  Level1.push_back(&WilsonNf2);
-  Level1.push_back(&Waction);
-  ActionSet<LatticeGaugeField> FullSet;
-  FullSet.push_back(Level1);
-
-
-  // Create integrator
-  typedef MinimumNorm2<LatticeGaugeField>  IntegratorAlgorithm;// change here to change the algorithm
-  //  typedef LeapFrog  IntegratorAlgorithm;// change here to change the algorithm
-  IntegratorParameters MDpar(20);
-  IntegratorAlgorithm MDynamics(&Fine,MDpar, FullSet);
-
-  // Create HMC
-  NerscHmcCheckpointer<LatticeGaugeField> Checkpoint(std::string("ckpoint_lat"),std::string("ckpoint_rng"),1);
-  HMCparameters HMCpar;
-  HybridMonteCarlo<LatticeGaugeField,IntegratorAlgorithm>  HMC(HMCpar, MDynamics,sRNG,pRNG,U);
-  HMC.AddObservable(&Checkpoint);
-
-  // Run it
-  HMC.evolve();
+  TheHMC.BuildTheAction(argc,argv);
 
 }
+
+
diff --git a/tests/Test_rhmc_EOWilson1p1.cc b/tests/Test_rhmc_EOWilson1p1.cc
index 65d0e324..9ae84860 100644
--- a/tests/Test_rhmc_EOWilson1p1.cc
+++ b/tests/Test_rhmc_EOWilson1p1.cc
@@ -1,67 +1,69 @@
 #include "Grid.h"
 
-
 using namespace std;
 using namespace Grid;
 using namespace Grid::QCD;
 
+namespace Grid { 
+  namespace QCD { 
+
+
+class HmcRunner : public NerscHmcRunner {
+public:
+
+  void BuildTheAction (int argc, char **argv)
+
+  {
+    typedef WilsonImplR ImplPolicy;
+    typedef WilsonFermionR FermionAction;
+    typedef typename FermionAction::FermionField FermionField;
+
+    UGrid   = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd,vComplex::Nsimd()),GridDefaultMpi());
+    UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid);
+  
+    FGrid   = UGrid;
+    FrbGrid = UrbGrid;
+
+    // temporarily need a gauge field
+    LatticeGaugeField  U(UGrid);
+
+    // Gauge action
+    WilsonGaugeActionR Waction(5.6);
+
+    Real mass=-0.77;
+    FermionAction FermOp(U,*FGrid,*FrbGrid,mass);
+
+    // 1+1 flavour
+    OneFlavourRationalParams Params(1.0e-4,64.0,1000,1.0e-6);
+    OneFlavourEvenOddRationalPseudoFermionAction<WilsonImplR> WilsonNf1a(FermOp,Params);
+    OneFlavourEvenOddRationalPseudoFermionAction<WilsonImplR> WilsonNf1b(FermOp,Params);
+
+    //Collect actions
+    ActionLevel<LatticeGaugeField> Level1;
+    Level1.push_back(&WilsonNf1a);
+    Level1.push_back(&WilsonNf1b);
+    Level1.push_back(&Waction);
+    
+    TheAction.push_back(Level1);
+
+    Run(argc,argv);
+  };
+
+};
+
+}}
+
 int main (int argc, char ** argv)
 {
   Grid_init(&argc,&argv);
 
-  std::vector<int> latt_size   = GridDefaultLatt();
-  std::vector<int> simd_layout = GridDefaultSimd(4,vComplex::Nsimd());
-  std::vector<int> mpi_layout  = GridDefaultMpi();
+  int threads = GridThread::GetThreads();
+  std::cout<<GridLogMessage << "Grid is setup to use "<<threads<<" threads"<<std::endl;
+
+  HmcRunner TheHMC;
   
-  GridCartesian            Fine(latt_size,simd_layout,mpi_layout);
-  GridRedBlackCartesian  RBFine(latt_size,simd_layout,mpi_layout);
-
-
-  std::vector<int> seeds({6,7,8,80});
-  GridParallelRNG  pRNG(&Fine);
-  pRNG.SeedFixedIntegers(seeds);
-
-  std::vector<int> seedsS({1,2,3,4});
-  GridSerialRNG    sRNG;
-  sRNG.SeedFixedIntegers(seedsS);
-
-  LatticeLorentzColourMatrix     U(&Fine);
-
-  SU3::HotConfiguration(pRNG, U);
-
-  // simplify template declaration? Strip the lorentz from the second template
-  WilsonGaugeActionR Waction(5.6);
-
-  Real mass=-0.77;
-  WilsonFermionR FermOp(U,Fine,RBFine,mass);
-
-  // 1+1 flavour
-  OneFlavourRationalParams Params(1.0e-4,64.0,1000,1.0e-6);
-  OneFlavourEvenOddRationalPseudoFermionAction<WilsonImplR> WilsonNf1a(FermOp,Params);
-  OneFlavourEvenOddRationalPseudoFermionAction<WilsonImplR> WilsonNf1b(FermOp,Params);
-  
-  //Collect actions
-  ActionLevel<LatticeGaugeField> Level1;
-  Level1.push_back(&WilsonNf1a);
-  Level1.push_back(&WilsonNf1b);
-  Level1.push_back(&Waction);
-
-  ActionSet<LatticeGaugeField> FullSet;
-  FullSet.push_back(Level1);
-
-  // Create integrator
-  typedef MinimumNorm2<LatticeGaugeField>  IntegratorAlgorithm;// change here to change the algorithm
-
-  IntegratorParameters MDpar(20);
-  IntegratorAlgorithm MDynamics(&Fine,MDpar, FullSet);
-
-  // Create HMC
-  NerscHmcCheckpointer<LatticeGaugeField> Checkpoint(std::string("ckpoint_lat"),std::string("ckpoint_rng"),1);
-  HMCparameters HMCpar;
-  HybridMonteCarlo<LatticeGaugeField,IntegratorAlgorithm>  HMC(HMCpar, MDynamics,sRNG,pRNG,U);
-  HMC.AddObservable(&Checkpoint);
-
-  // Run it
-  HMC.evolve();
+  TheHMC.BuildTheAction(argc,argv);
 
 }
+
+
diff --git a/tests/Test_rhmc_EOWilsonRatio.cc b/tests/Test_rhmc_EOWilsonRatio.cc
index 0a24874e..d993e894 100644
--- a/tests/Test_rhmc_EOWilsonRatio.cc
+++ b/tests/Test_rhmc_EOWilsonRatio.cc
@@ -1,70 +1,72 @@
 #include "Grid.h"
 
-
 using namespace std;
 using namespace Grid;
 using namespace Grid::QCD;
 
+namespace Grid { 
+  namespace QCD { 
+
+
+class HmcRunner : public NerscHmcRunner {
+public:
+
+  void BuildTheAction (int argc, char **argv)
+
+  {
+    typedef WilsonImplR ImplPolicy;
+    typedef WilsonFermionR FermionAction;
+    typedef typename FermionAction::FermionField FermionField;
+
+    UGrid   = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd,vComplex::Nsimd()),GridDefaultMpi());
+    UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid);
+  
+    FGrid   = UGrid;
+    FrbGrid = UrbGrid;
+
+    // temporarily need a gauge field
+    LatticeGaugeField  U(UGrid);
+
+    // Gauge action
+    WilsonGaugeActionR Waction(5.6);
+
+    RealD mass=-0.77;
+    RealD pv  =0.0;
+    FermionAction DenOp(U,*FGrid,*FrbGrid,mass);
+    FermionAction NumOp(U,*FGrid,*FrbGrid,pv);
+  
+
+  // erange,maxiter,resid,npoly
+    OneFlavourRationalParams Params(1.0e-2,64.0,1000,1.0e-6,6);
+    OneFlavourEvenOddRatioRationalPseudoFermionAction<WilsonImplR> WilsonNf1a(NumOp,DenOp,Params);
+    OneFlavourEvenOddRatioRationalPseudoFermionAction<WilsonImplR> WilsonNf1b(NumOp,DenOp,Params);
+  
+    //Collect actions
+    ActionLevel<LatticeGaugeField> Level1;
+    Level1.push_back(&WilsonNf1a);
+    Level1.push_back(&WilsonNf1b);
+    Level1.push_back(&Waction);
+
+    TheAction.push_back(Level1);
+
+    Run(argc,argv);
+  };
+
+};
+
+}}
+
 int main (int argc, char ** argv)
 {
   Grid_init(&argc,&argv);
 
-  std::vector<int> latt_size   = GridDefaultLatt();
-  std::vector<int> simd_layout = GridDefaultSimd(4,vComplex::Nsimd());
-  std::vector<int> mpi_layout  = GridDefaultMpi();
+  int threads = GridThread::GetThreads();
+  std::cout<<GridLogMessage << "Grid is setup to use "<<threads<<" threads"<<std::endl;
+
+  HmcRunner TheHMC;
   
-  GridCartesian            Fine(latt_size,simd_layout,mpi_layout);
-  GridRedBlackCartesian  RBFine(latt_size,simd_layout,mpi_layout);
-
-
-  std::vector<int> seeds({6,7,8,80});
-  GridParallelRNG  pRNG(&Fine);
-  pRNG.SeedFixedIntegers(seeds);
-
-  std::vector<int> seedsS({1,2,3,4});
-  GridSerialRNG    sRNG;
-  sRNG.SeedFixedIntegers(seedsS);
-
-  LatticeLorentzColourMatrix     U(&Fine);
-
-  SU3::HotConfiguration(pRNG, U);
-
-  // simplify template declaration? Strip the lorentz from the second template
-  WilsonGaugeActionR Waction(5.6);
-
-  Real mass=-0.77;
-  Real pv  =0.0;
-  WilsonFermionR DenOp(U,Fine,RBFine,mass);
-  WilsonFermionR NumOp(U,Fine,RBFine,pv);
-  
-  // erange,maxiter,resid,npoly
-  OneFlavourRationalParams Params(1.0e-2,64.0,1000,1.0e-6,6);
-  OneFlavourEvenOddRatioRationalPseudoFermionAction<WilsonImplR> WilsonNf1a(NumOp,DenOp,Params);
-  OneFlavourEvenOddRatioRationalPseudoFermionAction<WilsonImplR> WilsonNf1b(NumOp,DenOp,Params);
-  
-  //Collect actions
-  ActionLevel<LatticeGaugeField> Level1;
-  Level1.push_back(&WilsonNf1a);
-  Level1.push_back(&WilsonNf1b);
-  Level1.push_back(&Waction);
-  ActionSet<LatticeGaugeField> FullSet;
-  FullSet.push_back(Level1);
-
-
-  // Create integrator
-  typedef MinimumNorm2<LatticeGaugeField>  IntegratorAlgorithm;// change here to change the algorithm
-  //  typedef LeapFrog  IntegratorAlgorithm;// change here to change the algorithm
-  IntegratorParameters MDpar(20);
-  IntegratorAlgorithm MDynamics(&Fine,MDpar, FullSet);
-
-  // Create HMC
-  NerscHmcCheckpointer<LatticeGaugeField> Checkpoint(std::string("ckpoint_lat"),std::string("ckpoint_rng"),1);
-  HMCparameters HMCpar;
-  HybridMonteCarlo<LatticeGaugeField,IntegratorAlgorithm>  HMC(HMCpar, MDynamics,sRNG,pRNG,U);
-  HMC.AddObservable(&Checkpoint);
-
-  // Run it
-  HMC.evolve();
-
+  TheHMC.BuildTheAction(argc,argv);
 
 }
+
+
diff --git a/tests/Test_rhmc_Wilson1p1.cc b/tests/Test_rhmc_Wilson1p1.cc
index b0b04f6e..d5ebd362 100644
--- a/tests/Test_rhmc_Wilson1p1.cc
+++ b/tests/Test_rhmc_Wilson1p1.cc
@@ -1,68 +1,68 @@
 #include "Grid.h"
 
-
 using namespace std;
 using namespace Grid;
 using namespace Grid::QCD;
 
+namespace Grid { 
+  namespace QCD { 
+
+
+class HmcRunner : public NerscHmcRunner {
+public:
+
+  void BuildTheAction (int argc, char **argv)
+
+  {
+    typedef WilsonImplR ImplPolicy;
+    typedef WilsonFermionR FermionAction;
+    typedef typename FermionAction::FermionField FermionField;
+
+    UGrid   = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd,vComplex::Nsimd()),GridDefaultMpi());
+    UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid);
+  
+    FGrid   = UGrid;
+    FrbGrid = UrbGrid;
+
+    // temporarily need a gauge field
+    LatticeGaugeField  U(UGrid);
+
+    // Gauge action
+    WilsonGaugeActionR Waction(5.6);
+
+    Real mass=-0.77;
+    FermionAction FermOp(U,*FGrid,*FrbGrid,mass);
+
+    // 1+1 flavour
+    OneFlavourRationalParams Params(1.0e-4,64.0,1000,1.0e-6);
+    OneFlavourRationalPseudoFermionAction<WilsonImplR> WilsonNf1a(FermOp,Params);
+    OneFlavourRationalPseudoFermionAction<WilsonImplR> WilsonNf1b(FermOp,Params);
+
+    //Collect actions
+    ActionLevel<LatticeGaugeField> Level1;
+    Level1.push_back(&WilsonNf1a);
+    Level1.push_back(&WilsonNf1b);
+    Level1.push_back(&Waction);
+    
+    TheAction.push_back(Level1);
+
+    Run(argc,argv);
+  };
+
+};
+
+}}
+
 int main (int argc, char ** argv)
 {
   Grid_init(&argc,&argv);
 
-  std::vector<int> latt_size   = GridDefaultLatt();
-  std::vector<int> simd_layout = GridDefaultSimd(4,vComplex::Nsimd());
-  std::vector<int> mpi_layout  = GridDefaultMpi();
+  int threads = GridThread::GetThreads();
+  std::cout<<GridLogMessage << "Grid is setup to use "<<threads<<" threads"<<std::endl;
+
+  HmcRunner TheHMC;
   
-  GridCartesian            Fine(latt_size,simd_layout,mpi_layout);
-  GridRedBlackCartesian  RBFine(latt_size,simd_layout,mpi_layout);
-
-
-  std::vector<int> seeds({6,7,8,80});
-  GridParallelRNG  pRNG(&Fine);
-  pRNG.SeedFixedIntegers(seeds);
-
-  std::vector<int> seedsS({1,2,3,4});
-  GridSerialRNG    sRNG;
-  sRNG.SeedFixedIntegers(seedsS);
-
-  LatticeLorentzColourMatrix     U(&Fine);
-
-  SU3::HotConfiguration(pRNG, U);
-
-  // simplify template declaration? Strip the lorentz from the second template
-  WilsonGaugeActionR Waction(5.6);
-
-  Real mass=-0.77;
-  WilsonFermionR FermOp(U,Fine,RBFine,mass);
-
-  // 1+1 flavour
-  OneFlavourRationalParams Params(1.0e-4,64.0,1000,1.0e-6);
-  OneFlavourRationalPseudoFermionAction<WilsonImplR> WilsonNf1a(FermOp,Params);
-  OneFlavourRationalPseudoFermionAction<WilsonImplR> WilsonNf1b(FermOp,Params);
-  
-  //Collect actions
-  ActionLevel<LatticeGaugeField> Level1;
-  Level1.push_back(&WilsonNf1a);
-  Level1.push_back(&WilsonNf1b);
-  Level1.push_back(&Waction);
-
-  ActionSet<LatticeGaugeField> FullSet;
-  FullSet.push_back(Level1);
-
-  // Create integrator
-  typedef MinimumNorm2<LatticeGaugeField>  IntegratorAlgorithm;// change here to change the algorithm
-
-  IntegratorParameters MDpar(20);
-  IntegratorAlgorithm MDynamics(&Fine,MDpar, FullSet);
-
-
-  // Create HMC
-  NerscHmcCheckpointer<LatticeGaugeField> Checkpoint(std::string("ckpoint_lat"),std::string("ckpoint_rng"),1);
-  HMCparameters HMCpar;
-  HybridMonteCarlo<LatticeGaugeField,IntegratorAlgorithm>  HMC(HMCpar, MDynamics,sRNG,pRNG,U);
-  HMC.AddObservable(&Checkpoint);
-
-  // Run it
-  HMC.evolve();
+  TheHMC.BuildTheAction(argc,argv);
 
 }
+
diff --git a/tests/Test_rhmc_WilsonRatio.cc b/tests/Test_rhmc_WilsonRatio.cc
index e38243c2..ad38d4ac 100644
--- a/tests/Test_rhmc_WilsonRatio.cc
+++ b/tests/Test_rhmc_WilsonRatio.cc
@@ -1,69 +1,71 @@
 #include "Grid.h"
 
-
 using namespace std;
 using namespace Grid;
 using namespace Grid::QCD;
 
+namespace Grid { 
+  namespace QCD { 
+
+
+class HmcRunner : public NerscHmcRunner {
+public:
+
+  void BuildTheAction (int argc, char **argv)
+
+  {
+    typedef WilsonImplR ImplPolicy;
+    typedef WilsonFermionR FermionAction;
+    typedef typename FermionAction::FermionField FermionField;
+
+    UGrid   = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd,vComplex::Nsimd()),GridDefaultMpi());
+    UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid);
+  
+    FGrid   = UGrid;
+    FrbGrid = UrbGrid;
+
+    // temporarily need a gauge field
+    LatticeGaugeField  U(UGrid);
+
+    // Gauge action
+    WilsonGaugeActionR Waction(5.6);
+
+    RealD mass=-0.77;
+    RealD pv  =0.0;
+    FermionAction DenOp(U,*FGrid,*FrbGrid,mass);
+    FermionAction NumOp(U,*FGrid,*FrbGrid,pv);
+  
+
+  // erange,maxiter,resid,npoly
+    OneFlavourRationalParams Params(1.0e-2,64.0,1000,1.0e-6,6);
+    OneFlavourRatioRationalPseudoFermionAction<WilsonImplR> WilsonNf1a(NumOp,DenOp,Params);
+    OneFlavourRatioRationalPseudoFermionAction<WilsonImplR> WilsonNf1b(NumOp,DenOp,Params);
+  
+    //Collect actions
+    ActionLevel<LatticeGaugeField> Level1;
+    Level1.push_back(&WilsonNf1a);
+    Level1.push_back(&WilsonNf1b);
+    Level1.push_back(&Waction);
+
+    TheAction.push_back(Level1);
+
+    Run(argc,argv);
+  };
+
+};
+
+}}
+
 int main (int argc, char ** argv)
 {
   Grid_init(&argc,&argv);
 
-  std::vector<int> latt_size   = GridDefaultLatt();
-  std::vector<int> simd_layout = GridDefaultSimd(4,vComplex::Nsimd());
-  std::vector<int> mpi_layout  = GridDefaultMpi();
+  int threads = GridThread::GetThreads();
+  std::cout<<GridLogMessage << "Grid is setup to use "<<threads<<" threads"<<std::endl;
+
+  HmcRunner TheHMC;
   
-  GridCartesian            Fine(latt_size,simd_layout,mpi_layout);
-  GridRedBlackCartesian  RBFine(latt_size,simd_layout,mpi_layout);
-
-
-  std::vector<int> seeds({6,7,8,80});
-  GridParallelRNG  pRNG(&Fine);
-  pRNG.SeedFixedIntegers(seeds);
-
-  std::vector<int> seedsS({1,2,3,4});
-  GridSerialRNG    sRNG;
-  sRNG.SeedFixedIntegers(seedsS);
-
-  LatticeLorentzColourMatrix     U(&Fine);
-
-  SU3::HotConfiguration(pRNG, U);
-
-  // simplify template declaration? Strip the lorentz from the second template
-  WilsonGaugeActionR Waction(5.6);
-
-  Real mass=-0.77;
-  Real pv  =0.0;
-  WilsonFermionR DenOp(U,Fine,RBFine,mass);
-  WilsonFermionR NumOp(U,Fine,RBFine,pv);
-  
-  // erange,maxiter,resid,npoly
-  OneFlavourRationalParams Params(1.0e-2,64.0,1000,1.0e-6,6);
-  OneFlavourRatioRationalPseudoFermionAction<WilsonImplR> WilsonNf1a(NumOp,DenOp,Params);
-  OneFlavourRatioRationalPseudoFermionAction<WilsonImplR> WilsonNf1b(NumOp,DenOp,Params);
-  
-  //Collect actions
-  ActionLevel<LatticeGaugeField> Level1;
-  Level1.push_back(&WilsonNf1a);
-  Level1.push_back(&WilsonNf1b);
-  Level1.push_back(&Waction);
-  ActionSet<LatticeGaugeField> FullSet;
-  FullSet.push_back(Level1);
-
-
-  // Create integrator
-  typedef MinimumNorm2<LatticeGaugeField>  IntegratorAlgorithm;// change here to change the algorithm
-  //  typedef LeapFrog  IntegratorAlgorithm;// change here to change the algorithm
-  IntegratorParameters MDpar(20);
-  IntegratorAlgorithm MDynamics(&Fine,MDpar, FullSet);
-
-  // Create HMC
-  NerscHmcCheckpointer<LatticeGaugeField> Checkpoint(std::string("ckpoint_lat"),std::string("ckpoint_rng"),1);
-  HMCparameters HMCpar;
-  HybridMonteCarlo<LatticeGaugeField,IntegratorAlgorithm>  HMC(HMCpar, MDynamics,sRNG,pRNG,U);
-  HMC.AddObservable(&Checkpoint);
-
-  // Run it
-  HMC.evolve();
+  TheHMC.BuildTheAction(argc,argv);
 
 }
+