diff --git a/lib/qcd/QCD.h b/lib/qcd/QCD.h
index a318141d..a61102a3 100644
--- a/lib/qcd/QCD.h
+++ b/lib/qcd/QCD.h
@@ -252,47 +252,47 @@ namespace QCD {
     //////////////////////////////////////////////////////////////////////////////
 
     //spin
-    template<class vobj> auto peekSpin(const vobj &rhs,int i) -> decltype(peekIndex<SpinIndex>(rhs,0))
+    template<class vobj> auto peekSpin(const vobj &rhs,int i) -> decltype(PeekIndex<SpinIndex>(rhs,0))
     {
-      return peekIndex<SpinIndex>(rhs,i);
+      return PeekIndex<SpinIndex>(rhs,i);
     }
-    template<class vobj> auto peekSpin(const vobj &rhs,int i,int j) -> decltype(peekIndex<SpinIndex>(rhs,0,0))
+    template<class vobj> auto peekSpin(const vobj &rhs,int i,int j) -> decltype(PeekIndex<SpinIndex>(rhs,0,0))
     {
-      return peekIndex<SpinIndex>(rhs,i,j);
+      return PeekIndex<SpinIndex>(rhs,i,j);
     }
-    template<class vobj> auto peekSpin(const Lattice<vobj> &rhs,int i) -> decltype(peekIndex<SpinIndex>(rhs,0))
+    template<class vobj> auto peekSpin(const Lattice<vobj> &rhs,int i) -> decltype(PeekIndex<SpinIndex>(rhs,0))
     {
-      return peekIndex<SpinIndex>(rhs,i);
+      return PeekIndex<SpinIndex>(rhs,i);
     }
-    template<class vobj> auto peekSpin(const Lattice<vobj> &rhs,int i,int j) -> decltype(peekIndex<SpinIndex>(rhs,0,0))
+    template<class vobj> auto peekSpin(const Lattice<vobj> &rhs,int i,int j) -> decltype(PeekIndex<SpinIndex>(rhs,0,0))
     {
-      return peekIndex<SpinIndex>(rhs,i,j);
+      return PeekIndex<SpinIndex>(rhs,i,j);
     }
     //colour
-    template<class vobj> auto peekColour(const vobj &rhs,int i) -> decltype(peekIndex<ColourIndex>(rhs,0))
+    template<class vobj> auto peekColour(const vobj &rhs,int i) -> decltype(PeekIndex<ColourIndex>(rhs,0))
     {
-      return peekIndex<ColourIndex>(rhs,i);
+      return PeekIndex<ColourIndex>(rhs,i);
     }
-    template<class vobj> auto peekColour(const vobj &rhs,int i,int j) -> decltype(peekIndex<ColourIndex>(rhs,0,0))
+    template<class vobj> auto peekColour(const vobj &rhs,int i,int j) -> decltype(PeekIndex<ColourIndex>(rhs,0,0))
     {
-      return peekIndex<ColourIndex>(rhs,i,j);
+      return PeekIndex<ColourIndex>(rhs,i,j);
     }
-    template<class vobj> auto peekColour(const Lattice<vobj> &rhs,int i) -> decltype(peekIndex<ColourIndex>(rhs,0))
+    template<class vobj> auto peekColour(const Lattice<vobj> &rhs,int i) -> decltype(PeekIndex<ColourIndex>(rhs,0))
     {
-      return peekIndex<ColourIndex>(rhs,i);
+      return PeekIndex<ColourIndex>(rhs,i);
     }
-    template<class vobj> auto peekColour(const Lattice<vobj> &rhs,int i,int j) -> decltype(peekIndex<ColourIndex>(rhs,0,0))
+    template<class vobj> auto peekColour(const Lattice<vobj> &rhs,int i,int j) -> decltype(PeekIndex<ColourIndex>(rhs,0,0))
     {
-      return peekIndex<ColourIndex>(rhs,i,j);
+      return PeekIndex<ColourIndex>(rhs,i,j);
     }
     //lorentz
-    template<class vobj> auto peekLorentz(const vobj &rhs,int i) -> decltype(peekIndex<LorentzIndex>(rhs,0))
+    template<class vobj> auto peekLorentz(const vobj &rhs,int i) -> decltype(PeekIndex<LorentzIndex>(rhs,0))
     {
-      return peekIndex<LorentzIndex>(rhs,i);
+      return PeekIndex<LorentzIndex>(rhs,i);
     }
-    template<class vobj> auto peekLorentz(const Lattice<vobj> &rhs,int i) -> decltype(peekIndex<LorentzIndex>(rhs,0))
+    template<class vobj> auto peekLorentz(const Lattice<vobj> &rhs,int i) -> decltype(PeekIndex<LorentzIndex>(rhs,0))
     {
-      return peekIndex<LorentzIndex>(rhs,i);
+      return PeekIndex<LorentzIndex>(rhs,i);
     }
 
     //////////////////////////////////////////////
diff --git a/lib/qcd/action/ActionBase.h b/lib/qcd/action/ActionBase.h
index e3b13f0e..51531750 100644
--- a/lib/qcd/action/ActionBase.h
+++ b/lib/qcd/action/ActionBase.h
@@ -7,10 +7,10 @@ template<class GaugeField>
 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  init(const GaugeField &U, GridParallelRNG& pRNG) = 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 & ) {}                ; 
   // Boundary conditions?
   // Heatbath?
