diff --git a/lib/qcd/action/Actions.h b/lib/qcd/action/Actions.h
index 5f83aa55..d5888c13 100644
--- a/lib/qcd/action/Actions.h
+++ b/lib/qcd/action/Actions.h
@@ -101,11 +101,10 @@ typedef SymanzikGaugeAction<ConjugateGimplF>        ConjugateSymanzikGaugeAction
 typedef SymanzikGaugeAction<ConjugateGimplD>        ConjugateSymanzikGaugeActionD;
 
 
-  typedef ScalarAction<ScalarImplR>                 ScalarActionR;
-  typedef ScalarAction<ScalarImplF>                 ScalarActionF;
-  typedef ScalarAction<ScalarImplD>                 ScalarActionD;
-  
-  
+typedef ScalarAction<ScalarImplR>                 ScalarActionR;
+typedef ScalarAction<ScalarImplF>                 ScalarActionF;
+typedef ScalarAction<ScalarImplD>                 ScalarActionD;
+
 }}
 
 ////////////////////////////////////////////////////////////////////////////////////////////////////
diff --git a/lib/qcd/action/gauge/GaugeImplTypes.h b/lib/qcd/action/gauge/GaugeImplTypes.h
index 4cfb4780..1d716be1 100644
--- a/lib/qcd/action/gauge/GaugeImplTypes.h
+++ b/lib/qcd/action/gauge/GaugeImplTypes.h
@@ -81,7 +81,7 @@ public:
   ///////////////////////////////////////////////////////////
   // Move these to another class
   // HMC auxiliary functions
-  static inline void generate_momenta(Field& P, GridParallelRNG& pRNG){
+  static inline void generate_momenta(Field &P, GridParallelRNG &pRNG) {
     // specific for SU gauge fields
     LinkField Pmu(P._grid);
     Pmu = zero;
@@ -91,9 +91,7 @@ public:
     }
   }
 
-   static inline Field projectForce(Field& P){
-    return Ta(P);
-   }
+  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++) {
diff --git a/lib/qcd/action/scalar/ScalarImpl.h b/lib/qcd/action/scalar/ScalarImpl.h
index a259c34b..ee2d2fb8 100644
--- a/lib/qcd/action/scalar/ScalarImpl.h
+++ b/lib/qcd/action/scalar/ScalarImpl.h
@@ -45,6 +45,49 @@ namespace Grid {
     }
     
   };
