Hadrons: moving Hadrons to root directory, build system improvements

Grid/qcd/hmc/GenericHMCrunner.h

@@ -0,0 +1,219 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid
+
+Source file: ./lib/qcd/hmc/GenericHmcRunner.h
+
+Copyright (C) 2015
+
+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 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_GENERIC_HMC_RUNNER
+#define GRID_GENERIC_HMC_RUNNER
+
+#include <unordered_map>
+
+namespace Grid {
+namespace QCD {
+
+
+// very ugly here but possibly resolved if we had a base Reader class
+template < class ReaderClass >
+class HMCRunnerBase {
+public:
+  virtual void Run() = 0;
+  virtual void initialize(ReaderClass& ) = 0;
+};
+
+
+template <class Implementation,
+          template <typename, typename, typename> class Integrator,
+          class RepresentationsPolicy = NoHirep, class ReaderClass = XmlReader>
+class HMCWrapperTemplate: public HMCRunnerBase<ReaderClass> {
+ public:
+  INHERIT_FIELD_TYPES(Implementation);
+  typedef Implementation ImplPolicy;  // visible from outside
+  template <typename S = NoSmearing<Implementation> >
+  using IntegratorType = Integrator<Implementation, S, RepresentationsPolicy>;
+
+  HMCparameters Parameters;
+  std::string ParameterFile;
+  HMCResourceManager<Implementation> Resources;
+
+  // The set of actions (keep here for lower level users, for now)
+  ActionSet<Field, RepresentationsPolicy> TheAction;
+
+  HMCWrapperTemplate() = default;
+
+  HMCWrapperTemplate(HMCparameters Par){
+    Parameters = Par;
+  }
+
+  void initialize(ReaderClass & TheReader){
+    std::cout  << "Initialization of the HMC" << std::endl;
+    Resources.initialize(TheReader);
+
+    // eventually add smearing
+
+    Resources.GetActionSet(TheAction);    
+  }
+
+
+  void ReadCommandLine(int argc, char **argv) {
+    std::string arg;
+
+    if (GridCmdOptionExists(argv, argv + argc, "--StartingType")) {
+      arg = GridCmdOptionPayload(argv, argv + argc, "--StartingType");
+
+      if (arg != "HotStart" && arg != "ColdStart" && arg != "TepidStart" &&
+          arg != "CheckpointStart") {
+        std::cout << GridLogError << "Unrecognized option in --StartingType\n";
+        std::cout
+            << GridLogError
+            << "Valid [HotStart, ColdStart, TepidStart, CheckpointStart]\n";
+        exit(1);
+      }
+      Parameters.StartingType = arg;
+    }
+
+    if (GridCmdOptionExists(argv, argv + argc, "--StartingTrajectory")) {
+      arg = GridCmdOptionPayload(argv, argv + argc, "--StartingTrajectory");
+      std::vector<int> ivec(0);
+      GridCmdOptionIntVector(arg, ivec);
+      Parameters.StartTrajectory = ivec[0];
+    }
+
+    if (GridCmdOptionExists(argv, argv + argc, "--Trajectories")) {
+      arg = GridCmdOptionPayload(argv, argv + argc, "--Trajectories");
+      std::vector<int> ivec(0);
+      GridCmdOptionIntVector(arg, ivec);
+      Parameters.Trajectories = ivec[0];
+    }
+
+    if (GridCmdOptionExists(argv, argv + argc, "--Thermalizations")) {
+      arg = GridCmdOptionPayload(argv, argv + argc, "--Thermalizations");
+      std::vector<int> ivec(0);
+      GridCmdOptionIntVector(arg, ivec);
+      Parameters.NoMetropolisUntil = ivec[0];
+    }
+    if (GridCmdOptionExists(argv, argv + argc, "--ParameterFile")) {
+      arg = GridCmdOptionPayload(argv, argv + argc, "--ParameterFile");
+      ParameterFile = arg;
+    }
+  }
+
+
+  template <class SmearingPolicy>
+  void Run(SmearingPolicy &S) {
+    Runner(S);
+  }
+
+  void Run(){
+    NoSmearing<Implementation> S;
+    Runner(S);
+  }
+
+  //////////////////////////////////////////////////////////////////
+
+ private:
+  template <class SmearingPolicy>
+  void Runner(SmearingPolicy &Smearing) {
+    auto UGrid = Resources.GetCartesian();
+    Resources.AddRNGs();
+    Field U(UGrid);
+
+    // Can move this outside?
+    typedef IntegratorType<SmearingPolicy> TheIntegrator;
+    TheIntegrator MDynamics(UGrid, Parameters.MD, TheAction, Smearing);
+
+    if (Parameters.StartingType == "HotStart") {
+      // Hot start
+      Resources.SeedFixedIntegers();
+      Implementation::HotConfiguration(Resources.GetParallelRNG(), U);
+    } else if (Parameters.StartingType == "ColdStart") {
+      // Cold start
+      Resources.SeedFixedIntegers();
+      Implementation::ColdConfiguration(Resources.GetParallelRNG(), U);
+    } else if (Parameters.StartingType == "TepidStart") {
+      // Tepid start
+      Resources.SeedFixedIntegers();
+      Implementation::TepidConfiguration(Resources.GetParallelRNG(), U);
+    } else if (Parameters.StartingType == "CheckpointStart") {
+      // CheckpointRestart
+      Resources.GetCheckPointer()->CheckpointRestore(Parameters.StartTrajectory, U,
+                                   Resources.GetSerialRNG(),
+                                   Resources.GetParallelRNG());
+    }
+
+    Smearing.set_Field(U);
+
+    HybridMonteCarlo<TheIntegrator> HMC(Parameters, MDynamics,
+                                        Resources.GetSerialRNG(),
+                                        Resources.GetParallelRNG(), 
+                                        Resources.GetObservables(), U);
+
+    // Run it
+    HMC.evolve();
+  }
+};
+
+// These are for gauge fields, default integrator MinimumNorm2
+template <template <typename, typename, typename> class Integrator>
+using GenericHMCRunner = HMCWrapperTemplate<PeriodicGimplR, Integrator>;
+template <template <typename, typename, typename> class Integrator>
+using GenericHMCRunnerF = HMCWrapperTemplate<PeriodicGimplF, Integrator>;
+template <template <typename, typename, typename> class Integrator>
+using GenericHMCRunnerD = HMCWrapperTemplate<PeriodicGimplD, Integrator>;
+
+
+// These are for gauge fields, default integrator MinimumNorm2
+template <template <typename, typename, typename> class Integrator>
+using ConjugateHMCRunner = HMCWrapperTemplate<ConjugateGimplR, Integrator>;
+template <template <typename, typename, typename> class Integrator>
+using ConjugateHMCRunnerF = HMCWrapperTemplate<ConjugateGimplF, Integrator>;
+template <template <typename, typename, typename> class Integrator>
+using ConjugateHMCRunnerD = HMCWrapperTemplate<ConjugateGimplD, Integrator>;
+
+
+
+template <class RepresentationsPolicy,
+          template <typename, typename, typename> class Integrator>
+using GenericHMCRunnerHirep =
+    HMCWrapperTemplate<PeriodicGimplR, Integrator, RepresentationsPolicy>;
+
+template <class Implementation, class RepresentationsPolicy, 
+          template <typename, typename, typename> class Integrator>
+using GenericHMCRunnerTemplate = HMCWrapperTemplate<Implementation, Integrator, RepresentationsPolicy>;
+
+typedef HMCWrapperTemplate<ScalarImplR, MinimumNorm2, ScalarFields>
+    ScalarGenericHMCRunner;
+
+typedef HMCWrapperTemplate<ScalarAdjImplR, MinimumNorm2, ScalarMatrixFields>
+    ScalarAdjGenericHMCRunner;
+
+template <int Colours> 
+using ScalarNxNAdjGenericHMCRunner = HMCWrapperTemplate < ScalarNxNAdjImplR<Colours>, ForceGradient, ScalarNxNMatrixFields<Colours> >;
+
+}  // namespace QCD
+}  // namespace Grid
+
+#endif  // GRID_GENERIC_HMC_RUNNER

Grid/qcd/hmc/HMC.h

