diff --git a/lib/qcd/action/ActionBase.h b/lib/qcd/action/ActionBase.h
index 51531750..05edce86 100644
--- a/lib/qcd/action/ActionBase.h
+++ b/lib/qcd/action/ActionBase.h
@@ -10,7 +10,6 @@ class Action {
   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 c6ddd4a3..d5eb56e8 100644
--- a/lib/qcd/action/gauge/WilsonGaugeAction.h
+++ b/lib/qcd/action/gauge/WilsonGaugeAction.h
@@ -17,17 +17,28 @@ namespace Grid{
       virtual void init(const GaugeField &U, GridParallelRNG& pRNG) {};
       
       virtual RealD S(const GaugeField &U) {
-	return WilsonLoops<MatrixField,GaugeField>::sumPlaquette(U);
+	RealD plaq = WilsonLoops<MatrixField,GaugeField>::avgPlaquette(U);
+	std::cout << "Plaq : "<<plaq << "\n";
+	double vol = U._grid->gSites();
+	return beta*(1.0 -plaq)*(Nd*(Nd-1.0))*vol*0.5;
       };
       virtual void deriv(const GaugeField &U,GaugeField & dSdU) {
-	//FIXME loop on directions
+
+	//not optimal implementation FIXME
+	//extend Ta to include Lorentz indexes
+	RealD factor = 0.5*beta/RealD(Nc);
+	std::cout << "Debug force Wilson  "<< factor << "\n";
 	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);
+	MatrixField Umu(U._grid);
+	for (int mu=0; mu < Nd; mu++){
+	  Umu = PeekIndex<LorentzIndex>(U,mu);
+	  // Staple in direction mu
+	  WilsonLoops<MatrixField,GaugeField>::Staple(dSdU_mu,U,mu);
+	  std::cout << "Staple norm ["<<mu<<"] = "<< norm2(dSdU_mu) <<"\n";
+	  dSdU_mu = Ta(Umu*adj(dSdU_mu))*factor;
+	  std::cout << "Force norm ["<<mu<<"] = "<< norm2(dSdU_mu) << "  : Umu "<< norm2(Umu)<<"\n";
+	  pokeLorentz(dSdU, dSdU_mu, mu);
+	}
       };
     };
     