+
+  template <class S, unsigned int N>
+  class ScalarMatrixImplTypes {
+  public:
+    typedef S Simd;
+    
+    template <typename vtype>
+    using iImplField = iScalar<iScalar<iMatrix<vtype, N> > >;
+    
+    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 (TensorRemove(- sum(trace(U*U))*0.5).real());
+    }
+    
+    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;
diff --git a/lib/qcd/action/scalar/ScalarInteractionAction.h b/lib/qcd/action/scalar/ScalarInteractionAction.h
new file mode 100644
index 00000000..bd54a010
--- /dev/null
+++ b/lib/qcd/action/scalar/ScalarInteractionAction.h
@@ -0,0 +1,84 @@
+/*************************************************************************************
+
+  Grid physics library, www.github.com/paboyle/Grid
+
+  Source file: ./lib/qcd/action/gauge/WilsonGaugeAction.h
+
+  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
+  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 SCALAR_ACTION_H
+#define SCALAR_ACTION_H
+
+namespace Grid {
+  // FIXME drop the QCD namespace everywhere here
+  
+  template <class Impl>
+  class ScalarInteractionAction : 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 LogParameters(){
+      std::stringstream sstream;
+      sstream << GridLogMessage << "[ScalarAction] lambda      : " << lambda      << std::endl;
+      sstream << GridLogMessage << "[ScalarAction] mass_square : " << mass_square << std::endl;
+      return sstream.str();
+      
+    }
+    
+    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
+
+#endif // SCALAR_ACTION_H
diff --git a/lib/qcd/hmc/GenericHMCrunner.h b/lib/qcd/hmc/GenericHMCrunner.h
index 044af52a..085b5a35 100644
--- a/lib/qcd/hmc/GenericHMCrunner.h
+++ b/lib/qcd/hmc/GenericHMCrunner.h
@@ -94,8 +94,8 @@ class HMCWrapperTemplate: public HMCRunnerBase<ReaderClass> {
       Parameters.StartingType = arg;
     }
 
-    if (GridCmdOptionExists(argv, argv + argc, "--StartTrajectory")) {
-      arg = GridCmdOptionPayload(argv, argv + argc, "--StartTrajectory");
+    if (GridCmdOptionExists(argv, argv + argc, "--StartingTrajectory")) {
+      arg = GridCmdOptionPayload(argv, argv + argc, "--StartingTrajectory");
       std::vector<int> ivec(0);
       GridCmdOptionIntVector(arg, ivec);
       Parameters.StartTrajectory = ivec[0];
diff --git a/lib/qcd/hmc/HMC.h b/lib/qcd/hmc/HMC.h
index 98f4f6ad..be469f82 100644
--- a/lib/qcd/hmc/HMC.h
+++ b/lib/qcd/hmc/HMC.h
@@ -40,6 +40,7 @@ directory
 #define HMC_INCLUDED
 
 #include <string>
+#include <list>
 
 namespace Grid {
 namespace QCD {
@@ -158,7 +159,7 @@ class HybridMonteCarlo {
   GridParallelRNG &pRNG; 
 
   Field &Ucur;
-
+  
   IntegratorType &TheIntegrator;
 	ObsListType Observables;
 
@@ -195,7 +196,7 @@ class HybridMonteCarlo {
   /////////////////////////////////////////////////////////
   // Evolution
   /////////////////////////////////////////////////////////
-  RealD evolve_step(Field &U) {
+  RealD evolve_hmc_step(Field &U) {
     TheIntegrator.refresh(U, pRNG);  // set U and initialize P and phi's
 
     RealD H0 = TheIntegrator.S(U);  // initial state action
@@ -211,11 +212,14 @@ class HybridMonteCarlo {
 
     std::cout.precision(17);
     std::cout << GridLogMessage << "Total H after trajectory  = " << H1
-              << "  dH = " << H1 - H0 << "\n";
+	      << "  dH = " << H1 - H0 << "\n";
     std::cout.precision(current_precision);
-
+    
     return (H1 - H0);
   }
+  
+
+  
 
  public:
   /////////////////////////////////////////
@@ -224,7 +228,7 @@ class HybridMonteCarlo {
   HybridMonteCarlo(HMCparameters _Pams, IntegratorType &_Int,
                    GridSerialRNG &_sRNG, GridParallelRNG &_pRNG, 
                    ObsListType _Obs, Field &_U)
-      : Params(_Pams), TheIntegrator(_Int), sRNG(_sRNG), pRNG(_pRNG), Observables(_Obs), Ucur(_U) {}
+    : Params(_Pams), TheIntegrator(_Int), sRNG(_sRNG), pRNG(_pRNG), Observables(_Obs), Ucur(_U) {}
   ~HybridMonteCarlo(){};
 
   void evolve(void) {
@@ -241,13 +245,13 @@ class HybridMonteCarlo {
       std::cout << GridLogMessage << "-- # Trajectory = " << traj << "\n";
       if (traj < Params.StartTrajectory + Params.NoMetropolisUntil) {
       	std::cout << GridLogMessage << "-- Thermalization" << std::endl;
-    	}
-
-    	double t0=usecond();
+      }
+      
+      double t0=usecond();
       Ucopy = Ucur;
 
-      DeltaH = evolve_step(Ucopy);
-
+      DeltaH = evolve_hmc_step(Ucopy);
+      // Metropolis-Hastings test
       bool accept = true;
       if (traj >= Params.StartTrajectory + Params.NoMetropolisUntil) {
         accept = metropolis_test(DeltaH);
@@ -255,9 +259,11 @@ class HybridMonteCarlo {
       	std::cout << GridLogMessage << "Skipping Metropolis test" << std::endl;
       }
 
-      if (accept) {
-        Ucur = Ucopy;
-      }
+      if (accept)
+        Ucur = Ucopy; 
+      
+     
+      
       double t1=usecond();
       std::cout << GridLogMessage << "Total time for trajectory (s): " << (t1-t0)/1e6 << std::endl;
 
diff --git a/lib/qcd/hmc/HMCResourceManager.h b/lib/qcd/hmc/HMCResourceManager.h
index b535ed03..aaa49c8c 100644
--- a/lib/qcd/hmc/HMCResourceManager.h
+++ b/lib/qcd/hmc/HMCResourceManager.h
@@ -75,6 +75,13 @@ class HMCResourceManager {
   bool have_RNG;
   bool have_CheckPointer;
 
+  void output_vector_string(const std::vector<std::string> &vs){
+    for (auto &i: vs)
+      std::cout << i << " ";
+    std::cout << std::endl;
+  }
+
+
  public:
   HMCResourceManager() : have_RNG(false), have_CheckPointer(false) {}
 
@@ -93,7 +100,11 @@ class HMCResourceManager {
     std::string cp_type;
     read(Read,"name", cp_type);
     std::cout << "Registered types " << std::endl;
-    std::cout << CPfactory.getBuilderList() << std::endl;
+    // NOTE: operator << is not overloaded for std::vector<string> 
+    // so it complains here
+    //std::cout << CPfactory.getBuilderList() << std::endl;
+    output_vector_string(CPfactory.getBuilderList());
+
 
     CP = CPfactory.create(cp_type, Read);
     CP->print_parameters();
@@ -110,7 +121,7 @@ class HMCResourceManager {
       std::string obs_type;
       read(Read,"name", obs_type);
       std::cout << "Registered types " << std::endl;
-      std::cout << ObsFactory.getBuilderList() << std::endl;
+      output_vector_string(ObsFactory.getBuilderList() );
 
       ObservablesList.emplace_back(ObsFactory.create(obs_type, Read));
       ObservablesList[ObservablesList.size() - 1]->print_parameters();
@@ -271,8 +282,8 @@ private:
     do{
       auto &ActionFactory = HMC_ActionModuleFactory<gauge_string, typename ImplementationPolicy::Field, ReaderClass>::getInstance(); 
       std::string action_type;
-      Read.readDefault("name", action_type);
-      std::cout << ActionFactory.getBuilderList() << std::endl;  // temporary
+      Read.readDefault("name", action_type); 
+      output_vector_string(ActionFactory.getBuilderList() );
       ActionsList.emplace(m, ActionFactory.create(action_type, Read));
     } while (Read.nextElement("Action"));
     ActionsList.find(m)->second->print_parameters();    
diff --git a/lib/qcd/hmc/integrators/Integrator.h b/lib/qcd/hmc/integrators/Integrator.h
index 6a72919e..8033bd41 100644
--- a/lib/qcd/hmc/integrators/Integrator.h
+++ b/lib/qcd/hmc/integrators/Integrator.h
@@ -39,8 +39,8 @@ namespace QCD {
 
 class IntegratorParameters: Serializable {
 public:
-	GRID_SERIALIZABLE_CLASS_MEMBERS(IntegratorParameters,
-		std::string, name,      // name of the integrator
+        GRID_SERIALIZABLE_CLASS_MEMBERS(IntegratorParameters,
+                std::string, name,      // name of the integrator
     unsigned int, MDsteps,  // number of outer steps
     RealD, trajL,           // trajectory length
   )
@@ -55,7 +55,7 @@ public:
   template <class ReaderClass, typename std::enable_if<isReader<ReaderClass>::value, int >::type = 0 >
   IntegratorParameters(ReaderClass & Reader){
     std::cout << "Reading integrator\n";
-  	read(Reader, "Integrator", *this);
+        read(Reader, "Integrator", *this);
   }
 
   void print_parameters() const {
@@ -177,16 +177,16 @@ class Integrator {
   }
 
   void print_actions(){
-  	std::cout << GridLogMessage << ":::::::::::::::::::::::::::::::::::::::::" << std::endl;
-  	std::cout << GridLogMessage << "[Integrator] Action summary: "<<std::endl;
-  	for (int level = 0; level < as.size(); ++level) {
-  		std::cout << GridLogMessage << "[Integrator] ---- Level: "<< level << std::endl;
-  		for (int actionID = 0; actionID < as[level].actions.size(); ++actionID) {
-  			std::cout << GridLogMessage << "["<< as[level].actions.at(actionID)->action_name() << "] ID: " << actionID << std::endl;
-  			std::cout << as[level].actions.at(actionID)->LogParameters();
-  		}
-  	}
-  	std::cout << GridLogMessage << ":::::::::::::::::::::::::::::::::::::::::"<< std::endl;
+        std::cout << GridLogMessage << ":::::::::::::::::::::::::::::::::::::::::" << std::endl;
+        std::cout << GridLogMessage << "[Integrator] Action summary: "<<std::endl;
+        for (int level = 0; level < as.size(); ++level) {
+                std::cout << GridLogMessage << "[Integrator] ---- Level: "<< level << std::endl;
+                for (int actionID = 0; actionID < as[level].actions.size(); ++actionID) {
+                        std::cout << GridLogMessage << "["<< as[level].actions.at(actionID)->action_name() << "] ID: " << actionID << std::endl;
+                        std::cout << as[level].actions.at(actionID)->LogParameters();
+                }
+        }
+        std::cout << GridLogMessage << ":::::::::::::::::::::::::::::::::::::::::"<< std::endl;
 
   }
 
@@ -209,7 +209,7 @@ class Integrator {
     assert(P._grid == U._grid);
     std::cout << GridLogIntegrator << "Integrator refresh\n";
     FieldImplementation::generate_momenta(P, pRNG);
- 
+     
     // Update the smeared fields, can be implemented as observer
     // necessary to keep the fields updated even after a reject
     // of the Metropolis
diff --git a/lib/qcd/modules/ActionModules.h b/lib/qcd/modules/ActionModules.h
index cf7f4c07..fbfbbf33 100644
--- a/lib/qcd/modules/ActionModules.h
+++ b/lib/qcd/modules/ActionModules.h
@@ -510,28 +510,6 @@ class HMC_ActionModuleFactory
 };
 
 
-/*
-template <char const *str, class ReaderClass >
-class HMC_ScalarActionModuleFactory
-    : public Factory < HMC_ScalarActionModBase , Reader<ReaderClass> > {
- public:
-  typedef Reader<ReaderClass> TheReader; 
-  // use SINGLETON FUNCTOR MACRO HERE
-  HMC_ScalarActionModuleFactory(const HMC_ScalarActionModuleFactory& e) = delete;
-  void operator=(const HMC_ScalarActionModuleFactory& e) = delete;
-  static HMC_ScalarActionModuleFactory& getInstance(void) {
-    static HMC_ScalarActionModuleFactory e;
-    return e;
-  }
-
- private:
-  HMC_ScalarActionModuleFactory(void) = default;
-    std::string obj_type() const {
-        return std::string(str);
-  }
-};
-*/
-
 extern char gauge_string[];
 } // Grid
 