@@ -0,0 +1,251 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid
+
+Source file: ./lib/qcd/hmc/HMC.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 */
+//--------------------------------------------------------------------
+/*! @file HMC.h
+ * @brief Classes for Hybrid Monte Carlo update
+ *
+ * @author Guido Cossu
+ */
+//--------------------------------------------------------------------
+#ifndef HMC_INCLUDED
+#define HMC_INCLUDED
+
+#include <string>
+#include <list>
+
+
+
+#include <Grid/qcd/hmc/integrators/Integrator.h>
+#include <Grid/qcd/hmc/integrators/Integrator_algorithm.h>
+
+namespace Grid {
+namespace QCD {
+
+struct HMCparameters: Serializable {
+	GRID_SERIALIZABLE_CLASS_MEMBERS(HMCparameters,
+                                  Integer, StartTrajectory,
+                                  Integer, Trajectories, /* @brief Number of sweeps in this run */
+                                  bool, MetropolisTest,
+                                  Integer, NoMetropolisUntil,
+                                  std::string, StartingType,
+                                  IntegratorParameters, MD)
+
+  HMCparameters() {
+    ////////////////////////////// Default values
+    MetropolisTest    = true;
+    NoMetropolisUntil = 10;
+    StartTrajectory   = 0;
+    Trajectories      = 10;
+    StartingType      = "HotStart";
+    /////////////////////////////////
+  }
+
+  template <class ReaderClass >
+  HMCparameters(Reader<ReaderClass> & TheReader){
+  	initialize(TheReader);
+  }
+
+  template < class ReaderClass > 
+  void initialize(Reader<ReaderClass> &TheReader){
+  	std::cout << GridLogMessage << "Reading HMC\n";
+  	read(TheReader, "HMC", *this);
+  }
+
+
+  void print_parameters() const {
+    std::cout << GridLogMessage << "[HMC parameters] Trajectories            : " << Trajectories << "\n";
+    std::cout << GridLogMessage << "[HMC parameters] Start trajectory        : " << StartTrajectory << "\n";
+    std::cout << GridLogMessage << "[HMC parameters] Metropolis test (on/off): " << std::boolalpha << MetropolisTest << "\n";
+    std::cout << GridLogMessage << "[HMC parameters] Thermalization trajs    : " << NoMetropolisUntil << "\n";
+    std::cout << GridLogMessage << "[HMC parameters] Starting type           : " << StartingType << "\n";
+    MD.print_parameters();
+  }
+  
+};
+	
+template <class IntegratorType>
+class HybridMonteCarlo {
+ private:
+  const HMCparameters Params;
+
+  typedef typename IntegratorType::Field Field;
+  typedef std::vector< HmcObservable<Field> * > ObsListType;
+  
+  	//pass these from the resource manager
+  GridSerialRNG &sRNG;   
+  GridParallelRNG &pRNG; 
+
+  Field &Ucur;
+  
+  IntegratorType &TheIntegrator;
+	ObsListType Observables;
+
+  /////////////////////////////////////////////////////////
+  // Metropolis step
+  /////////////////////////////////////////////////////////
+  bool metropolis_test(const RealD DeltaH) {
+    RealD rn_test;
+
+    RealD prob = std::exp(-DeltaH);
+
+    random(sRNG, rn_test);
+
+    std::cout << GridLogMessage
+              << "--------------------------------------------------\n";
+    std::cout << GridLogMessage << "exp(-dH) = " << prob
+              << "  Random = " << rn_test << "\n";
+    std::cout << GridLogMessage
+              << "Acc. Probability = " << ((prob < 1.0) ? prob : 1.0) << "\n";
+
+    if ((prob > 1.0) || (rn_test <= prob)) {  // accepted
+      std::cout << GridLogMessage << "Metropolis_test -- ACCEPTED\n";
+      std::cout << GridLogMessage
+                << "--------------------------------------------------\n";
+      return true;
+    } else {  // rejected
+      std::cout << GridLogMessage << "Metropolis_test -- REJECTED\n";
+      std::cout << GridLogMessage
+                << "--------------------------------------------------\n";
+      return false;
+    }
+  }
+
+  /////////////////////////////////////////////////////////
+  // Evolution
+  /////////////////////////////////////////////////////////
+  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
+
+    std::streamsize current_precision = std::cout.precision();
+    std::cout.precision(15);
+    std::cout << GridLogMessage << "Total H before trajectory = " << H0 << "\n";
+    std::cout.precision(current_precision);
+
+    TheIntegrator.integrate(U);
+
+    RealD H1 = TheIntegrator.S(U);  // updated state action
+
+    ///////////////////////////////////////////////////////////
+    if(0){
+      std::cout << "------------------------- Reversibility test" << std::endl;
+      TheIntegrator.reverse_momenta();
+      TheIntegrator.integrate(U);
+
+      H1 = TheIntegrator.S(U);  // updated state action
+      std::cout << "--------------------------------------------" << std::endl;
+    }
+    ///////////////////////////////////////////////////////////
+
+
+    std::cout.precision(15);
+    std::cout << GridLogMessage << "Total H after trajectory  = " << H1
+	      << "  dH = " << H1 - H0 << "\n";
+    std::cout.precision(current_precision);
+    
+    return (H1 - H0);
+  }
+  
+
+  
+
+ public:
+  /////////////////////////////////////////
+  // Constructor
+  /////////////////////////////////////////
+  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) {}
+  ~HybridMonteCarlo(){};
+
+  void evolve(void) {
+    Real DeltaH;
+
+    Field Ucopy(Ucur._grid);
+
+    Params.print_parameters();
+    TheIntegrator.print_actions();
+
+    // Actual updates (evolve a copy Ucopy then copy back eventually)
+    unsigned int FinalTrajectory = Params.Trajectories + Params.NoMetropolisUntil + Params.StartTrajectory;
+    for (int traj = Params.StartTrajectory; traj < FinalTrajectory; ++traj) {
+      std::cout << GridLogMessage << "-- # Trajectory = " << traj << "\n";
+      if (traj < Params.StartTrajectory + Params.NoMetropolisUntil) {
+      	std::cout << GridLogMessage << "-- Thermalization" << std::endl;
+      }
+      
+      double t0=usecond();
+      Ucopy = Ucur;
+
+      DeltaH = evolve_hmc_step(Ucopy);
+      // Metropolis-Hastings test
+      bool accept = true;
+      if (traj >= Params.StartTrajectory + Params.NoMetropolisUntil) {
+        accept = metropolis_test(DeltaH);
+      } else {
+      	std::cout << GridLogMessage << "Skipping Metropolis test" << std::endl;
+      }
+
+      if (accept)
+        Ucur = Ucopy; 
+      
+     
+      
+      double t1=usecond();
+      std::cout << GridLogMessage << "Total time for trajectory (s): " << (t1-t0)/1e6 << std::endl;
+
+
+      for (int obs = 0; obs < Observables.size(); obs++) {
+      	std::cout << GridLogDebug << "Observables # " << obs << std::endl;
+      	std::cout << GridLogDebug << "Observables total " << Observables.size() << std::endl;
+      	std::cout << GridLogDebug << "Observables pointer " << Observables[obs] << std::endl;
+        Observables[obs]->TrajectoryComplete(traj + 1, Ucur, sRNG, pRNG);
+      }
+      std::cout << GridLogMessage << ":::::::::::::::::::::::::::::::::::::::::::" << std::endl;
+    }
+  }
+
+};
+
+
+}  // QCD
+}  // Grid
+
+
+// april 11 2017 merge, Guido, commenting out
+//#include <Grid/parallelIO/NerscIO.h>
+//#include <Grid/qcd/hmc/NerscCheckpointer.h>
+//#include <Grid/qcd/hmc/HmcRunner.h>
+
+#endif

Grid/qcd/hmc/HMCModules.h

@@ -0,0 +1,111 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid
+
+Source file: ./lib/qcd/hmc/GenericHmcRunner.h
+
+Copyright (C) 2015
+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 GRID_HMC_MODULES
+#define GRID_HMC_MODULES
+
+
+#include "HMC_GridModules.h"
+
+namespace Grid {
+namespace QCD {
+
+////////////////////////////////////////////////////////////////////
+struct RNGModuleParameters: Serializable {
+  GRID_SERIALIZABLE_CLASS_MEMBERS(RNGModuleParameters,
+  std::string, serial_seeds,
+  std::string, parallel_seeds,);
+
+  std::vector<int> getSerialSeeds(){return strToVec<int>(serial_seeds);}
+  std::vector<int> getParallelSeeds(){return strToVec<int>(parallel_seeds);}
+
+  RNGModuleParameters(): serial_seeds("1"), parallel_seeds("1"){}
+
+  template <class ReaderClass >
+  RNGModuleParameters(Reader<ReaderClass>& Reader){
+    read(Reader, "RandomNumberGenerator", *this); 
+  }
+  
+};
+
+// Random number generators module
+class RNGModule{
+   GridSerialRNG sRNG_;
+   std::unique_ptr<GridParallelRNG> pRNG_;
+   RNGModuleParameters Params_;
+
+public:
+
+  RNGModule(){};
+
+  void set_pRNG(GridParallelRNG* pRNG){
+    pRNG_.reset(pRNG);
+  }
+
+  void set_RNGSeeds(RNGModuleParameters& Params) {
+    Params_ = Params;
+  }
+
+  GridSerialRNG& get_sRNG() { return sRNG_; }
+  GridParallelRNG& get_pRNG() { return *pRNG_.get(); }
+
+  void seed() {
+    auto SerialSeeds   = Params_.getSerialSeeds();
+    auto ParallelSeeds = Params_.getParallelSeeds();
+    if (SerialSeeds.size() == 0 && ParallelSeeds.size() == 0) {
+      std::cout << GridLogError << "Seeds not initialized" << std::endl;
+      exit(1);
+    }
+    sRNG_.SeedFixedIntegers(SerialSeeds);
+    pRNG_->SeedFixedIntegers(ParallelSeeds);
+  }
+};
+
+
+/*
+///////////////////////////////////////////////////////////////////
+/// Smearing module
+template <class ImplementationPolicy>
+class SmearingModule{
+   virtual void get_smearing();
+};
+
+template <class ImplementationPolicy>
+class StoutSmearingModule: public SmearingModule<ImplementationPolicy>{
+   SmearedConfiguration<ImplementationPolicy> SmearingPolicy;
+};
+
+*/
+
+
+
+}  // namespace QCD
+}  // namespace Grid
+
+#endif  // GRID_HMC_MODULES

Grid/qcd/hmc/HMCResourceManager.h