diff --git a/lib/qcd/action/gauge/WilsonGaugeAction.h b/lib/qcd/action/gauge/WilsonGaugeAction.h
index 7726db80..c6ddd4a3 100644
--- a/lib/qcd/action/gauge/WilsonGaugeAction.h
+++ b/lib/qcd/action/gauge/WilsonGaugeAction.h
@@ -14,8 +14,7 @@ namespace Grid{
     public:
     WilsonGaugeAction(RealD b):beta(b){};
       
-      virtual void init(const GaugeField &U) {//FIXME
-      };
+      virtual void init(const GaugeField &U, GridParallelRNG& pRNG) {};
       
       virtual RealD S(const GaugeField &U) {
 	return WilsonLoops<MatrixField,GaugeField>::sumPlaquette(U);
diff --git a/lib/qcd/hmc/HMC.h b/lib/qcd/hmc/HMC.h
index e4d8740c..e783bee0 100644
--- a/lib/qcd/hmc/HMC.h
+++ b/lib/qcd/hmc/HMC.h
@@ -1,6 +1,6 @@
 //--------------------------------------------------------------------
 /*! @file HMC.h
- * @brief Declaration of classes for HybridMonteCarlo update
+ * @brief Declaration of classes for Hybrid Monte Carlo update
  *
  * @author Guido Cossu
  */
@@ -30,7 +30,7 @@ namespace Grid{
       const HMCparameters Params;
       GridSerialRNG sRNG;
       GridParallelRNG pRNG;
-      std::unique_ptr< Integrator<IntegType> > MD;
+      Integrator<IntegType>& MD;
       
 
       bool metropolis_test(const RealD DeltaH){
@@ -51,13 +51,13 @@ namespace Grid{
 	}
       }
 
-      RealD evolve_step(LatticeColourMatrix& Uin){
-     
-	MD->init(Uin,pRNG);     // set U and initialize P and phi's 
+      RealD evolve_step(LatticeLorentzColourMatrix& U){
+
+	MD->init(U,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);
+	MD->integrate(U,0);
       
 	RealD H1 = MD->S();     // updated state            
 	std::cout<<"Total H_after = "<< H1 << "\n";
@@ -68,7 +68,10 @@ namespace Grid{
       
       
     public:
-      HybridMonteCarlo(Integrator<IntegType>& MolDyn, GridBase* grid):MD(&MolDyn),pRNG(grid){
+    HybridMonteCarlo(HMCparameters Pms, 
+		     Integrator<IntegType>& MolDyn, 
+		     GridBase* grid):
+      Params(Pms),MD(MolDyn),pRNG(grid){
 	//FIXME
 
 	// initialize RNGs
@@ -79,7 +82,7 @@ namespace Grid{
 
 
 
-      void evolve(LatticeColourMatrix& Uin){
+      void evolve(LatticeLorentzColourMatrix& Uin){
 	Real DeltaH;
 	Real timer;
       
@@ -90,20 +93,21 @@ namespace Grid{
 	  DeltaH = evolve_step(Uin);
 	
 	  std::cout<< "[Timing] Trajectory time (s) : "<< timer/1000.0 << "\n";
-	  std::cout<< "dH = "<< DeltaH << "\n";
+	  std::cout<< " dH = "<< DeltaH << "\n";
 	
 	  // Update matrix
-	  Uin = MD->get_U();  //accept every time
+	  //Uin = MD->get_U();  //accept every time
 	}
 
-	// Actual updates
+	// Actual updates (evolve a copy Ucopy then copy back eventually)
+	LatticeLorentzColourMatrix Ucopy(Uin._grid);
 	for(int iter=Params.StartingConfig; 
 	    iter < Params.Nsweeps+Params.StartingConfig; ++iter){
 	  std::cout << "-- # Sweep = "<< iter <<  "\n";
+	  Ucopy = Uin;
+	  DeltaH = evolve_step(Ucopy);
 		
-	  DeltaH = evolve_step(Uin);
-		
-	  if(metropolis_test(DeltaH))  Uin = MD->get_U();
+	  if(metropolis_test(DeltaH))  Uin = Ucopy;
 	  // need sync?	
 	}
 
diff --git a/lib/qcd/hmc/integrators/Integrator_base.h b/lib/qcd/hmc/integrators/Integrator_base.h
index 61abf118..3f610958 100644
--- a/lib/qcd/hmc/integrators/Integrator_base.h
+++ b/lib/qcd/hmc/integrators/Integrator_base.h
@@ -9,6 +9,8 @@
 #ifndef INTEGRATOR_INCLUDED
 #define INTEGRATOR_INCLUDED
 
+#include <memory>
+
 class Observer;
 
 
@@ -22,59 +24,188 @@ namespace Grid{
     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 LeapFrog;
+
+
+    struct IntegratorParameters{
+      int Nexp;
+      int MDsteps;  // number of outer steps
+      RealD trajL;  // trajectory length 
+      RealD stepsize;
+
+      IntegratorParameters(int Nexp_,
+			   int MDsteps_, 
+			   RealD trajL_):
+	Nexp(Nexp_),MDsteps(MDsteps_),trajL(trajL_),stepsize(trajL/MDsteps){};
     };
+
+
+    namespace MDutils{
+      void generate_momenta(LatticeLorentzColourMatrix&,GridParallelRNG&);
+      void generate_momenta_su3(LatticeLorentzColourMatrix&,GridParallelRNG&);
+    }
+
     
-    class IntegratorLeapFrog{
-      void step (LatticeColourMatrix&, LatticeColourMatrix&,
-		 int, std::vector<int>&);
-    };
     
 
 
     template< class IntegratorPolicy >
     class Integrator{
     private:
-      int Nexp;
-      int MDsteps;  // number of outer steps
-      RealD trajL;  // trajectory length 
-      RealD stepsize;
-
-      const std::vector<int> Nrel; // relative steps per level
+      IntegratorParameters Params;
       const ActionSet as;
-      ObserverList observers;
-      //      LatticeColourMatrix* const U;  // is shared among all actions - or use a singleton...
-      LatticeColourMatrix P;
+      const std::vector<int> Nrel; //relative step size per level
+      //ObserverList observers; // not yet
+      std::unique_ptr<LatticeLorentzColourMatrix> P;
       
       IntegratorPolicy TheIntegrator;// contains parameters too
-      void update_P(int lv,double ep);
-      void update_U(double ep);
+
+
+      void update_P(LatticeLorentzColourMatrix&U, int level,double ep){
+	for(int a=0; a<as[level].size(); ++a){
+	  LatticeLorentzColourMatrix force(U._grid);
+	  as[level].at(a)->deriv(U,force);
+	  *P -= force*ep;
+	}
+      }
+
+
+      void update_U(LatticeLorentzColourMatrix&U, double ep){
+	//rewrite exponential to deal with the lorentz index?
+	LatticeColourMatrix Umu(U._grid);
+	LatticeColourMatrix Pmu(U._grid);
+	for (int mu = 0; mu < Nd; mu++){
+	  Umu=peekLorentz(U, mu);
+	  Pmu=peekLorentz(*P, mu);
+	  Umu = expMat(Pmu, Complex(ep, 0.0))*Umu;
+	}
+
+      }
       
       void register_observers();
       void notify_observers();
-      void integrator_step(int level ,std::vector<Integer>& clock);
+      
+      friend void IntegratorPolicy::step (LatticeLorentzColourMatrix& U, 
+				     int level, std::vector<int>& clock,
+				     Integrator<LeapFrog>* Integ);
+    public:
+      Integrator(IntegratorParameters Par,
+		 ActionSet& Aset, std::vector<int> Nrel_):
+	Params(Par),as(Aset),Nrel(Nrel_){
+	assert(as.size() == Nrel.size());
+      };
+      
+      ~Integrator(){}
+
+
+      //Initialization of momenta and actions
+      void init(LatticeLorentzColourMatrix& U,
+		GridParallelRNG& pRNG){
+	std::cout<< "Integrator init\n";
+	if (!P)
+	  P = new LatticeLorentzColourMatrix(U._grid);
+	MDutils::generate_momenta(*P,pRNG);
+	
+	for(int level=0; level< as.size(); ++level){
+	  for(int actionID=0; actionID<as.at(level).size(); ++actionID){
+	    as[level].at(actionID)->init(U, pRNG);
+	  }
+	}
+      }
 
       
-    public:
-      Integrator(int Nexp_, int MDsteps_, RealD trajL_,
-		 ActionSet& Aset, ObserverList obs):as(Aset), observers(obs){};
-      ~Integrator(){}
-      void init(LatticeLorentzColourMatrix&,
-		GridParallelRNG&);
-      double S();
-      void integrate(int level);
-      LatticeColourMatrix get_U();
+
+      RealD S(LatticeLorentzColourMatrix& U){
+	// Momenta
+	LatticeComplex Hloc = - trace((*P)*adj(*P));
+	Complex Hsum = sum(Hloc);
+	
+	RealD H = Hsum.real();
+
+	// Actions
+	for(int level=0; level<as.size(); ++level)
+	  for(int actionID=0; actionID<as.at(level).size(); ++actionID)
+	    H += as[level].at(actionID)->S(U);
+	
+	return H;
+      }
+      
+      
+
+      void integrate(LatticeLorentzColourMatrix& U, int level){
+	std::vector<int> clock;
+	clock.resize(as.size(),0);
+	for(int step=0; step< Params.MDsteps; ++step)   // MD step
+	  TheIntegrator.step(U,0,clock, *(this));
+      }
     };
     
-    namespace MDutils{
-      void generate_momenta(LatticeLorentzColourMatrix&,GridParallelRNG&);
-      void generate_momenta_su3(LatticeLorentzColourMatrix&,GridParallelRNG&);
-    }
+
+    class MinimumNorm2{
+      const double lambda = 0.1931833275037836;
+    public:
+      void step (LatticeLorentzColourMatrix& U, int level, std::vector<int>& clock);
+      
+    };
+    
+    class LeapFrog{
+    public:
+      void step (LatticeLorentzColourMatrix& U, 
+		 int level, std::vector<int>& clock,
+		 Integrator<LeapFrog>* Integ){
+	// cl  : current level
+	// fl  : final level
+	// eps : current step size
+	
+	int fl = Integ->as.size() -1;
+	double eps = Integ->Params.stepsize;
+	
+	// Get current level step size
+	for(int l=0; l<=level; ++l) eps/= Integ->Nrel[l];
+	
+	int fin = 1;
+	for(int l=0; l<=level; ++l) fin*= Integ->Nrel[l];
+	fin = 2*Integ->Params.MDsteps*fin - 1;
+	
+	for(int e=0; e<Integ->Nrel[level]; ++e){
+	  
+	  if(clock[level] == 0){    // initial half step
+	    Integ->update_P(U, level,eps/2);
+	    ++clock[level];
+	    for(int l=0; l<level;++l) std::cout<<"   ";
+	    std::cout<<"P "<< 0.5*clock[level] <<std::endl;
+	  }
+	  if(level == fl){          // lowest level
+	    Integ->update_U(U, eps);
+	    for(int l=0; l<level;++l) std::cout<<"   ";
+	    std::cout<<"U "<< 0.5*(clock[level]+1) <<std::endl;
+	  }else{                 // recursive function call
+	    step(U, level+1,clock, Integ);
+	  }
+	  if(clock[level] == fin){  // final half step
+	    Integ->update_P(U, level,eps/2);
+	    
+	    ++clock[level];
+	    for(int l=0; l<level;++l) std::cout<<"   ";
+	    std::cout<<"P "<< 0.5*clock[level] <<std::endl;
+	  }else{                  // bulk step
+	    Integ->update_P(U, level,eps);
+	    
+	    clock[level]+=2;
+	    for(int l=0; l<level;++l) std::cout<<"   ";
+	    std::cout<<"P "<< 0.5*clock[level] <<std::endl;
+	  }
+	}
+
+
+
+
+      }
+    };
+
+
+
+
     
   }
 }
diff --git a/tests/Test_hmc.cc b/tests/Test_hmc.cc
index 2a3d5d74..cbae27ab 100644
--- a/tests/Test_hmc.cc
+++ b/tests/Test_hmc.cc
@@ -22,15 +22,23 @@ int main (int argc, char ** argv)
   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();
   
+  // simplify template?
   WilsonGaugeAction<LatticeLorentzColourMatrix, LatticeColourMatrix> Waction(6.0);
-  //Collect actions
-  ActionLevel Level;
-  Level.push_back(&Waction);
-  
 
-  //Integrator<IntegratorLeapFrog(12,10,1.0);
+  //Collect actions
+  ActionLevel Level1;
+  Level1.push_back(&Waction);
+  ActionSet FullSet;
+  FullSet.push_back(Level1);
+
+  // Create integrator
+  IntegratorParameters MDpar(12,10,1.0);
+  std::vector<int> rel ={1};
+  Integrator<LeapFrog> MDleapfrog(MDpar, FullSet,rel);
+
+  // Create HMC
+  HMCparameters HMCpar;
+  HybridMonteCarlo<LeapFrog>  HMCrun(HMCpar, MDleapfrog, &Fine);
+
 }