diff --git a/lib/qcd/representations/hmc_types.h b/lib/qcd/representations/hmc_types.h
index 66f178c4..3701c9b2 100644
--- a/lib/qcd/representations/hmc_types.h
+++ b/lib/qcd/representations/hmc_types.h
@@ -62,6 +62,7 @@ class Representations {
 
 typedef Representations<FundamentalRepresentation> NoHirep;
 typedef Representations<EmptyRep<typename ScalarImplR::Field> > ScalarFields;
+  //typedef Representations<EmptyRep<typename ScalarMatrixImplR::Field> > ScalarMatrixFields;
 
 // Helper classes to access the elements
 // Strips the first N parameters from the tuple
diff --git a/tests/hmc/Test_hmc_EODWFRatio_Binary.cc b/tests/hmc/Test_hmc_EODWFRatio_Binary.cc
index 49d4e94d..1acd05f8 100644
--- a/tests/hmc/Test_hmc_EODWFRatio_Binary.cc
+++ b/tests/hmc/Test_hmc_EODWFRatio_Binary.cc
@@ -2,12 +2,11 @@
 
 Grid physics library, www.github.com/paboyle/Grid
 
-Source file: ./tests/Test_hmc_WilsonFermionGauge.cc
+Source file: ./tests/Test_hmc_EODWFRatio.cc
 
-Copyright (C) 2015
+Copyright (C) 2015-2016
 
 Author: Peter Boyle <pabobyle@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
@@ -30,31 +29,135 @@ directory
 /*  END LEGAL */
 #include <Grid/Grid.h>
 
-using namespace std;
-using namespace Grid;
-using namespace Grid::QCD;
 
-namespace Grid {
-namespace QCD {
 
-  class HMCRunnerParameters : Serializable {
-  public:
-    GRID_SERIALIZABLE_CLASS_MEMBERS(HMCRunnerParameters,
-      double, beta,
-      double, mass,
-      int, MaxCGIterations,
-      double, StoppingCondition,
-      bool, smearedAction,
-      int, SaveInterval,
-      std::string, format,
-      std::string, conf_prefix,
-      std::string, rng_prefix,
-      double, rho,
-      int, SmearingLevels,
-      );
+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 GenericHMCRunner<MinimumNorm2> HMCWrapper;  // Uses the default minimum norm
+  typedef WilsonImplR FermionImplPolicy;
+  typedef DomainWallFermionR FermionAction;
+  typedef typename FermionAction::FermionField FermionField;
+
+
+  //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
+  HMCWrapper TheHMC;
+
+  // Grid from the command line
+  TheHMC.Resources.AddFourDimGrid("gauge");
+  // Possibile to create the module by hand 
+  // hardcoding parameters or using a Reader
+
+
+  // Checkpointer definition
+  CheckpointerParameters CPparams;  
+  CPparams.config_prefix = "ckpoint_EODWF_lat";
+  CPparams.rng_prefix = "ckpoint_EODWF_rng";
+  CPparams.saveInterval = 5;
+  CPparams.format = "IEEE64BIG";
+  
+  TheHMC.Resources.LoadBinaryCheckpointer(CPparams);
+
+  RNGModuleParameters RNGpar;
+  RNGpar.SerialSeed = {1,2,3,4,5};
+  RNGpar.ParallelSeed = {6,7,8,9,10};
+  TheHMC.Resources.SetRNGSeeds(RNGpar);
+
+  // Construct observables
+  // here there is too much indirection 
+  PlaquetteObsParameters PlPar;
+  PlPar.output_prefix = "Plaquette";
+  PlaquetteMod<HMCWrapper::ImplPolicy> PlaqModule(PlPar);
+  TheHMC.Resources.AddObservable(&PlaqModule);
+  //////////////////////////////////////////////
+
+  /////////////////////////////////////////////////////////////
+  // Collect actions, here use more encapsulation
+  // need wrappers of the fermionic classes 
+  // that have a complex construction
+  // standard
+  RealD beta = 5.6 ;
+  WilsonGaugeActionR Waction(beta);
+    
+  auto GridPtr   = TheHMC.Resources.GetCartesian();
+  auto GridRBPtr = TheHMC.Resources.GetRBCartesian();
+  auto FGrid     = SpaceTimeGrid::makeFiveDimGrid(Ls,GridPtr);
+  auto FrbGrid   = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,GridPtr);
+
+
+  // temporarily need a gauge field
+  LatticeGaugeField U(GridPtr);
+
+  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);
+
+  double StoppingCondition = 1.0e-8;
+  double MaxCGIterations = 10000;
+  ConjugateGradient<FermionField>  CG(StoppingCondition,MaxCGIterations);
+  TwoFlavourEvenOddRatioPseudoFermionAction<ImplPolicy> Nf2(NumOp, DenOp,CG,CG);
+
+    // Set smearing (true/false), default: false
+  Nf2.is_smeared = true;
+
+  // Collect actions
+  ActionLevel<HMCWrapper::Field> Level1(1);
+  Level1.push_back(&Nf2);
+
+  ActionLevel<HMCWrapper::Field> Level2(4);
+  Level2.push_back(&Waction);
+
+  TheHMC.TheAction.push_back(Level1);
+  TheHMC.TheAction.push_back(Level2);
+  /////////////////////////////////////////////////////////////
+
+  /*
+    double rho = 0.1;  // smearing parameter
+    int Nsmear = 2;    // number of smearing levels
+    Smear_Stout<HMCWrapper::ImplPolicy> Stout(rho);
+    SmearedConfiguration<HMCWrapper::ImplPolicy> SmearingPolicy(
+        UGrid, Nsmear, Stout);
+  */
+
+  // HMC parameters are serialisable 
+  TheHMC.Parameters.MD.MDsteps = 20;
+  TheHMC.Parameters.MD.trajL   = 1.0;
+
+  TheHMC.ReadCommandLine(argc, argv); // these can be parameters from file
+
+  // Reset performance counters 
+  NumOp.ZeroCounters();
+  DenOp.ZeroCounters();
+  TheHMC.Run();  // no smearing
+  // TheHMC.Run(SmearingPolicy); // for smearing
+
+  std::cout << GridLogMessage << "Numerator report, Pauli-Villars term         : " << std::endl;
+  NumOp.Report();
+  std::cout << GridLogMessage << "Denominator report, Dw(m) term (includes CG) : " << std::endl;
+  DenOp.Report();
+
+
+
+
+  Grid_finalize();
+
+
+
+
+} // main
+
+
 
-    HMCRunnerParameters() {}
-  };
 
 // Derive from the BinaryHmcRunner (templated for gauge fields)
 class HmcRunner : public BinaryHmcRunner {
@@ -156,25 +259,3 @@ public:
 };
 }
 }
