diff --git a/lib/Grid.h b/lib/Grid.h
index c2f6e00f..47b3100c 100644
--- a/lib/Grid.h
+++ b/lib/Grid.h
@@ -89,7 +89,6 @@ Author: paboyle <paboyle@ph.ed.ac.uk>
 #include <Grid/qcd/hmc/HMCModules.h>
 #include <Grid/qcd/modules/mods.h>
 #include <Grid/qcd/hmc/HMCResourceManager.h>
-#include <Grid/qcd/hmc/HmcRunner.h>
 #include <Grid/qcd/hmc/GenericHMCrunner.h>
 #include <Grid/qcd/hmc/HMCRunnerModule.h>
 
diff --git a/lib/qcd/action/ActionBase.h b/lib/qcd/action/ActionBase.h
index dbabfbe2..8d853d45 100644
--- a/lib/qcd/action/ActionBase.h
+++ b/lib/qcd/action/ActionBase.h
@@ -4,10 +4,11 @@ Grid physics library, www.github.com/paboyle/Grid
 
 Source file: ./lib/qcd/action/ActionBase.h
 
-Copyright (C) 2015
+Copyright (C) 2015-2016
 
 Author: Peter Boyle <paboyle@ph.ed.ac.uk>
 Author: neo <cossu@post.kek.jp>
+Author: Guido Cossu <guido.cossu@ed.ac.uk>
 
 This program is free software; you can redistribute it and/or modify
 it under the terms of the GNU General Public License as published by
@@ -27,97 +28,29 @@ See the full license in the file "LICENSE" in the top level distribution
 directory
 *************************************************************************************/
 /*  END LEGAL */
-#ifndef QCD_ACTION_BASE
-#define QCD_ACTION_BASE
+
+#ifndef ACTION_BASE_H
+#define ACTION_BASE_H
+
 namespace Grid {
 namespace QCD {
 
-template <class GaugeField>
-class Action {
+template <class GaugeField >
+class Action 
+{
+
  public:
   bool is_smeared = false;
   // Heatbath?
   virtual void refresh(const GaugeField& U, GridParallelRNG& pRNG) = 0; // refresh pseudofermions
-  virtual RealD S(const GaugeField& U) = 0;                       // evaluate the action
-  virtual void deriv(const GaugeField& U, GaugeField& dSdU) = 0;  // evaluate the action derivative
-  virtual std::string action_name()    = 0;                       // return the action name
-  virtual std::string LogParameters()  = 0;                       // prints action parameters
+  virtual RealD S(const GaugeField& U) = 0;                             // evaluate the action
+  virtual void deriv(const GaugeField& U, GaugeField& dSdU) = 0;        // evaluate the action derivative
+  virtual std::string action_name()    = 0;                             // return the action name
+  virtual std::string LogParameters()  = 0;                             // prints action parameters
   virtual ~Action(){}
 };
 
-// Indexing of tuple types
-template <class T, class Tuple>
-struct Index;
-
-template <class T, class... Types>
-struct Index<T, std::tuple<T, Types...>> {
-  static const std::size_t value = 0;
-};
-
-template <class T, class U, class... Types>
-struct Index<T, std::tuple<U, Types...>> {
-  static const std::size_t value = 1 + Index<T, std::tuple<Types...>>::value;
-};
-
-template <class Field, class Repr = NoHirep >
-struct ActionLevel {
- public:
-  unsigned int multiplier;
-
-  // Fundamental repr actions separated because of the smearing
-  typedef Action<Field>* ActPtr;
-
-  // construct a tuple of vectors of the actions for the corresponding higher
-  // representation fields
-  typedef typename AccessTypes<Action, Repr>::VectorCollection action_collection;
-  action_collection actions_hirep;
-  typedef typename  AccessTypes<Action, Repr>::FieldTypeCollection action_hirep_types;
-
-  std::vector<ActPtr>& actions;
-
-  explicit ActionLevel(unsigned int mul = 1) : actions(std::get<0>(actions_hirep)), multiplier(mul) {
-    // initialize the hirep vectors to zero.
-    // apply(this->resize, actions_hirep, 0); //need a working resize
-    assert(mul >= 1);
-  }
-
-  // void push_back(ActPtr ptr) { actions.push_back(ptr); }
-
-  template < class GenField >
-  void push_back(Action<GenField>* ptr) {
-    // insert only in the correct vector
-    std::get< Index < GenField, action_hirep_types>::value >(actions_hirep).push_back(ptr);
-  };
-
- 
-
-  template < class ActPtr>
-  static void resize(ActPtr ap, unsigned int n){
-    ap->resize(n);
-
-  }
-
-  //template <std::size_t I>
-  //auto getRepresentation(Repr& R)->decltype(std::get<I>(R).U)  {return std::get<I>(R).U;}
-
-  // Loop on tuple for a callable function
-  template <std::size_t I = 1, typename Callable, typename ...Args>
-  inline typename std::enable_if<I == std::tuple_size<action_collection>::value, void>::type apply(
-      Callable, Repr& R,Args&...) const {}
-
-  template <std::size_t I = 1, typename Callable, typename ...Args>
-  inline typename std::enable_if<I < std::tuple_size<action_collection>::value, void>::type apply(
-      Callable fn, Repr& R, Args&... arguments) const {
-    fn(std::get<I>(actions_hirep), std::get<I>(R.rep), arguments...);
-    apply<I + 1>(fn, R, arguments...);
-  }  
-
-};
-
-
-template <class GaugeField, class R>
-using ActionSet = std::vector<ActionLevel<GaugeField, R> >;
-
 }
 }
-#endif
+
+#endif // ACTION_BASE_H
diff --git a/lib/qcd/action/ActionParams.h b/lib/qcd/action/ActionParams.h
index 48ed23cf..b972df38 100644
--- a/lib/qcd/action/ActionParams.h
+++ b/lib/qcd/action/ActionParams.h
@@ -1,68 +1,80 @@
-    /*************************************************************************************
+/*************************************************************************************
 
-    Grid physics library, www.github.com/paboyle/Grid 
+Grid physics library, www.github.com/paboyle/Grid
 
-    Source file: ./lib/qcd/action/ActionParams.h
+Source file: ./lib/qcd/action/ActionParams.h
 
-    Copyright (C) 2015
+Copyright (C) 2015
 
 Author: Peter Boyle <paboyle@ph.ed.ac.uk>
 Author: paboyle <paboyle@ph.ed.ac.uk>
+Author: Guido Cossu <guido.cossu@ed.ac.uk>
 
-    This program is free software; you can redistribute it and/or modify
-    it under the terms of the GNU General Public License as published by
-    the Free Software Foundation; either version 2 of the License, or
-    (at your option) any later version.
+This program is free software; you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation; either version 2 of the License, or
+(at your option) any later version.
 
-    This program is distributed in the hope that it will be useful,
-    but WITHOUT ANY WARRANTY; without even the implied warranty of
-    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-    GNU General Public License for more details.
+This program is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
 
-    You should have received a copy of the GNU General Public License along
-    with this program; if not, write to the Free Software Foundation, Inc.,
-    51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+You should have received a copy of the GNU General Public License along
+with this program; if not, write to the Free Software Foundation, Inc.,
+51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+
+See the full license in the file "LICENSE" in the top level distribution
+directory
+*************************************************************************************/
+/*  END LEGAL */
 
-    See the full license in the file "LICENSE" in the top level distribution directory
-    *************************************************************************************/
-    /*  END LEGAL */
 #ifndef GRID_QCD_ACTION_PARAMS_H
 #define GRID_QCD_ACTION_PARAMS_H
 
 namespace Grid {
 namespace QCD {
 
-    // These can move into a params header and be given MacroMagic serialisation
-    struct GparityWilsonImplParams {
-      bool overlapCommsCompute;
-      std::vector<int> twists; 
-      GparityWilsonImplParams () : twists(Nd,0), overlapCommsCompute(false) {};
+// These can move into a params header and be given MacroMagic serialisation
+struct GparityWilsonImplParams {
+  bool overlapCommsCompute;
+  std::vector<int> twists;
+  GparityWilsonImplParams() : twists(Nd, 0), overlapCommsCompute(false){};
+};
 
-    };
+struct WilsonImplParams {
+  bool overlapCommsCompute;
+  WilsonImplParams() : overlapCommsCompute(false){};
+};
 
-    struct WilsonImplParams {
-      bool overlapCommsCompute;
-      WilsonImplParams() : overlapCommsCompute(false) {};
-    };
-
-  struct OneFlavourRationalParams : Serializable {
-      GRID_SERIALIZABLE_CLASS_MEMBERS(OneFlavourRationalParams,
-				      RealD,  lo,
-				      RealD,  hi,
-				      int, MaxIter,   
-				      RealD, tolerance, 
-				      int,    degree,
-				      int, precision);
-
-      // MaxIter and tolerance, vectors??
-  public:
-      OneFlavourRationalParams (RealD _lo   = 0.0,  RealD _hi      = 0.0 ,
-				int _maxit  = 1000, RealD tol      = 1.0e-8,
-				int _degree = 10,   int _precision = 64) :
-        lo(_lo), hi(_hi), MaxIter(_maxit), tolerance(tol), degree(_degree), precision(_precision)
-      {};
-    };
-
-}}
+struct OneFlavourRationalParams : Serializable {
+  GRID_SERIALIZABLE_CLASS_MEMBERS(OneFlavourRationalParams, 
+				  RealD, lo, 
+				  RealD, hi, 
+				  int,   MaxIter, 
+				  RealD, tolerance, 
+				  int,   degree, 
+				  int,   precision);
+  
+  // MaxIter and tolerance, vectors??
+  
+  // constructor 
+  OneFlavourRationalParams(	RealD _lo      = 0.0, 
+				RealD _hi      = 1.0, 
+				int _maxit     = 1000,
+				RealD tol      = 1.0e-8, 
+                           	int _degree    = 10,
+				int _precision = 64)
+    : lo(_lo),
+      hi(_hi),
+      MaxIter(_maxit),
+      tolerance(tol),
+      degree(_degree),
+      precision(_precision){};
+};
+  
+  
+}
+}
 
 #endif