@@ -0,0 +1,328 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid
+
+Source file: ./lib/qcd/hmc/GenericHmcRunner.h
+
+Copyright (C) 2015
+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 HMC_RESOURCE_MANAGER_H
+#define HMC_RESOURCE_MANAGER_H
+
+#include <unordered_map>
+
+// One function per Checkpointer, use a macro to simplify
+#define RegisterLoadCheckPointerFunction(NAME)                           \
+  void Load##NAME##Checkpointer(const CheckpointerParameters& Params_) { \
+    if (!have_CheckPointer) {                                            \
+      std::cout << GridLogDebug << "Loading Checkpointer " << #NAME      \
+                << std::endl;                                            \
+      CP = std::unique_ptr<CheckpointerBaseModule>(                      \
+        new NAME##CPModule<ImplementationPolicy>(Params_));              \
+      have_CheckPointer = true;                                          \
+    } else {                                                             \
+      std::cout << GridLogError << "Checkpointer already loaded "        \
+                << std::endl;                                            \
+      exit(1);                                                           \
+    }                                                                    \
+  }
+
+#define RegisterLoadCheckPointerMetadataFunction(NAME)                   \
+  template < class Metadata >                                            \
+  void Load##NAME##Checkpointer(const CheckpointerParameters& Params_, const Metadata& M_) { \
+    if (!have_CheckPointer) {                                            \
+      std::cout << GridLogDebug << "Loading Metadata Checkpointer " << #NAME      \
+                << std::endl;                                            \
+      CP = std::unique_ptr<CheckpointerBaseModule>(                      \
+        new NAME##CPModule<ImplementationPolicy, Metadata >(Params_, M_));   \
+      have_CheckPointer = true;                                          \
+    } else {                                                             \
+      std::cout << GridLogError << "Checkpointer already loaded "        \
+                << std::endl;                                            \
+      exit(1);                                                           \
+    }                                                                    \
+  }
+
+namespace Grid {
+namespace QCD {
+
+// HMC Resource manager
+template <class ImplementationPolicy>
+class HMCResourceManager {
+  typedef HMCModuleBase< QCD::BaseHmcCheckpointer<ImplementationPolicy> > CheckpointerBaseModule;
+  typedef HMCModuleBase< QCD::HmcObservable<typename ImplementationPolicy::Field> > ObservableBaseModule;
+  typedef ActionModuleBase< QCD::Action<typename ImplementationPolicy::Field>, GridModule > ActionBaseModule;
+
+  // Named storage for grid pairs (std + red-black)
+  std::unordered_map<std::string, GridModule> Grids;
+  RNGModule RNGs;
+
+  // SmearingModule<ImplementationPolicy> Smearing;
+  std::unique_ptr<CheckpointerBaseModule> CP;
+
+  // A vector of HmcObservable modules
+  std::vector<std::unique_ptr<ObservableBaseModule> > ObservablesList;
+
+
+  // A vector of HmcObservable modules
+  std::multimap<int, std::unique_ptr<ActionBaseModule> > ActionsList;
+  std::vector<int> multipliers;
+
+  bool have_RNG;
+  bool have_CheckPointer;
+
+  // NOTE: operator << is not overloaded for std::vector<string> 
+  // so this function is necessary
+  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) {}
+
+  template <class ReaderClass, class vector_type = vComplex >
+  void initialize(ReaderClass &Read){
+    // assumes we are starting from the main node
+
+    // Geometry
+    GridModuleParameters GridPar(Read);
+    GridFourDimModule<vector_type> GridMod( GridPar) ;
+    AddGrid("gauge", GridMod);
+
+    // Checkpointer
+    auto &CPfactory = HMC_CPModuleFactory<cp_string, ImplementationPolicy, ReaderClass >::getInstance();
+    Read.push("Checkpointer");
+    std::string cp_type;
+    read(Read,"name", cp_type);
+    std::cout << "Registered types " << std::endl;
+    output_vector_string(CPfactory.getBuilderList());
+
+
+    CP = CPfactory.create(cp_type, Read);
+    CP->print_parameters();
+    Read.pop();    
+    have_CheckPointer = true;  
+
+    RNGModuleParameters RNGpar(Read);
+    SetRNGSeeds(RNGpar);
+
+    // Observables
+    auto &ObsFactory = HMC_ObservablesModuleFactory<observable_string, typename ImplementationPolicy::Field, ReaderClass>::getInstance(); 
+    Read.push(observable_string);// here must check if existing...
+    do {
+      std::string obs_type;
+      read(Read,"name", obs_type);
+      std::cout << "Registered types " << std::endl;
+      output_vector_string(ObsFactory.getBuilderList() );
+
+      ObservablesList.emplace_back(ObsFactory.create(obs_type, Read));
+      ObservablesList[ObservablesList.size() - 1]->print_parameters();
+    } while (Read.nextElement(observable_string));
+    Read.pop();
+
+    // Loop on levels
+    if(!Read.push("Actions")){
+      std::cout << "Actions not found" << std::endl; 
+      exit(1);
+    }
+
+    if(!Read.push("Level")){// push must check if the node exist
+         std::cout << "Level not found" << std::endl; 
+      exit(1);
+    }
+    do
+    {
+      fill_ActionsLevel(Read); 
+    }
+    while(Read.push("Level"));
+
+    Read.pop();
+  }
+
+
+ 
+  template <class RepresentationPolicy>
+  void GetActionSet(ActionSet<typename ImplementationPolicy::Field, RepresentationPolicy>& Aset){
+    Aset.resize(multipliers.size());
+ 
+    for(auto it = ActionsList.begin(); it != ActionsList.end(); it++){
+      (*it).second->acquireResource(Grids["gauge"]);
+      Aset[(*it).first-1].push_back((*it).second->getPtr());
+    }
+  }
+
+
+
+  //////////////////////////////////////////////////////////////
+  // Grids
+  //////////////////////////////////////////////////////////////
+
+  void AddGrid(const std::string s, GridModule& M) {
+    // Check for name clashes
+    auto search = Grids.find(s);
+    if (search != Grids.end()) {
+      std::cout << GridLogError << "Grid with name \"" << search->first
+                << "\" already present. Terminating\n";
+      exit(1);
+    }
+    Grids[s] = std::move(M);
+    std::cout << GridLogMessage << "::::::::::::::::::::::::::::::::::::::::" <<std::endl;
+    std::cout << GridLogMessage << "HMCResourceManager:" << std::endl;
+    std::cout << GridLogMessage << "Created grid set with name '" << s << "' and decomposition for the full cartesian " << std::endl;
+    Grids[s].show_full_decomposition();
+    std::cout << GridLogMessage << "::::::::::::::::::::::::::::::::::::::::" <<std::endl;
+  }
+
+  // Add a named grid set, 4d shortcut
+  void AddFourDimGrid(const std::string s) {
+    GridFourDimModule<vComplex> Mod;
+    AddGrid(s, Mod);
+  }
+
+  // Add a named grid set, 4d shortcut + tweak simd lanes
+  void AddFourDimGrid(const std::string s, const std::vector<int> simd_decomposition) {
+    GridFourDimModule<vComplex> Mod(simd_decomposition);
+    AddGrid(s, Mod);
+  }
+
+
+  GridCartesian* GetCartesian(std::string s = "") {
+    if (s.empty()) s = Grids.begin()->first;
+    std::cout << GridLogDebug << "Getting cartesian grid from: " << s
+              << std::endl;
+    return Grids[s].get_full();
+  }
+
+  GridRedBlackCartesian* GetRBCartesian(std::string s = "") {
+    if (s.empty()) s = Grids.begin()->first;
+    std::cout << GridLogDebug << "Getting rb-cartesian grid from: " << s
+              << std::endl;
+    return Grids[s].get_rb();
+  }
+
+  //////////////////////////////////////////////////////
+  // Random number generators
+  //////////////////////////////////////////////////////
+
+  void AddRNGs(std::string s = "") {
+    // Couple the RNGs to the GridModule tagged by s
+    // the default is the first grid registered
+    assert(Grids.size() > 0 && !have_RNG);
+    if (s.empty()) s = Grids.begin()->first;
+    std::cout << GridLogDebug << "Adding RNG to grid: " << s << std::endl;
+    RNGs.set_pRNG(new GridParallelRNG(GetCartesian(s)));
+    have_RNG = true;
+  }
+
+  void SetRNGSeeds(RNGModuleParameters& Params) { RNGs.set_RNGSeeds(Params); }
+
+  GridSerialRNG& GetSerialRNG() { return RNGs.get_sRNG(); }
+
+  GridParallelRNG& GetParallelRNG() {
+    assert(have_RNG);
+    return RNGs.get_pRNG();
+  }
+
+  void SeedFixedIntegers() {
+    assert(have_RNG);
+    RNGs.seed();
+  }
+
+  //////////////////////////////////////////////////////
+  // Checkpointers
+  //////////////////////////////////////////////////////
+
+  BaseHmcCheckpointer<ImplementationPolicy>* GetCheckPointer() {
+    if (have_CheckPointer)
+      return CP->getPtr();
+    else {
+      std::cout << GridLogError << "Error: no checkpointer defined"
+                << std::endl;
+      exit(1);
+    }
+  }
+
+  RegisterLoadCheckPointerFunction(Binary);
+  RegisterLoadCheckPointerFunction(Nersc);
+  #ifdef HAVE_LIME
+  RegisterLoadCheckPointerFunction(ILDG);
+  RegisterLoadCheckPointerMetadataFunction(Scidac);
+  #endif
+
+  ////////////////////////////////////////////////////////
+  // Observables
+  ////////////////////////////////////////////////////////
+
+  template<class T, class... Types>
+  void AddObservable(Types&&... Args){
+    ObservablesList.push_back(std::unique_ptr<T>(new T(std::forward<Types>(Args)...)));
+    ObservablesList.back()->print_parameters();
+  }
+
+  std::vector<HmcObservable<typename ImplementationPolicy::Field>* > GetObservables(){
+    std::vector<HmcObservable<typename ImplementationPolicy::Field>* > out;
+    for (auto &i : ObservablesList){
+      out.push_back(i->getPtr());
+    }
+
+    // Add the checkpointer to the observables
+    out.push_back(GetCheckPointer());
+    return out;
+  }
+
+
+
+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_ActionModuleFactory<gauge_string, typename ImplementationPolicy::Field, ReaderClass>::getInstance(); 
+      std::string action_type;
+      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();    
+    Read.pop();
+
+  }
+
+
+
+};
+}
+}
+
+#endif  // HMC_RESOURCE_MANAGER_H

Grid/qcd/hmc/HMCRunnerModule.h

@@ -0,0 +1,137 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid
+
+Source file: ./lib/qcd/hmc/GenericHmcRunner.h
+
+Copyright (C) 2015
+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 HMC_RUNNER_MODULE
+#define HMC_RUNNER_MODULE
+
+namespace Grid {
+
+// the reader class is necessary here for the automatic initialization of the resources
+// if we had a virtual reader would have been unecessary
+template <class HMCType, class ReaderClass >
+class HMCModule
+    : public Parametrized< QCD::HMCparameters >,
+      public HMCModuleBase< QCD::HMCRunnerBase<ReaderClass> > {
+ public:
+  typedef HMCModuleBase< QCD::HMCRunnerBase<ReaderClass> > Base;
+  typedef typename Base::Product Product;
+
+  std::unique_ptr<HMCType> HMCPtr;
+
+  HMCModule(QCD::HMCparameters Par) : Parametrized<QCD::HMCparameters>(Par) {}
+
+  template <class ReaderCl>
+  HMCModule(Reader<ReaderCl>& R) : Parametrized<QCD::HMCparameters>(R, "HMC"){};
+
+  Product* getPtr() {
+    if (!HMCPtr) initialize();
+ 
+    return HMCPtr.get();
+  }
+
+ private:
+  virtual void initialize() = 0;
+};
+
+// Factory
+template <char const *str, class ReaderClass >
+class HMCRunnerModuleFactory
+    : public Factory < HMCModuleBase< QCD::HMCRunnerBase<ReaderClass> > ,	Reader<ReaderClass> > {
+ public:
+ 	typedef Reader<ReaderClass> TheReader; 
+ 	// use SINGLETON FUNCTOR MACRO HERE
+  HMCRunnerModuleFactory(const HMCRunnerModuleFactory& e) = delete;
+  void operator=(const HMCRunnerModuleFactory& e) = delete;
+  static HMCRunnerModuleFactory& getInstance(void) {
+    static HMCRunnerModuleFactory e;
+    return e;
+  }
+
+ private:
+  HMCRunnerModuleFactory(void) = default;
+  std::string obj_type() const {
+  	return std::string(str);
+  }
+};
+
+
+
+
+
+///////////////
+// macro for these
+
+template < class ImplementationPolicy, class RepresentationPolicy, class ReaderClass >
+class HMCLeapFrog: public HMCModule< QCD::GenericHMCRunnerTemplate<ImplementationPolicy, RepresentationPolicy, QCD::LeapFrog>, ReaderClass >{
+  typedef HMCModule< QCD::GenericHMCRunnerTemplate<ImplementationPolicy, RepresentationPolicy, QCD::LeapFrog>, ReaderClass  > HMCBaseMod;
+  using HMCBaseMod::HMCBaseMod;
+
+  // aquire resource
+  virtual void initialize(){
+    this->HMCPtr.reset(new QCD::GenericHMCRunnerTemplate<ImplementationPolicy, RepresentationPolicy, QCD::LeapFrog>(this->Par_) );
+  }
+};
+
+template < class ImplementationPolicy, class RepresentationPolicy, class ReaderClass >
+class HMCMinimumNorm2: public HMCModule< QCD::GenericHMCRunnerTemplate<ImplementationPolicy, RepresentationPolicy, QCD::MinimumNorm2>, ReaderClass  >{
+  typedef HMCModule< QCD::GenericHMCRunnerTemplate<ImplementationPolicy, RepresentationPolicy, QCD::MinimumNorm2>, ReaderClass  > HMCBaseMod;
+  using HMCBaseMod::HMCBaseMod;
+
+  // aquire resource
+  virtual void initialize(){
+    this->HMCPtr.reset(new QCD::GenericHMCRunnerTemplate<ImplementationPolicy, RepresentationPolicy, QCD::MinimumNorm2>(this->Par_));
+  }
+};
+
+
+template < class ImplementationPolicy, class RepresentationPolicy, class ReaderClass >
+class HMCForceGradient: public HMCModule< QCD::GenericHMCRunnerTemplate<ImplementationPolicy, RepresentationPolicy, QCD::ForceGradient>, ReaderClass  >{
+  typedef HMCModule< QCD::GenericHMCRunnerTemplate<ImplementationPolicy, RepresentationPolicy, QCD::ForceGradient>, ReaderClass   > HMCBaseMod;
+  using HMCBaseMod::HMCBaseMod;
+
+  // aquire resource
+  virtual void initialize(){
+    this->HMCPtr.reset(new QCD::GenericHMCRunnerTemplate<ImplementationPolicy, RepresentationPolicy, QCD::ForceGradient>(this->Par_) );
+  }
+};
+
+extern char hmc_string[];
+
+
+
+
+
+
+//////////////////////////////////////////////////////////////
+
+
+
+}
+
+#endif

