More progress in the HMC construction

lib/Make.inc

@@ -1,4 +1,4 @@
 
-HFILES=./Cshift.h ./simd/Grid_avx.h ./simd/Grid_vector_types.h ./simd/Grid_sse4.h ./simd/Grid_avx512.h ./simd/Grid_empty.h ./simd/Grid_vector_unops.h ./simd/Grid_neon.h ./simd/Grid_qpx.h ./Tensors.h ./Algorithms.h ./communicator/Communicator_base.h ./lattice/Lattice_rng.h ./lattice/Lattice_reduction.h ./lattice/Lattice_transfer.h ./lattice/Lattice_unary.h ./lattice/Lattice_peekpoke.h ./lattice/Lattice_coordinate.h ./lattice/Lattice_comparison.h ./lattice/Lattice_overload.h ./lattice/Lattice_reality.h ./lattice/Lattice_local.h ./lattice/Lattice_conformable.h ./lattice/Lattice_where.h ./lattice/Lattice_comparison_utils.h ./lattice/Lattice_arith.h ./lattice/Lattice_base.h ./lattice/Lattice_ET.h ./lattice/Lattice_transpose.h ./lattice/Lattice_trace.h ./Stencil.h ./tensors/Tensor_arith_sub.h ./tensors/Tensor_exp.h ./tensors/Tensor_arith_mul.h ./tensors/Tensor_class.h ./tensors/Tensor_logical.h ./tensors/Tensor_transpose.h ./tensors/Tensor_arith_mac.h ./tensors/Tensor_arith_scalar.h ./tensors/Tensor_reality.h ./tensors/Tensor_trace.h ./tensors/Tensor_index.h ./tensors/Tensor_arith_add.h ./tensors/Tensor_outer.h ./tensors/Tensor_inner.h ./tensors/Tensor_traits.h ./tensors/Tensor_Ta.h ./tensors/Tensor_unary.h ./tensors/Tensor_determinant.h ./tensors/Tensor_arith.h ./tensors/Tensor_extract_merge.h ./Communicator.h ./Cartesian.h ./parallelIO/NerscIO.h ./qcd/QCD.h ./qcd/hmc/integrators/Integrator_base.h ./qcd/hmc/HMC.h ./qcd/utils/SpaceTimeGrid.h ./qcd/utils/SUn.h ./qcd/utils/LinalgUtils.h ./qcd/utils/CovariantCshift.h ./qcd/utils/WilsonLoops.h ./qcd/action/gauge/WilsonGaugeAction.h ./qcd/action/gauge/GaugeActionBase.h ./qcd/action/Actions.h ./qcd/action/fermion/CayleyFermion5D.h ./qcd/action/fermion/ScaledShamirFermion.h ./qcd/action/fermion/MobiusFermion.h ./qcd/action/fermion/OverlapWilsonContfracTanhFermion.h ./qcd/action/fermion/PartialFractionFermion5D.h ./qcd/action/fermion/ShamirZolotarevFermion.h ./qcd/action/fermion/FermionOperator.h ./qcd/action/fermion/WilsonFermion5D.h ./qcd/action/fermion/WilsonCompressor.h ./qcd/action/fermion/OverlapWilsonPartialFractionZolotarevFermion.h ./qcd/action/fermion/WilsonKernels.h ./qcd/action/fermion/DomainWallFermion.h ./qcd/action/fermion/OverlapWilsonPartialFractionTanhFermion.h ./qcd/action/fermion/OverlapWilsonContfracZolotarevFermion.h ./qcd/action/fermion/MobiusZolotarevFermion.h ./qcd/action/fermion/g5HermitianLinop.h ./qcd/action/fermion/OverlapWilsonCayleyTanhFermion.h ./qcd/action/fermion/WilsonFermion.h ./qcd/action/fermion/ContinuedFractionFermion5D.h ./qcd/action/fermion/OverlapWilsonCayleyZolotarevFermion.h ./qcd/spin/TwoSpinor.h ./qcd/spin/Dirac.h ./cshift/Cshift_common.h ./cshift/Cshift_none.h ./cshift/Cshift_mpi.h ./Simd.h ./GridConfig.h ./cartesian/Cartesian_base.h ./cartesian/Cartesian_red_black.h ./cartesian/Cartesian_full.h ./AlignedAllocator.h ./Lattice.h ./Old/Tensor_poke.h ./Old/Tensor_peek.h ./Threads.h ./Grid.h ./algorithms/Preconditioner.h ./algorithms/iterative/ConjugateResidual.h ./algorithms/iterative/PrecGeneralisedConjugateResidual.h ./algorithms/iterative/ConjugateGradientMultiShift.h ./algorithms/iterative/SchurRedBlack.h ./algorithms/iterative/NormalEquations.h ./algorithms/iterative/ConjugateGradient.h ./algorithms/iterative/AdefGeneric.h ./algorithms/approx/Chebyshev.h ./algorithms/approx/Zolotarev.h ./algorithms/approx/MultiShiftFunction.h ./algorithms/approx/bigfloat.h ./algorithms/approx/bigfloat_double.h ./algorithms/approx/Remez.h ./algorithms/LinearOperator.h ./algorithms/SparseMatrix.h ./algorithms/CoarsenedMatrix.h ./stencil/Lebesgue.h
+HFILES=./Cshift.h ./simd/Grid_avx.h ./simd/Grid_vector_types.h ./simd/Grid_sse4.h ./simd/Grid_avx512.h ./simd/Grid_empty.h ./simd/Grid_vector_unops.h ./simd/Grid_neon.h ./simd/Grid_qpx.h ./Tensors.h ./Algorithms.h ./communicator/Communicator_base.h ./lattice/Lattice_rng.h ./lattice/Lattice_reduction.h ./lattice/Lattice_transfer.h ./lattice/Lattice_unary.h ./lattice/Lattice_peekpoke.h ./lattice/Lattice_coordinate.h ./lattice/Lattice_comparison.h ./lattice/Lattice_overload.h ./lattice/Lattice_reality.h ./lattice/Lattice_local.h ./lattice/Lattice_conformable.h ./lattice/Lattice_where.h ./lattice/Lattice_comparison_utils.h ./lattice/Lattice_arith.h ./lattice/Lattice_base.h ./lattice/Lattice_ET.h ./lattice/Lattice_transpose.h ./lattice/Lattice_trace.h ./Stencil.h ./tensors/Tensor_arith_sub.h ./tensors/Tensor_exp.h ./tensors/Tensor_arith_mul.h ./tensors/Tensor_class.h ./tensors/Tensor_logical.h ./tensors/Tensor_transpose.h ./tensors/Tensor_arith_mac.h ./tensors/Tensor_arith_scalar.h ./tensors/Tensor_reality.h ./tensors/Tensor_trace.h ./tensors/Tensor_index.h ./tensors/Tensor_arith_add.h ./tensors/Tensor_outer.h ./tensors/Tensor_inner.h ./tensors/Tensor_traits.h ./tensors/Tensor_Ta.h ./tensors/Tensor_unary.h ./tensors/Tensor_determinant.h ./tensors/Tensor_arith.h ./tensors/Tensor_extract_merge.h ./Communicator.h ./Cartesian.h ./parallelIO/NerscIO.h ./qcd/QCD.h ./qcd/hmc/integrators/Integrator_base.h ./qcd/hmc/HMC.h ./qcd/utils/SpaceTimeGrid.h ./qcd/utils/SUn.h ./qcd/utils/LinalgUtils.h ./qcd/utils/CovariantCshift.h ./qcd/utils/WilsonLoops.h ./qcd/action/ActionBase.h ./qcd/action/gauge/WilsonGaugeAction.h ./qcd/action/Actions.h ./qcd/action/fermion/CayleyFermion5D.h ./qcd/action/fermion/ScaledShamirFermion.h ./qcd/action/fermion/MobiusFermion.h ./qcd/action/fermion/OverlapWilsonContfracTanhFermion.h ./qcd/action/fermion/PartialFractionFermion5D.h ./qcd/action/fermion/ShamirZolotarevFermion.h ./qcd/action/fermion/FermionOperator.h ./qcd/action/fermion/WilsonFermion5D.h ./qcd/action/fermion/WilsonCompressor.h ./qcd/action/fermion/OverlapWilsonPartialFractionZolotarevFermion.h ./qcd/action/fermion/WilsonKernels.h ./qcd/action/fermion/DomainWallFermion.h ./qcd/action/fermion/OverlapWilsonPartialFractionTanhFermion.h ./qcd/action/fermion/OverlapWilsonContfracZolotarevFermion.h ./qcd/action/fermion/MobiusZolotarevFermion.h ./qcd/action/fermion/g5HermitianLinop.h ./qcd/action/fermion/OverlapWilsonCayleyTanhFermion.h ./qcd/action/fermion/WilsonFermion.h ./qcd/action/fermion/ContinuedFractionFermion5D.h ./qcd/action/fermion/OverlapWilsonCayleyZolotarevFermion.h ./qcd/spin/TwoSpinor.h ./qcd/spin/Dirac.h ./cshift/Cshift_common.h ./cshift/Cshift_none.h ./cshift/Cshift_mpi.h ./Simd.h ./GridConfig.h ./cartesian/Cartesian_base.h ./cartesian/Cartesian_red_black.h ./cartesian/Cartesian_full.h ./AlignedAllocator.h ./Lattice.h ./Old/Tensor_poke.h ./Old/Tensor_peek.h ./Threads.h ./Grid.h ./algorithms/Preconditioner.h ./algorithms/iterative/ConjugateResidual.h ./algorithms/iterative/PrecGeneralisedConjugateResidual.h ./algorithms/iterative/ConjugateGradientMultiShift.h ./algorithms/iterative/SchurRedBlack.h ./algorithms/iterative/NormalEquations.h ./algorithms/iterative/ConjugateGradient.h ./algorithms/iterative/AdefGeneric.h ./algorithms/approx/Chebyshev.h ./algorithms/approx/Zolotarev.h ./algorithms/approx/MultiShiftFunction.h ./algorithms/approx/bigfloat.h ./algorithms/approx/bigfloat_double.h ./algorithms/approx/Remez.h ./algorithms/LinearOperator.h ./algorithms/SparseMatrix.h ./algorithms/CoarsenedMatrix.h ./stencil/Lebesgue.h
 
 CCFILES=./qcd/hmc/integrators/Integrator_base.cc ./qcd/hmc/HMC.cc ./qcd/utils/SpaceTimeGrid.cc ./qcd/action/fermion/WilsonKernels.cc ./qcd/action/fermion/PartialFractionFermion5D.cc ./qcd/action/fermion/CayleyFermion5D.cc ./qcd/action/fermion/WilsonKernelsHand.cc ./qcd/action/fermion/WilsonFermion.cc ./qcd/action/fermion/ContinuedFractionFermion5D.cc ./qcd/action/fermion/WilsonFermion5D.cc ./qcd/spin/Dirac.cc ./GridInit.cc ./algorithms/approx/MultiShiftFunction.cc ./algorithms/approx/Remez.cc ./algorithms/approx/Zolotarev.cc ./stencil/Lebesgue.cc ./stencil/Stencil_common.cc