diff --git a/lib/qcd/action/ActionSet.h b/lib/qcd/action/ActionSet.h
new file mode 100644
index 00000000..4ed6a582
--- /dev/null
+++ b/lib/qcd/action/ActionSet.h
@@ -0,0 +1,116 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid
+
+Source file: ./lib/qcd/action/ActionSet.h
+
+Copyright (C) 2015
+
+Author: Peter Boyle <paboyle@ph.ed.ac.uk>
+Author: neo <cossu@post.kek.jp>
+
+This program is free software; you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation; either version 2 of the License, or
+(at your option) any later version.
+
+This program is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License along
+with this program; if not, write to the Free Software Foundation, Inc.,
+51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+
+See the full license in the file "LICENSE" in the top level distribution
+directory
+*************************************************************************************/
+/*  END LEGAL */
+#ifndef ACTION_SET_H
+#define ACTION_SET_H
+
+namespace Grid {
+
+// Should drop this namespace here
+namespace QCD {
+
+//////////////////////////////////
+// Indexing of tuple types
+//////////////////////////////////
+
+template <class T, class Tuple>
+struct Index;
+
+template <class T, class... Types>
+struct Index<T, std::tuple<T, Types...>> {
+  static const std::size_t value = 0;
+};
+
+template <class T, class U, class... Types>
+struct Index<T, std::tuple<U, Types...>> {
+  static const std::size_t value = 1 + Index<T, std::tuple<Types...>>::value;
+};
+
+
+////////////////////////////////////////////
+// Action Level
+// Action collection 
+// in a integration level
+// (for multilevel integration schemes)
+////////////////////////////////////////////
+
+template <class Field, class Repr = NoHirep >
+struct ActionLevel {
+ public:
+  unsigned int multiplier;
+
+  // Fundamental repr actions separated because of the smearing
+  typedef Action<Field>* ActPtr;
+
+  // construct a tuple of vectors of the actions for the corresponding higher
+  // representation fields
+  typedef typename AccessTypes<Action, Repr>::VectorCollection action_collection;
+  typedef typename  AccessTypes<Action, Repr>::FieldTypeCollection action_hirep_types;
+
+  action_collection actions_hirep;
+  std::vector<ActPtr>& actions;
+
+  explicit ActionLevel(unsigned int mul = 1) : 
+  actions(std::get<0>(actions_hirep)), multiplier(mul) {
+    // initialize the hirep vectors to zero.
+    // apply(this->resize, actions_hirep, 0); //need a working resize
+    assert(mul >= 1);
+  }
+
+  template < class GenField >
+  void push_back(Action<GenField>* ptr) {
+    // insert only in the correct vector
+    std::get< Index < GenField, action_hirep_types>::value >(actions_hirep).push_back(ptr);
+  };
+
+  template <class ActPtr>
+  static void resize(ActPtr ap, unsigned int n) {
+    ap->resize(n);
+  }
+
+  // Loop on tuple for a callable function
+  template <std::size_t I = 1, typename Callable, typename ...Args>
+  inline typename std::enable_if<I == std::tuple_size<action_collection>::value, void>::type apply(Callable, Repr& R,Args&...) const {}
+
+  template <std::size_t I = 1, typename Callable, typename ...Args>
+  inline typename std::enable_if<I < std::tuple_size<action_collection>::value, void>::type apply(Callable fn, Repr& R, Args&... arguments) const {
+    fn(std::get<I>(actions_hirep), std::get<I>(R.rep), arguments...);
+    apply<I + 1>(fn, R, arguments...);
+  }  
+
+};
+
+// Define the ActionSet
+template <class GaugeField, class R>
+using ActionSet = std::vector<ActionLevel<GaugeField, R> >;
+
+} // QCD
+} // Grid
+
+#endif  // ACTION_SET_H
diff --git a/lib/qcd/action/Actions.h b/lib/qcd/action/Actions.h
index f073427a..5f83aa55 100644
--- a/lib/qcd/action/Actions.h
+++ b/lib/qcd/action/Actions.h
@@ -1,4 +1,4 @@
-    /*************************************************************************************
+/*************************************************************************************
 
     Grid physics library, www.github.com/paboyle/Grid 
 
@@ -29,7 +29,8 @@ Author: paboyle <paboyle@ph.ed.ac.uk>
 
     See the full license in the file "LICENSE" in the top level distribution directory
     *************************************************************************************/