Grid/qcd/hmc/HMC_GridModules.h

@@ -0,0 +1,177 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid
+
+Source file: ./lib/qcd/hmc/HMC_GridModules.h
+
+Copyright (C) 2015
+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 HMC_GRID_MODULES
+#define HMC_GRID_MODULES
+
+namespace Grid {
+
+// Resources
+// Modules for grids
+
+// Introduce another namespace HMCModules?
+
+class GridModuleParameters: Serializable{
+public:
+  GRID_SERIALIZABLE_CLASS_MEMBERS(GridModuleParameters,
+  std::string, lattice,
+  std::string, mpi);
+
+  std::vector<int> getLattice() const {return strToVec<int>(lattice);}
+  std::vector<int> getMpi()     const {return strToVec<int>(mpi);}
+
+
+  void check() const {
+    if (getLattice().size() != getMpi().size() ) {
+      std::cout << GridLogError
+                << "Error in GridModuleParameters: lattice and mpi dimensions "
+                   "do not match"
+                << std::endl;
+      exit(1);
+    }
+  }
+
+  template <class ReaderClass>
+  GridModuleParameters(Reader<ReaderClass>& Reader, std::string n = "LatticeGrid"):name(n) {
+    read(Reader, name, *this);
+    check();
+  }
+
+  // Save on file
+  template< class WriterClass>
+  void save(Writer<WriterClass>& Writer){
+    check();
+    write(Writer, name, *this);
+  }
+private:
+    std::string name;
+};
+
+// Lower level class
+class GridModule {
+ public:
+  GridCartesian* get_full() {
+    std::cout << GridLogDebug << "Getting cartesian in module"<< std::endl;
+    return grid_.get(); }
+  GridRedBlackCartesian* get_rb() {
+    std::cout << GridLogDebug << "Getting rb-cartesian in module"<< std::endl;
+    return rbgrid_.get(); }
+
+  void set_full(GridCartesian* grid) { grid_.reset(grid); }
+  void set_rb(GridRedBlackCartesian* rbgrid) { rbgrid_.reset(rbgrid); }
+  void show_full_decomposition(){ grid_->show_decomposition(); }
+  void show_rb_decomposition(){ rbgrid_->show_decomposition(); }
+
+ protected:
+  std::unique_ptr<GridCartesian> grid_;
+  std::unique_ptr<GridRedBlackCartesian> rbgrid_;
+
+};
+
+////////////////////////////////////
+// Classes for the user
+////////////////////////////////////
+// Note: the space time grid should be out of the QCD namespace
+template <class vector_type>
+class GridFourDimModule : public GridModule
+{
+public:
+  GridFourDimModule()
+  {
+    using namespace QCD;
+    set_full(SpaceTimeGrid::makeFourDimGrid(
+        GridDefaultLatt(), 
+        GridDefaultSimd(4, vector_type::Nsimd()),
+        GridDefaultMpi()));
+    set_rb(SpaceTimeGrid::makeFourDimRedBlackGrid(grid_.get()));
+  }
+
+  GridFourDimModule(const std::vector<int> tweak_simd)
+  {
+    using namespace QCD;
+    if (tweak_simd.size() != 4)
+    {
+      std::cout << GridLogError
+                << "Error in GridFourDimModule: SIMD size different from 4" 
+                << std::endl;
+      exit(1);
+    }
+
+    // Checks that the product agrees with the expectation
+    int simd_sum = 1;
+    for (auto &n : tweak_simd)
+      simd_sum *= n;
+    std::cout << GridLogDebug << "TweakSIMD: " << tweak_simd << "  Sum: " << simd_sum << std::endl;
+
+    if (simd_sum == vector_type::Nsimd())
+    {
+      set_full(SpaceTimeGrid::makeFourDimGrid(
+          GridDefaultLatt(), 
+          tweak_simd, 
+          GridDefaultMpi()));
+      set_rb(SpaceTimeGrid::makeFourDimRedBlackGrid(grid_.get()));
+    }
+    else
+    {
+      std::cout << GridLogError 
+                << "Error in GridFourDimModule: SIMD lanes must sum to " 
+                << vector_type::Nsimd() 
+                << std::endl;
+    }
+  }
+
+  GridFourDimModule(const GridModuleParameters Params)
+  {
+    using namespace QCD;
+    std::vector<int> lattice_v = Params.getLattice();
+    std::vector<int> mpi_v = Params.getMpi();
+    if (lattice_v.size() == 4)
+    {
+      set_full(SpaceTimeGrid::makeFourDimGrid(
+          lattice_v, 
+          GridDefaultSimd(4, vector_type::Nsimd()),
+          mpi_v));
+      set_rb(SpaceTimeGrid::makeFourDimRedBlackGrid(grid_.get()));
+    }
+    else
+    {
+      std::cout << GridLogError
+                << "Error in GridFourDimModule: lattice dimension different from 4"
+                << std::endl;
+      exit(1);
+    }
+  }
+};
+
+typedef GridFourDimModule<vComplex> GridDefaultFourDimModule;
+
+
+}  // namespace Grid
+
+#endif  // HMC_GRID_MODULES

Grid/qcd/hmc/HMC_aggregate.h

@@ -0,0 +1,52 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid
+
+Source file: ./lib/qcd/hmc/HMC.h
+
+Copyright (C) 2015
+
+Author: Peter Boyle <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_AGGREGATE_INCLUDED
+#define HMC_AGGREGATE_INCLUDED
+
+#include <string>
+
+#include <Grid/qcd/observables/hmc_observable.h>
+#include <Grid/qcd/hmc/HMC.h>
+
+
+// annoying location; should move this ?
+#include <Grid/parallelIO/IldgIOtypes.h>
+#include <Grid/parallelIO/IldgIO.h>
+#include <Grid/parallelIO/NerscIO.h>
+
+#include <Grid/qcd/hmc/checkpointers/CheckPointers.h>
+#include <Grid/qcd/hmc/HMCModules.h>
+#include <Grid/qcd/modules/mods.h>
+#include <Grid/qcd/hmc/HMCResourceManager.h>
+#include <Grid/qcd/hmc/GenericHMCrunner.h>
+#include <Grid/qcd/hmc/HMCRunnerModule.h>
+
+#endif

Grid/qcd/hmc/UsingHMC.md

@@ -0,0 +1,107 @@
+Using HMC in Grid version 0.5.1
+
+These are the instructions to use the Generalised HMC on Grid version 0.5.1.
+Disclaimer: GRID is still under active development so any information here can be changed in future releases.
+
+
+Command line options
+===================
+(relevant file GenericHMCrunner.h)
+The initial configuration can be changed at the command line using 
+--StartType <your choice>
+valid choices, one among these
+HotStart, ColdStart, TepidStart, CheckpointStart
+default: HotStart
+
+example
+./My_hmc_exec  --StartType HotStart
+
+The CheckpointStart option uses the prefix for the configurations and rng seed files defined in your executable and the initial configuration is specified by
+--StartTrajectory <integer>
+default: 0
+
+The number of trajectories for a specific run are specified at command line by
+--Trajectories <integer>
+default: 1
+
+The number of thermalization steps (i.e. steps when the Metropolis acceptance check is turned off) is specified by
+--Thermalizations <integer>
+default: 10
+
+
+Any other parameter is defined in the source for the executable.
+
+HMC controls
+===========
+
+The lines 
+
+  std::vector<int> SerSeed({1, 2, 3, 4, 5});
+  std::vector<int> ParSeed({6, 7, 8, 9, 10});
+
+define the seeds for the serial and the parallel RNG.
+
+The line 
+
+  TheHMC.MDparameters.set(20, 1.0);// MDsteps, traj length
+
+declares the number of molecular dynamics steps and the total trajectory length.
+
+
+Actions
+======
+
+Action names are defined in the file
+lib/qcd/Actions.h
+
+Gauge actions list:
+
+WilsonGaugeActionR;
+WilsonGaugeActionF;
+WilsonGaugeActionD;
+PlaqPlusRectangleActionR;
+PlaqPlusRectangleActionF;
+PlaqPlusRectangleActionD;
+IwasakiGaugeActionR;
+IwasakiGaugeActionF;
+IwasakiGaugeActionD;
+SymanzikGaugeActionR;
+SymanzikGaugeActionF;
+SymanzikGaugeActionD;
+
+
+ConjugateWilsonGaugeActionR;
+ConjugateWilsonGaugeActionF;
+ConjugateWilsonGaugeActionD;
+ConjugatePlaqPlusRectangleActionR;
+ConjugatePlaqPlusRectangleActionF;
+ConjugatePlaqPlusRectangleActionD;
+ConjugateIwasakiGaugeActionR;
+ConjugateIwasakiGaugeActionF;
+ConjugateIwasakiGaugeActionD;
+ConjugateSymanzikGaugeActionR;
+ConjugateSymanzikGaugeActionF;
+ConjugateSymanzikGaugeActionD;
+
+
+ScalarActionR;
+ScalarActionF;
+ScalarActionD;
+
+
+each of these action accept one single parameter at creation time (beta).
+Example for creating a Symanzik action with beta=4.0
+
+	SymanzikGaugeActionR(4.0)
+
+The suffixes R,F,D in the action names refer to the Real
+(the precision is defined at compile time by the --enable-precision flag in the configure),
+Float and Double, that force the precision of the action to be 32, 64 bit respectively.
+
+
+
+
+
+
+
+