lib/qcd/QCD.h

@@ -303,35 +303,35 @@ namespace QCD {
 		      const Lattice<decltype(peekIndex<ColourIndex>(lhs._odata[0],0))> & rhs,
 		      int i)
     {
-      pokeIndex<ColourIndex>(lhs,rhs,i);
+      PokeIndex<ColourIndex>(lhs,rhs,i);
     }
     template<class vobj> 
       void pokeColour(Lattice<vobj> &lhs,
 		      const Lattice<decltype(peekIndex<ColourIndex>(lhs._odata[0],0,0))> & rhs,
 		      int i,int j)
     {
-      pokeIndex<ColourIndex>(lhs,rhs,i,j);
+      PokeIndex<ColourIndex>(lhs,rhs,i,j);
     }
     template<class vobj> 
       void pokeSpin(Lattice<vobj> &lhs,
 		      const Lattice<decltype(peekIndex<SpinIndex>(lhs._odata[0],0))> & rhs,
 		      int i)
     {
-      pokeIndex<SpinIndex>(lhs,rhs,i);
+      PokeIndex<SpinIndex>(lhs,rhs,i);
     }
     template<class vobj> 
       void pokeSpin(Lattice<vobj> &lhs,
 		      const Lattice<decltype(peekIndex<SpinIndex>(lhs._odata[0],0,0))> & rhs,
 		      int i,int j)
     {
-      pokeIndex<SpinIndex>(lhs,rhs,i,j);
+      PokeIndex<SpinIndex>(lhs,rhs,i,j);
     }
     template<class vobj> 
       void pokeLorentz(Lattice<vobj> &lhs,
 		      const Lattice<decltype(peekIndex<LorentzIndex>(lhs._odata[0],0))> & rhs,
 		      int i)
     {
-      pokeIndex<LorentzIndex>(lhs,rhs,i);
+      PokeIndex<LorentzIndex>(lhs,rhs,i);
     }
 
     //////////////////////////////////////////////