-    /*  END LEGAL */
+/*  END LEGAL */
+
 #ifndef GRID_QCD_ACTIONS_H
 #define GRID_QCD_ACTIONS_H
 
@@ -41,12 +42,14 @@ Author: paboyle <paboyle@ph.ed.ac.uk>
 // Abstract base interface
 ////////////////////////////////////////////
 #include <Grid/qcd/action/ActionBase.h>
+#include <Grid/qcd/action/ActionSet.h>
 #include <Grid/qcd/action/ActionParams.h>
 
 ////////////////////////////////////////////
 // Utility functions
 ////////////////////////////////////////////
-#include <Grid/qcd/action/gauge/GaugeImpl.h>
+//#include <Grid/qcd/action/gauge/GaugeImplTypes.h>
+#include <Grid/qcd/action/gauge/GaugeImplementations.h>
 #include <Grid/qcd/utils/WilsonLoops.h>
 
 #include <Grid/qcd/action/fermion/WilsonCompressor.h>     //used by all wilson type fermions
@@ -63,7 +66,7 @@ Author: paboyle <paboyle@ph.ed.ac.uk>
 ////////////////////////////////////////////
 // Scalar Actions
 ////////////////////////////////////////////
-#include <Grid/qcd/action/scalar/scalarImpl.h>
+#include <Grid/qcd/action/scalar/ScalarImpl.h>
 #include <Grid/qcd/action/scalar/ScalarAction.h>
 
 namespace Grid {
diff --git a/lib/qcd/action/gauge/GaugeImplTypes.h b/lib/qcd/action/gauge/GaugeImplTypes.h
new file mode 100644
index 00000000..4cfb4780
--- /dev/null
+++ b/lib/qcd/action/gauge/GaugeImplTypes.h
@@ -0,0 +1,144 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid
+
+Source file: ./lib/qcd/action/gauge/GaugeImpl.h
+
+Copyright (C) 2015
+
+Author: paboyle <paboyle@ph.ed.ac.uk>
+
+This program is free software; you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation; either version 2 of the License, or
+(at your option) any later version.
+
+This program is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License along
+with this program; if not, write to the Free Software Foundation, Inc.,
+51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+
+See the full license in the file "LICENSE" in the top level distribution
+directory
+*************************************************************************************/
+/*  END LEGAL */
+#ifndef GRID_GAUGE_IMPL_TYPES_H
+#define GRID_GAUGE_IMPL_TYPES_H
+
+namespace Grid {
+namespace QCD {
+
+////////////////////////////////////////////////////////////////////////
+// Implementation dependent gauge types
+////////////////////////////////////////////////////////////////////////
+
+#define INHERIT_GIMPL_TYPES(GImpl)                  \
+  typedef typename GImpl::Simd Simd;                \
+  typedef typename GImpl::LinkField GaugeLinkField; \
+  typedef typename GImpl::Field GaugeField;         \
+  typedef typename GImpl::SiteField SiteGaugeField; \
+  typedef typename GImpl::SiteLink SiteGaugeLink;
+
+#define INHERIT_FIELD_TYPES(Impl)             \
+  typedef typename Impl::Simd Simd;           \
+  typedef typename Impl::SiteField SiteField; \
+  typedef typename Impl::Field Field;
+
+// hardcodes the exponential approximation in the template
+template <class S, int Nrepresentation = Nc, int Nexp = 12 > class GaugeImplTypes {
+public:
+  typedef S Simd;
+
+  template <typename vtype> using iImplGaugeLink  = iScalar<iScalar<iMatrix<vtype, Nrepresentation>>>;
+  template <typename vtype> using iImplGaugeField = iVector<iScalar<iMatrix<vtype, Nrepresentation>>, Nd>;
+
+  typedef iImplGaugeLink<Simd>  SiteLink;
+  typedef iImplGaugeField<Simd> SiteField;
+
+  typedef Lattice<SiteLink>  LinkField; 
+  typedef Lattice<SiteField> Field;
+
+  // Guido: we can probably separate the types from the HMC functions
+  // this will create 2 kind of implementations
+  // probably confusing the users
+  // Now keeping only one class
+
+
+  // Move this elsewhere? FIXME
+  static inline void AddLink(Field &U, LinkField &W,
+                                  int mu) { // U[mu] += W
+    PARALLEL_FOR_LOOP
+    for (auto ss = 0; ss < U._grid->oSites(); ss++) {
+      U._odata[ss]._internal[mu] =
+          U._odata[ss]._internal[mu] + W._odata[ss]._internal;
+    }
+  }
+
+  ///////////////////////////////////////////////////////////
+  // Move these to another class
+  // HMC auxiliary functions
+  static inline void generate_momenta(Field& P, GridParallelRNG& pRNG){
+    // specific for SU gauge fields
+    LinkField Pmu(P._grid);
+    Pmu = zero;
+    for (int mu = 0; mu < Nd; mu++) {
+      SU<Nrepresentation>::GaussianFundamentalLieAlgebraMatrix(pRNG, Pmu);
+      PokeIndex<LorentzIndex>(P, Pmu, mu);
+    }
+  }
+
+   static inline Field projectForce(Field& P){
+    return Ta(P);
+   }
+  
+  static inline void update_field(Field& P, Field& U, double ep){
+    for (int mu = 0; mu < Nd; mu++) {
+      auto Umu = PeekIndex<LorentzIndex>(U, mu);
+      auto Pmu = PeekIndex<LorentzIndex>(P, mu);
+      Umu = expMat(Pmu, ep, Nexp) * Umu;
+      PokeIndex<LorentzIndex>(U, ProjectOnGroup(Umu), mu);
+    }
+  }
+
+  static inline RealD FieldSquareNorm(Field& U){
+    LatticeComplex Hloc(U._grid);
+    Hloc = zero;
+    for (int mu = 0; mu < Nd; mu++) {
+      auto Umu = PeekIndex<LorentzIndex>(U, mu);
+      Hloc += trace(Umu * Umu);
+    }
+    Complex Hsum = sum(Hloc);
+    return Hsum.real();
+  }
+
+  static inline void HotConfiguration(GridParallelRNG &pRNG, Field &U) {
+    SU<Nc>::HotConfiguration(pRNG, U);
+  }
+
+  static inline void TepidConfiguration(GridParallelRNG &pRNG, Field &U) {
+    SU<Nc>::TepidConfiguration(pRNG, U);
+  }
+
+  static inline void ColdConfiguration(GridParallelRNG &pRNG, Field &U) {
+    SU<Nc>::ColdConfiguration(pRNG, U);
+  }
+};
+
+
+typedef GaugeImplTypes<vComplex, Nc> GimplTypesR;
+typedef GaugeImplTypes<vComplexF, Nc> GimplTypesF;
+typedef GaugeImplTypes<vComplexD, Nc> GimplTypesD;
+
+typedef GaugeImplTypes<vComplex, SU<Nc>::AdjointDimension> GimplAdjointTypesR;
+typedef GaugeImplTypes<vComplexF, SU<Nc>::AdjointDimension> GimplAdjointTypesF;
+typedef GaugeImplTypes<vComplexD, SU<Nc>::AdjointDimension> GimplAdjointTypesD;
+
+
+} // QCD
+} // Grid
+
+#endif // GRID_GAUGE_IMPL_TYPES_H
diff --git a/lib/qcd/action/gauge/GaugeImpl.h b/lib/qcd/action/gauge/GaugeImplementations.h
similarity index 60%
rename from lib/qcd/action/gauge/GaugeImpl.h
rename to lib/qcd/action/gauge/GaugeImplementations.h
index 08fbfa68..e06a813d 100644
--- a/lib/qcd/action/gauge/GaugeImpl.h
+++ b/lib/qcd/action/gauge/GaugeImplementations.h
@@ -26,107 +26,14 @@ See the full license in the file "LICENSE" in the top level distribution
 directory
 *************************************************************************************/
 /*  END LEGAL */
-#ifndef GRID_QCD_GAUGE_IMPL_H
-#define GRID_QCD_GAUGE_IMPL_H
+#ifndef GRID_QCD_GAUGE_IMPLEMENTATIONS_H
+#define GRID_QCD_GAUGE_IMPLEMENTATIONS_H
+
+#include "GaugeImplTypes.h"
 
 namespace Grid {
 namespace QCD {
 
-////////////////////////////////////////////////////////////////////////
-// Implementation dependent gauge types
-////////////////////////////////////////////////////////////////////////
-
-template <class Gimpl> class WilsonLoops;
-
-// 
-#define INHERIT_GIMPL_TYPES(GImpl)             \
-  typedef typename GImpl::Simd Simd;           \
-  typedef typename GImpl::LinkField GaugeLinkField; \
-  typedef typename GImpl::Field GaugeField;         \
-  typedef typename GImpl::SiteField SiteGaugeField; \
-  typedef typename GImpl::SiteLink SiteGaugeLink;
-
-#define INHERIT_FIELD_TYPES(Impl) \
-  typedef typename Impl::Simd Simd;           \
-  typedef typename Impl::SiteField SiteField; \
-  typedef typename Impl::Field Field; 
-
-// hard codes the exponential approximation in the template
-template <class S, int Nrepresentation = Nc, int Nexp = 12 > class GaugeImplTypes {
-public:
-  typedef S Simd;
-
-  template <typename vtype>
-  using iImplGaugeLink  = iScalar<iScalar<iMatrix<vtype, Nrepresentation>>>;
-  template <typename vtype>
-  using iImplGaugeField = iVector<iScalar<iMatrix<vtype, Nrepresentation>>, Nd>;
-
-  typedef iImplGaugeLink<Simd> SiteLink;
-  typedef iImplGaugeField<Simd> SiteField;
-
-  typedef Lattice<SiteLink>  LinkField; 
-  typedef Lattice<SiteField> Field;
-
-  // Move this elsewhere? FIXME
-  static inline void AddLink(Field &U, LinkField &W,
-                                  int mu) { // U[mu] += W
-    PARALLEL_FOR_LOOP
-    for (auto ss = 0; ss < U._grid->oSites(); ss++) {
-      U._odata[ss]._internal[mu] =
-          U._odata[ss]._internal[mu] + W._odata[ss]._internal;
-    }
-  }
-
-  ///////////////////////////////////////////////////////////
-  // Move these to another class
-  // HMC auxiliary functions
-  static inline void generate_momenta(Field& P, GridParallelRNG& pRNG){
-    // specific for SU gauge fields
-    LinkField Pmu(P._grid);
-    Pmu = zero;
-    for (int mu = 0; mu < Nd; mu++) {
-      SU<Nrepresentation>::GaussianFundamentalLieAlgebraMatrix(pRNG, Pmu);
-      PokeIndex<LorentzIndex>(P, Pmu, mu);
-    }
-  }
-
-   static inline Field projectForce(Field& P){
-    return Ta(P);
-   }
-  
-  static inline void update_field(Field& P, Field& U, double ep){
-    for (int mu = 0; mu < Nd; mu++) {
-      auto Umu = PeekIndex<LorentzIndex>(U, mu);
-      auto Pmu = PeekIndex<LorentzIndex>(P, mu);
-      Umu = expMat(Pmu, ep, Nexp) * Umu;
-      PokeIndex<LorentzIndex>(U, ProjectOnGroup(Umu), mu);
-    }
-  }
-
-  static inline RealD FieldSquareNorm(Field& U){
-    LatticeComplex Hloc(U._grid);
-    Hloc = zero;
-    for (int mu = 0; mu < Nd; mu++) {
-      auto Umu = PeekIndex<LorentzIndex>(U, mu);
-      Hloc += trace(Umu * Umu);
-    }
-    Complex Hsum = sum(Hloc);
-    return Hsum.real();
-  }
-
-  static inline void HotConfiguration(GridParallelRNG &pRNG, Field &U) {
-    SU<Nc>::HotConfiguration(pRNG, U);
-  }
-
-  static inline void TepidConfiguration(GridParallelRNG &pRNG, Field &U) {
-    SU<Nc>::TepidConfiguration(pRNG, U);
-  }
-
-  static inline void ColdConfiguration(GridParallelRNG &pRNG, Field &U) {
-    SU<Nc>::ColdConfiguration(pRNG, U);
-  }
-};
-
 // Composition with smeared link, bc's etc.. probably need multiple inheritance
 // Variable precision "S" and variable Nc
 template <class GimplTypes> class PeriodicGaugeImpl : public GimplTypes {
@@ -222,14 +129,6 @@ public:
   static inline bool isPeriodicGaugeField(void) { return false; }
 };
 
-typedef GaugeImplTypes<vComplex, Nc> GimplTypesR;
-typedef GaugeImplTypes<vComplexF, Nc> GimplTypesF;
-typedef GaugeImplTypes<vComplexD, Nc> GimplTypesD;
-
-typedef GaugeImplTypes<vComplex, SU<Nc>::AdjointDimension> GimplAdjointTypesR;
-typedef GaugeImplTypes<vComplexF, SU<Nc>::AdjointDimension> GimplAdjointTypesF;
-typedef GaugeImplTypes<vComplexD, SU<Nc>::AdjointDimension> GimplAdjointTypesD;
-
 typedef PeriodicGaugeImpl<GimplTypesR> PeriodicGimplR; // Real.. whichever prec
 typedef PeriodicGaugeImpl<GimplTypesF> PeriodicGimplF; // Float
 typedef PeriodicGaugeImpl<GimplTypesD> PeriodicGimplD; // Double
@@ -241,6 +140,8 @@ typedef PeriodicGaugeImpl<GimplAdjointTypesD> PeriodicGimplAdjD; // Double
 typedef ConjugateGaugeImpl<GimplTypesR> ConjugateGimplR; // Real.. whichever prec
 typedef ConjugateGaugeImpl<GimplTypesF> ConjugateGimplF; // Float
 typedef ConjugateGaugeImpl<GimplTypesD> ConjugateGimplD; // Double
+
+
 }
 }
 
diff --git a/lib/qcd/action/scalar/ScalarAction.h b/lib/qcd/action/scalar/ScalarAction.h
index b3e14109..a584cad6 100644
--- a/lib/qcd/action/scalar/ScalarAction.h
+++ b/lib/qcd/action/scalar/ScalarAction.h
@@ -29,51 +29,48 @@ Author: paboyle <paboyle@ph.ed.ac.uk>
 directory
   *************************************************************************************/
 /*  END LEGAL */
+
 #ifndef SCALAR_ACTION_H
 #define SCALAR_ACTION_H
 
 namespace Grid {
-namespace QCD {
-
-////////////////////////////////////////////////////////////////////////
-// Wilson Gauge Action .. should I template the Nc etc..
-////////////////////////////////////////////////////////////////////////
-
-template <class Impl>
-class ScalarAction : public Action<typename Impl::Field> {
- public:
-  INHERIT_FIELD_TYPES(Impl);
-
- private:
-  RealD mass_square;
-  RealD lambda;
-
- public:
-  ScalarAction(RealD ms, RealD l) : mass_square(ms), lambda(l){};
-
-  virtual std::string action_name(){return "ScalarAction";}
+  // FIXME drop the QCD namespace everywhere here
   
-  virtual void refresh(const Field &U,
-                       GridParallelRNG &pRNG){};  // noop as no pseudoferms
-
-  virtual RealD S(const Field &p) {
-    return (mass_square * 0.5 + Nd) * ScalarObs<Impl>::sumphisquared(p) +
-           (lambda / 24.) * ScalarObs<Impl>::sumphifourth(p) +
-           ScalarObs<Impl>::sumphider(p);
+  template <class Impl>
+  class ScalarAction : public QCD::Action<typename Impl::Field> {
+  public:
+    INHERIT_FIELD_TYPES(Impl);
+    
+  private:
+    RealD mass_square;
+    RealD lambda;
+    
+  public:
+    ScalarAction(RealD ms, RealD l) : mass_square(ms), lambda(l){};
+    
+    virtual std::string action_name(){return "ScalarAction";}
+    
+    virtual void refresh(const Field &U,
+			 GridParallelRNG &pRNG){};  // noop as no pseudoferms
+    
+    virtual RealD S(const Field &p) {
+      return (mass_square * 0.5 + QCD::Nd) * ScalarObs<Impl>::sumphisquared(p) +
+	(lambda / 24.) * ScalarObs<Impl>::sumphifourth(p) +
+	ScalarObs<Impl>::sumphider(p);
+    };
+    
+    virtual void deriv(const Field &p,
+		       Field &force) {
+      Field tmp(p._grid);
+      Field p2(p._grid);
+      ScalarObs<Impl>::phisquared(p2, p);
+      tmp = -(Cshift(p, 0, -1) + Cshift(p, 0, 1));
+      for (int mu = 1; mu < QCD::Nd; mu++) tmp -= Cshift(p, mu, -1) + Cshift(p, mu, 1);
+      
+      force=+(mass_square + 2. * QCD::Nd) * p + (lambda / 6.) * p2 * p + tmp;
+    };
   };
+  
+} // Grid
 
-  virtual void deriv(const Field &p,
-                     Field &force) {
-    Field tmp(p._grid);
-    Field p2(p._grid);
-    ScalarObs<Impl>::phisquared(p2, p);
-    tmp = -(Cshift(p, 0, -1) + Cshift(p, 0, 1));
-    for (int mu = 1; mu < Nd; mu++) tmp -= Cshift(p, mu, -1) + Cshift(p, mu, 1);
-
-    force=+(mass_square + 2. * Nd) * p + (lambda / 6.) * p2 * p + tmp;
-  };
-};
-}
-}
-
-#endif
+#endif // SCALAR_ACTION_H
diff --git a/lib/qcd/action/scalar/ScalarImpl.h b/lib/qcd/action/scalar/ScalarImpl.h
new file mode 100644
index 00000000..a259c34b
--- /dev/null
+++ b/lib/qcd/action/scalar/ScalarImpl.h
@@ -0,0 +1,57 @@
+#ifndef SCALAR_IMPL
+#define SCALAR_IMPL
+
+
+namespace Grid {
+  //namespace QCD {
+
+  template <class S>
+  class ScalarImplTypes {
+  public:
+    typedef S Simd;
+    
+    template <typename vtype>
+    using iImplField = iScalar<iScalar<iScalar<vtype> > >;
+    
+    typedef iImplField<Simd> SiteField;
+    
+    
+    typedef Lattice<SiteField> Field;
+    
+    static inline void generate_momenta(Field& P, GridParallelRNG& pRNG){
+      gaussian(pRNG, P);
+    }
+    
+    static inline Field projectForce(Field& P){return P;}
+    
+    static inline void update_field(Field& P, Field& U, double ep){
+      U += P*ep;
+    }
+    
+    static inline RealD FieldSquareNorm(Field& U){
+      return (- sum(trace(U*U))/2.0);
+    }
+    
+    static inline void HotConfiguration(GridParallelRNG &pRNG, Field &U) {
+      gaussian(pRNG, U);
+    }
+    
+    static inline void TepidConfiguration(GridParallelRNG &pRNG, Field &U) {
+      gaussian(pRNG, U);
+    }
+    
+    static inline void ColdConfiguration(GridParallelRNG &pRNG, Field &U) {
+      U = 1.0;
+    }
+    
+  };
+  
+  
+  typedef ScalarImplTypes<vReal> ScalarImplR;
+  typedef ScalarImplTypes<vRealF> ScalarImplF;
+  typedef ScalarImplTypes<vRealD> ScalarImplD;
+  
+  //} 
+} 
+
+#endif
diff --git a/lib/qcd/action/scalar/scalarImpl.h b/lib/qcd/action/scalar/scalarImpl.h
deleted file mode 100644
index f3a7932d..00000000
--- a/lib/qcd/action/scalar/scalarImpl.h
+++ /dev/null
@@ -1,56 +0,0 @@
-#ifndef SCALAR_IMPL
-#define SCALAR_IMPL
-
-
-namespace Grid {
-namespace QCD {
-
-template <class S>
-class ScalarImplTypes {
- public:
-  typedef S Simd;
-
-  template <typename vtype>
-  using iImplField = iScalar<iScalar<iScalar<vtype> > >;
-
-  typedef iImplField<Simd> SiteField;
-
-  
-  typedef Lattice<SiteField> Field;
-
-  static inline void generate_momenta(Field& P, GridParallelRNG& pRNG){
-  	gaussian(pRNG, P);
-  }
-
-  static inline Field projectForce(Field& P){return P;}
-
-  static inline void update_field(Field& P, Field& U, double ep){
-  	U += P*ep;
-  }
-
-  static inline RealD FieldSquareNorm(Field& U){
-  	return (- sum(trace(U*U))/2.0);
-  }
-
-  static inline void HotConfiguration(GridParallelRNG &pRNG, Field &U) {
-    gaussian(pRNG, U);
-  }
-
-  static inline void TepidConfiguration(GridParallelRNG &pRNG, Field &U) {
-    gaussian(pRNG, U);
-  }
-
-  static inline void ColdConfiguration(GridParallelRNG &pRNG, Field &U) {
-    U = 1.0;
-  }
-
-};
-
-
-typedef ScalarImplTypes<vReal> ScalarImplR;
-typedef ScalarImplTypes<vRealF> ScalarImplF;
-typedef ScalarImplTypes<vRealD> ScalarImplD;
-
-}}
-
-#endif
\ No newline at end of file
diff --git a/lib/qcd/hmc/HmcRunner.h b/lib/qcd/hmc/HmcRunner.h
deleted file mode 100644
index c2fb5078..00000000
--- a/lib/qcd/hmc/HmcRunner.h
+++ /dev/null
@@ -1,191 +0,0 @@
-/*************************************************************************************
-
-Grid physics library, www.github.com/paboyle/Grid
-
-Source file: ./lib/qcd/hmc/HmcRunner.h
-
-Copyright (C) 2015
-
-Author: paboyle <paboyle@ph.ed.ac.uk>
-
-This program is free software; you can redistribute it and/or modify
-it under the terms of the GNU General Public License as published by
-the Free Software Foundation; either version 2 of the License, or
-(at your option) any later version.
-
-This program is distributed in the hope that it will be useful,
-but WITHOUT ANY WARRANTY; without even the implied warranty of
-MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-GNU General Public License for more details.
-
-You should have received a copy of the GNU General Public License along
-with this program; if not, write to the Free Software Foundation, Inc.,
-51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
-
-See the full license in the file "LICENSE" in the top level distribution
-directory
-*************************************************************************************/
-/*  END LEGAL */
-#ifndef HMC_RUNNER
-#define HMC_RUNNER
-
-namespace Grid {
-namespace QCD {
-
-// Class for HMC specific for gauge theories
-template <class Gimpl, class RepresentationsPolicy = NoHirep >
-class NerscHmcRunnerTemplate {
- public:
-  INHERIT_GIMPL_TYPES(Gimpl);
-
-  enum StartType_t { ColdStart, HotStart, TepidStart, CheckpointStart };
-
-  ActionSet<GaugeField, RepresentationsPolicy> TheAction;
-
-  GridCartesian *UGrid;
-  GridCartesian *FGrid;
-  GridRedBlackCartesian *UrbGrid;
-  GridRedBlackCartesian *FrbGrid;
-
-  virtual void BuildTheAction(int argc, char **argv) = 0;  // necessary?
-
-  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 {
-        std::cout << GridLogError << "Unrecognized option in --StartType\n";
-        std::cout << GridLogError << "Valid [HotStart, ColdStart, TepidStart, CheckpointStart]\n";
-        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];
-    }
-
-    int NumThermalizations = 10;
-    if (GridCmdOptionExists(argv, argv + argc, "--Thermalizations")) {
-      arg = GridCmdOptionPayload(argv, argv + argc, "--Thermalizations");
-      std::vector<int> ivec(0);
-      GridCmdOptionIntVector(arg, ivec);
-      NumThermalizations = ivec[0];
-    }
-
-    GridSerialRNG sRNG;
-    GridParallelRNG pRNG(UGrid);
-    LatticeGaugeField U(UGrid);  // change this to an extended field (smearing class)?
-
-    std::vector<int> SerSeed({1, 2, 3, 4, 5});
-    std::vector<int> ParSeed({6, 7, 8, 9, 10});
-
-    // Create integrator, including the smearing policy
-    // Smearing policy, only defined for Nc=3
-    /*
-    std::cout << GridLogDebug << " Creating the Stout class\n";
-    double rho = 0.1;  // smearing parameter, now hardcoded
-    int Nsmear = 1;    // number of smearing levels
-    Smear_Stout<Gimpl> Stout(rho);
-    std::cout << GridLogDebug << " Creating the SmearedConfiguration class\n";
-    //SmearedConfiguration<Gimpl> SmearingPolicy(UGrid, Nsmear, Stout);
-    std::cout << GridLogDebug << " done\n";
-    */
-    //////////////
-    NoSmearing<Gimpl> SmearingPolicy;
-    typedef MinimumNorm2<Gimpl, NoSmearing<Gimpl>, RepresentationsPolicy >
-        IntegratorType;  // change here to change the algorithm
-    IntegratorParameters MDpar(40, 1.0);
-    IntegratorType MDynamics(UGrid, MDpar, TheAction, SmearingPolicy);
-
-    // Checkpoint strategy
-    NerscHmcCheckpointer<Gimpl> Checkpoint(std::string("ckpoint_lat"),
-                                           std::string("ckpoint_rng"), 1);
-    PlaquetteLogger<Gimpl> PlaqLog(std::string("plaq"));
-
-    HMCparameters HMCpar;
-    HMCpar.StartTrajectory = StartTraj;
-    HMCpar.Trajectories = NumTraj;
-    HMCpar.NoMetropolisUntil = NumThermalizations;
-
-    if (StartType == HotStart) {
-      // Hot start
-      HMCpar.MetropolisTest = true;
-      sRNG.SeedFixedIntegers(SerSeed);
-      pRNG.SeedFixedIntegers(ParSeed);
-      SU<Nc>::HotConfiguration(pRNG, U);
-    } else if (StartType == ColdStart) {
-      // Cold start
-      HMCpar.MetropolisTest = true;
-      sRNG.SeedFixedIntegers(SerSeed);
-      pRNG.SeedFixedIntegers(ParSeed);
-      SU<Nc>::ColdConfiguration(pRNG, U);
-    } else if (StartType == TepidStart) {
-      // Tepid start
-      HMCpar.MetropolisTest = true;
-      sRNG.SeedFixedIntegers(SerSeed);
-      pRNG.SeedFixedIntegers(ParSeed);
-      SU<Nc>::TepidConfiguration(pRNG, U);
-    } else if (StartType == CheckpointStart) {
-      HMCpar.MetropolisTest = true;
-      // CheckpointRestart
-      Checkpoint.CheckpointRestore(StartTraj, U, sRNG, pRNG);
-    }
-
-    // Attach the gauge field to the smearing Policy and create the fill the
-    // smeared set
-    // notice that the unit configuration is singular in this procedure
-    std::cout << GridLogMessage << "Filling the smeared set\n";
-    SmearingPolicy.set_Field(U);
-
-    HybridMonteCarlo<IntegratorType> HMC(HMCpar, MDynamics, sRNG, pRNG, U);
-    HMC.AddObservable(&Checkpoint);
-    HMC.AddObservable(&PlaqLog);
-
-    // Run it
-    HMC.evolve();
-  }
-};
-
-typedef NerscHmcRunnerTemplate<PeriodicGimplR> NerscHmcRunner;
-typedef NerscHmcRunnerTemplate<PeriodicGimplF> NerscHmcRunnerF;
-typedef NerscHmcRunnerTemplate<PeriodicGimplD> NerscHmcRunnerD;
-
-typedef NerscHmcRunnerTemplate<PeriodicGimplR> PeriodicNerscHmcRunner;
-typedef NerscHmcRunnerTemplate<PeriodicGimplF> PeriodicNerscHmcRunnerF;
-typedef NerscHmcRunnerTemplate<PeriodicGimplD> PeriodicNerscHmcRunnerD;
-
-typedef NerscHmcRunnerTemplate<ConjugateGimplR> ConjugateNerscHmcRunner;
-typedef NerscHmcRunnerTemplate<ConjugateGimplF> ConjugateNerscHmcRunnerF;
-typedef NerscHmcRunnerTemplate<ConjugateGimplD> ConjugateNerscHmcRunnerD;
-
-template <class RepresentationsPolicy>
-using NerscHmcRunnerHirep = NerscHmcRunnerTemplate<PeriodicGimplR, RepresentationsPolicy>;
-
-
-
-}
-}
-#endif
diff --git a/lib/qcd/modules/Registration.h b/lib/qcd/modules/Registration.h
new file mode 100644
index 00000000..e82cb252
--- /dev/null
+++ b/lib/qcd/modules/Registration.h
@@ -0,0 +1,120 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid
+
+Source file: ./lib/qcd/modules/Registration.h
+
+Copyright (C) 2016
+
+Author: Guido Cossu <guido.cossu@ed.ac.uk>
+
+This program is free software; you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation; either version 2 of the License, or
+(at your option) any later version.
+
+This program is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License along
+with this program; if not, write to the Free Software Foundation, Inc.,
+51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+
+See the full license in the file "LICENSE" in the top level distribution
+directory
+*************************************************************************************/
+/*  END LEGAL */
+
+#ifndef MODULES_REGISTRATION_H
+#define MODULES_REGISTRATION_H
+
+///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+// Actions
+///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+
+
+static Registrar<QCD::WilsonGMod,            HMC_LGTActionModuleFactory<gauge_string, Serialiser > > __WGmodXMLInit("Wilson"); 
+static Registrar<QCD::SymanzikGMod,          HMC_LGTActionModuleFactory<gauge_string, Serialiser> > __SymGmodXMLInit("Symanzik"); 
+static Registrar<QCD::IwasakiGMod,           HMC_LGTActionModuleFactory<gauge_string, Serialiser> > __IwGmodXMLInit("Iwasaki"); 
+static Registrar<QCD::DBW2GMod,              HMC_LGTActionModuleFactory<gauge_string, Serialiser> > __DBW2GmodXMLInit("DBW2"); 
+static Registrar<QCD::RBCGMod,               HMC_LGTActionModuleFactory<gauge_string, Serialiser> > __RBCGmodXMLInit("RBC"); 
+static Registrar<QCD::PlaqPlusRectangleGMod, HMC_LGTActionModuleFactory<gauge_string, Serialiser> > __PPRectGmodXMLInit("PlaqPlusRect"); 
+
+
+// FIXME more general implementation
+static Registrar<QCD::TwoFlavourFModule<FermionImplementationPolicy> ,      
+	HMC_LGTActionModuleFactory<gauge_string, Serialiser> > __TwoFlavourFmodXMLInit("TwoFlavours"); 
+static Registrar<QCD::TwoFlavourRatioFModule<FermionImplementationPolicy> , 
+	HMC_LGTActionModuleFactory<gauge_string, Serialiser> > __TwoFlavourRatioFmodXMLInit("TwoFlavoursRatio"); 
+static Registrar<QCD::TwoFlavourEOFModule<FermionImplementationPolicy> ,    
+	HMC_LGTActionModuleFactory<gauge_string, Serialiser> > __TwoFlavourEOFmodXMLInit("TwoFlavoursEvenOdd"); 
+static Registrar<QCD::TwoFlavourRatioEOFModule<FermionImplementationPolicy>,
+	HMC_LGTActionModuleFactory<gauge_string, Serialiser> > __TwoFlavourRatioEOFmodXMLInit("TwoFlavoursEvenOddRatio"); 
+static Registrar<QCD::OneFlavourFModule<FermionImplementationPolicy> ,      
+	HMC_LGTActionModuleFactory<gauge_string, Serialiser> > __OneFlavourFmodXMLInit("OneFlavour"); 
+static Registrar<QCD::OneFlavourEOFModule<FermionImplementationPolicy> ,    
+	HMC_LGTActionModuleFactory<gauge_string, Serialiser> > __OneFlavourEOFmodXMLInit("OneFlavourEvenOdd"); 
+static Registrar<QCD::OneFlavourRatioFModule<FermionImplementationPolicy> , 
+	HMC_LGTActionModuleFactory<gauge_string, Serialiser> > __OneFlavourRatioFmodXMLInit("OneFlavourRatio"); 
+static Registrar<QCD::OneFlavourRatioEOFModule<FermionImplementationPolicy>,
+	HMC_LGTActionModuleFactory<gauge_string, Serialiser> > __OneFlavourRatioEOFmodXMLInit("OneFlavourEvenOddRatio"); 
+
+
+
+///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+// Solvers
+///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+
+
+// Now a specific registration with a fermion field
+// here must instantiate CG and CR for every new fermion field type (macro!!)
+
+static Registrar< ConjugateGradientModule<QCD::WilsonFermionR::FermionField>,   
+                  HMC_SolverModuleFactory<solver_string, QCD::WilsonFermionR::FermionField, Serialiser> > __CGWFmodXMLInit("ConjugateGradientWF"); 
+static Registrar< ConjugateResidualModule<QCD::WilsonFermionR::FermionField>,   
+                  HMC_SolverModuleFactory<solver_string, QCD::WilsonFermionR::FermionField, Serialiser> > __CRWFmodXMLInit("ConjugateResidualWF"); 
+
+
+///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+// Fermion operators
+///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+
+static Registrar< QCD::WilsonFermionModule<FermionImplementationPolicy>,   
+                  HMC_FermionOperatorModuleFactory<fermionop_string, FermionImplementationPolicy, Serialiser> > __WilsonFOPmodXMLInit("Wilson"); 
+static Registrar< QCD::MobiusFermionModule<FermionImplementationPolicy>,   
+                  HMC_FermionOperatorModuleFactory<fermionop_string, FermionImplementationPolicy, Serialiser> > __MobiusFOPmodXMLInit("Mobius");
+static Registrar< QCD::DomainWallFermionModule<FermionImplementationPolicy>,   
+                  HMC_FermionOperatorModuleFactory<fermionop_string, FermionImplementationPolicy, Serialiser> > __DWFOPmodXMLInit("DomainWall");
+
+
+///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+// Observables
+///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+
+static Registrar<QCD::PlaquetteMod<ImplementationPolicy>, HMC_ObservablesModuleFactory<observable_string, Serialiser> > __OBSPLmodXMLInit("Plaquette"); 
+
+///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+// Checkpointers
+///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+
+static Registrar<QCD::BinaryCPModule<ImplementationPolicy>, HMC_CPModuleFactory<cp_string, ImplementationPolicy, Serialiser> > __CPBinarymodXMLInit("Binary");
+static Registrar<QCD::NerscCPModule<ImplementationPolicy> , HMC_CPModuleFactory<cp_string, ImplementationPolicy, Serialiser> > __CPNerscmodXMLInit("Nersc");
+
+#ifdef HAVE_LIME
+static Registrar<QCD::ILDGCPModule<ImplementationPolicy>  , HMC_CPModuleFactory<cp_string, ImplementationPolicy, Serialiser> > __CPILDGmodXMLInit("ILDG");
+#endif
+
+
+
+///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+// Integrators
+///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+static Registrar< HMCLeapFrog<ImplementationPolicy, RepresentationPolicy, Serialiser>      , HMCRunnerModuleFactory<hmc_string, Serialiser> > __HMCLFmodXMLInit("LeapFrog");
+static Registrar< HMCMinimumNorm2<ImplementationPolicy, RepresentationPolicy, Serialiser>  , HMCRunnerModuleFactory<hmc_string, Serialiser> > __HMCMN2modXMLInit("MinimumNorm2");
+static Registrar< HMCForceGradient<ImplementationPolicy, RepresentationPolicy, Serialiser> , HMCRunnerModuleFactory<hmc_string, Serialiser> > __HMCFGmodXMLInit("ForceGradient");
+
+typedef HMCRunnerModuleFactory<hmc_string, Serialiser > HMCModuleFactory;
+
+#endif
\ No newline at end of file
diff --git a/lib/qcd/representations/hmc_types.h b/lib/qcd/representations/hmc_types.h
index 65ba7062..66f178c4 100644
--- a/lib/qcd/representations/hmc_types.h
+++ b/lib/qcd/representations/hmc_types.h
@@ -4,7 +4,7 @@
 #include <Grid/qcd/representations/adjoint.h>
 #include <Grid/qcd/representations/two_index.h>
 #include <Grid/qcd/representations/fundamental.h>