Grid/qcd/hmc/checkpointers/BaseCheckpointer.h

@@ -0,0 +1,97 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid
+
+Source file: ./lib/qcd/hmc/BaseCheckpointer.h
+
+Copyright (C) 2015
+
+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 BASE_CHECKPOINTER
+#define BASE_CHECKPOINTER
+
+namespace Grid {
+namespace QCD {
+
+class CheckpointerParameters : Serializable {
+ public:
+  GRID_SERIALIZABLE_CLASS_MEMBERS(CheckpointerParameters, 
+  	std::string, config_prefix, 
+  	std::string, rng_prefix, 
+  	int, saveInterval, 
+  	std::string, format, );
+
+  CheckpointerParameters(std::string cf = "cfg", std::string rn = "rng",
+   		      int savemodulo = 1, const std::string &f = "IEEE64BIG")
+      : config_prefix(cf),
+        rng_prefix(rn),
+        saveInterval(savemodulo),
+        format(f){};
+
+
+  template <class ReaderClass >
+  CheckpointerParameters(Reader<ReaderClass> &Reader) {
+    read(Reader, "Checkpointer", *this);
+  }
+ 
+
+};
+
+//////////////////////////////////////////////////////////////////////////////
+// Base class for checkpointers
+template <class Impl>
+class BaseHmcCheckpointer : public HmcObservable<typename Impl::Field> {
+ public:
+  void build_filenames(int traj, CheckpointerParameters &Params,
+                       std::string &conf_file, std::string &rng_file) {
+    {
+      std::ostringstream os;
+      os << Params.rng_prefix << "." << traj;
+      rng_file = os.str();
+    }
+
+    {
+      std::ostringstream os;
+      os << Params.config_prefix << "." << traj;
+      conf_file = os.str();
+    }
+ 	} 
+
+  void check_filename(const std::string &filename){
+    std::ifstream f(filename.c_str());
+    if(!f.good()){
+      std::cout << GridLogError << "Filename " << filename << " not found. Aborting. " << std::endl;
+      abort();
+    };
+  }
+
+  virtual void initialize(const CheckpointerParameters &Params) = 0;
+
+  virtual void CheckpointRestore(int traj, typename Impl::Field &U,
+                                 GridSerialRNG &sRNG,
+                                 GridParallelRNG &pRNG) = 0;
+
+};  // class BaseHmcCheckpointer
+///////////////////////////////////////////////////////////////////////////////
+}
+}
+#endif

Grid/qcd/hmc/checkpointers/BinaryCheckpointer.h

@@ -0,0 +1,116 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid
+
+Source file: ./lib/qcd/hmc/BinaryCheckpointer.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 BINARY_CHECKPOINTER
+#define BINARY_CHECKPOINTER
+
+#include <iostream>
+#include <sstream>
+#include <string>
+
+namespace Grid {
+namespace QCD {
+
+// Simple checkpointer, only binary file
+template <class Impl>
+class BinaryHmcCheckpointer : public BaseHmcCheckpointer<Impl> {
+ private:
+  CheckpointerParameters Params;
+
+ public:
+  INHERIT_FIELD_TYPES(Impl);  // Gets the Field type, a Lattice object
+
+  // Extract types from the Field
+  typedef typename Field::vector_object vobj;
+  typedef typename vobj::scalar_object sobj;
+  typedef typename getPrecision<sobj>::real_scalar_type sobj_stype;
+  typedef typename sobj::DoublePrecision sobj_double;
+
+  BinaryHmcCheckpointer(const CheckpointerParameters &Params_) {
+    initialize(Params_);
+  }
+
+  void initialize(const CheckpointerParameters &Params_) { Params = Params_; }
+
+  void truncate(std::string file) {
+    std::ofstream fout(file, std::ios::out);
+    fout.close();
+  }
+
+  void TrajectoryComplete(int traj, Field &U, GridSerialRNG &sRNG, GridParallelRNG &pRNG) {
+
+    if ((traj % Params.saveInterval) == 0) {
+      std::string config, rng;
+      this->build_filenames(traj, Params, config, rng);
+
+      uint32_t nersc_csum;
+      uint32_t scidac_csuma;
+      uint32_t scidac_csumb;
+      
+      BinarySimpleUnmunger<sobj_double, sobj> munge;
+      truncate(rng);
+      BinaryIO::writeRNG(sRNG, pRNG, rng, 0,nersc_csum,scidac_csuma,scidac_csumb);
+      truncate(config);
+
+      BinaryIO::writeLatticeObject<vobj, sobj_double>(U, config, munge, 0, Params.format,
+						      nersc_csum,scidac_csuma,scidac_csumb);
+
+      std::cout << GridLogMessage << "Written Binary Configuration " << config
+                << " checksum " << std::hex 
+		<< nersc_csum   <<"/"
+		<< scidac_csuma   <<"/"
+		<< scidac_csumb 
+		<< std::dec << std::endl;
+    }
+
+  };
+
+  void CheckpointRestore(int traj, Field &U, GridSerialRNG &sRNG, GridParallelRNG &pRNG) {
+    std::string config, rng;
+    this->build_filenames(traj, Params, config, rng);
+    this->check_filename(rng);
+    this->check_filename(config);
+
+
+    BinarySimpleMunger<sobj_double, sobj> munge;
+
+    uint32_t nersc_csum;
+    uint32_t scidac_csuma;
+    uint32_t scidac_csumb;
+    BinaryIO::readRNG(sRNG, pRNG, rng, 0,nersc_csum,scidac_csuma,scidac_csumb);
+    BinaryIO::readLatticeObject<vobj, sobj_double>(U, config, munge, 0, Params.format,
+						   nersc_csum,scidac_csuma,scidac_csumb);
+    
+    std::cout << GridLogMessage << "Read Binary Configuration " << config
+              << " checksums " << std::hex << nersc_csum<<"/"<<scidac_csuma<<"/"<<scidac_csumb 
+	      << std::dec << std::endl;
+  };
+};
+}
+}
+#endif

Grid/qcd/hmc/checkpointers/CheckPointerModules.h

@@ -0,0 +1,172 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid
+
+Source file: ./lib/qcd/action/gauge/WilsonGaugeAction.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 CP_MODULES_H
+#define CP_MODULES_H
+
+
+// FIXME  Reorganize QCD namespace 
+
+
+namespace Grid {
+
+////////////////////////////////////////////////////////////////////////
+// Checkpoint module, owns the Checkpointer
+////////////////////////////////////////////////////////////////////////
+
+template <class ImplementationPolicy>
+class CheckPointerModule: public Parametrized<QCD::CheckpointerParameters>, public HMCModuleBase< QCD::BaseHmcCheckpointer<ImplementationPolicy> >  {
+ public:
+ 	std::unique_ptr<QCD::BaseHmcCheckpointer<ImplementationPolicy> > CheckPointPtr;
+ 	typedef QCD::CheckpointerParameters APar;
+  typedef HMCModuleBase< QCD::BaseHmcCheckpointer<ImplementationPolicy> > Base;
+  typedef typename Base::Product Product;
+
+  CheckPointerModule(APar Par): Parametrized<APar>(Par) {}
+  template <class ReaderClass>
+  CheckPointerModule(Reader<ReaderClass>& Reader) : Parametrized<APar>(Reader){};
+
+  virtual void print_parameters(){
+  	std::cout << this->Par_ << std::endl;
+  }
+
+  Product* getPtr() {
+    if (!CheckPointPtr) initialize();
+
+    return CheckPointPtr.get();
+  }
+
+ private:
+  virtual void initialize() = 0;
+
+};
+
+
+
+template <char const *str, class ImplementationPolicy, class ReaderClass >
+class HMC_CPModuleFactory
+    : public Factory < HMCModuleBase< QCD::BaseHmcCheckpointer<ImplementationPolicy> > ,	Reader<ReaderClass> > {
+ public:
+ 	typedef Reader<ReaderClass> TheReader; 
+ 	// use SINGLETON FUNCTOR MACRO HERE
+  HMC_CPModuleFactory(const HMC_CPModuleFactory& e) = delete;
+  void operator=(const HMC_CPModuleFactory& e) = delete;
+  static HMC_CPModuleFactory& getInstance(void) {
+    static HMC_CPModuleFactory e;
+    return e;
+  }
+
+ private:
+  HMC_CPModuleFactory(void) = default;
+  std::string obj_type() const {
+  	return std::string(str);
+  }
+};
+
+
+
+/////////////////////////////////////////////////////////////////////
+// Concrete classes
+/////////////////////////////////////////////////////////////////////
+namespace QCD{
+
+template<class ImplementationPolicy>
+class BinaryCPModule: public CheckPointerModule< ImplementationPolicy> {
+  typedef CheckPointerModule< ImplementationPolicy> CPBase;
+  using CPBase::CPBase; // for constructors
+
+  // acquire resource
+  virtual void initialize(){
+    this->CheckPointPtr.reset(new BinaryHmcCheckpointer<ImplementationPolicy>(this->Par_));
+  }
+
+};
+
+
+template<class ImplementationPolicy>
+class NerscCPModule: public CheckPointerModule< ImplementationPolicy> {
+  typedef CheckPointerModule< ImplementationPolicy> CPBase;
+  using CPBase::CPBase; // for constructors inheritance
+
+  // acquire resource
+  virtual void initialize(){
+     this->CheckPointPtr.reset(new NerscHmcCheckpointer<ImplementationPolicy>(this->Par_));
+  }
+
+};
+
+
+#ifdef HAVE_LIME
+  
+template<class ImplementationPolicy>
+class ILDGCPModule: public CheckPointerModule< ImplementationPolicy> {
+  typedef CheckPointerModule< ImplementationPolicy> CPBase;
+  using CPBase::CPBase; // for constructors
+
+  // acquire resource
+  virtual void initialize(){
+     this->CheckPointPtr.reset(new ILDGHmcCheckpointer<ImplementationPolicy>(this->Par_));
+  }
+
+};
+
+template<class ImplementationPolicy, class Metadata>
+class ScidacCPModule: public CheckPointerModule< ImplementationPolicy> {
+  typedef CheckPointerModule< ImplementationPolicy> CPBase;
+  Metadata M;
+
+  // acquire resource
+  virtual void initialize(){
+     this->CheckPointPtr.reset(new ScidacHmcCheckpointer<ImplementationPolicy, Metadata>(this->Par_, M));
+  }
+public:
+  ScidacCPModule(typename CPBase::APar Par, Metadata M_):M(M_), CPBase(Par) {}
+  template <class ReaderClass>
+  ScidacCPModule(Reader<ReaderClass>& Reader) : Parametrized<typename CPBase::APar>(Reader), M(Reader){};
+};
+#endif
+
+
+}// QCD temporarily here
+
+
+extern char cp_string[];
+
+/*
+// use macros?
+static Registrar<QCD::BinaryCPModule<QCD::PeriodicGimplR>, HMC_CPModuleFactory<cp_string, QCD::PeriodicGimplR, XmlReader> > __CPBinarymodXMLInit("Binary");
+static Registrar<QCD::NerscCPModule<QCD::PeriodicGimplR> , HMC_CPModuleFactory<cp_string, QCD::PeriodicGimplR, XmlReader> > __CPNerscmodXMLInit("Nersc");
+
+#ifdef HAVE_LIME
+static Registrar<QCD::ILDGCPModule<QCD::PeriodicGimplR>  , HMC_CPModuleFactory<cp_string, QCD::PeriodicGimplR, XmlReader> > __CPILDGmodXMLInit("ILDG");
+#endif
+*/
+
+}// Grid
+#endif //CP_MODULES_H