lib/qcd/action/ActionBase.h

@@ -1,19 +1,20 @@
-#ifndef QCD_GAUGE_ACTION_BASE
+#ifndef QCD_ACTION_BASE
-#define QCD_GAUGE_ACTION_BASE
+#define QCD_ACTION_BASE
 namespace Grid {
 namespace QCD{
 
 template<class GaugeField>
-class GaugeActionBase { // derive this from TermInAction?
+class Action { 
 
  public:
+  virtual void  init(const GaugeField &U)                      = 0;
   virtual RealD S(const GaugeField &U)                        = 0;  // evaluate the action
   virtual void  deriv(const GaugeField &U,GaugeField & dSdU ) = 0;  // evaluate the action derivative
   virtual void  staple(const GaugeField &U,GaugeField & dSdU ) = 0;  // evaluate the action derivative
-  virtual void refresh(const GaugeField & ) {}; 
+  //virtual void  refresh(const GaugeField & ) {}                ; 
   // Boundary conditions?
   // Heatbath?
-  virtual ~GaugeActionBase() {};
+  virtual ~Action() {};
 };
 
 }}

lib/qcd/action/Actions.h

@@ -15,8 +15,17 @@
 ////////////////////////////////////////////
 // Abstract base interface
 ////////////////////////////////////////////
+#include <qcd/action/ActionBase.h>
+
 #include <qcd/action/fermion/FermionOperator.h>
 
+
+////////////////////////////////////////////
+// Gauge Actions
+////////////////////////////////////////////
+#include <qcd/action/gauge/WilsonGaugeAction.h>
+
+
 ////////////////////////////////////////////
 // Utility functions
 ////////////////////////////////////////////

lib/qcd/action/gauge/WilsonGaugeAction.h

@@ -1,24 +1,38 @@
 #ifndef QCD_WILSON_GAUGE_ACTION_H
 #define QCD_WILSON_GAUGE_ACTION_H
 
-////////////////////////////////////////////////////////////////////////
+namespace Grid{
-// Wilson Gauge Action .. should I template the Nc etc..
+  namespace QCD{
-////////////////////////////////////////////////////////////////////////
+    
-template<class GaugeField,class MatrixField>
+    ////////////////////////////////////////////////////////////////////////
-class WilsonGaugeAction : public GaugeActionBase<GaugeField> {
+    // Wilson Gauge Action .. should I template the Nc etc..
- public:
+    ////////////////////////////////////////////////////////////////////////
-  virtual RealD S(const GaugeField &U) {
+    template<class GaugeField,class MatrixField>
-    return WilsonLoops<MatrixField,GaugeField>::sumPlaquette(U);
+      class WilsonGaugeAction : public Action<GaugeField> {
-  };
+    private:
-  virtual RealD deriv(GaugeField &U,GaugeField & dSdU ) {
+      RealD beta;
-    WilsonLoops<MatrixField,GaugeField>::Staple(dSdU,U,mu);
+    public:
-  };
+    WilsonGaugeAction(RealD b):beta(b){};
-  virtual void  staple(const MatrixField &stap,GaugeField & U,int mu ) {
+      
-    WilsonLoops<MatrixField,GaugeField>::Staple(stap,U,mu);
+      virtual void init(const GaugeField &U) {//FIXME
-  };
+      };
-};
+      
-
+      virtual RealD S(const GaugeField &U) {
-
+	return WilsonLoops<MatrixField,GaugeField>::sumPlaquette(U);
-#endif
+      };
+      virtual void deriv(const GaugeField &U,GaugeField & dSdU) {
+	//FIXME loop on directions
+	MatrixField dSdU_mu(U._grid);
+	WilsonLoops<MatrixField,GaugeField>::Staple(dSdU_mu,U,0);
+      };
+      virtual void  staple(const GaugeField &stap,GaugeField & U) {
+	//FIXME loop on directions
+	MatrixField stap_mu(U._grid);
+	WilsonLoops<MatrixField,GaugeField>::Staple(stap_mu,U,0);
+      };
+    };
+    
+  }
+}
 
 #endif

lib/qcd/hmc/HMC.cc

@@ -18,81 +18,7 @@ namespace Grid{
       }
 
 
-    /////////////////////////////////////////////////////////////////
+   
-
-    HybridMonteCarlo::HybridMonteCarlo(GridSerialRNG& R):RNG(R){
-	//FIXME
-      }
-
-
-    void HybridMonteCarlo::evolve(LatticeColourMatrix& Uin)const{
-      Real DeltaH;
-      Real timer;
-      
-      // Thermalizations
-      for(int iter=1; iter <= Params.ThermalizationSteps; ++iter){
-	std::cout << "-- # Thermalization step = "<< iter <<  "\n";
-	
-	DeltaH = evolve_step(Uin);
-	
-	std::cout<< "[Timing] Trajectory time (s) : "<< timer/1000.0 << "\n";
-	std::cout<< "dH = "<< DeltaH << "\n";
-	
-	// Update matrix
-	//Uin = md_->get_U();  //accept every time
-      }
-
-      // Actual updates
-      for(int iter=Params.StartingConfig; 
-	  iter < Params.Nsweeps+Params.StartingConfig; ++iter){
-	std::cout << "-- # Sweep = "<< iter <<  "\n";
-		
-	DeltaH = evolve_step(Uin);
-		
-	if(metropolis_test(DeltaH))  {};//Uin = md_->get_U();
-	// need sync?	
-      }
-
-      
-    }
-
-
-    RealD HybridMonteCarlo::evolve_step(LatticeColourMatrix& Uin)const{
-      /*
-      md_->init(Uin,RNG);     // set U and initialize P and phi's 
-      RealD H0 = md_->S();     // current state            
-      std::cout<<"Total H_before = "<< H0 << "\n";
-      
-      md_->integrator();
-      
-      RealD H1 = md_->calc_H();     // updated state            
-      std::cout<<"Total H_after = "<< H1 << "\n";
-      
-      return (H1-H0);
-      */
-      return 0;
-    }
-
-
-
-
-    bool HybridMonteCarlo::metropolis_test(const RealD DeltaH)const{
-      RealD rn_test;
-      RealD prob = std::exp(-DeltaH);
-      random(RNG,rn_test);
-      
-      std::cout<< "--------------------------------------------\n";
-      std::cout<< "dH = "<<DeltaH << "  Random = "<< rn_test 
-	       << "\nAcc. Probability = " << ((prob<1.0)? prob: 1.0)<< "   ";
-      
-      if((prob >1.0) || (rn_test <= prob)){       // accepted
-	std::cout <<"-- ACCEPTED\n";
-	return true;
-      } else {                               // rejected
-	std::cout <<"-- REJECTED\n";
-	return false;
-      }
-    }
 
   }
 }