-#include <Grid/qcd/action/scalar/scalarImpl.h>
+#include <Grid/qcd/action/scalar/ScalarImpl.h>
 
 #include <tuple>
 #include <utility>
diff --git a/lib/qcd/utils/ScalarObjs.h b/lib/qcd/utils/ScalarObjs.h
index 9915a89b..a284af40 100644
--- a/lib/qcd/utils/ScalarObjs.h
+++ b/lib/qcd/utils/ScalarObjs.h
@@ -6,10 +6,7 @@
 
     Copyright (C) 2015
 
-Author: Azusa Yamaguchi <ayamaguc@staffmail.ed.ac.uk>
-Author: Peter Boyle <paboyle@ph.ed.ac.uk>
 Author: neo <cossu@post.kek.jp>
-Author: paboyle <paboyle@ph.ed.ac.uk>
 
     This program is free software; you can redistribute it and/or modify
     it under the terms of the GNU General Public License as published by
@@ -32,7 +29,9 @@ directory
 #ifndef SCALAR_OBJS_H
 #define SCALAR_OBJS_H
 namespace Grid {
-namespace QCD {
+
+  // FIXME drop the QCD namespace in Nd
+  
 
 // Scalar field obs
 template <class Impl>
@@ -62,7 +61,7 @@ class ScalarObs {
   static void phider(typename Impl::Field &fsq,
                      const typename Impl::Field &f) {
     fsq = Cshift(f, 0, -1) * f;
-    for (int mu = 1; mu < Nd; mu++) fsq += Cshift(f, mu, -1) * f;
+    for (int mu = 1; mu < QCD::Nd; mu++) fsq += Cshift(f, mu, -1) * f;
   }
 
   //////////////////////////////////////////////////
@@ -72,7 +71,7 @@ class ScalarObs {
   static RealD sumphider(const typename Impl::Field &f) {
     typename Impl::Field tmp(f._grid);
     tmp = Cshift(f, 0, -1) * f;
-    for (int mu = 1; mu < Nd; mu++) {
+    for (int mu = 1; mu < QCD::Nd; mu++) {
       tmp += Cshift(f, mu, -1) * f;
     }
     return -sum(trace(tmp));
@@ -90,7 +89,8 @@ class ScalarObs {
     return sum(trace(tmp));
   }
 };
-}
+
+
 }
 
-#endif
\ No newline at end of file
+#endif
diff --git a/tests/hmc/Test_hmc_Factories.cc b/tests/hmc/Test_hmc_Factories.cc
index 50c47ed2..de571ab7 100644
--- a/tests/hmc/Test_hmc_Factories.cc
+++ b/tests/hmc/Test_hmc_Factories.cc
@@ -31,36 +31,35 @@ directory
 namespace Grid{
 
 // Put this section in a separate header
-
 // ifdefs ?? Local makefile suggestion , policy as make parameter
-typedef QCD::PeriodicGimplR ImplementationPolicy;
-typedef QCD::NoHirep RepresentationPolicy;
-typedef QCD::WilsonFermionR FermionImplementation;
+typedef QCD::PeriodicGimplR   ImplementationPolicy;
+typedef QCD::WilsonImplR      FermionImplementationPolicy;
+typedef QCD::NoHirep          RepresentationPolicy;
+typedef Grid::XmlReader       Serialiser;
 
-///////////////////////////////////////////////////////////////////////
-// Put all registrations in an header file
-static Registrar< HMCLeapFrog<ImplementationPolicy, RepresentationPolicy, XmlReader>      , HMCRunnerModuleFactory<hmc_string, XmlReader> > __HMCLFmodXMLInit("LeapFrog");
-static Registrar< HMCMinimumNorm2<ImplementationPolicy, RepresentationPolicy, XmlReader>  , HMCRunnerModuleFactory<hmc_string, XmlReader> > __HMCMN2modXMLInit("MinimumNorm2");
-static Registrar< HMCForceGradient<ImplementationPolicy, RepresentationPolicy, XmlReader> , HMCRunnerModuleFactory<hmc_string, XmlReader> > __HMCFGmodXMLInit("ForceGradient");
-}
+// Register all object names
+#include "Grid/qcd/modules/Registration.h"
+
+} // Grid
 
 int main(int argc, char **argv) {
   using namespace Grid;
   using namespace Grid::QCD;
 
   Grid_init(&argc, &argv);
+
   int threads = GridThread::GetThreads();
   // here make a routine to print all the relevant information on the run
   std::cout << GridLogMessage << "Grid is setup to use " << threads << " threads" << std::endl;
 
    // Typedefs to simplify notation 
-  typedef Grid::XmlReader InputFileReader; 
+  //typedef XmlReader InputFileReader; 
 
   // Reader, file should come from command line
-  InputFileReader Reader("input.wilson_gauge.params.xml");
+  Serialiser Reader("input.wilson_gauge.params.xml");
 
   // Test HMC factory (put in an external file)
-  auto &HMCfactory = HMCRunnerModuleFactory<hmc_string, InputFileReader >::getInstance();
+  auto &HMCfactory = HMCModuleFactory::getInstance();
   // Simplify this step (IntergratorName field?)
   HMCparameters HMCpar(Reader);
   
@@ -71,5 +70,7 @@ int main(int argc, char **argv) {
   HMCmodule->getPtr()->Run();
 
   Grid_finalize();
+  return 0;
 
 } // main
+