Grid/qcd/hmc/checkpointers/CheckPointers.h

@@ -0,0 +1,41 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid
+
+Source file: ./lib/qcd/action/gauge/WilsonGaugeAction.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 CHECKPOINTERS_H
+#define CHECKPOINTERS_H
+
+#include <Grid/qcd/hmc/checkpointers/BaseCheckpointer.h>
+#include <Grid/qcd/hmc/checkpointers/NerscCheckpointer.h>
+#include <Grid/qcd/hmc/checkpointers/BinaryCheckpointer.h>
+#include <Grid/qcd/hmc/checkpointers/ILDGCheckpointer.h>
+#include <Grid/qcd/hmc/checkpointers/ScidacCheckpointer.h>
+//#include <Grid/qcd/hmc/checkpointers/CheckPointerModules.h>
+
+
+#endif // CHECKPOINTERS_H

Grid/qcd/hmc/checkpointers/ILDGCheckpointer.h

@@ -0,0 +1,124 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid
+
+Source file: ./lib/qcd/hmc/ILDGCheckpointer.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 ILDG_CHECKPOINTER
+#define ILDG_CHECKPOINTER
+
+#ifdef HAVE_LIME
+
+#include <iostream>
+#include <sstream>
+#include <string>
+
+namespace Grid {
+namespace QCD {
+
+// Only for Gauge fields
+template <class Implementation>
+class ILDGHmcCheckpointer : public BaseHmcCheckpointer<Implementation> {
+ private:
+  CheckpointerParameters Params;
+
+ public:
+  INHERIT_GIMPL_TYPES(Implementation);
+
+  ILDGHmcCheckpointer(const CheckpointerParameters &Params_) { initialize(Params_); }
+
+  void initialize(const CheckpointerParameters &Params_) {
+    Params = Params_;
+
+    // check here that the format is valid
+    int ieee32big = (Params.format == std::string("IEEE32BIG"));
+    int ieee32    = (Params.format == std::string("IEEE32"));
+    int ieee64big = (Params.format == std::string("IEEE64BIG"));
+    int ieee64    = (Params.format == std::string("IEEE64"));
+
+    if (!(ieee64big || ieee32 || ieee32big || ieee64)) {
+      std::cout << GridLogError << "Unrecognized file format " << Params.format
+                << std::endl;
+      std::cout << GridLogError
+                << "Allowed: IEEE32BIG | IEEE32 | IEEE64BIG | IEEE64"
+                << std::endl;
+
+      exit(1);
+    }
+  }
+
+  void TrajectoryComplete(int traj, GaugeField &U, GridSerialRNG &sRNG,
+                          GridParallelRNG &pRNG) {
+    if ((traj % Params.saveInterval) == 0) {
+      std::string config, rng;
+      this->build_filenames(traj, Params, config, rng);
+      GridBase *grid = U._grid;
+      uint32_t nersc_csum,scidac_csuma,scidac_csumb;
+      BinaryIO::writeRNG(sRNG, pRNG, rng, 0,nersc_csum,scidac_csuma,scidac_csumb);
+      IldgWriter _IldgWriter(grid->IsBoss());
+      _IldgWriter.open(config);
+      _IldgWriter.writeConfiguration(U, traj, config, config);
+      _IldgWriter.close();
+
+      std::cout << GridLogMessage << "Written ILDG Configuration on " << config
+                << " checksum " << std::hex 
+		<< nersc_csum<<"/"
+		<< scidac_csuma<<"/"
+		<< scidac_csumb
+		<< std::dec << std::endl;
+    }
+  };
+
+  void CheckpointRestore(int traj, GaugeField &U, GridSerialRNG &sRNG,
+                         GridParallelRNG &pRNG) {
+    std::string config, rng;
+    this->build_filenames(traj, Params, config, rng);
+    this->check_filename(rng);
+    this->check_filename(config);
+
+    
+
+    uint32_t nersc_csum,scidac_csuma,scidac_csumb;
+    BinaryIO::readRNG(sRNG, pRNG, rng, 0,nersc_csum,scidac_csuma,scidac_csumb);
+
+    FieldMetaData header;
+    IldgReader _IldgReader;
+    _IldgReader.open(config);
+    _IldgReader.readConfiguration(U,header);  // format from the header
+    _IldgReader.close();
+
+    std::cout << GridLogMessage << "Read ILDG Configuration from " << config
+              << " checksum " << std::hex 
+	      << nersc_csum<<"/"
+	      << scidac_csuma<<"/"
+	      << scidac_csumb
+	      << std::dec << std::endl;
+  };
+};
+}
+}
+
+#endif  // HAVE_LIME
+#endif  // ILDG_CHECKPOINTER

Grid/qcd/hmc/checkpointers/NerscCheckpointer.h

@@ -0,0 +1,83 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid
+
+Source file: ./lib/qcd/hmc/NerscCheckpointer.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 NERSC_CHECKPOINTER
+#define NERSC_CHECKPOINTER
+
+#include <iostream>
+#include <sstream>
+#include <string>
+
+namespace Grid {
+namespace QCD {
+
+// Only for Gauge fields
+template <class Gimpl>
+class NerscHmcCheckpointer : public BaseHmcCheckpointer<Gimpl> {
+ private:
+  CheckpointerParameters Params;
+
+ public:
+  INHERIT_GIMPL_TYPES(Gimpl);  // only for gauge configurations
+
+  NerscHmcCheckpointer(const CheckpointerParameters &Params_) { initialize(Params_); }
+
+  void initialize(const CheckpointerParameters &Params_) {
+    Params = Params_;
+    Params.format = "IEEE64BIG";  // fixed, overwrite any other choice
+  }
+
+  void TrajectoryComplete(int traj, GaugeField &U, GridSerialRNG &sRNG,
+                          GridParallelRNG &pRNG) {
+    if ((traj % Params.saveInterval) == 0) {
+      std::string config, rng;
+      this->build_filenames(traj, Params, config, rng);
+
+      int precision32 = 1;
+      int tworow = 0;
+      NerscIO::writeRNGState(sRNG, pRNG, rng);
+      NerscIO::writeConfiguration(U, config, tworow, precision32);
+    }
+  };
+
+  void CheckpointRestore(int traj, GaugeField &U, GridSerialRNG &sRNG,
+                         GridParallelRNG &pRNG) {
+    std::string config, rng;
+    this->build_filenames(traj, Params, config, rng);
+    this->check_filename(rng);
+    this->check_filename(config);
+
+
+    FieldMetaData header;
+    NerscIO::readRNGState(sRNG, pRNG, header, rng);
+    NerscIO::readConfiguration(U, header, config);
+  };
+};
+}
+}
+#endif

Grid/qcd/hmc/checkpointers/ScidacCheckpointer.h

@@ -0,0 +1,122 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid
+
+Source file: ./lib/qcd/hmc/ScidacCheckpointer.h
+
+Copyright (C) 2018
+
+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 SCIDAC_CHECKPOINTER
+#define SCIDAC_CHECKPOINTER
+
+#ifdef HAVE_LIME
+
+#include <iostream>
+#include <sstream>
+#include <string>
+
+namespace Grid {
+namespace QCD {
+
+// For generic fields
+template <class Implementation, class Metadata>
+class ScidacHmcCheckpointer : public BaseHmcCheckpointer<Implementation> {
+ private:
+  CheckpointerParameters Params;
+  Metadata MData;
+
+  typedef typename Implementation::Field Field;
+
+ public:
+  //INHERIT_GIMPL_TYPES(Implementation);
+
+  ScidacHmcCheckpointer(const CheckpointerParameters &Params_) { initialize(Params_); }
+  ScidacHmcCheckpointer(const CheckpointerParameters &Params_, const Metadata& M_):MData(M_) { initialize(Params_); }
+
+  void initialize(const CheckpointerParameters &Params_) {
+    Params = Params_;
+
+    // check here that the format is valid
+    int ieee32big = (Params.format == std::string("IEEE32BIG"));
+    int ieee32    = (Params.format == std::string("IEEE32"));
+    int ieee64big = (Params.format == std::string("IEEE64BIG"));
+    int ieee64    = (Params.format == std::string("IEEE64"));
+
+    if (!(ieee64big || ieee32 || ieee32big || ieee64)) {
+      std::cout << GridLogError << "Unrecognized file format " << Params.format
+                << std::endl;
+      std::cout << GridLogError
+                << "Allowed: IEEE32BIG | IEEE32 | IEEE64BIG | IEEE64"
+                << std::endl;
+
+      exit(1);
+    }
+  }
+
+  void TrajectoryComplete(int traj, Field &U, GridSerialRNG &sRNG,
+                          GridParallelRNG &pRNG) {
+    if ((traj % Params.saveInterval) == 0) {
+      std::string config, rng;
+      this->build_filenames(traj, Params, config, rng);
+      GridBase *grid = U._grid;
+      uint32_t nersc_csum,scidac_csuma,scidac_csumb;
+      BinaryIO::writeRNG(sRNG, pRNG, rng, 0,nersc_csum,scidac_csuma,scidac_csumb);
+      ScidacWriter _ScidacWriter(grid->IsBoss());
+      _ScidacWriter.open(config);
+      _ScidacWriter.writeScidacFieldRecord(U, MData);
+      _ScidacWriter.close();
+
+      std::cout << GridLogMessage << "Written Scidac Configuration on " << config << std::endl;
+    }
+  };
+
+  void CheckpointRestore(int traj, Field &U, GridSerialRNG &sRNG,
+                         GridParallelRNG &pRNG) {
+    std::string config, rng;
+    this->build_filenames(traj, Params, config, rng);
+    this->check_filename(rng);
+    this->check_filename(config);
+
+
+    uint32_t nersc_csum,scidac_csuma,scidac_csumb;
+    BinaryIO::readRNG(sRNG, pRNG, rng, 0,nersc_csum,scidac_csuma,scidac_csumb);
+
+    Metadata md_content;
+    ScidacReader _ScidacReader;
+    _ScidacReader.open(config);
+    _ScidacReader.readScidacFieldRecord(U,md_content);  // format from the header
+    _ScidacReader.close();
+
+    std::cout << GridLogMessage << "Read Scidac Configuration from " << config
+              << " checksum " << std::hex 
+	      << nersc_csum<<"/"
+	      << scidac_csuma<<"/"
+	      << scidac_csumb
+	      << std::dec << std::endl;
+  };
+};
+}
+}
+
+#endif  // HAVE_LIME
+#endif  // ILDG_CHECKPOINTER