lib/qcd/hmc/HMC.h

@@ -25,21 +25,91 @@ namespace Grid{
       HMCparameters();
     };
     
- 
+    template <class IntegType> 
     class HybridMonteCarlo{
       const HMCparameters Params;
-      GridSerialRNG& RNG; 
+      GridSerialRNG sRNG;
-      // FIXME need the integrator
+      GridParallelRNG pRNG;
+      std::unique_ptr< Integrator<IntegType> > MD;
       
-      bool metropolis_test(const RealD DeltaH)const;
+
-      RealD evolve_step(LatticeColourMatrix&)const;
+      bool metropolis_test(const RealD DeltaH){
+	RealD rn_test;
+	RealD prob = std::exp(-DeltaH);
+	random(sRNG,rn_test);
+      
+	std::cout<< "--------------------------------------------\n";
+	std::cout<< "dH = "<<DeltaH << "  Random = "<< rn_test 
+		 << "\nAcc. Probability = " << ((prob<1.0)? prob: 1.0)<< "   ";
+      
+	if((prob >1.0) || (rn_test <= prob)){       // accepted
+	  std::cout <<"-- ACCEPTED\n";
+	  return true;
+	} else {                               // rejected
+	  std::cout <<"-- REJECTED\n";
+	  return false;
+	}
+      }
+
+      RealD evolve_step(LatticeColourMatrix& Uin){
+     
+	MD->init(Uin,pRNG);     // set U and initialize P and phi's 
+	RealD H0 = MD->S();     // current state            
+	std::cout<<"Total H_before = "<< H0 << "\n";
+      
+	MD->integrate(0);
+      
+	RealD H1 = MD->S();     // updated state            
+	std::cout<<"Total H_after = "<< H1 << "\n";
+      
+	return (H1-H0);
+      }
+
       
       
     public:
-      HybridMonteCarlo(GridSerialRNG&);
+      HybridMonteCarlo(Integrator<IntegType>& MolDyn, GridBase* grid):MD(&MolDyn),pRNG(grid){
+	//FIXME
+
+	// initialize RNGs
+	sRNG.SeedRandomDevice();
+	pRNG.SeedRandomDevice();
+      }
       ~HybridMonteCarlo(){};
 
-      void evolve(LatticeColourMatrix& Uin)const;
+
+
+      void evolve(LatticeColourMatrix& Uin){
+	Real DeltaH;
+	Real timer;
+      
+	// Thermalizations
+	for(int iter=1; iter <= Params.ThermalizationSteps; ++iter){
+	  std::cout << "-- # Thermalization step = "<< iter <<  "\n";
+	
+	  DeltaH = evolve_step(Uin);
+	
+	  std::cout<< "[Timing] Trajectory time (s) : "<< timer/1000.0 << "\n";
+	  std::cout<< "dH = "<< DeltaH << "\n";
+	
+	  // Update matrix
+	  Uin = MD->get_U();  //accept every time
+	}
+
+	// Actual updates
+	for(int iter=Params.StartingConfig; 
+	    iter < Params.Nsweeps+Params.StartingConfig; ++iter){
+	  std::cout << "-- # Sweep = "<< iter <<  "\n";
+		
+	  DeltaH = evolve_step(Uin);
+		
+	  if(metropolis_test(DeltaH))  Uin = MD->get_U();
+	  // need sync?	
+	}
+
+      
+      }
+
 
     };
     