-
-int main(int argc, char **argv) {
-  Grid_init(&argc, &argv);
-
-  int threads = GridThread::GetThreads();
-  std::cout << GridLogMessage << "Grid is setup to use " << threads
-            << " threads" << std::endl;
-
-  HmcRunner TheHMC;
-
-
-  // Seeds for the random number generators
-  std::vector<int> SerSeed({1, 2, 3, 4, 5});
-  std::vector<int> ParSeed({6, 7, 8, 9, 10});
-  TheHMC.RNGSeeds(SerSeed, ParSeed);
-
-  TheHMC.MDparameters.set(40, 1.0);// MDsteps, traj length
-
-  TheHMC.BuildTheAction(argc, argv);
-
-  Grid_finalize();
-}
diff --git a/tests/hmc/Test_hmc_ScalarAction.cc b/tests/hmc/Test_hmc_ScalarAction.cc
index d40d43ef..437854df 100644
--- a/tests/hmc/Test_hmc_ScalarAction.cc
+++ b/tests/hmc/Test_hmc_ScalarAction.cc
@@ -37,8 +37,9 @@ int main(int argc, char **argv) {
   std::cout << GridLogMessage << "Grid is setup to use " << threads << " threads" << std::endl;
 
    // Typedefs to simplify notation
-  typedef ScalarGenericHMCRunner HMCWrapper;  // Uses the default minimum norm
-
+  typedef ScalarGenericHMCRunner HMCWrapper;  // Uses the default minimum norm, real scalar fields
+  //typedef Representations<EmptyRep<typename ScalarMatrixImplTypes<vComplex, 3>::Field> > ScalarMatrixFields;
+  //typedef HMCWrapperTemplate<ScalarMatrixImplTypes<vComplex, 3>, MinimumNorm2, ScalarMatrixFields> HMCWrapper;
   //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
   HMCWrapper TheHMC;
 
@@ -76,6 +77,9 @@ int main(int argc, char **argv) {
 
   // Real Scalar action
   ScalarActionR Saction(0.11,0.);
+  
+  //typedef ScalarAction<ScalarMatrixImplTypes<vComplex, 3>> ScalarMatrixActionR;
+  //ScalarMatrixActionR Saction(0.11,0.);
 
   // Collect actions
   ActionLevel<ScalarActionR::Field, ScalarFields> Level1(1);
diff --git a/tests/hmc/Test_hmc_WilsonFermionGauge.cc b/tests/hmc/Test_hmc_WilsonFermionGauge.cc
deleted file mode 100644
index 351d1e68..00000000
--- a/tests/hmc/Test_hmc_WilsonFermionGauge.cc
+++ /dev/null
@@ -1,99 +0,0 @@
-/*************************************************************************************
-
-Grid physics library, www.github.com/paboyle/Grid
-
-Source file: ./tests/Test_hmc_WilsonFermionGauge.cc
-
-Copyright (C) 2015
-
-Author: Peter Boyle <pabobyle@ph.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
-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 */
-#include <Grid/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);
-
-    // Set smearing (true/false), default: false
-    Nf2.is_smeared = false;
-
-    // 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);
-
-  int threads = GridThread::GetThreads();
-  std::cout << GridLogMessage << "Grid is setup to use " << threads
-            << " threads" << std::endl;
-
-  HmcRunner TheHMC;
-
-  TheHMC.BuildTheAction(argc, argv);
-}
diff --git a/tests/hmc/Test_hmc_WilsonGauge.cc b/tests/hmc/Test_hmc_WilsonGauge.cc
deleted file mode 100644
index d2779ad2..00000000
--- a/tests/hmc/Test_hmc_WilsonGauge.cc
+++ /dev/null
@@ -1,86 +0,0 @@
-    /*************************************************************************************
-
-    Grid physics library, www.github.com/paboyle/Grid 
-
-    Source file: ./tests/Test_hmc_WilsonGauge.cc
-
-    Copyright (C) 2015
-
-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
-    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 */
-#include <Grid/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);
-
-  int threads = GridThread::GetThreads();
-  std::cout<<GridLogMessage << "Grid is setup to use "<<threads<<" threads"<<std::endl;
-
-  HmcRunner TheHMC;
-  
-  TheHMC.BuildTheAction(argc,argv);
-
-}
-
diff --git a/tests/hmc/Test_hmc_WilsonGauge_Binary.cc b/tests/hmc/Test_hmc_WilsonGauge_Binary.cc
index 28bbafc3..f37c0f03 100644
--- a/tests/hmc/Test_hmc_WilsonGauge_Binary.cc
+++ b/tests/hmc/Test_hmc_WilsonGauge_Binary.cc
@@ -53,7 +53,7 @@ int main(int argc, char **argv) {
   CheckpointerParameters CPparams;  
   CPparams.config_prefix = "ckpoint_lat";
   CPparams.rng_prefix = "ckpoint_rng";
-  CPparams.saveInterval = 5;
+  CPparams.saveInterval = 20;
   CPparams.format = "IEEE64BIG";
   
   TheHMC.Resources.LoadBinaryCheckpointer(CPparams);
@@ -76,7 +76,7 @@ int main(int argc, char **argv) {
   // need wrappers of the fermionic classes 
   // that have a complex construction
   // standard
-  RealD beta = 5.6 ;
+  RealD beta = 6.0 ;
   WilsonGaugeActionR Waction(beta);
   
   ActionLevel<HMCWrapper::Field> Level1(1);
@@ -86,8 +86,8 @@ int main(int argc, char **argv) {
   /////////////////////////////////////////////////////////////
 
   // HMC parameters are serialisable 
-  TheHMC.Parameters.MD.MDsteps = 20;
-  TheHMC.Parameters.MD.trajL   = 1.0;
+  TheHMC.Parameters.MD.MDsteps = 40;
+  TheHMC.Parameters.MD.trajL   = 2.0;
 
   TheHMC.ReadCommandLine(argc, argv); // these can be parameters from file
   TheHMC.Run();  // no smearing