Grid/qcd/hmc/integrators/Integrator.h

@@ -0,0 +1,325 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid
+
+Source file: ./lib/qcd/hmc/integrators/Integrator.h
+
+Copyright (C) 2015
+
+Author: Azusa Yamaguchi <ayamaguc@staffmail.ed.ac.uk>
+Author: Peter Boyle <paboyle@ph.ed.ac.uk>
+Author: Guido Cossu <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 INTEGRATOR_INCLUDED
+#define INTEGRATOR_INCLUDED
+
+#include <memory>
+
+namespace Grid {
+namespace QCD {
+
+class IntegratorParameters: Serializable {
+public:
+        GRID_SERIALIZABLE_CLASS_MEMBERS(IntegratorParameters,
+                std::string, name,      // name of the integrator
+    unsigned int, MDsteps,  // number of outer steps
+    RealD, trajL,           // trajectory length
+  )
+
+  IntegratorParameters(int MDsteps_ = 10, RealD trajL_ = 1.0)
+      : MDsteps(MDsteps_),
+        trajL(trajL_){
+        // empty body constructor
+  };
+
+
+  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);
+  }
+
+  void print_parameters() const {
+    std::cout << GridLogMessage << "[Integrator] Type               : " << name << std::endl;
+    std::cout << GridLogMessage << "[Integrator] Trajectory length  : " << trajL << std::endl;
+    std::cout << GridLogMessage << "[Integrator] Number of MD steps : " << MDsteps << std::endl;
+    std::cout << GridLogMessage << "[Integrator] Step size          : " << trajL/MDsteps << std::endl;
+  }
+};
+
+/*! @brief Class for Molecular Dynamics management */
+template <class FieldImplementation, class SmearingPolicy, class RepresentationPolicy>
+class Integrator {
+ protected:
+  typedef typename FieldImplementation::Field MomentaField;  //for readability
+  typedef typename FieldImplementation::Field Field;
+
+  int levels;  // number of integration levels
+  double t_U;  // Track time passing on each level and for U and for P
+  std::vector<double> t_P;  
+
+  MomentaField P;
+  SmearingPolicy& Smearer;
+  RepresentationPolicy Representations;
+  IntegratorParameters Params;
+
+  const ActionSet<Field, RepresentationPolicy> as;
+
+  void update_P(Field& U, int level, double ep) {
+    t_P[level] += ep;
+    update_P(P, U, level, ep);
+
+    std::cout << GridLogIntegrator << "[" << level << "] P "
+              << " dt " << ep << " : t_P " << t_P[level] << std::endl;
+  }
+
+  // to be used by the actionlevel class to iterate
+  // over the representations
+  struct _updateP {
+    template <class FieldType, class GF, class Repr>
+    void operator()(std::vector<Action<FieldType>*> repr_set, Repr& Rep,
+                    GF& Mom, GF& U, double ep) {
+      for (int a = 0; a < repr_set.size(); ++a) {
+        FieldType forceR(U._grid);
+        // Implement smearing only for the fundamental representation now
+        repr_set.at(a)->deriv(Rep.U, forceR);
+        GF force = Rep.RtoFundamentalProject(forceR);  // Ta for the fundamental rep
+        Real force_abs = std::sqrt(norm2(force)/(U._grid->gSites()));
+        std::cout << GridLogIntegrator << "Hirep Force average: " << force_abs << std::endl;
+        Mom -= force * ep ;
+      }
+    }
+  } update_P_hireps{};
+
+  void update_P(MomentaField& Mom, Field& U, int level, double ep) {
+    // input U actually not used in the fundamental case
+    // Fundamental updates, include smearing
+
+    for (int a = 0; a < as[level].actions.size(); ++a) {
+      double start_full = usecond();
+      Field force(U._grid);
+      conformable(U._grid, Mom._grid);
+
+      Field& Us = Smearer.get_U(as[level].actions.at(a)->is_smeared);
+      double start_force = usecond();
+      as[level].actions.at(a)->deriv(Us, force);  // deriv should NOT include Ta
+
+      std::cout << GridLogIntegrator << "Smearing (on/off): " << as[level].actions.at(a)->is_smeared << std::endl;
+      if (as[level].actions.at(a)->is_smeared) Smearer.smeared_force(force);
+      force = FieldImplementation::projectForce(force); // Ta for gauge fields
+      double end_force = usecond();
+      Real force_abs = std::sqrt(norm2(force)/U._grid->gSites());
+      std::cout << GridLogIntegrator << "["<<level<<"]["<<a<<"] Force average: " << force_abs << std::endl;
+      Mom -= force * ep; 
+      double end_full = usecond();
+      double time_full  = (end_full - start_full) / 1e3;
+      double time_force = (end_force - start_force) / 1e3;
+      std::cout << GridLogIntegrator << "["<<level<<"]["<<a<<"] P update elapsed time: " << time_full << " ms (force: " << time_force << " ms)"  << std::endl;
+    }
+
+    // Force from the other representations
+    as[level].apply(update_P_hireps, Representations, Mom, U, ep);
+  }
+
+  void update_U(Field& U, double ep) {
+    update_U(P, U, ep);
+
+    t_U += ep;
+    int fl = levels - 1;
+    std::cout << GridLogIntegrator << "   " << "[" << fl << "] U " << " dt " << ep << " : t_U " << t_U << std::endl;
+  }
+  
+  void update_U(MomentaField& Mom, Field& U, double ep) {
+    // exponential of Mom*U in the gauge fields case
+    FieldImplementation::update_field(Mom, U, ep);
+
+    // Update the smeared fields, can be implemented as observer
+    Smearer.set_Field(U);
+
+    // Update the higher representations fields
+    Representations.update(U);  // void functions if fundamental representation
+  }
+
+  virtual void step(Field& U, int level, int first, int last) = 0;
+
+ public:
+  Integrator(GridBase* grid, IntegratorParameters Par,
+             ActionSet<Field, RepresentationPolicy>& Aset,
+             SmearingPolicy& Sm)
+      : Params(Par),
+        as(Aset),
+        P(grid),
+        levels(Aset.size()),
+        Smearer(Sm),
+        Representations(grid) {
+    t_P.resize(levels, 0.0);
+    t_U = 0.0;
+    // initialization of smearer delegated outside of Integrator
+  };
+
+  virtual ~Integrator() {}
+
+  virtual std::string integrator_name() = 0;
+
+  void print_parameters(){
+    std::cout << GridLogMessage << "[Integrator] Name : "<< integrator_name() << std::endl;
+    Params.print_parameters();
+  }
+
+  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;
+
+  }
+
+  void reverse_momenta(){
+    P *= -1.0;
+  }
+
+  // to be used by the actionlevel class to iterate
+  // over the representations
+  struct _refresh {
+    template <class FieldType, class Repr>
+    void operator()(std::vector<Action<FieldType>*> repr_set, Repr& Rep,
+                    GridParallelRNG& pRNG) {
+      for (int a = 0; a < repr_set.size(); ++a){
+        repr_set.at(a)->refresh(Rep.U, pRNG);
+      
+      std::cout << GridLogDebug << "Hirep refreshing pseudofermions" << std::endl;
+    }
+    }
+  } refresh_hireps{};
+
+  // Initialization of momenta and actions
+  void refresh(Field& U, GridParallelRNG& pRNG) {
+    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
+    Smearer.set_Field(U);
+    // Set the (eventual) representations gauge fields
+    Representations.update(U);
+
+    // The Smearer is attached to a pointer of the gauge field
+    // automatically gets the correct field
+    // whether or not has been accepted in the previous sweep
+    for (int level = 0; level < as.size(); ++level) {
+      for (int actionID = 0; actionID < as[level].actions.size(); ++actionID) {
+        // get gauge field from the SmearingPolicy and
+        // based on the boolean is_smeared in actionID
+        Field& Us =
+            Smearer.get_U(as[level].actions.at(actionID)->is_smeared);
+        as[level].actions.at(actionID)->refresh(Us, pRNG);
+      }
+
+      // Refresh the higher representation actions
+      as[level].apply(refresh_hireps, Representations, pRNG);
+    }
+  }
+
+  // to be used by the actionlevel class to iterate
+  // over the representations
+  struct _S {
+    template <class FieldType, class Repr>
+    void operator()(std::vector<Action<FieldType>*> repr_set, Repr& Rep,
+                    int level, RealD& H) {
+      
+      for (int a = 0; a < repr_set.size(); ++a) {
+        RealD Hterm = repr_set.at(a)->S(Rep.U);
+        std::cout << GridLogMessage << "S Level " << level << " term " << a
+                  << " H Hirep = " << Hterm << std::endl;
+        H += Hterm;
+
+      }
+    }
+  } S_hireps{};
+
+  // Calculate action
+  RealD S(Field& U) {  // here also U not used
+
+    RealD H = - FieldImplementation::FieldSquareNorm(P); // - trace (P*P)
+    RealD Hterm;
+    std::cout << GridLogMessage << "Momentum action H_p = " << H << "\n";
+
+    // Actions
+    for (int level = 0; level < as.size(); ++level) {
+      for (int actionID = 0; actionID < as[level].actions.size(); ++actionID) {
+        // get gauge field from the SmearingPolicy and
+        // based on the boolean is_smeared in actionID
+        Field& Us =
+            Smearer.get_U(as[level].actions.at(actionID)->is_smeared);
+        Hterm = as[level].actions.at(actionID)->S(Us);
+        std::cout << GridLogMessage << "S Level " << level << " term "
+                  << actionID << " H = " << Hterm << std::endl;
+        H += Hterm;
+      }
+      as[level].apply(S_hireps, Representations, level, H);
+    }
+
+    return H;
+  }
+
+  void integrate(Field& U) {
+    // reset the clocks
+    t_U = 0;
+    for (int level = 0; level < as.size(); ++level) {
+      t_P[level] = 0;
+    }
+
+    for (int step = 0; step < Params.MDsteps; ++step) {  // MD step
+      int first_step = (step == 0);
+      int last_step = (step == Params.MDsteps - 1);
+      this->step(U, 0, first_step, last_step);
+    }
+
+    // Check the clocks all match on all levels
+    for (int level = 0; level < as.size(); ++level) {
+      assert(fabs(t_U - t_P[level]) < 1.0e-6);  // must be the same
+      std::cout << GridLogIntegrator << " times[" << level
+                << "]= " << t_P[level] << " " << t_U << std::endl;
+    }
+
+    // and that we indeed got to the end of the trajectory
+    assert(fabs(t_U - Params.trajL) < 1.0e-6);
+
+  }
+
+
+  
+
+};
+}
+}
+#endif  // INTEGRATOR_INCLUDED
+