lib/qcd/hmc/integrators/Integrator_base.cc

@@ -3,18 +3,22 @@
   @brief utilities for MD including funcs to generate initial HMC momentum
  */
 
-
 #include <Grid.h>
 
-static const double sq3i = 1.0/sqrt(3.0);
-
 namespace Grid{
   namespace QCD{
 
-    
+    void MDutils::generate_momenta(LatticeLorentzColourMatrix& P,GridParallelRNG& pRNG){
+      // for future support of different groups
+      MDutils::generate_momenta_su3(P, pRNG);
+    }
 
-    void MDutils::generate_momenta_su3(LatticeColourMatrix& P,GridParallelRNG& pRNG){
+    void MDutils::generate_momenta_su3(LatticeLorentzColourMatrix& P,GridParallelRNG& pRNG){
-      SU3::GaussianLieAlgebraMatrix(pRNG, P);
+      LatticeColourMatrix Pmu(P._grid);
+      for(int mu=0;mu<Nd;mu++){
+	SU3::GaussianLieAlgebraMatrix(pRNG, Pmu);
+	pokeLorentz(P, Pmu, mu);
+      }
     }
 
 

lib/qcd/hmc/integrators/Integrator_base.h

@@ -9,40 +9,71 @@
 #ifndef INTEGRATOR_INCLUDED
 #define INTEGRATOR_INCLUDED
 
-class Action;
-class RandNum;
 class Observer;
 
-typedef std::vector<Action*> ActionLevel;
-typedef std::vector<ActionLevel> ActionSet;
-typedef std::vector<Observer*> ObserverList;
 
 /*! @brief Abstract base class for Molecular Dynamics management */
 
 namespace Grid{
   namespace QCD{
+
+    typedef Action<LatticeLorentzColourMatrix>*  ActPtr; // now force the same size as the rest of the code
+    typedef std::vector<ActPtr> ActionLevel;
+    typedef std::vector<ActionLevel> ActionSet;
+    typedef std::vector<Observer*> ObserverList;
     
+    
+    class Integrator2MN{
+      const double lambda = 0.1931833275037836;
+      void step (LatticeColourMatrix&, LatticeColourMatrix&, 
+		 int, std::vector<int>&);
+      
+    };
+    
+    class IntegratorLeapFrog{
+      void step (LatticeColourMatrix&, LatticeColourMatrix&,
+		 int, std::vector<int>&);
+    };
+    
+
+
+    template< class IntegratorPolicy >
     class Integrator{
     private:
-      virtual void update_P(int lv,double ep) = 0;
+      int Nexp;
-      virtual void update_U(double ep) = 0;
+      int MDsteps;  // number of outer steps
+      RealD trajL;  // trajectory length 
+      RealD stepsize;
+
+      const std::vector<int> Nrel; // relative steps per level
+      const ActionSet as;
+      ObserverList observers;
+      //      LatticeColourMatrix* const U;  // is shared among all actions - or use a singleton...
+      LatticeColourMatrix P;
       
-      virtual void register_observers() = 0;
+      IntegratorPolicy TheIntegrator;// contains parameters too
-      virtual void notify_observers() = 0;
+      void update_P(int lv,double ep);
+      void update_U(double ep);
+      
+      void register_observers();
+      void notify_observers();
+      void integrator_step(int level ,std::vector<Integer>& clock);
+
       
     public:
-      virtual ~Integrator(){}
+      Integrator(int Nexp_, int MDsteps_, RealD trajL_,
-      virtual void init(const LatticeColourMatrix&,
+		 ActionSet& Aset, ObserverList obs):as(Aset), observers(obs){};
-			const GridParallelRNG& RNG)=0;
+      ~Integrator(){}
-      virtual double S()const =0;
+      void init(LatticeLorentzColourMatrix&,
-      virtual void integrate(int level) =0;
+		GridParallelRNG&);
-      virtual const LatticeColourMatrix get_U() const =0;
+      double S();
-      
+      void integrate(int level);
-      void generate_momenta(LatticeColourMatrix& P,const RandNum& rand);
+      LatticeColourMatrix get_U();
     };
     
     namespace MDutils{
-      void generate_momenta_su3(LatticeColourMatrix& P,GridParallelRNG& RNG);
+      void generate_momenta(LatticeLorentzColourMatrix&,GridParallelRNG&);
+      void generate_momenta_su3(LatticeLorentzColourMatrix&,GridParallelRNG&);
     }
     
   }

lib/qcd/utils/SUn.h

@@ -574,14 +574,14 @@ Note that in step D setting B ~ X - A and using B in place of A in step E will g
     LatticeMatrix Umu(out._grid);
     for(int mu=0;mu<Nd;mu++){
       LieRandomize(pRNG,Umu,0.01);
-      PokeLorentz(out,Umu,mu);
+      pokeLorentz(out,Umu,mu);
     }
   }
   static void ColdConfiguration(GridParallelRNG &pRNG,LatticeGaugeField &out){
     LatticeMatrix Umu(out._grid);
     Umu=1.0;
     for(int mu=0;mu<Nd;mu++){
-      PokeLorentz(out,Umu,mu);
+      pokeLorentz(out,Umu,mu);
     }
   }
 

tests/Make.inc

@@ -1,5 +1,5 @@
 
-bin_PROGRAMS = Test_cayley_cg Test_cayley_coarsen_support Test_cayley_even_odd Test_cayley_ldop_cr Test_cf_coarsen_support Test_cf_cr_unprec Test_contfrac_cg Test_contfrac_even_odd Test_cshift Test_cshift_red_black Test_dwf_cg_prec Test_dwf_cg_schur Test_dwf_cg_unprec Test_dwf_cr_unprec Test_dwf_even_odd Test_dwf_fpgcr Test_dwf_hdcr Test_gamma Test_GaugeAction Test_lie_generators Test_main Test_multishift_sqrt Test_nersc_io Test_quenched_update Test_remez Test_rng Test_rng_fixed Test_simd Test_stencil Test_wilson_cg_prec Test_wilson_cg_schur Test_wilson_cg_unprec Test_wilson_cr_unprec Test_wilson_even_odd
+bin_PROGRAMS = Test_cayley_cg Test_cayley_coarsen_support Test_cayley_even_odd Test_cayley_ldop_cr Test_cf_coarsen_support Test_cf_cr_unprec Test_contfrac_cg Test_contfrac_even_odd Test_cshift Test_cshift_red_black Test_dwf_cg_prec Test_dwf_cg_schur Test_dwf_cg_unprec Test_dwf_cr_unprec Test_dwf_even_odd Test_dwf_fpgcr Test_dwf_hdcr Test_gamma Test_GaugeAction Test_hmc Test_lie_generators Test_main Test_multishift_sqrt Test_nersc_io Test_quenched_update Test_remez Test_rng Test_rng_fixed Test_simd Test_stencil Test_wilson_cg_prec Test_wilson_cg_schur Test_wilson_cg_unprec Test_wilson_cr_unprec Test_wilson_even_odd
 
 
 Test_cayley_cg_SOURCES=Test_cayley_cg.cc
@@ -78,6 +78,10 @@ Test_GaugeAction_SOURCES=Test_GaugeAction.cc
 Test_GaugeAction_LDADD=-lGrid
 
 
+Test_hmc_SOURCES=Test_hmc.cc
+Test_hmc_LDADD=-lGrid
+
+
 Test_lie_generators_SOURCES=Test_lie_generators.cc
 Test_lie_generators_LDADD=-lGrid
 

tests/Test_hmc.cc

@@ -0,0 +1,36 @@
+#include "Grid.h"
+
+
+using namespace std;
+using namespace Grid;
+using namespace Grid::QCD;
+
+int main (int argc, char ** argv)
+{
+  Grid_init(&argc,&argv);
+
+  std::vector<int> latt_size   = GridDefaultLatt();
+  std::vector<int> simd_layout = GridDefaultSimd(4,vComplex::Nsimd());
+  std::vector<int> mpi_layout  = GridDefaultMpi();
+
+  latt_size.resize(4);
+
+  latt_size[0] = 8;
+  latt_size[1] = 8;
+  latt_size[2] = 8;
+  latt_size[3] = 8;
+  double volume = latt_size[0]*latt_size[1]*latt_size[2]*latt_size[3];
+  
+  GridCartesian           Fine(latt_size,simd_layout,mpi_layout);
+  GridParallelRNG       FineRNG(&Fine);
+  GridSerialRNG       SerialRNG;
+  FineRNG.SeedRandomDevice();
+  
+  WilsonGaugeAction<LatticeLorentzColourMatrix, LatticeColourMatrix> Waction(6.0);
+  //Collect actions
+  ActionLevel Level;
+  Level.push_back(&Waction);
+  
+
+  //Integrator<IntegratorLeapFrog(12,10,1.0);
+}