Adding fermions actions support in the factories

lib/lattice/Lattice_rng.h

@@ -278,7 +278,7 @@ namespace Grid {
     GridBase *_grid;
     unsigned int _vol;
 
-    int generator_idx(int os,int is){
+    int generator_idx(int os,int is) {
       return is*_grid->oSites()+os;
     }
 

lib/qcd/action/ActionBase.h

@@ -37,11 +37,11 @@ class Action {
  public:
   bool is_smeared = false;
   // Heatbath?
-  virtual void refresh(const GaugeField& U, const GridParallelRNG& pRNG) = 0; // refresh pseudofermions
+  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 std::string LogParameters()  = 0;                       // prints action parameters
   virtual ~Action(){}
 };
 

lib/qcd/action/gauge/PlaqPlusRectangleAction.h

@@ -52,6 +52,14 @@ namespace Grid{
       
       virtual void refresh(const GaugeField &U, GridParallelRNG& pRNG) {}; // noop as no pseudoferms
       
+      virtual std::string LogParameters(){
+      	std::stringstream sstream;
+      	sstream << GridLogMessage << "["<<action_name() <<"] c_plaq: " << c_plaq << std::endl;
+      	sstream << GridLogMessage << "["<<action_name() <<"] c_rect: " << c_rect << std::endl;
+      	return sstream.str();
+      }
+
+
       virtual RealD S(const GaugeField &U) {
 	RealD vol = U._grid->gSites();
 

lib/qcd/action/gauge/WilsonGaugeAction.h

@@ -56,7 +56,7 @@ class WilsonGaugeAction : public Action<typename Gimpl::GaugeField> {
   }
 
   virtual void refresh(const GaugeField &U,
-                       const GridParallelRNG &pRNG){};  // noop as no pseudoferms
+                       GridParallelRNG &pRNG){};  // noop as no pseudoferms
 
   virtual RealD S(const GaugeField &U) {
     RealD plaq = WilsonLoops<Gimpl>::avgPlaquette(U);

lib/qcd/action/pseudofermion/TwoFlavour.h

@@ -64,6 +64,12 @@ class TwoFlavourPseudoFermionAction : public Action<typename Impl::GaugeField> {
 
 
   virtual std::string action_name(){return "TwoFlavourPseudoFermionAction";}
+
+  virtual std::string LogParameters(){
+    std::stringstream sstream;
+    sstream << GridLogMessage << "["<<action_name()<<"] has no parameters" << std::endl;
+    return sstream.str();
+  }  
   
   //////////////////////////////////////////////////////////////////////////////////////
   // Push the gauge field in to the dops. Assume any BC's and smearing already
@@ -85,7 +91,7 @@ class TwoFlavourPseudoFermionAction : public Action<typename Impl::GaugeField> {
     //         the partition function.
     //
 
-		RealD scale = std::sqrt(0.5);
+    RealD scale = std::sqrt(0.5);
 
     FermionField eta(FermOp.FermionGrid());
 

lib/qcd/hmc/HMCResourceManager.h

@@ -48,7 +48,6 @@ with this program; if not, write to the Free Software Foundation, Inc.,
     }                                                                    \
   }
 
-
 namespace Grid {
 namespace QCD {
 
@@ -131,31 +130,7 @@ class HMCResourceManager {
   }
 
 
-  // this private
+ 
-  template <class ReaderClass >
-  void fill_ActionsLevel(ReaderClass &Read){
-    // Actions set
-    int m;
-    Read.readDefault("multiplier",m);
-    multipliers.push_back(m);
-    std::cout << "Level : " << multipliers.size()  << " with multiplier : " << m << std::endl; 
-    // here gauge
-    Read.push("Action");
-    do{
-      auto &ActionFactory = HMC_LGTActionModuleFactory<gauge_string, ReaderClass>::getInstance(); 
-      std::string action_type;
-      Read.readDefault("name", action_type);
-      std::cout << "Registered types " << std::endl;
-      std::cout << ActionFactory.getBuilderList() << std::endl;  
-
-      ActionsList.emplace(m, ActionFactory.create(action_type, Read));
-
-    } while (Read.nextElement("Action"));
-
-      ActionsList.find(m)->second->print_parameters();
-
-  }
-
   template <class RepresentationPolicy>
   void GetActionSet(ActionSet<typename ImplementationPolicy::Field, RepresentationPolicy>& Aset){
     Aset.resize(multipliers.size());
@@ -271,6 +246,35 @@ class HMCResourceManager {
   }
 
 
+
+private:
+   // this private
+  template <class ReaderClass >
+  void fill_ActionsLevel(ReaderClass &Read){
+    // Actions set
+    int m;
+    Read.readDefault("multiplier",m);
+    multipliers.push_back(m);
+    std::cout << "Level : " << multipliers.size()  << " with multiplier : " << m << std::endl; 
+    // here gauge
+    Read.push("Action");
+    do{
+      auto &ActionFactory = HMC_LGTActionModuleFactory<gauge_string, ReaderClass>::getInstance(); 
+      std::string action_type;
+      Read.readDefault("name", action_type);
+      std::cout << "Registered types " << std::endl;
+      std::cout << ActionFactory.getBuilderList() << std::endl;  
+
+      ActionsList.emplace(m, ActionFactory.create(action_type, Read));
+
+    } while (Read.nextElement("Action"));
+
+      ActionsList.find(m)->second->print_parameters();
+
+  }
+
+
+
 };
 }
 }

lib/qcd/modules/ActionModules.h

@@ -0,0 +1,246 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid
+
+Source file: ./lib/qcd/modules/ActionModules.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 ACTION_MODULES_H
+#define ACTION_MODULES_H
+
+/*
+Define loadable, serializable modules
+for the HMC execution
+*/
+
+namespace Grid {
+
+//////////////////////////////////////////////
+//		Actions
+//////////////////////////////////////////////
+
+template <class ActionType, class APar>
+class ActionModule
+    : public Parametrized<APar>,
+      public HMCModuleBase<QCD::Action<typename ActionType::GaugeField> > {
+ public:
+  typedef HMCModuleBase< QCD::Action<typename ActionType::GaugeField> > Base;
+  typedef typename Base::Product Product;
+
+  std::unique_ptr<ActionType> ActionPtr;
+
+  ActionModule(APar Par) : Parametrized<APar>(Par) {}
+
+  template <class ReaderClass>
+  ActionModule(Reader<ReaderClass>& Reader) : Parametrized<APar>(Reader){};
+
+  virtual void print_parameters(){
+    std::cout << this->Par_ << std::endl;
+  }
+
+  Product* getPtr() {
+    if (!ActionPtr) initialize();
+
+    return ActionPtr.get();
+  }
+
+ private:
+  virtual void initialize() = 0;
+};
+
+
+////////////////
+// Modules
+////////////////
+
+namespace QCD{
+
+class PlaqPlusRectangleGaugeActionParameters : Serializable {
+ public:
+  GRID_SERIALIZABLE_CLASS_MEMBERS(PlaqPlusRectangleGaugeActionParameters, 
+    RealD, c_plaq,
+    RealD, c_rect);
+
+};
+
+class RBCGaugeActionParameters : Serializable {
+ public:
+  GRID_SERIALIZABLE_CLASS_MEMBERS(RBCGaugeActionParameters, 
+    RealD, beta,
+    RealD, c1);
+
+};
+
+class BetaGaugeActionParameters : Serializable {
+ public:
+  GRID_SERIALIZABLE_CLASS_MEMBERS(BetaGaugeActionParameters, 
+    RealD, beta);
+};
+
+
+
+
+template <class Impl >
+class WilsonGModule: public ActionModule<WilsonGaugeAction<Impl>, BetaGaugeActionParameters> {
+  typedef ActionModule<WilsonGaugeAction<Impl>, BetaGaugeActionParameters> ActionBase;
+  using ActionBase::ActionBase; // for constructors
+
+  // acquire resource
+  virtual void initialize(){
+    this->ActionPtr.reset(new WilsonGaugeAction<Impl>(this->Par_.beta));
+  }
+
+};
+
+template <class Impl >
+class PlaqPlusRectangleGModule: public ActionModule<PlaqPlusRectangleAction<Impl>, PlaqPlusRectangleGaugeActionParameters> {
+  typedef ActionModule<PlaqPlusRectangleAction<Impl>, PlaqPlusRectangleGaugeActionParameters> ActionBase;
+  using ActionBase::ActionBase; // for constructors
+
+  // acquire resource
+  virtual void initialize(){
+    this->ActionPtr.reset(new PlaqPlusRectangleAction<Impl>(this->Par_.c_plaq, this->Par_.c_rect));
+  }
+
+};
+
+template <class Impl >
+class RBCGModule: public ActionModule<RBCGaugeAction<Impl>, RBCGaugeActionParameters> {
+  typedef ActionModule<RBCGaugeAction<Impl>, RBCGaugeActionParameters> ActionBase;
+  using ActionBase::ActionBase; // for constructors
+
+  // acquire resource
+  virtual void initialize(){
+    this->ActionPtr.reset(new RBCGaugeAction<Impl>(this->Par_.beta, this->Par_.c1));
+  }
+
+};
+
+
+
+
+template <class Impl >
+class SymanzikGModule: public ActionModule<SymanzikGaugeAction<Impl>, BetaGaugeActionParameters> {
+  typedef ActionModule<SymanzikGaugeAction<Impl>, BetaGaugeActionParameters> ActionBase;
+  using ActionBase::ActionBase; // for constructors
+
+  // acquire resource
+  virtual void initialize(){
+    this->ActionPtr.reset(new SymanzikGaugeAction<Impl>(this->Par_.beta));
+  }
+
+};
+
+template <class Impl >
+class IwasakiGModule: public ActionModule<IwasakiGaugeAction<Impl>, BetaGaugeActionParameters> {
+  typedef ActionModule<IwasakiGaugeAction<Impl>, BetaGaugeActionParameters> ActionBase;
+  using ActionBase::ActionBase; // for constructors
+
+  // acquire resource
+  virtual void initialize(){
+    this->ActionPtr.reset(new IwasakiGaugeAction<Impl>(this->Par_.beta));
+  }
+
+};
+
+
+template <class Impl >
+class DBW2GModule: public ActionModule<DBW2GaugeAction<Impl>, BetaGaugeActionParameters> {
+  typedef ActionModule<DBW2GaugeAction<Impl>, BetaGaugeActionParameters> ActionBase;
+  using ActionBase::ActionBase; // for constructors
+
+  // acquire resource
+  virtual void initialize(){
+    this->ActionPtr.reset(new DBW2GaugeAction<Impl>(this->Par_.beta));
+  }
+
+};
+
+
+
+
+
+typedef WilsonGModule<PeriodicGimplR> WilsonGMod;
+typedef SymanzikGModule<PeriodicGimplR> SymanzikGMod;
+typedef IwasakiGModule<PeriodicGimplR> IwasakiGMod;
+typedef DBW2GModule<PeriodicGimplR> DBW2GMod;
+typedef RBCGModule<PeriodicGimplR> RBCGMod;
+typedef PlaqPlusRectangleGModule<PeriodicGimplR> PlaqPlusRectangleGMod;
+
+}// QCD temporarily here
+
+
+
+
+
+
+
+////////////////////////////////////////
+// Factories specialisations
+////////////////////////////////////////
+
+
+
+// use the same classed defined by Antonin, does not make sense to rewrite
+// Factory is perfectly fine
+// Registar must be changed because I do not want to use the ModuleFactory
+
+// explicit ref to LatticeGaugeField must be changed of put in the factory
+typedef HMCModuleBase< QCD::Action< QCD::LatticeGaugeField > > HMC_LGTActionModBase;
+
+template <char const *str, class ReaderClass >
+class HMC_LGTActionModuleFactory
+    : public Factory < HMC_LGTActionModBase ,	Reader<ReaderClass> > {
+ public:
+ 	typedef Reader<ReaderClass> TheReader; 
+ 	// use SINGLETON FUNCTOR MACRO HERE
+  HMC_LGTActionModuleFactory(const HMC_LGTActionModuleFactory& e) = delete;
+  void operator=(const HMC_LGTActionModuleFactory& e) = delete;
+  static HMC_LGTActionModuleFactory& getInstance(void) {
+    static HMC_LGTActionModuleFactory e;
+    return e;
+  }
+
+ private:
+  HMC_LGTActionModuleFactory(void) = default;
+    std::string obj_type() const {
+  	return std::string(str);
+  }
+};
+
+extern char gauge_string[];
+static Registrar<QCD::WilsonGMod,            HMC_LGTActionModuleFactory<gauge_string, XmlReader> > __WGmodXMLInit("Wilson"); 
+static Registrar<QCD::SymanzikGMod,          HMC_LGTActionModuleFactory<gauge_string, XmlReader> > __SymGmodXMLInit("Symanzik"); 
+static Registrar<QCD::IwasakiGMod,           HMC_LGTActionModuleFactory<gauge_string, XmlReader> > __IwGmodXMLInit("Iwasaki"); 
+static Registrar<QCD::DBW2GMod,              HMC_LGTActionModuleFactory<gauge_string, XmlReader> > __DBW2GmodXMLInit("DBW2"); 
+static Registrar<QCD::RBCGMod,               HMC_LGTActionModuleFactory<gauge_string, XmlReader> > __RBCGmodXMLInit("RBC"); 
+static Registrar<QCD::PlaqPlusRectangleGMod, HMC_LGTActionModuleFactory<gauge_string, XmlReader> > __PPRectGmodXMLInit("PlaqPlusRect"); 
+
+// add here the registration for other implementations and readers
+
+
+}
+
+
+#endif //HMC_MODULES_H

lib/qcd/modules/Modules.cc

@@ -34,4 +34,8 @@ char cp_string[]    = "CheckPointer";
 char hmc_string[]   = "HMC";
 char observable_string[] = "Observable";
 
+
+
+
+
 }

lib/qcd/modules/Modules.h

@@ -102,43 +102,6 @@ class HMCModuleBase {
   virtual void print_parameters(){};  // default to nothing
 };
 
-
-//////////////////////////////////////////////
-//		Actions
-//////////////////////////////////////////////
-
-template <class ActionType, class APar>
-class ActionModule
-    : public Parametrized<APar>,
-      public HMCModuleBase<QCD::Action<typename ActionType::GaugeField> > {
- public:
-  typedef HMCModuleBase< QCD::Action<typename ActionType::GaugeField> > Base;
-  typedef typename Base::Product Product;
-
-  std::unique_ptr<ActionType> ActionPtr;
-
-  ActionModule(APar Par) : Parametrized<APar>(Par) {}
-
-  template <class ReaderClass>
-  ActionModule(Reader<ReaderClass>& Reader) : Parametrized<APar>(Reader){};
-
-  virtual void print_parameters(){
-    std::cout << this->Par_ << std::endl;
-  }
-
-  Product* getPtr() {
-    if (!ActionPtr) initialize();
-
-    return ActionPtr.get();
-  }
-
- private:
-  virtual void initialize() = 0;
-};
-
-
-
-
 /////////////////////////////
 // Observables
 /////////////////////////////
@@ -174,48 +137,12 @@ class ObservableModule
 
 
 
-
-
-
-
-
 ////////////////
 // Modules
 ////////////////
 
 namespace QCD{
 
-class WilsonGaugeActionParameters : Serializable {
- public:
-  GRID_SERIALIZABLE_CLASS_MEMBERS(WilsonGaugeActionParameters, 
-    RealD, beta);
-};
-
-
-
-
-template <class Impl >
-class WilsonGModule: public ActionModule<WilsonGaugeAction<Impl>, WilsonGaugeActionParameters> {
-  typedef ActionModule<WilsonGaugeAction<Impl>, WilsonGaugeActionParameters> ActionBase;
-  using ActionBase::ActionBase; // for constructors
-
-  // acquire resource
-  virtual void initialize(){
-    this->ActionPtr.reset(new WilsonGaugeAction<Impl>(this->Par_.beta));
-  }
-
-};
-
-typedef WilsonGModule<PeriodicGimplR> WilsonGMod;
-
-
-
-
-
-
-
-
-
 //// Observables module
 class PlaquetteObsParameters : Serializable {
  public:
@@ -235,18 +162,6 @@ class PlaquetteMod: public ObservableModule<PlaquetteLogger<Impl>, PlaquetteObsP
 };
 
 
-
-
-
-
-
-
-
-
-
-
-
-
 }// QCD temporarily here
 
 
@@ -260,34 +175,6 @@ class PlaquetteMod: public ObservableModule<PlaquetteLogger<Impl>, PlaquetteObsP
 ////////////////////////////////////////
 
 
-
-// use the same classed defined by Antonin, does not make sense to rewrite
-// Factory is perfectly fine
-// Registar must be changed because I do not want to use the ModuleFactory
-
-// explicit ref to LatticeGaugeField must be changed of put in the factory
-typedef HMCModuleBase< QCD::Action< QCD::LatticeGaugeField > > HMC_LGTActionModBase;
-
-template <char const *str, class ReaderClass >
-class HMC_LGTActionModuleFactory
-    : public Factory < HMC_LGTActionModBase ,	Reader<ReaderClass> > {
- public:
- 	typedef Reader<ReaderClass> TheReader; 
- 	// use SINGLETON FUNCTOR MACRO HERE
-  HMC_LGTActionModuleFactory(const HMC_LGTActionModuleFactory& e) = delete;
-  void operator=(const HMC_LGTActionModuleFactory& e) = delete;
-  static HMC_LGTActionModuleFactory& getInstance(void) {
-    static HMC_LGTActionModuleFactory e;
-    return e;
-  }
-
- private:
-  HMC_LGTActionModuleFactory(void) = default;
-    std::string obj_type() const {
-  	return std::string(str);
-  }
-};
-
 // explicit ref to LatticeGaugeField must be changed of put in the factory
 typedef HMCModuleBase< QCD::HmcObservable<QCD::LatticeGaugeField> > HMC_ObsModBase;
 
@@ -313,11 +200,6 @@ class HMC_ObservablesModuleFactory
 
 
 
-
-
-
-
-
 template <class T, class TheFactory>
 class Registrar {
  public:
@@ -330,12 +212,6 @@ class Registrar {
 
 
 
-extern char gauge_string[];
-static Registrar<QCD::WilsonGMod, HMC_LGTActionModuleFactory<gauge_string, XmlReader> > __WGmodXMLInit("Wilson"); 
-// add here the registration for other implementations and readers
-
-
-
 extern char observable_string[];
 static Registrar<QCD::PlaquetteMod<QCD::PeriodicGimplR>, HMC_ObservablesModuleFactory<observable_string, XmlReader> > __OBSPLmodXMLInit("Plaquette"); 
 

lib/qcd/modules/mods.h

@@ -34,5 +34,7 @@ directory
 #include <Grid/qcd/modules/Factory.h>
 #include <Grid/qcd/modules/Modules.h>
 
+#include <Grid/qcd/modules/ActionModules.h>
+
 
 #endif //MODS_H

tests/hmc/Test_hmc_WilsonFermionGauge_Binary.cc

@@ -30,132 +30,110 @@ directory
 /*  END LEGAL */
 #include <Grid/Grid.h>
 
-using namespace std;
+int main(int argc, char **argv) {
-using namespace Grid;
+  using namespace Grid;
-using namespace Grid::QCD;
+  using namespace Grid::QCD;
 
-namespace Grid {
+  Grid_init(&argc, &argv);
-namespace QCD {
+  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;
 
-class HMCRunnerParameters : Serializable {
+   // Typedefs to simplify notation
- public:
+  typedef GenericHMCRunner<MinimumNorm2> HMCWrapper;  // Uses the default minimum norm
-  GRID_SERIALIZABLE_CLASS_MEMBERS(HMCRunnerParameters,
+  typedef WilsonImplR FermionImplPolicy;
-                                  double, beta,
+  typedef WilsonFermionR FermionAction;
-                                  double, mass,
+  typedef typename FermionAction::FermionField FermionField;
-                                  int, MaxCGIterations,
-                                  double, StoppingCondition,
-                                  bool, smearedAction,
-                                  int, SaveInterval,
-                                  std::string, format,
-                                  std::string, conf_prefix,
-                                  std::string, rng_prefix,
-                                  double, rho,
-                                  int, SmearingLevels,
-                                  );
 
-  HMCRunnerParameters() {}
-};
 
-// Derive from the BinaryHmcRunner (templated for gauge fields)
+  //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
-class HmcRunner : public BinaryHmcRunner {
+  HMCWrapper TheHMC;
- public:
-  void BuildTheAction(int argc, char **argv)
 
-  {
+  // Grid from the command line
-    typedef WilsonImplR ImplPolicy;
+  TheHMC.Resources.AddFourDimGrid("gauge");
-    typedef WilsonFermionR FermionAction;
+  // Possibile to create the module by hand 
-    typedef typename FermionAction::FermionField FermionField;
+  // hardcoding parameters or using a Reader
 
-    UGrid = SpaceTimeGrid::makeFourDimGrid(
-        GridDefaultLatt(), GridDefaultSimd(Nd, vComplex::Nsimd()),
-        GridDefaultMpi());
-    UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid);
 
-    FGrid = UGrid;
+  // Checkpointer definition
-    FrbGrid = UrbGrid;
+  CheckpointerParameters CPparams;  
+  CPparams.config_prefix = "ckpoint_lat";
+  CPparams.rng_prefix = "ckpoint_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);
+  
     // temporarily need a gauge field
-    LatticeGaugeField U(UGrid);
+  auto GridPtr = TheHMC.Resources.GetCartesian();
+  auto GridRBPtr = TheHMC.Resources.GetRBCartesian();
 
-    // Gauge action
+  LatticeGaugeField U(GridPtr);
-    WilsonGaugeActionR Waction(5.6);
 
-    Real mass = -0.77;
+  Real mass = -0.77;
 
-    // Can we define an overloaded operator that does not need U and initialises
+  // Can we define an overloaded operator that does not need U and initialises
-    // it with zeroes?
+  // it with zeroes?
-    FermionAction FermOp(U, *FGrid, *FrbGrid, mass);
+  FermionAction FermOp(U, *GridPtr, *GridRBPtr, mass);
 
-    ConjugateGradient<FermionField> CG(1.0e-8, 2000);
+  ConjugateGradient<FermionField> CG(1.0e-8, 2000);
 
-    TwoFlavourPseudoFermionAction<ImplPolicy> Nf2(FermOp, CG, CG);
+  TwoFlavourPseudoFermionAction<FermionImplPolicy> Nf2(FermOp, CG, CG);
 
     // Set smearing (true/false), default: false
-    Nf2.is_smeared = false;
+  Nf2.is_smeared = false;
+
 
     // Collect actions
-    ActionLevel<Field> Level1(1);
+  ActionLevel<HMCWrapper::Field> Level1(1);
-    Level1.push_back(&Nf2);
+  Level1.push_back(&Nf2);
 
-    ActionLevel<Field> Level2(4);
+  ActionLevel<HMCWrapper::Field> Level2(4);
-    Level2.push_back(&Waction);
+  Level2.push_back(&Waction);
 
-    TheAction.push_back(Level1);
+  TheHMC.TheAction.push_back(Level1);
-    TheAction.push_back(Level2);
+  TheHMC.TheAction.push_back(Level2);
+  /////////////////////////////////////////////////////////////
 
-    // Add observables
+  /*
-    int SaveInterval = 1;
-    std::string format = std::string("IEEE64BIG");
-    std::string conf_prefix = std::string("ckpoint_lat");
-    std::string rng_prefix = std::string("ckpoint_rng");
-    BinaryHmcCheckpointer<BinaryHmcRunner::ImplPolicy> Checkpoint(
-        conf_prefix, rng_prefix, SaveInterval, format);
-    // Can implement also a specific function in the hmcrunner
-    // AddCheckpoint (...) that takes the same parameters + a string/tag
-    // defining the type of the checkpointer
-    // with tags can be implemented by overloading and no ifs
-    // Then force all checkpoint to have few common functions
-    // return an object that is then passed to the Run function
-
-    PlaquetteLogger<BinaryHmcRunner::ImplPolicy> PlaqLog(
-        std::string("Plaquette"));
-    ObservablesList.push_back(&PlaqLog);
-    ObservablesList.push_back(&Checkpoint);
-
-    // Smearing section, omit if not needed
-    /*
     double rho = 0.1;  // smearing parameter
     int Nsmear = 2;    // number of smearing levels
-    Smear_Stout<BinaryHmcRunner::ImplPolicy> Stout(rho);
+    Smear_Stout<HMCWrapper::ImplPolicy> Stout(rho);
-    SmearedConfiguration<BinaryHmcRunner::ImplPolicy> SmearingPolicy(
+    SmearedConfiguration<HMCWrapper::ImplPolicy> SmearingPolicy(
         UGrid, Nsmear, Stout);
-        */
+  */
-    ///////////////////
 
-    // Run(argc, argv, Checkpoint, SmearingPolicy);
+  // HMC parameters are serialisable 
-    Run(argc, argv, Checkpoint);  // no smearing
+  TheHMC.Parameters.MD.MDsteps = 20;
-  };
+  TheHMC.Parameters.MD.trajL   = 1.0;
-};
-}
-}
 
-int main(int argc, char **argv) {
+  TheHMC.ReadCommandLine(argc, argv); // these can be parameters from file
-  Grid_init(&argc, &argv);
+  TheHMC.Run();  // no smearing
-
+  // TheHMC.Run(SmearingPolicy); // for smearing
-  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(10, 1.0);// MDsteps, traj length
-
-  TheHMC.BuildTheAction(argc, argv);
 
   Grid_finalize();
-}
+
+} // main
+
+

tests/hmc/Test_hmc_WilsonGauge_Binary.cc

@@ -78,7 +78,7 @@ int main(int argc, char **argv) {
   // standard
   RealD beta = 5.6 ;
   WilsonGaugeActionR Waction(beta);
-
+  
   ActionLevel<HMCWrapper::Field> Level1(1);
   Level1.push_back(&Waction);
   //Level1.push_back(WGMod.getPtr());
@@ -89,7 +89,7 @@ int main(int argc, char **argv) {
   TheHMC.Parameters.MD.MDsteps = 20;
   TheHMC.Parameters.MD.trajL   = 1.0;
 
-  TheHMC.ReadCommandLine(argc, argv); // these must be parameters from file
+  TheHMC.ReadCommandLine(argc, argv); // these can be parameters from file
   TheHMC.Run();  // no smearing
 
   Grid_finalize();