diff --git a/lib/qcd/hmc/HMC.cc b/lib/qcd/hmc/HMC.cc
index 859d6646..343fb32c 100644
--- a/lib/qcd/hmc/HMC.cc
+++ b/lib/qcd/hmc/HMC.cc
@@ -9,7 +9,7 @@ namespace Grid{
 	////////////////////////////// Default values
 	Nsweeps             = 10;
 	TotalSweeps         = 10;
-	ThermalizationSteps = 0;
+	ThermalizationSteps = 1;
 	StartingConfig      = 0;
 	SaveInterval        = 1;
 	Filename_prefix     = "Conf_";
diff --git a/lib/qcd/hmc/HMC.h b/lib/qcd/hmc/HMC.h
index e783bee0..83a66a48 100644
--- a/lib/qcd/hmc/HMC.h
+++ b/lib/qcd/hmc/HMC.h
@@ -53,13 +53,13 @@ namespace Grid{
 
       RealD evolve_step(LatticeLorentzColourMatrix& U){
 
-	MD->init(U,pRNG);     // set U and initialize P and phi's 
-	RealD H0 = MD->S();     // current state            
+	MD.init(U,pRNG);     // set U and initialize P and phi's 
+	RealD H0 = MD.S(U);     // current state            
 	std::cout<<"Total H_before = "<< H0 << "\n";
       
-	MD->integrate(U,0);
+	MD.integrate(U,0);
       
-	RealD H1 = MD->S();     // updated state            
+	RealD H1 = MD.S(U);     // updated state            
 	std::cout<<"Total H_after = "<< H1 << "\n";
       
 	return (H1-H0);
@@ -104,10 +104,11 @@ namespace Grid{
 	for(int iter=Params.StartingConfig; 
 	    iter < Params.Nsweeps+Params.StartingConfig; ++iter){
 	  std::cout << "-- # Sweep = "<< iter <<  "\n";
+	  
 	  Ucopy = Uin;
 	  DeltaH = evolve_step(Ucopy);
 		
-	  if(metropolis_test(DeltaH))  Uin = Ucopy;
+	  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 3f610958..acb420f7 100644
--- a/lib/qcd/hmc/integrators/Integrator_base.h
+++ b/lib/qcd/hmc/integrators/Integrator_base.h
@@ -25,7 +25,7 @@ namespace Grid{
     typedef std::vector<Observer*> ObserverList;
     
     class LeapFrog;
-
+   
 
     struct IntegratorParameters{
       int Nexp;
@@ -77,7 +77,9 @@ namespace Grid{
 	for (int mu = 0; mu < Nd; mu++){
 	  Umu=peekLorentz(U, mu);
 	  Pmu=peekLorentz(*P, mu);
+	  std::cout << "U norm ["<<mu<<"] = "<< norm2(Umu) << "  : Pmu "<< norm2(Pmu)<<"\n";
 	  Umu = expMat(Pmu, Complex(ep, 0.0))*Umu;
+	  pokeLorentz(U, Umu, mu);
 	}
 
       }
@@ -102,10 +104,15 @@ namespace Grid{
       void init(LatticeLorentzColourMatrix& U,
 		GridParallelRNG& pRNG){
 	std::cout<< "Integrator init\n";
-	if (!P)
-	  P = new LatticeLorentzColourMatrix(U._grid);
+	if (!P){
+	  std::unique_ptr<LatticeLorentzColourMatrix> Pnew(new LatticeLorentzColourMatrix(U._grid));
+	  P = std::move(Pnew);
+
+	}
 	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);
@@ -116,12 +123,19 @@ namespace Grid{
       
 
       RealD S(LatticeLorentzColourMatrix& U){
+	LatticeComplex Hloc(U._grid);
+	Hloc = zero;
 	// Momenta
-	LatticeComplex Hloc = - trace((*P)*adj(*P));
+	for (int mu=0; mu <Nd; mu++){
+	  LatticeColourMatrix Pmu = peekLorentz(*P, mu);
+	  Hloc -= trace(Pmu*adj(Pmu));
+	}
 	Complex Hsum = sum(Hloc);
 	
 	RealD H = Hsum.real();
 
+	std::cout << "H_p = "<< H << "\n";
+
 	// Actions
 	for(int level=0; level<as.size(); ++level)
 	  for(int actionID=0; actionID<as.at(level).size(); ++actionID)
@@ -136,7 +150,7 @@ namespace Grid{
 	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));
+	  TheIntegrator.step(U,0,clock, (this));
       }
     };
     
diff --git a/lib/simd/Grid_vector_unops.h b/lib/simd/Grid_vector_unops.h
index 37d1f6de..e4942937 100644
--- a/lib/simd/Grid_vector_unops.h
+++ b/lib/simd/Grid_vector_unops.h
@@ -94,6 +94,11 @@ namespace Grid {
   inline Grid_simd<S,V> rsqrt(const Grid_simd<S,V> &r) {
     return SimdApply(RSqrtRealFunctor<S>(),r);
   }
+  template < class Scalar > 
+  inline Scalar rsqrt(const Scalar &r) {
+    return (RSqrtRealFunctor<Scalar>(),r);
+  }
+
   template < class S, class V > 
   inline Grid_simd<S,V> cos(const Grid_simd<S,V> &r) {
     return SimdApply(CosRealFunctor<S>(),r);
diff --git a/lib/tensors/Tensor_Ta.h b/lib/tensors/Tensor_Ta.h
index e49ede35..cba07362 100644
--- a/lib/tensors/Tensor_Ta.h
+++ b/lib/tensors/Tensor_Ta.h
@@ -62,14 +62,16 @@ namespace Grid {
     {
       // need a check for the group type?
       iMatrix<vtype,N> ret(arg);
-      RealD nrm;
+      vtype nrm;
       vtype inner;
       for(int c1=0;c1<N;c1++){
 	zeroit(inner);	
 	for(int c2=0;c2<N;c2++)
 	  inner += innerProduct(ret._internal[c1][c2],ret._internal[c1][c2]);
 
-	nrm = 1.0/sqrt(Reduce(toReal(inner)));
+	//nrm = 1.0/sqrt(Reduce(toReal(inner)));
+	nrm = rsqrt(inner);
+
 	for(int c2=0;c2<N;c2++)
 	  ret._internal[c1][c2]*= nrm;
       
diff --git a/lib/tensors/Tensor_exp.h b/lib/tensors/Tensor_exp.h
index f8437116..6a514330 100644
--- a/lib/tensors/Tensor_exp.h
+++ b/lib/tensors/Tensor_exp.h
@@ -28,7 +28,9 @@ namespace Grid {
 	temp *= alpha/ComplexD(i);
 	temp = unit + temp*arg;
       }
+      
       return ProjectOnGroup(temp);//maybe not strictly necessary
+      
     }
 
 
diff --git a/tests/Test_hmc.cc b/tests/Test_hmc.cc
index cbae27ab..b292c57d 100644
--- a/tests/Test_hmc.cc
+++ b/tests/Test_hmc.cc
@@ -22,7 +22,13 @@ 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  pRNG(&Fine);
+  pRNG.SeedRandomDevice();
+  LatticeLorentzColourMatrix     U(&Fine);
+
+  SU3::ColdConfiguration(pRNG, U);
+ 
+
   // simplify template?
   WilsonGaugeAction<LatticeLorentzColourMatrix, LatticeColourMatrix> Waction(6.0);
 
@@ -33,12 +39,14 @@ int main (int argc, char ** argv)
   FullSet.push_back(Level1);
 
   // Create integrator
-  IntegratorParameters MDpar(12,10,1.0);
+  IntegratorParameters MDpar(12,50,1.0);
   std::vector<int> rel ={1};
   Integrator<LeapFrog> MDleapfrog(MDpar, FullSet,rel);
 
   // Create HMC
   HMCparameters HMCpar;
-  HybridMonteCarlo<LeapFrog>  HMCrun(HMCpar, MDleapfrog, &Fine);
+  HybridMonteCarlo<LeapFrog>  HMC(HMCpar, MDleapfrog, &Fine);
+
+  HMC.evolve(U);
 
 }
diff --git a/tests/Test_main.cc b/tests/Test_main.cc
index 295b1981..a98fbeb7 100644
--- a/tests/Test_main.cc
+++ b/tests/Test_main.cc
@@ -220,15 +220,19 @@ int main (int argc, char ** argv)
     std::cout << cm << std::endl;
 
 
-    cm = Exponentiate(cm, 1.0, 12);
+    cm = Exponentiate(cm, 2.0, 12);
     std::cout << cm << "  " << std::endl;
     Complex det = Determinant(cm);
     std::cout << "determinant: " << det <<  std::endl;
+    std::cout << "norm: " << norm2(cm) <<  std::endl;
 
     cm = ProjectOnGroup(cm);
     std::cout << cm << "  " << std::endl;
+    std::cout << "norm: " << norm2(cm) <<  std::endl;
     cm = ProjectOnGroup(cm);
     std::cout << cm << "  " << std::endl;
+    std::cout << "norm: " << norm2(cm) <<  std::endl;
+
 
     //    det = Determinant(cm);
     //    std::cout << "determinant: " << det <<  std::endl;
@@ -245,7 +249,24 @@ int main (int argc, char ** argv)
     trscMat = trace(scMat); // Trace
 
     // Exponentiate test
+    std::vector<int> mysite {0,0,0,0};
+    random(FineRNG,cMat);
+    cMat = Ta(cMat);
+    peekSite(cm, cMat, mysite);
+    std::cout << cm << "  " << std::endl;
+    cm = Exponentiate(cm, 1.0, 12);
+    std::cout << cm << "  " << std::endl;
+    std::cout << "norm: " << norm2(cm) <<  std::endl;
+
+
+    std::cout << "norm cMmat : " << norm2(cMat) <<  std::endl;
     cMat = expMat(cMat, ComplexD(1.0, 0.0));
+    std::cout << "norm expMat: " << norm2(cMat) <<  std::endl;
+    peekSite(cm, cMat, mysite);
+    std::cout << cm << "  " << std::endl;
+    std::cout << "determinant: " << Determinant(cm) <<  std::endl;
+    std::cout << "norm: " << norm2(cm) <<  std::endl;
+
 
     // LatticeComplex trlcMat(&Fine);
     // trlcMat = trace(lcMat); // Trace involving iVector - now generates error