Grid/qcd/hmc/integrators/Integrator_algorithm.h

@@ -0,0 +1,295 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid
+
+Source file: ./lib/qcd/hmc/integrators/Integrator_algorithm.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 */
+//--------------------------------------------------------------------
+
+
+/*! @file Integrator_algorithm.h
+ * @brief Declaration of classes for the Molecular Dynamics algorithms
+ *
+ */
+//--------------------------------------------------------------------
+
+#ifndef INTEGRATOR_ALG_INCLUDED
+#define INTEGRATOR_ALG_INCLUDED
+
+namespace Grid {
+namespace QCD {
+
+/* PAB:
+ *
+ * Recursive leapfrog; explanation of nested stepping
+ *
+ * Nested 1:4; units in dt for top level integrator
+ *
+ * CHROMA                           IroIro
+ *   0        1                      0
+ *  P 1/2                           P 1/2
+ *          P 1/16                                  P1/16
+ *                 U 1/8                                   U1/8
+ *          P 1/8                                   P1/8
+ *                 U 1/8                                   U1/8
+ *          P 1/8                                   P1/8
+ *                 U 1/8                                 U1/8
+ *          P 1/8                                   P1/8
+ *                 U 1/8                                 U1/8
+ *          P 1/16                                  P1/8
+ *  P 1                             P 1
+ *          P 1/16                    * skipped --- avoids revaluating force
+ *                 U 1/8                                 U1/8
+ *          P 1/8                                   P1/8
+ *                 U 1/8                                 U1/8
+ *          P 1/8                                   P1/8
+ *                 U 1/8                                 U1/8
+ *          P 1/8                                   P1/8
+ *                 U 1/8                                 U1/8
+ *          P 1/16                                  P1/8
+ *  P 1                             P 1
+ *          P 1/16                    * skipped
+ *                 U 1/8                                 U1/8
+ *          P 1/8                                   P1/8
+ *                 U 1/8                                 U1/8
+ *          P 1/8                                   P1/8
+ *                 U 1/8                                 U1/8
+ *          P 1/8                                   P1/8
+ *                 U 1/8                                 U1/8
+ *          P 1/16                    * skipped
+ *  P 1                             P 1
+ *          P 1/16                                  P1/8
+ *                 U 1/8                                 U1/8
+ *          P 1/8                                   P1/8
+ *                 U 1/8                                 U1/8
+ *          P 1/8                                   P1/8
+ *                 U 1/8                                 U1/8
+ *          P 1/8                                   P1/8
+ *                 U 1/8                                 U1/8
+ *          P 1/16                                  P1/16
+ *  P 1/2                            P 1/2
+ */
+
+template <class FieldImplementation, class SmearingPolicy,
+          class RepresentationPolicy =
+              Representations<FundamentalRepresentation> >
+class LeapFrog : public Integrator<FieldImplementation, SmearingPolicy,
+                                   RepresentationPolicy> {
+ public:
+  typedef LeapFrog<FieldImplementation, SmearingPolicy, RepresentationPolicy>
+      Algorithm;
+  INHERIT_FIELD_TYPES(FieldImplementation);
+
+  std::string integrator_name(){return "LeapFrog";}
+
+  LeapFrog(GridBase* grid, IntegratorParameters Par,
+           ActionSet<Field, RepresentationPolicy>& Aset, SmearingPolicy& Sm)
+      : Integrator<FieldImplementation, SmearingPolicy, RepresentationPolicy>(
+            grid, Par, Aset, Sm){};
+
+  void step(Field& U, int level, int _first, int _last) {
+    int fl = this->as.size() - 1;
+    // level  : current level
+    // fl     : final level
+    // eps    : current step size
+
+    // Get current level step size
+    RealD eps = this->Params.trajL/this->Params.MDsteps;
+    for (int l = 0; l <= level; ++l) eps /= this->as[l].multiplier;
+
+    int multiplier = this->as[level].multiplier;
+    for (int e = 0; e < multiplier; ++e) {
+      int first_step = _first && (e == 0);
+      int last_step = _last && (e == multiplier - 1);
+
+      if (first_step) {  // initial half step
+        this->update_P(U, level, eps / 2.0);
+      }
+
+      if (level == fl) {  // lowest level
+        this->update_U(U, eps);
+      } else {  // recursive function call
+        this->step(U, level + 1, first_step, last_step);
+      }
+
+      int mm = last_step ? 1 : 2;
+      this->update_P(U, level, mm * eps / 2.0);
+    }
+  }
+};
+
+template <class FieldImplementation, class SmearingPolicy,
+          class RepresentationPolicy =
+              Representations<FundamentalRepresentation> >
+class MinimumNorm2 : public Integrator<FieldImplementation, SmearingPolicy,
+                                       RepresentationPolicy> {
+ private:
+  const RealD lambda = 0.1931833275037836;
+
+ public:
+  INHERIT_FIELD_TYPES(FieldImplementation);
+
+  MinimumNorm2(GridBase* grid, IntegratorParameters Par,
+               ActionSet<Field, RepresentationPolicy>& Aset, SmearingPolicy& Sm)
+      : Integrator<FieldImplementation, SmearingPolicy, RepresentationPolicy>(
+            grid, Par, Aset, Sm){};
+
+  std::string integrator_name(){return "MininumNorm2";}
+
+  void step(Field& U, int level, int _first, int _last) {
+    // level  : current level
+    // fl     : final level
+    // eps    : current step size
+
+    int fl = this->as.size() - 1;
+
+    RealD eps = this->Params.trajL/this->Params.MDsteps * 2.0;
+    for (int l = 0; l <= level; ++l) eps /= 2.0 * this->as[l].multiplier;
+
+    // Nesting:  2xupdate_U of size eps/2
+    // Next level is eps/2/multiplier
+
+    int multiplier = this->as[level].multiplier;
+    for (int e = 0; e < multiplier; ++e) {  // steps per step
+
+      int first_step = _first && (e == 0);
+      int last_step = _last && (e == multiplier - 1);
+
+      if (first_step) {  // initial half step
+        this->update_P(U, level, lambda * eps);
+      }
+
+      if (level == fl) {  // lowest level
+        this->update_U(U, 0.5 * eps);
+      } else {  // recursive function call
+        this->step(U, level + 1, first_step, 0);
+      }
+
+      this->update_P(U, level, (1.0 - 2.0 * lambda) * eps);
+
+      if (level == fl) {  // lowest level
+        this->update_U(U, 0.5 * eps);
+      } else {  // recursive function call
+        this->step(U, level + 1, 0, last_step);
+      }
+
+      int mm = (last_step) ? 1 : 2;
+      this->update_P(U, level, lambda * eps * mm);
+    }
+  }
+};
+
+template <class FieldImplementation, class SmearingPolicy,
+          class RepresentationPolicy =
+              Representations<FundamentalRepresentation> >
+class ForceGradient : public Integrator<FieldImplementation, SmearingPolicy,
+                                        RepresentationPolicy> {
+ private:
+  const RealD lambda = 1.0 / 6.0;
+  ;
+  const RealD chi = 1.0 / 72.0;
+  const RealD xi = 0.0;
+  const RealD theta = 0.0;
+
+ public:
+  INHERIT_FIELD_TYPES(FieldImplementation);
+
+  // Looks like dH scales as dt^4. tested wilson/wilson 2 level.
+  ForceGradient(GridBase* grid, IntegratorParameters Par,
+                ActionSet<Field, RepresentationPolicy>& Aset,
+                SmearingPolicy& Sm)
+      : Integrator<FieldImplementation, SmearingPolicy, RepresentationPolicy>(
+            grid, Par, Aset, Sm){};
+
+  std::string integrator_name(){return "ForceGradient";}
+  
+  void FG_update_P(Field& U, int level, double fg_dt, double ep) {
+    Field Ufg(U._grid);
+    Field Pfg(U._grid);
+    Ufg = U;
+    Pfg = zero;
+    std::cout << GridLogIntegrator << "FG update " << fg_dt << " " << ep
+              << std::endl;
+    // prepare_fg; no prediction/result cache for now
+    // could relax CG stopping conditions for the
+    // derivatives in the small step since the force gets multiplied by
+    // a tiny dt^2 term relative to main force.
+    //
+    // Presently 4 force evals, and should have 3, so 1.33x too expensive.
+    // could reduce this with sloppy CG to perhaps 1.15x too expensive
+    // even without prediction.
+    this->update_P(Pfg, Ufg, level, 1.0);
+    this->update_U(Pfg, Ufg, fg_dt);
+    this->update_P(Ufg, level, ep);
+  }
+
+  void step(Field& U, int level, int _first, int _last) {
+    RealD eps = this->Params.trajL/this->Params.MDsteps * 2.0;
+    for (int l = 0; l <= level; ++l) eps /= 2.0 * this->as[l].multiplier;
+
+    RealD Chi = chi * eps * eps * eps;
+
+    int fl = this->as.size() - 1;
+
+    int multiplier = this->as[level].multiplier;
+
+    for (int e = 0; e < multiplier; ++e) {  // steps per step
+
+      int first_step = _first && (e == 0);
+      int last_step = _last && (e == multiplier - 1);
+
+      if (first_step) {  // initial half step
+        this->update_P(U, level, lambda * eps);
+      }
+
+      if (level == fl) {  // lowest level
+        this->update_U(U, 0.5 * eps);
+      } else {  // recursive function call
+        this->step(U, level + 1, first_step, 0);
+      }
+
+      this->FG_update_P(U, level, 2 * Chi / ((1.0 - 2.0 * lambda) * eps),
+                        (1.0 - 2.0 * lambda) * eps);
+
+      if (level == fl) {  // lowest level
+        this->update_U(U, 0.5 * eps);
+      } else {  // recursive function call
+        this->step(U, level + 1, 0, last_step);
+      }
+
+      int mm = (last_step) ? 1 : 2;
+      this->update_P(U, level, lambda * eps * mm);
+    }
+  }
+};
+
+
+
+
+}
+}
+
+#endif  // INTEGRATOR_INCLUDED