diff --git a/Grid/Grid.h b/Grid/Grid.h
index 70859c8154..4bf573b77d 100644
--- a/Grid/Grid.h
+++ b/Grid/Grid.h
@@ -40,11 +40,11 @@ Author: paboyle <paboyle@ph.ed.ac.uk>
 
 #include <Grid/GridCore.h>
 #include <Grid/GridQCDcore.h>
-#include <Grid/qcd/action/Action.h>
+#include <Grid/qcd/action/Action.h>// now contains #include <Grid/parallelIO/MetaData.h>
 #include <Grid/qcd/utils/GaugeFix.h>
 #include <Grid/qcd/utils/CovariantSmearing.h>
 #include <Grid/qcd/smearing/Smearing.h>
-#include <Grid/parallelIO/MetaData.h>
+//#include <Grid/parallelIO/MetaData.h>//moved to ActionCore.h 
 #include <Grid/qcd/hmc/HMC_aggregate.h>
 
 #endif
diff --git a/Grid/GridCore.h b/Grid/GridCore.h
index 90d01f2e46..cdb1aeb755 100644
--- a/Grid/GridCore.h
+++ b/Grid/GridCore.h
@@ -42,12 +42,12 @@ Author: paboyle <paboyle@ph.ed.ac.uk>
 #include <Grid/DisableWarnings.h>
 #include <Grid/Namespace.h>
 #include <Grid/GridStd.h>
-#include <Grid/threads/Pragmas.h>
 #include <Grid/perfmon/Timer.h>
-//#include <Grid/perfmon/PerfCount.h>
-#include <Grid/util/Util.h>
 #include <Grid/log/Log.h>
 #include <Grid/perfmon/Tracing.h>
+#include <Grid/threads/Pragmas.h>
+//#include <Grid/perfmon/PerfCount.h>
+#include <Grid/util/Util.h>
 #include <Grid/allocator/Allocator.h>
 #include <Grid/simd/Simd.h>
 #include <Grid/threads/ThreadReduction.h>
diff --git a/Grid/Grid_Eigen_Dense.h b/Grid/Grid_Eigen_Dense.h
index 8bd1d11357..0d2249ec25 100644
--- a/Grid/Grid_Eigen_Dense.h
+++ b/Grid/Grid_Eigen_Dense.h
@@ -39,6 +39,7 @@
 #endif
 
 /* HIP save and restore compile environment*/
+
 #ifdef GRID_HIP
 #pragma push
 #pragma push_macro("__HIP_DEVICE_COMPILE__")
@@ -63,11 +64,12 @@
 #endif
 
 /*HIP restore*/
+/*
 #ifdef __HIP__REDEFINE__
 #pragma pop_macro("__HIP_DEVICE_COMPILE__")
 #pragma pop
 #endif
-
+*/
 #if defined __GNUC__
 #pragma GCC diagnostic pop
 #endif
diff --git a/Grid/cartesian/Cartesian_base.h b/Grid/cartesian/Cartesian_base.h
index 1a5049ca21..c32de62b24 100644
--- a/Grid/cartesian/Cartesian_base.h
+++ b/Grid/cartesian/Cartesian_base.h
@@ -92,6 +92,7 @@ class GridBase : public CartesianCommunicator , public GridThread {
   ////////////////////////////////////////////////////////////////
   virtual int CheckerBoarded(int dim) =0;
   virtual int CheckerBoard(const Coordinate &site)=0;
+  virtual int CheckerDim(void){ return 0; };
   virtual int CheckerBoardDestination(int source_cb,int shift,int dim)=0;
   virtual int CheckerBoardShift(int source_cb,int dim,int shift,int osite)=0;
   virtual int CheckerBoardShiftForCB(int source_cb,int dim,int shift,int cb)=0;
diff --git a/Grid/cartesian/Cartesian_red_black.h b/Grid/cartesian/Cartesian_red_black.h
index f218c2b56d..3099907a8c 100644
--- a/Grid/cartesian/Cartesian_red_black.h
+++ b/Grid/cartesian/Cartesian_red_black.h
@@ -60,6 +60,7 @@ class GridRedBlackCartesian : public GridBase
   //  int              _checker_dim;
   std::vector<int> _checker_board;
 
+  virtual int CheckerDim(void){ return _checker_dim; };
   virtual int isCheckerBoarded(void) const { return 1; };
   virtual int CheckerBoarded(int dim){
     if( dim==_checker_dim) return 1;
diff --git a/Grid/cshift/Cshift_common.h b/Grid/cshift/Cshift_common.h
index fa3f27a56f..19eb5ba7ed 100644
--- a/Grid/cshift/Cshift_common.h
+++ b/Grid/cshift/Cshift_common.h
@@ -52,6 +52,7 @@ inline std::pair<int,int> *MapCshiftTable(void)
 template<class vobj> void 
 Gather_plane_simple (const Lattice<vobj> &rhs,deviceVector<vobj> &buffer,int dimension,int plane,int cbmask, int off=0)
 {
+  GRID_TRACE("Gather_plane_simple");
   int rd = rhs.Grid()->_rdimensions[dimension];
 
   if ( !rhs.Grid()->CheckerBoarded(dimension) ) {
@@ -105,6 +106,7 @@ Gather_plane_extract(const Lattice<vobj> &rhs,
 		     ExtractPointerArray<typename vobj::scalar_object> pointers,
 		     int dimension,int plane,int cbmask)
 {
+  GRID_TRACE("Gather_plane_extract");
   int rd = rhs.Grid()->_rdimensions[dimension];
 
   if ( !rhs.Grid()->CheckerBoarded(dimension) ) {
@@ -160,6 +162,7 @@ Gather_plane_extract(const Lattice<vobj> &rhs,
 //////////////////////////////////////////////////////
 template<class vobj> void Scatter_plane_simple (Lattice<vobj> &rhs,deviceVector<vobj> &buffer, int dimension,int plane,int cbmask)
 {
+  GRID_TRACE("Scatter_plane_simple");
   int rd = rhs.Grid()->_rdimensions[dimension];
 
   if ( !rhs.Grid()->CheckerBoarded(dimension) ) {
@@ -214,6 +217,7 @@ template<class vobj> void Scatter_plane_simple (Lattice<vobj> &rhs,deviceVector<
 //////////////////////////////////////////////////////
 template<class vobj> void Scatter_plane_merge(Lattice<vobj> &rhs,ExtractPointerArray<typename vobj::scalar_object> pointers,int dimension,int plane,int cbmask)
 {
+  GRID_TRACE("Scatter_plane_merge");
   int rd = rhs.Grid()->_rdimensions[dimension];
 
   if ( !rhs.Grid()->CheckerBoarded(dimension) ) {
@@ -263,6 +267,7 @@ template<class vobj> void Scatter_plane_merge(Lattice<vobj> &rhs,ExtractPointerA
 
 template<class vobj> void Copy_plane(Lattice<vobj>& lhs,const Lattice<vobj> &rhs, int dimension,int lplane,int rplane,int cbmask)
 {
+  GRID_TRACE("Copy_plane");
   int rd = rhs.Grid()->_rdimensions[dimension];
 
   if ( !rhs.Grid()->CheckerBoarded(dimension) ) {
@@ -311,6 +316,7 @@ template<class vobj> void Copy_plane(Lattice<vobj>& lhs,const Lattice<vobj> &rhs
 
 template<class vobj> void Copy_plane_permute(Lattice<vobj>& lhs,const Lattice<vobj> &rhs, int dimension,int lplane,int rplane,int cbmask,int permute_type)
 {
+  GRID_TRACE("Copy_plane_permute");
   int rd = rhs.Grid()->_rdimensions[dimension];
 
   if ( !rhs.Grid()->CheckerBoarded(dimension) ) {
@@ -358,6 +364,7 @@ template<class vobj> void Copy_plane_permute(Lattice<vobj>& lhs,const Lattice<vo
 //////////////////////////////////////////////////////
 template<class vobj> void Cshift_local(Lattice<vobj>& ret,const Lattice<vobj> &rhs,int dimension,int shift)
 {
+  GRID_TRACE("Cshift_local");
   int sshift[2];
 
   sshift[0] = rhs.Grid()->CheckerBoardShiftForCB(rhs.Checkerboard(),dimension,shift,Even);
diff --git a/Grid/cshift/Cshift_mpi.h b/Grid/cshift/Cshift_mpi.h
index 252478ade2..4688ff25ad 100644
--- a/Grid/cshift/Cshift_mpi.h
+++ b/Grid/cshift/Cshift_mpi.h
@@ -38,6 +38,7 @@ extern uint64_t checksum_index;
 const int Cshift_verbose=0;
 template<class vobj> Lattice<vobj> Cshift(const Lattice<vobj> &rhs,int dimension,int shift)
 {
+  GRID_TRACE("Cshift");
   typedef typename vobj::vector_type vector_type;
   typedef typename vobj::scalar_type scalar_type;
 
@@ -138,6 +139,7 @@ template<class vobj> void Cshift_simple(Lattice<vobj>& ret,const Lattice<vobj> &
 }
 template<class vobj> void Cshift_comms(Lattice<vobj>& ret,const Lattice<vobj> &rhs,int dimension,int shift)
 {
+  GRID_TRACE("Cshift_comms");
   int sshift[2];
 
   sshift[0] = rhs.Grid()->CheckerBoardShiftForCB(rhs.Checkerboard(),dimension,shift,Even);
@@ -156,6 +158,7 @@ template<class vobj> void Cshift_comms(Lattice<vobj>& ret,const Lattice<vobj> &r
 
 template<class vobj> void Cshift_comms_simd(Lattice<vobj>& ret,const Lattice<vobj> &rhs,int dimension,int shift)
 {
+  GRID_TRACE("Cshift_comms_simd");
   int sshift[2];
 
   sshift[0] = rhs.Grid()->CheckerBoardShiftForCB(rhs.Checkerboard(),dimension,shift,Even);
@@ -173,6 +176,7 @@ template<class vobj> void Cshift_comms_simd(Lattice<vobj>& ret,const Lattice<vob
 }
 template<class vobj> void Cshift_comms(Lattice<vobj> &ret,const Lattice<vobj> &rhs,int dimension,int shift,int cbmask)
 {
+  GRID_TRACE("Cshift_comms_cb");
   typedef typename vobj::vector_type vector_type;
   typedef typename vobj::scalar_type scalar_type;
 
@@ -305,6 +309,7 @@ template<class vobj> void Cshift_comms(Lattice<vobj> &ret,const Lattice<vobj> &r
 
 template<class vobj> void  Cshift_comms_simd(Lattice<vobj> &ret,const Lattice<vobj> &rhs,int dimension,int shift,int cbmask)
 {
+  GRID_TRACE("Cshift_comms_simd_cb");
   GridBase *grid=rhs.Grid();
   const int Nsimd = grid->Nsimd();
   typedef typename vobj::vector_type vector_type;
diff --git a/Grid/lattice/Lattice_trace.h b/Grid/lattice/Lattice_trace.h
index 9a0c6f9647..507a6804d3 100644
--- a/Grid/lattice/Lattice_trace.h
+++ b/Grid/lattice/Lattice_trace.h
@@ -97,6 +97,7 @@ Lattice<iScalar<iScalar<iScalar<Vec> > > > Determinant(const Lattice<iScalar<iSc
 template<int N>
 Lattice<iScalar<iScalar<iMatrix<vComplexD, N> > > > Inverse(const Lattice<iScalar<iScalar<iMatrix<vComplexD, N> > > > &Umu)
 {
+#if 0
   GridBase *grid=Umu.Grid();
   auto lvol = grid->lSites();
   Lattice<iScalar<iScalar<iMatrix<vComplexD, N> > > > ret(grid);
@@ -121,9 +122,163 @@ Lattice<iScalar<iScalar<iMatrix<vComplexD, N> > > > Inverse(const Lattice<iScala
       }}
     pokeLocalSite(Ui,ret_v,lcoor);
   });
+#else
+  GridBase *grid=Umu.Grid();
+  auto osites = grid->oSites();
+  const int Nsimd=grid->Nsimd();
+  Lattice<iScalar<iScalar<iMatrix<vComplexD, N> > > > ret(grid);
+  autoView(Umu_v,Umu,CpuRead);
+  autoView(ret_v,ret,CpuWrite);
+  thread_for(site,osites,{
+    Eigen::MatrixXcd EigenU = Eigen::MatrixXcd::Zero(N,N);
+
+    iScalar<iScalar<iMatrix<ComplexD, N> > > Us;
+    iScalar<iScalar<iMatrix<ComplexD, N> > > Ui;
+
+    for(int lane=0;lane<Nsimd;lane++){
+      Us = extractLane(lane,Umu_v[site]);
+      for(int i=0;i<N;i++){
+	for(int j=0;j<N;j++){
+	  EigenU(i,j) = Us()()(i,j);
+	}}
+      Eigen::MatrixXcd EigenUinv = EigenU.inverse();
+      for(int i=0;i<N;i++){
+	for(int j=0;j<N;j++){
+	  Ui()()(i,j) = EigenUinv(i,j);
+	}}
+      insertLane(lane,ret_v[site],Ui);
+    }
+  });
+#endif  
   return ret;
 }
 
+/* Helper functions for inversion of real matrix on GPU based on Nobu's code*/
+template<class type1,class type2,int N>
+accelerator_inline void LUdcmp( iMatrix<type1,N> &LU, iVector<type2,N> &P)
+{
+
+  const RealD TINY=1.0e-40;
+  
+  type1 vv[N], tmp, _max;
+  for(int i=0; i<N; i++) vv[i]=0.0;
+  
+  for(int i=0; i<N; i++){
+    _max=0.0;
+    for(int j=0; j<N; j++) if( (tmp=abs(LU(i,j))) > _max ) _max = tmp;
+    assert( abs(_max) > TINY );
+    vv[i] = abs(1.0/_max);
+  }
+
+  int imax;
+  for(int k=0; k<N; k++){
+    _max=0.0;
+    for(int i=k; i<N; i++){
+      tmp = vv[i] * abs( LU(i,k) );
+      if(tmp>_max) {
+	_max = tmp;
+	imax = i;
+      }
+    }
+    if(k!=imax){
+      for(int j=0; j<N; j++){
+	tmp = LU(imax,j);
+        LU(imax,j) = LU(k,j);
+	LU(k,j) = tmp;
+      }
+      vv[imax] = vv[k];
+    }
+    P(k)=imax;
+    
+    for(int i=k+1; i<N; i++){
+      LU(i,k) = LU(i,k) / LU(k,k);
+      tmp = LU(i,k);
+      for(int j=k+1; j<N; j++) LU(i,j) = LU(i,j) - tmp * LU(k,j);
+    } // end i
+  } // end k
+};
+  
+template<class type1,class type2,int N>
+accelerator_inline void solve( iVector<type1,N> &x, const iMatrix<type1,N> LU, const iVector<type2,N> P){
+  
+  type1 sum = 0.0;
+
+  int ii=0;
+  for(int i=0; i<N; i++){
+    int ip = P(i);
+    sum = x(ip);
+    x(ip) = x(i);
+    if(ii!=0)for(int j=ii-1;j<i;j++) sum = sum - LU(i,j)*x(j);
+    else if (abs(sum)>0.0) ii=i+1;
+    x(i) = sum;
+  }
+  for(int i=N-1; i>=0; i--){
+    sum = x(i);
+    for(int j=i+1; j<N; j++) sum = sum - LU(i,j)*x(j);
+    x(i) = sum/LU(i,i);
+  }
+};
+
+template<int N>
+Lattice<iScalar<iScalar<iMatrix<vComplexD, N> > > > Inverse_RealPart(const Lattice<iScalar<iScalar<iMatrix<vComplexD, N> > > > &Umu)
+{
+  GridBase *grid=Umu.Grid();
+  auto osites = grid->oSites();
+  const int Nsimd=grid->Nsimd();
+  Lattice<iScalar<iScalar<iMatrix<vComplexD, N> > > > ret(grid);
+#if 0 // CPU version
+  autoView(Umu_v,Umu,CpuRead);
+  autoView(ret_v,ret,CpuWrite);
+  thread_for(site,osites,{
+    Eigen::MatrixXd EigenU = Eigen::MatrixXd::Zero(N,N);
+
+    iScalar<iScalar<iMatrix<ComplexD, N> > > Us;
+    iScalar<iScalar<iMatrix<ComplexD, N> > > Ui;
+
+    for(int lane=0;lane<Nsimd;lane++){
+      Us = extractLane(lane,Umu_v[site]);
+      for(int i=0;i<N;i++){
+        for(int j=0;j<N;j++){
+          EigenU(i,j) = real(Us()()(i,j));
+        }}
+      Eigen::MatrixXd EigenUinv = EigenU.inverse();
+      for(int i=0;i<N;i++){
+        for(int j=0;j<N;j++){
+          Ui()()(i,j) = EigenUinv(i,j);
+        }}
+      insertLane(lane,ret_v[site],Ui);
+    }
+  });
+#else //GPU version
+  autoView(Umu_v,Umu,AcceleratorRead);
+  autoView(ret_v,ret,AcceleratorWrite);
+  accelerator_for(ss,grid->oSites(),vComplex::Nsimd(),{
+      iMatrix<RealD, N>  LU;
+      iVector<Integer, N> P;
+      iVector<RealD, N> e;
+      // scalar layout won't coalesce
+#ifdef GRID_SIMT
+      {
+	int blane=acceleratorSIMTlane(Nsimd); // buffer lane
+#else
+      for(int blane=0;blane<Nsimd;blane++) {
+#endif
+	
+	for(int i=0;i<N;i++){
+	  for(int j=0;j<N;j++){
+	    LU(i,j) = getlane(toReal(TensorRemove(Umu_v(ss)()()(i,j))),blane);
+	  }}
+	LUdcmp(LU,P);
+	for(int j=0; j<N; j++){
+	  for(int i=0; i<N; i++) e(i) = (i==j);
+	  solve(e,LU,P);
+	  for(int i=0; i<N; i++) putlane(ret_v[ss]()()(i,j),(ComplexD) e(i),blane);
+	}
+      }
+    });
+#endif
+ return ret;
+}
 
 NAMESPACE_END(Grid);
 #endif
diff --git a/Grid/lattice/Lattice_transfer.h b/Grid/lattice/Lattice_transfer.h
index 20ca49eb77..fc52ac0e3a 100644
--- a/Grid/lattice/Lattice_transfer.h
+++ b/Grid/lattice/Lattice_transfer.h
@@ -42,50 +42,21 @@ inline void subdivides(GridBase *coarse,GridBase *fine)
     GRID_ASSERT((fine->_rdimensions[d] / coarse->_rdimensions[d])* coarse->_rdimensions[d]==fine->_rdimensions[d]); 
   }
 }
-
  
 ////////////////////////////////////////////////////////////////////////////////////////////
 // remove and insert a half checkerboard
 ////////////////////////////////////////////////////////////////////////////////////////////
+
 template<class vobj> inline void pickCheckerboard(int cb,Lattice<vobj> &half,const Lattice<vobj> &full)
 {
-  half.Checkerboard() = cb;
-
-  autoView( half_v, half, CpuWrite);
-  autoView( full_v, full, CpuRead);
-  thread_for(ss, full.Grid()->oSites(),{
-    int cbos;
-    Coordinate coor;
-    full.Grid()->oCoorFromOindex(coor,ss);
-    cbos=half.Grid()->CheckerBoard(coor);
-
-    if (cbos==cb) {
-      int ssh=half.Grid()->oIndex(coor);
-      half_v[ssh] = full_v[ss];
-    }
-  });
+  acceleratorPickCheckerboard(cb,half,full);
 }
 template<class vobj> inline void setCheckerboard(Lattice<vobj> &full,const Lattice<vobj> &half)
 {
-  int cb = half.Checkerboard();
-  autoView( half_v , half, CpuRead);
-  autoView( full_v , full, CpuWrite);
-  thread_for(ss,full.Grid()->oSites(),{
-
-    Coordinate coor;
-    int cbos;
-
-    full.Grid()->oCoorFromOindex(coor,ss);
-    cbos=half.Grid()->CheckerBoard(coor);
-      
-    if (cbos==cb) {
-      int ssh=half.Grid()->oIndex(coor);
-      full_v[ss]=half_v[ssh];
-    }
-  });
+  acceleratorSetCheckerboard(full,half);
 }
 
-template<class vobj> inline void acceleratorPickCheckerboard(int cb,Lattice<vobj> &half,const Lattice<vobj> &full, int checker_dim_half=0)
+template<class vobj> inline void acceleratorPickCheckerboard(int cb,Lattice<vobj> &half,const Lattice<vobj> &full, int dummy=0)
 {
   half.Checkerboard() = cb;
   autoView(half_v, half, AcceleratorWrite);
@@ -95,6 +66,7 @@ template<class vobj> inline void acceleratorPickCheckerboard(int cb,Lattice<vobj
   unsigned long ndim_half          = half.Grid()->_ndimension;
   Coordinate checker_dim_mask_half = half.Grid()->_checker_dim_mask;
   Coordinate ostride_half          = half.Grid()->_ostride;
+  int checker_dim_half             = half.Grid()->CheckerDim();
   accelerator_for(ss, full.Grid()->oSites(),full.Grid()->Nsimd(),{
     
     Coordinate coor;
@@ -119,7 +91,7 @@ template<class vobj> inline void acceleratorPickCheckerboard(int cb,Lattice<vobj
     }
   });
 }
-template<class vobj> inline void acceleratorSetCheckerboard(Lattice<vobj> &full,const Lattice<vobj> &half, int checker_dim_half=0)
+template<class vobj> inline void acceleratorSetCheckerboard(Lattice<vobj> &full,const Lattice<vobj> &half, int dummy=0)
 {
   int cb = half.Checkerboard();
   autoView(half_v , half, AcceleratorRead);
@@ -129,6 +101,7 @@ template<class vobj> inline void acceleratorSetCheckerboard(Lattice<vobj> &full,
   unsigned long ndim_half          = half.Grid()->_ndimension;
   Coordinate checker_dim_mask_half = half.Grid()->_checker_dim_mask;
   Coordinate ostride_half          = half.Grid()->_ostride;
+  int checker_dim_half             = half.Grid()->CheckerDim();
   accelerator_for(ss,full.Grid()->oSites(),full.Grid()->Nsimd(),{
 
     Coordinate coor;
diff --git a/Grid/lattice/Lattice_view.h b/Grid/lattice/Lattice_view.h
index 4affe542b6..98df3a104c 100644
--- a/Grid/lattice/Lattice_view.h
+++ b/Grid/lattice/Lattice_view.h
@@ -144,7 +144,7 @@ class ViewCloser
 
 #define autoView(l_v,l,mode)				\
 	  auto l_v = l.View(mode);			\
-	  ViewCloser<decltype(l_v)> _autoView##l_v(l_v,__FILE__,__LINE__,mode);
+	  ViewCloser<decltype(l_v)> _autoView##l_v(l_v,(const char *)__FILE__,(int)__LINE__,(int)mode);
 
 #else
 // Little autoscope assister
diff --git a/Grid/lattice/PaddedCell.h b/Grid/lattice/PaddedCell.h
index 23f6c13af6..43e75c6ab6 100644
--- a/Grid/lattice/PaddedCell.h
+++ b/Grid/lattice/PaddedCell.h
@@ -298,6 +298,7 @@ class PaddedCell {
       if( processors[d]==1 ) fll[d]=0;
     }
     localCopyRegion(in,out,fll,tll,local);
+    //DumpSliceNorm("PaddedCell Extract ",out);
     return out;
   }
   template<class vobj>
@@ -309,6 +310,7 @@ class PaddedCell {
     for(int d=0;d<dims;d++){
       tmp = Expand(d,tmp,cshift); // rvalue && assignment
     }
+    //    DumpSliceNorm("PaddedCell Exchange",tmp);
     return tmp;
   }
   template<class vobj>
@@ -320,6 +322,7 @@ class PaddedCell {
     for(int d=0;d<dims;d++){
       tmp = ExpandPeriodic(d,tmp); // rvalue && assignment
     }
+    //    DumpSliceNorm("PaddedCell ExchangePeriodic",tmp);
     return tmp;
   }
   // expand up one dim at a time
@@ -391,6 +394,7 @@ class PaddedCell {
 
     }
     std::cout << GridLogPerformance << "PaddedCell::Expand timings: cshift:" << tshift/1000 << "ms, insert-slice:" << tins/1000 << "ms" << std::endl;
+    //    DumpSliceNorm("PaddedCell Expand",padded);
     
     return padded;
   }
@@ -418,8 +422,10 @@ class PaddedCell {
       padded=in; // slightly different interface could avoid a copy operation
     } else {
       Face_exchange(in,padded,dim,depth);
+      //      DumpSliceNorm("PaddedCell ExpandPeriodic",padded);
       return padded;
     }
+    //    DumpSliceNorm("PaddedCell ExpandPeriodic",padded);
     return padded;
   }
   template<class vobj>
@@ -436,8 +442,8 @@ class PaddedCell {
     RealD t_comms=0.0;
     RealD t_copy=0.0;
     
-    //    std::cout << GridLogMessage << "dimension " <<dimension<<std::endl;
-    //    DumpSliceNorm(std::string("Face_exchange from"),from,dimension);
+    //    std::cout << GridLogMessage << "PDCF dimension " <<dimension<<" "<<depth<<std::endl;
+    //    DumpSliceNorm(std::string("PaddedCell::Face_exchange from "),from,dimension);
     GridBase *grid=from.Grid();
     GridBase *new_grid=to.Grid();
 
@@ -582,6 +588,8 @@ class PaddedCell {
     t_scatter+= usecond() - t;
     t_tot+=usecond();
 
+    //    DumpSliceNorm("PaddedCell::FaceExchange to ",to,dimension);
+
     std::cout << GridLogPerformance << "PaddedCell::Expand new timings: gather :" << t_gather/1000  << "ms"<<std::endl;
     std::cout << GridLogPerformance << "PaddedCell::Expand new timings: scatter:" << t_scatter/1000   << "ms"<<std::endl;
     std::cout << GridLogPerformance << "PaddedCell::Expand new timings: copy   :" << t_copy/1000      << "ms"<<std::endl;
diff --git a/Grid/parallelIO/BinaryIO.h b/Grid/parallelIO/BinaryIO.h
index 4df9fdf9e9..96190b6000 100644
--- a/Grid/parallelIO/BinaryIO.h
+++ b/Grid/parallelIO/BinaryIO.h
@@ -560,8 +560,11 @@ class BinaryIO {
     std::vector<sobj> scalardata(lsites); 
     std::vector<fobj>     iodata(lsites); // Munge, checksum, byte order in here
     
+    FlightRecorder::StepLog("readLatticeObject BEGIN");
+    grid->Barrier();
     IOobject(w,grid,iodata,file,offset,format,BINARYIO_READ|control,
 	     nersc_csum,scidac_csuma,scidac_csumb);
+    FlightRecorder::StepLog("readLatticeObject END");
 
     GridStopWatch timer; 
     timer.Start();
@@ -611,9 +614,12 @@ class BinaryIO {
     timer.Stop();
     while (attemptsLeft >= 0)
     {
+      FlightRecorder::StepLog("writeLatticeObject BEGIN");
       grid->Barrier();
       IOobject(w,grid,iodata,file,offset,format,BINARYIO_WRITE|control,
 	             nersc_csum,scidac_csuma,scidac_csumb);
+      FlightRecorder::StepLog("writeLatticeObject END");
+      grid->Barrier();
       if (checkWrite)
       {
         std::vector<fobj> ckiodata(lsites);
diff --git a/Grid/parallelIO/IldgIO.h b/Grid/parallelIO/IldgIO.h
index 8ed7df5b57..7b916ac16b 100644
--- a/Grid/parallelIO/IldgIO.h
+++ b/Grid/parallelIO/IldgIO.h
@@ -24,7 +24,8 @@ See the full license in the file "LICENSE" in the top level distribution
 directory
 *************************************************************************************/
 /*  END LEGAL */
-#pragma once
+#ifndef GRID_ILDG_IO_H
+#define GRID_ILDG_IO_H
 
 #ifdef HAVE_LIME
 #include <algorithm>
@@ -951,4 +952,4 @@ NAMESPACE_END(Grid);
 
 //HAVE_LIME
 #endif
-
+#endif
diff --git a/Grid/parallelIO/MetaData.h b/Grid/parallelIO/MetaData.h
index 1fc3fe31d5..1e38cd4879 100644
--- a/Grid/parallelIO/MetaData.h
+++ b/Grid/parallelIO/MetaData.h
@@ -27,6 +27,9 @@
 *************************************************************************************/
 /*  END LEGAL */
 
+#ifndef GRID_METADATA_IO_H
+#define GRID_METADATA_IO_H
+
 #include <algorithm>
 #include <iostream>
 #include <iomanip>
@@ -36,6 +39,7 @@
 #include <sys/utsname.h>
 #include <pwd.h>
 
+
 NAMESPACE_BEGIN(Grid);
 
 ///////////////////////////////////////////////////////
@@ -343,3 +347,4 @@ struct Gauge3x2unmunger{
 
 NAMESPACE_END(Grid);
 
+#endif
diff --git a/Grid/qcd/action/Action.h b/Grid/qcd/action/Action.h
index 737c1ff02b..8b8da5b21d 100644
--- a/Grid/qcd/action/Action.h
+++ b/Grid/qcd/action/Action.h
@@ -33,6 +33,7 @@ Author: paboyle <paboyle@ph.ed.ac.uk>
 #ifndef GRID_QCD_ACTION_H
 #define GRID_QCD_ACTION_H
 
+
 ////////////////////////////////////////////
 // Abstract base interface
 ////////////////////////////////////////////
diff --git a/Grid/qcd/action/ActionBase.h b/Grid/qcd/action/ActionBase.h
index 96afea0a5a..a765a9d1a9 100644
--- a/Grid/qcd/action/ActionBase.h
+++ b/Grid/qcd/action/ActionBase.h
@@ -32,8 +32,67 @@ directory
 #ifndef ACTION_BASE_H
 #define ACTION_BASE_H
 
+//#define TIMom
+#define PRINT_SNAPSHOTS
+#ifdef PRINT_SNAPSHOTS
+#include <algorithm>
+#include <cctype>
+#endif
+
+
 NAMESPACE_BEGIN(Grid);
 
+#ifdef PRINT_SNAPSHOTS
+template <class vobj> LatticeComplex LocalFieldSquareNorm(Lattice<vobj>& F){
+  LatticeComplex Hloc(F.Grid()), tmp(F.Grid());
+  LatticeColourMatrix Fmu(F.Grid());
+  Hloc = Zero();
+  for (int mu = 0; mu < Nd; mu++) {
+    Fmu = PeekIndex<LorentzIndex>(F, mu);
+    tmp = (-1.0/HMC_MOMENTUM_DENOMINATOR)*trace(Fmu * Fmu);
+    Hloc = Hloc + tmp;
+  }
+  return Hloc;
+}
+template <class vobj> void writeField(Lattice<vobj> &F, std::string filename, bool reduce=true, bool alg=true) {
+  
+  std::replace(filename.begin(),filename.end(), '/', '_');
+  std::replace(filename.begin(),filename.end(), '(', '_');
+  std::replace(filename.begin(),filename.end(), ')', ' '); // wanted to eliminate ')' from file name
+  filename.erase(std::remove_if(filename.begin(), filename.end(), [](unsigned char x) { return std::isspace(x); }), filename.end());
+
+#ifdef HAVE_LIME
+  emptyUserRecord record;
+  ScidacWriter WR(F.Grid()->IsBoss());
+  WR.open(filename);
+  if (reduce) {
+    auto F_reduced=alg? LocalFieldSquareNorm(F) : localNorm2(F);
+    WR.writeScidacFieldRecord(F_reduced,record,0);
+  }
+  else
+    WR.writeScidacFieldRecord(F,record,0);
+  WR.close();
+#else
+  assert( false && "writeField function needs c-lime support!");
+#endif
+}
+
+// vobj = vLorentzColourMatrixD
+template <class vobj> void writeConfig(Lattice<vobj> &F, std::string filename){
+
+  std::replace(filename.begin(),filename.end(), '/', '_');
+  std::replace(filename.begin(),filename.end(), '(', '_');
+  std::replace(filename.begin(),filename.end(), ')', ' '); // wanted to eliminate ')' from file name
+  filename.erase(std::remove_if(filename.begin(), filename.end(), [](unsigned char x) { return std::isspace(x); }), filename.end());
+  
+  // Use NERSC IO for consistency with our HMC runs
+  // Then, GaugeStats = PeriodicGaugeStatistics <- the default
+  int precision32 = 1;
+  int tworow = 0;
+  NerscIO::writeConfiguration(F, filename, tworow, precision32); 
+}
+#endif
+
 ///////////////////////////////////
 // Smart configuration base class
 ///////////////////////////////////
@@ -53,7 +112,7 @@ template <class GaugeField >
 class Action 
 {
 public:
-  bool is_smeared = false;
+  bool  is_smeared = false;
   RealD deriv_norm_sum;
   RealD deriv_max_sum;
   RealD Fdt_norm_sum;
@@ -116,9 +175,17 @@ class Action
   }
   virtual void deriv(ConfigurationBase<GaugeField>& U, GaugeField& dSdU)
   {
-    deriv(U.get_U(is_smeared),dSdU); 
+    // PRINT_SNAPSHOTS: For some actions, deriv is overridden => could be better if we put writeField in update_P
+    deriv(U.get_U(is_smeared),dSdU);
+#ifdef PRINT_SNAPSHOTS
+    writeField(dSdU, "F_"+action_name()+"_lat."+std::to_string(deriv_num));
+#endif
+    
     if ( is_smeared ) {
       U.smeared_force(dSdU);
+#ifdef PRINT_SNAPSHOTS
+      writeField(dSdU, "F_"+action_name()+"_smr."+std::to_string(deriv_num));
+#endif
     }
   }
   ///////////////////////////////
diff --git a/Grid/qcd/action/ActionCore.h b/Grid/qcd/action/ActionCore.h
index db916b0217..6c5f214741 100644
--- a/Grid/qcd/action/ActionCore.h
+++ b/Grid/qcd/action/ActionCore.h
@@ -30,7 +30,24 @@ directory
 #ifndef QCD_ACTION_CORE
 #define QCD_ACTION_CORE
 
-#include <Grid/qcd/action/gauge/GaugeImplementations.h>
+
+#include <Grid/qcd/action/gauge/GaugeImplementations.h> 
+
+////////////////////////////////////////////
+// I/O Interfaces for Action and HMC
+////////////////////////////////////////////
+#include <Grid/qcd/utils/WilsonLoops.h>    // also included from Grid/qcd/action/gauge/Gauge.h
+#include <Grid/parallelIO/MetaData.h>      // moved from Grid.h
+// moved from HMC_aggregate.h to allow IO op's in action
+// Action.h preceeds HMC_aggregate.h in Grid.h => safe to move here from HMC_aggregate.h
+#include <Grid/parallelIO/IldgIOtypes.h>
+#include <Grid/parallelIO/IldgIO.h>
+#include <Grid/parallelIO/NerscIO.h>
+#include <Grid/parallelIO/OpenQcdIO.h>
+#if !defined(GRID_COMMS_NONE)
+#include <Grid/parallelIO/OpenQcdIOChromaReference.h>
+#endif
+NAMESPACE_CHECK(Ildg);
 
 #include <Grid/qcd/action/ActionBase.h>
 NAMESPACE_CHECK(ActionBase);
diff --git a/Grid/qcd/action/gauge/GaugeImplTypes.h b/Grid/qcd/action/gauge/GaugeImplTypes.h
index 95b692905f..3137952d0f 100644
--- a/Grid/qcd/action/gauge/GaugeImplTypes.h
+++ b/Grid/qcd/action/gauge/GaugeImplTypes.h
@@ -125,7 +125,22 @@ template <class S, int Nrepresentation = Nc, int Nexp = 12, class Group = SU<Nc>
       PokeIndex<LorentzIndex>(P, Pmu, mu);
     }
   }
-    
+
+  //#ifdef TIMom
+  static inline void generate_momenta(Field &P, GridSerialRNG & sRNG, GridParallelRNG &pRNG, LatticeComplex W)
+  {
+    LinkField Pmu(P.Grid());
+    Pmu = Zero();
+
+    for (int mu = 0; mu < Nd; mu++) {
+      Group::GaussianFundamentalLieAlgebraMatrix(pRNG, Pmu);
+      RealD scale = ::sqrt(HMC_MOMENTUM_DENOMINATOR) ;
+      Pmu = scale*W*Pmu;
+      PokeIndex<LorentzIndex>(P, Pmu, mu);
+    }
+  }
+  //#endif
+  
   static inline Field projectForce(Field &P) {
       Field ret(P.Grid());
       Group::taProj(P, ret);
@@ -151,7 +166,7 @@ template <class S, int Nrepresentation = Nc, int Nexp = 12, class Group = SU<Nc>
    // diff += end - start;
    // std::cout << "Time to exponentiate matrix " << diff.count() << " s\n";
   }
-    
+
   static inline RealD FieldSquareNorm(Field& U){
     LatticeComplex Hloc(U.Grid());
     Hloc = Zero();
diff --git a/Grid/qcd/hmc/HMC.h b/Grid/qcd/hmc/HMC.h
index c56ae6ad8c..1087dfd2fe 100644
--- a/Grid/qcd/hmc/HMC.h
+++ b/Grid/qcd/hmc/HMC.h
@@ -263,10 +263,14 @@ class HybridMonteCarlo {
       if (traj < Params.StartTrajectory + Params.NoMetropolisUntil) {
       	std::cout << GridLogHMC << "-- Thermalization" << std::endl;
       }
+
+#ifndef PRINT_SNAPSHOTS
+      Grid_heartbeat();
+#endif
       
       double t0=usecond();
       Ucopy = Ucur;
-
+      
       DeltaH = evolve_hmc_step(Ucopy);
       // Metropolis-Hastings test
       bool accept = true;
@@ -282,6 +286,10 @@ class HybridMonteCarlo {
       double t1=usecond();
       std::cout << GridLogHMC << "Total time for trajectory (s): " << (t1-t0)/1e6 << std::endl;
 
+#ifndef PRINT_SNAPSHOTS
+      Grid_heartbeat_off();
+#endif
+
       TheIntegrator.print_timer();
       
       TheIntegrator.Smearer.set_Field(Ucur);
diff --git a/Grid/qcd/hmc/HMC_aggregate.h b/Grid/qcd/hmc/HMC_aggregate.h
index cb5109532b..f8edf9dc0c 100644
--- a/Grid/qcd/hmc/HMC_aggregate.h
+++ b/Grid/qcd/hmc/HMC_aggregate.h
@@ -35,16 +35,6 @@ directory
 #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/parallelIO/OpenQcdIO.h>
-#if !defined(GRID_COMMS_NONE)
-#include <Grid/parallelIO/OpenQcdIOChromaReference.h>
-#endif
-NAMESPACE_CHECK(Ildg);
-
 #include <Grid/qcd/hmc/checkpointers/CheckPointers.h>
 #include <Grid/qcd/hmc/HMCModules.h>
 #include <Grid/qcd/modules/mods.h>
diff --git a/Grid/qcd/hmc/integrators/Integrator.h b/Grid/qcd/hmc/integrators/Integrator.h
index f497f7498f..b0b8b3bf40 100644
--- a/Grid/qcd/hmc/integrators/Integrator.h
+++ b/Grid/qcd/hmc/integrators/Integrator.h
@@ -27,8 +27,9 @@ with this program; if not, write to the Free Software Foundation, Inc.,
 See the full license in the file "LICENSE" in the top level distribution
 directory
 *************************************************************************************/
-			   /*  END LEGAL */
-			   //--------------------------------------------------------------------
+                             /*  END LEGAL */
+//--------------------------------------------------------------------
+
 #ifndef INTEGRATOR_INCLUDED
 #define INTEGRATOR_INCLUDED
 
@@ -41,6 +42,8 @@ class IntegratorParameters: Serializable {
   GRID_SERIALIZABLE_CLASS_MEMBERS(IntegratorParameters,
 				  std::string, name,      // name of the integrator
 				  unsigned int, MDsteps,  // number of outer steps
+				  RealD, alpha,
+				  int, box_dim,
 				  RealD, trajL)           // trajectory length
 
   IntegratorParameters(int MDsteps_ = 10, RealD trajL_ = 1.0)
@@ -62,6 +65,7 @@ class IntegratorParameters: Serializable {
   }
 };
 
+
 /*! @brief Class for Molecular Dynamics management */
 template <class FieldImplementation_, class SmearingPolicy, class RepresentationPolicy>
 class Integrator {
@@ -94,11 +98,30 @@ class Integrator {
     static MomentumFilterNone<MomentaField> filter;
     return &filter;
   }
+#ifdef TIMom
+
+  // TODO: make mask as a member; put most of the code to the constructor
+  static void scaleMom(LatticeComplex& mask, RealD alpha, int box_dim ){
+    LatticeInteger coor(mask.Grid()),tmp(mask.Grid());
+    LatticeComplex ones(mask.Grid());
+    mask = Zero(); ones = ComplexD(1.0,0.0); tmp = Zero();
+    for(int i=0; i<Nd; i++) {
+      LatticeCoordinate(coor,i);
+      tmp = tmp + mod(coor,box_dim);
+    }
+    mask =  where( tmp==Integer(0), alpha*ones, ones); // or alpha*mask + (ones - mask) if mask is precomputed
+  }
+#endif
 
   void update_P(Field& U, int level, double ep) 
   {
     t_P[level] += ep;
     update_P(P, U, level, ep);
+#ifdef PRINT_SNAPSHOTS
+    std::cout << GridLogMessage << " Half Mom Squared: " << level << " "<< t_P[level] << " "
+	      << - FieldImplementation::FieldSquareNorm(P)/HMC_MOMENTUM_DENOMINATOR << std::endl;
+    //writeField(P,"Mom_Field_"+std::to_string(t_P[level]));
+#endif
     std::cout << GridLogIntegrator << "[" << level << "] P " << " dt " << ep << " : t_P " << t_P[level] << std::endl;
   }
 
@@ -206,15 +229,29 @@ class Integrator {
     MomentaField MomFiltered(Mom.Grid());
     MomFiltered = Mom;
     MomFilter->applyFilter(MomFiltered);
-
+#ifdef TIMom
+    LatticeComplex W(U.Grid()); W = Zero();
+    scaleMom(W,Params.alpha*Params.alpha, Params.box_dim);
+    for (int mu = 0; mu < Nd; mu++) {
+      typename FieldImplementation::LinkField Pmu(W.Grid());
+      Pmu = PeekIndex<LorentzIndex>(MomFiltered, mu);
+      Pmu = W*Pmu;
+      PokeIndex<LorentzIndex>(MomFiltered, Pmu, mu);
+    }
+#endif
     // exponential of Mom*U in the gauge fields case
     FieldImplementation::update_field(MomFiltered, U, ep);
-
+#ifdef PRINT_SNAPSHOTS
+    writeConfig(U,"U_lat_"+std::to_string(t_U+ep));
+#endif
     // 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
+#ifdef PRINT_SNAPSHOTS
+    writeConfig(Smearer.get_SmearedU(),"U_smr_"+std::to_string(t_U+ep));
+#endif
   }
 
   virtual void step(Field& U, int level, int first, int last) = 0;
@@ -402,8 +439,19 @@ class Integrator {
     std::cout << GridLogIntegrator << "Integrator refresh" << std::endl;
 
     std::cout << GridLogIntegrator << "Generating momentum" << std::endl;
+#ifndef TIMom
     FieldImplementation::generate_momenta(P, sRNG, pRNG);
-
+#else
+    std::cout << GridLogIntegrator << "Scale momentum" << std::endl;
+
+    LatticeComplex W(U.Grid()); W = Zero();
+    ////////////////////
+    // Setup the mask
+    ////////////////////
+    scaleMom(W,1/Params.alpha, Params.box_dim);
+    FieldImplementation::generate_momenta(P, sRNG, pRNG, W);
+#endif
+    
     // Update the smeared fields, can be implemented as observer
     // necessary to keep the fields updated even after a reject
     // of the Metropolis
@@ -459,7 +507,10 @@ class Integrator {
     std::cout << GridLogIntegrator << "Integrator action\n";
 
     RealD H = - FieldImplementation::FieldSquareNorm(P)/HMC_MOMENTUM_DENOMINATOR; // - trace (P*P)/denom
-
+#ifdef PRINT_SNAPSHOTS
+    std::cout << GridLogMessage << "Final Half Mom Squared: " << H << std::endl;
+#endif
+    
     RealD Hterm;
 
     // Actions
@@ -505,7 +556,10 @@ class Integrator {
     std::cout << GridLogIntegrator << "Integrator initial action\n";
 
     RealD H = - FieldImplementation::FieldSquareNorm(P)/HMC_MOMENTUM_DENOMINATOR; // - trace (P*P)/denom
-
+#ifdef PRINT_SNAPSHOTS
+    std::cout << GridLogMessage << "Initial Half Mom Squared: " << H << std::endl;
+#endif
+    
     RealD Hterm;
 
     // Actions
diff --git a/Grid/qcd/smearing/GaugeConfigurationMasked.h b/Grid/qcd/smearing/GaugeConfigurationMasked.h
index 73243be207..0d79c1775f 100644
--- a/Grid/qcd/smearing/GaugeConfigurationMasked.h
+++ b/Grid/qcd/smearing/GaugeConfigurationMasked.h
@@ -4,18 +4,23 @@
   @brief Declares the GaugeConfiguration class
 */
 #pragma once
+#ifndef GAUGE_CONFIG_MASK_H
 
 NAMESPACE_BEGIN(Grid);
 
-
+//#define DEBUG
 template<class T> void Dump(const Lattice<T> & lat,
 			    std::string s,
 			    Coordinate site = Coordinate({0,0,0,0}))
 {
   typename T::scalar_object tmp;
+  if (lat.Checkerboard() != lat.Grid()->CheckerBoard(site))
+    site = Coordinate({0,0,0,1});
   peekSite(tmp,lat,site);
   std::cout << " Dump "<<s<<" "<<tmp<<std::endl;
 }
+
+
 /*!
   @brief Smeared configuration masked container
   Modified for a multi-subset smearing (aka Luscher Flowed HMC)
@@ -25,24 +30,33 @@ class SmearedConfigurationMasked : public SmearedConfiguration<Gimpl>
 {
 public:
   INHERIT_GIMPL_TYPES(Gimpl);
-
+  
 private:
   // These live in base class
   //  const unsigned int smearingLevels;
   //  Smear_Stout<Gimpl> *StoutSmearing;
   //  std::vector<GaugeField> SmearedSet;
-  
-  std::vector<LatticeLorentzComplex> masks;
 
+  GridCartesian         * UGrid; // keep a copy of the grid 
+  GridRedBlackCartesian * UrbGrid; // keep a copy of the redblack grid for life of object
+  PaddedCell              Ghost;
+
+  std::vector<LatticeLorentzComplex> masks;
+  std::vector<int> cbs;
+  std::vector<GeneralLocalStencil> gStencils;
+  std::vector<GeneralLocalStencil> gStencils_smear;
+  
   typedef typename SU3Adjoint::AMatrix AdjMatrix;
   typedef typename SU3Adjoint::LatticeAdjMatrix  AdjMatrixField;
   typedef typename SU3Adjoint::LatticeAdjVector  AdjVectorField;
-
+  typedef typename SU3::vAlgebraMatrix vAlgebraMatrix;
+  
+  
   void BaseSmearDerivative(GaugeField& SigmaTerm,
 			   const GaugeField& iLambda,
 			   const GaugeField& U,
 			   int mmu, RealD rho)
-  {
+  {GRID_TRACE("BaseSmearDerivative");
     // Reference
     // Morningstar, Peardon, Phys.Rev.D69,054501(2004)
     // Equation 75
@@ -58,6 +72,9 @@ class SmearedConfigurationMasked : public SmearedConfiguration<Gimpl>
     GaugeLinkField sh_field(grid), temp_Sigma(grid);
     Real rho_munu, rho_numu;
 
+    RealD t = 0;
+    t-=usecond();
+    
     rho_munu = rho;
     rho_numu = rho;
     for(int mu = 0; mu < Nd; ++mu){
@@ -100,65 +117,975 @@ class SmearedConfigurationMasked : public SmearedConfiguration<Gimpl>
 	//	staple = Zero();
 	sh_field = Cshift(U_nu, mu, 1);
 
-	temp_Sigma = Zero();
+	temp_Sigma = Zero();
+
+	if ( mu == mmu )
+	  temp_Sigma = -rho_munu*adj(sh_field)*adj(U_mu)*iLambda_mu*U_nu;
+
+	if ( nu == mmu ) {
+	  temp_Sigma += rho_numu*adj(sh_field)*adj(U_mu)*iLambda_nu*U_nu;
+
+	  u_tmp = adj(U_nu)*iLambda_nu;
+	  sh_field = Cshift(u_tmp, mu, 1);
+	  temp_Sigma += -rho_numu*sh_field*adj(U_mu)*U_nu;
+	}
+	
+	sh_field = Cshift(temp_Sigma, nu, -1);
+	Gimpl::AddLink(SigmaTerm, sh_field, mu);
+
+      }
+    }
+
+    t+=usecond();
+    std::cout << GridLogMessage << " BaseSmearDerivative " << t/1e3 << " ms " << std::endl;
+    
+  }
+#ifdef DEBUG
+  //for debugging
+  void BaseSmear_cb(GaugeLinkField& Cup, const GaugeField& U,int mu,RealD rho) {
+    GRID_TRACE("BaseSmear_cb");
+    GridBase *grid = U.Grid();
+    GridBase *hgrid = Cup.Grid();
+    GaugeLinkField tmp_stpl(grid);
+    GaugeLinkField tmp_stpl_eo(hgrid);
+    WilsonLoops<Gimpl> WL;
+    int cb = Cup.Checkerboard();
+    RealD t = 0;
+
+    t-=usecond();
+    Cup = Zero();
+    for(int nu=0; nu<Nd; ++nu){
+      if (nu != mu) {
+        // get the staple in direction mu, nu
+        WL.Staple(tmp_stpl, U, mu, nu);  //nb staple conventions of IroIro and Grid differ by a dagger
+        pickCheckerboard(cb,tmp_stpl_eo,tmp_stpl); // ideally, compute tmp_stpl only on the current checkerboard
+        Cup += adj(tmp_stpl_eo*rho);
+      }
+    }
+    t+=usecond();
+    std::cout << GridLogMessage << " BaseSmear " << t/1e3 << " ms " << std::endl;
+  }
+#endif
+  void BaseSmear(GaugeLinkField& Cup, const GaugeField& U,int mu,RealD rho) {
+    GRID_TRACE("BaseSmear");
+    GridBase *grid = U.Grid();
+    GaugeLinkField tmp_stpl(grid);
+    WilsonLoops<Gimpl> WL;
+    RealD t = 0;
+
+    t-=usecond();
+    Cup = Zero();
+    for(int nu=0; nu<Nd; ++nu){
+      if (nu != mu) {
+	// get the staple in direction mu, nu
+	WL.Staple(tmp_stpl, U, mu, nu);  //nb staple conventions of IroIro and Grid differ by a dagger
+	Cup += adj(tmp_stpl*rho);
+      }
+    }
+    t+=usecond();
+    std::cout << GridLogMessage << " BaseSmear " << t/1e3 << " ms " << std::endl;
+  }
+  
+  // Assume: gU is extended gauge field with its boundary of size 1
+  void BaseSmear_ghost(GaugeLinkField& Cup, const GaugeField& gU,int mu,RealD rho) {
+    GRID_TRACE("BaseSmear_ghost");
+    GridBase *ggrid = gU.Grid();
+    GaugeLinkField tmp_stpl(ggrid);
+    int cb = Cup.Checkerboard();
+    RealD t = 0;
+
+    //std::cout << GridLogMessage << "DEBUG: BaseSmear_ghost " << ggrid->Nsimd()<<std::endl;//debug
+    t-=usecond();
+    autoView( tmp_stpl_v , tmp_stpl, AcceleratorWrite);
+    autoView( gU_v , gU, AcceleratorRead);
+    autoView( gStencil_v, gStencils_smear[mu], AcceleratorRead);
+    accelerator_for(ss, ggrid->oSites(), ggrid->Nsimd(), {
+	typedef decltype(coalescedRead(tmp_stpl_v[0])) LinkMat;
+
+	LinkMat tmp = Zero();
+	for(int nu=0; nu<Nd; ++nu){
+	  int inc = 6*(nu - (mu<=nu));
+	  if (nu != mu) {
+	    GeneralStencilEntry const* e = gStencil_v.GetEntry(0+inc,ss);
+	    auto U_nu_x = coalescedReadGeneralPermute(gU_v[e->_offset], e->_permute, Nd)(nu)();
+	    e = gStencil_v.GetEntry(1+inc,ss);
+	    auto U_mu_xpnu = coalescedReadGeneralPermute(gU_v[e->_offset], e->_permute, Nd)(mu)();
+	    e = gStencil_v.GetEntry(2+inc,ss);
+	    auto Udag_nu_xpmu = adj(coalescedReadGeneralPermute(gU_v[e->_offset], e->_permute, Nd))(nu)();
+      
+	    tmp()() = tmp()() + U_nu_x * U_mu_xpnu * Udag_nu_xpmu;
+
+	    e = gStencil_v.GetEntry(3+inc,ss);
+	    auto Udag_nu_xmnu = adj(coalescedReadGeneralPermute(gU_v[e->_offset], e->_permute, Nd))(nu)();
+	    e = gStencil_v.GetEntry(4+inc,ss);
+	    auto U_mu_xmnu = coalescedReadGeneralPermute(gU_v[e->_offset], e->_permute, Nd)(mu)();
+	    e = gStencil_v.GetEntry(5+inc,ss);
+	    auto U_nu_xpmu_mnu = coalescedReadGeneralPermute(gU_v[e->_offset], e->_permute, Nd)(nu)();
+
+	    tmp()() = tmp()() + Udag_nu_xmnu * U_mu_xmnu * U_nu_xpmu_mnu;
+            
+	  }
+	}
+	coalescedWrite(tmp_stpl_v[ss],rho*tmp);
+      });
+    pickCheckerboard(cb,Cup,Ghost.Extract(tmp_stpl));
+    
+    t+=usecond();
+    std::cout << GridLogMessage << " BaseSmear_ghost " << t/1e3 << " ms " << std::endl;
+  }
+
+  // Adjoint vector to GaugeField force
+  void InsertForce(GaugeField &Fdet,AdjVectorField &Fdet_nu,int nu)
+  {
+    Complex ci(0,1);
+    GaugeLinkField Fdet_pol(Fdet.Grid());
+    RealD t = 0;
+
+    t-=usecond();
+    Fdet_pol=Zero();
+    for(int e=0;e<8;e++){
+      ColourMatrix te;
+      SU3::generator(e, te);
+      auto tmp=peekColour(Fdet_nu,e);
+      Fdet_pol=Fdet_pol + ci*tmp*te; // but norm of te is different.. why?
+    }
+    pokeLorentz(Fdet, Fdet_pol, nu);
+
+    t+=usecond();
+    std::cout << GridLogPerformance << " InsertForce " << t/1e3 << " ms " << std::endl;
+  }
+
+  void Compute_MpInvJx_dNxxdSy(const GaugeLinkField &PlaqL,const GaugeLinkField &PlaqR, AdjMatrixField MpInvJx,AdjVectorField &Fdet2 )
+  {
+    GRID_TRACE("Compute_MpInvJx_dNxxdSy");
+    int cb = PlaqL.Checkerboard();
+    GridBase *grid = PlaqL.Grid();
+    const int Ngen = SU3Adjoint::Dimension;
+    Complex ci(0,1);
+    RealD t;
+
+    t=-usecond();
+#if 0 // Due to Aurora Issue
+    autoView(Fdet2_v,Fdet2,AcceleratorWrite);
+    autoView(PlaqL_v,PlaqL,AcceleratorRead);
+    autoView(PlaqR_v,PlaqR,AcceleratorRead);
+    autoView(MpInvJx_v,MpInvJx,AcceleratorRead);
+    const int nsimd = vAlgebraMatrix::Nsimd();
+    accelerator_for(ss,grid->oSites(),nsimd,{
+	typedef decltype(coalescedRead(MpInvJx_v[0])) adj_mat;
+	typedef decltype(coalescedRead(Fdet2_v[0]))   adj_vec;
+	
+	adj_mat Dbc;
+	adj_vec Fdet;
+	ColourMatrix ta;
+	
+	for(int a=0;a<Ngen;a++) {
+	  // Qlat Tb = 2i Tb^Grid
+	  SU3::generator(a, ta);
+	  ta = 2.0 * ci * ta;
+	  auto UtaU = adj(PlaqL_v(ss))*ta*PlaqR_v(ss);
+	  SU3::LieAlgebraProject(Dbc,UtaU);
+	  Fdet()()(a) = traceProduct(MpInvJx_v(ss),Dbc)()()();
+	}
+	coalescedWrite(Fdet2_v[ss],Fdet);
+      });
+#else
+    autoView(Fdet2_v,Fdet2,AcceleratorWrite);
+    autoView(PlaqL_v,PlaqL,AcceleratorRead);
+    autoView(PlaqR_v,PlaqR,AcceleratorRead);
+    autoView(MpInvJx_v,MpInvJx,AcceleratorRead);
+    const int nsimd = vAlgebraMatrix::Nsimd();
+    accelerator_for2d(ss,grid->oSites(),a,Ngen,nsimd,{
+        typedef decltype(coalescedRead(MpInvJx_v[0])) adj_mat;
+        typedef decltype(coalescedRead(Fdet2_v[0]))   adj_vec;
+
+        adj_mat Dbc;
+        ColourMatrix ta;
+
+          // Qlat Tb = 2i Tb^Grid
+	SU3::generator(a, ta);
+	ta = 2.0 * ci * ta;
+	auto UtaU = adj(PlaqL_v(ss))*ta*PlaqR_v(ss);
+	SU3::LieAlgebraProject(Dbc,UtaU);
+
+        coalescedWrite(Fdet2_v[ss]()()(a),traceProduct(MpInvJx_v(ss),Dbc)()()());
+      });
+    
+#endif
+    t+=usecond();
+    //std::cout << GridLogMessage << "DEBUG: Compute_MpInvJx_dNxxdSy " << nsimd<<" " <<norm2(PlaqL)<<" "<<norm2(PlaqR)<<" "<<norm2(MpInvJx)<<" "<<norm2(Fdet2)<<std::endl;//debug
+    std::cout << GridLogPerformance << " Compute_MpInvJx_dNxxdSy " << t/1e3 <<" ms"<<std::endl;
+  }
+
+  void ComputeNxy(const GaugeLinkField &PlaqL,const GaugeLinkField &PlaqR,AdjMatrixField &NxAd)
+  {
+    GRID_TRACE("ComputeNxy");
+    GridBase *grid = PlaqL.Grid();
+    RealD t = 0;
+
+    t-=usecond();
+    autoView(NxAd_v,NxAd,AcceleratorWrite);
+    autoView(PlaqL_v,PlaqL,AcceleratorRead);
+    autoView(PlaqR_v,PlaqR,AcceleratorRead);
+    const int nsimd = vAlgebraMatrix::Nsimd();
+    accelerator_for(ss,grid->oSites(),nsimd,{
+        typedef decltype(coalescedRead(NxAd_v[0]))  adj_mat;
+        adj_mat NxAd_site;
+	SU3::LieAlgebraProject(NxAd_site,PlaqL_v(ss),PlaqR_v(ss));
+        coalescedWrite(NxAd_v[ss],NxAd_site);
+      });
+    t+=usecond();
+    //    std::cout << GridLogMessage << "DEBUG: CNxy" << nsimd<<std::endl;//debug    
+    std::cout << GridLogPerformance << " ComputeNxy " << t/1e3 <<" ms"<<std::endl;
+  }
+
+  void ApplyMask(GaugeField &U,int smr)
+  {
+    LatticeComplex tmp(U.Grid());
+    GaugeLinkField Umu(U.Grid());
+    for(int mu=0;mu<Nd;mu++){
+      Umu=PeekIndex<LorentzIndex>(U,mu);
+      tmp=PeekIndex<LorentzIndex>(masks[smr],mu);
+      Umu=Umu*tmp;
+      PokeIndex<LorentzIndex>(U, Umu, mu);
+    }
+  }
+public:
+
+  void logDetJacobianForceLevel(const GaugeField &U, GaugeField &force ,int smr)
+  {
+    GRID_TRACE("logDetJacobianForceLevel");
+    GridBase* grid = U.Grid();
+    GridBase* hgrid = UrbGrid; // For now, assume masking is based on red-black checkerboarding
+    assert(grid==UGrid);
+    //GaugeField C(grid);
+    GaugeField Umsk(grid);
+    std::vector<GaugeLinkField> Umu(Nd,grid);
+    GaugeLinkField Cmu(hgrid); // U and staple; C contains factor of epsilon
+    GaugeLinkField Zx(hgrid);  // U times Staple, contains factor of epsilon
+    GaugeLinkField Utmp(grid);
+    GaugeLinkField Ueo(hgrid);
+    GaugeLinkField PlaqL(hgrid);
+    GaugeLinkField PlaqR(hgrid);
+    const int Ngen = SU3Adjoint::Dimension;
+    ColourMatrix Ident;
+    const int Nsimd = vAlgebraMatrix::Nsimd();
+
+    AdjVectorField  dJdXe_nMpInv(hgrid); 
+    AdjVectorField  dJdXe_nMpInv_y(hgrid); 
+    AdjMatrixField  MpAd(hgrid);    // Mprime luchang's notes
+    AdjMatrixField  MpAdInv(hgrid); // Mprime inverse
+    AdjMatrixField  NxxAd(hgrid);    // Nxx in adjoint space
+    AdjMatrixField  JxAd(hgrid);     
+    AdjMatrixField  ZxAd(hgrid);
+
+    RealD t0 = usecond();
+    Ident = ComplexD(1.0);
+    for(int d=0;d<Nd;d++){
+      Umu[d] = peekLorentz(U, d);
+    }
+    int mu= (smr/2) %Nd;
+
+    ////////////////////////////////////////////////////////////////////////////////
+    // Mask the gauge field
+    ////////////////////////////////////////////////////////////////////////////////
+    int cb = cbs[smr];
+    auto mask=PeekIndex<LorentzIndex>(masks[smr],mu); // the cb mask
+
+    Umsk = U;
+    ApplyMask(Umsk,smr);
+    Utmp = peekLorentz(Umsk,mu);
+    pickCheckerboard(cb,Ueo,Utmp);
+    
+    GaugeField gU(grid);
+    std::vector<GaugeLinkField> gUmu(Nd,grid);
+    {
+      GRID_TRACE("ExchangePeriodic");
+      gU = Ghost.ExchangePeriodic(U);
+      for(int d=0; d<Nd;d++)
+	gUmu[d] = peekLorentz(gU, d);
+    }
+    GridBase       *ggrid = gUmu[0].Grid();
+    
+    Cmu.Checkerboard() = cb;
+    Zx.Checkerboard() = cb;
+    Ueo.Checkerboard() = cb;
+    PlaqL.Checkerboard() = cb;
+    PlaqR.Checkerboard() = cb;
+    MpAd.Checkerboard() = cb;
+    MpAdInv.Checkerboard() = cb;
+    dJdXe_nMpInv.Checkerboard() = cb;
+    dJdXe_nMpInv_y.Checkerboard() = cb;
+    NxxAd.Checkerboard() = cb;
+    JxAd.Checkerboard() = cb;
+    ZxAd.Checkerboard() = cb;
+    
+    ////////////////////////////////////////////////////////////////////////////////
+    // Retrieve the eps/rho parameter(s) -- could allow all different but not so far
+    ////////////////////////////////////////////////////////////////////////////////
+    double rho=this->StoutSmearing->SmearRho[1];
+    int idx=0;
+    for(int mu=0;mu<4;mu++){
+      for(int nu=0;nu<4;nu++){
+	if ( mu!=nu) assert(this->StoutSmearing->SmearRho[idx]==rho);
+	else         assert(this->StoutSmearing->SmearRho[idx]==0.0);
+	idx++;
+      }}
+    
+    //////////////////////////////////////////////////////////////////
+    // Assemble the N matrix
+    //////////////////////////////////////////////////////////////////
+    // Computes ALL the staples -- could compute one only and do it here
+    RealD time;
+    time=-usecond();
+    BaseSmear_ghost(Cmu, gU, mu, rho);
+#ifdef DEBUG
+    GaugeLinkField Cmu2(hgrid);
+    Cmu2.Checkerboard() = cb;
+    BaseSmear_cb(Cmu2, U, mu, rho);
+    std::cout << GridLogMessage << " DEBUG: BaseSmear " <<smr<<" "<<mu<<" "<<cb<<" "<<" simd "<<AdjMatrix::Nsimd()<<" "<<norm2(Cmu-Cmu2)<<std::endl;
+#endif
+    
+    //////////////////////////////////////////////////////////////////
+    // Assemble Luscher exp diff map J matrix 
+    //////////////////////////////////////////////////////////////////
+    // Ta so Z lives in Lie algabra
+    {GRID_TRACE("Zx");
+    Zx  = Ta(Cmu * adj(Ueo));
+    }
+    time+=usecond();
+    std::cout << GridLogMessage << "Z took "<<time<< " us"<<std::endl;
+    
+    // Move Z to the Adjoint Rep == make_adjoint_representation
+    time=-usecond();
+    {GRID_TRACE("ZxAd");
+    SU3Adjoint::make_adjoint_rep(ZxAd, Zx);
+    }
+    time+=usecond();
+    std::cout << GridLogPerformance << "ZxAd took "<<time<< " us"<<std::endl;
+#ifdef DEBUG
+    AdjMatrixField  ZxAd2(hgrid); ZxAd2.Checkerboard() = cb;
+    LatticeComplex  cplx(hgrid);  cplx.Checkerboard() = cb;
+    AdjMatrix TRb;
+    Complex ci(0,1);
+    ColourMatrix   tb;
+
+    ZxAd2 = Zero();
+    for(int b=0;b<8;b++) {
+      // Adj group sets traceless antihermitian T's -- Guido, really????
+      SU3::generator(b, tb);         // Fund group sets traceless hermitian T's
+      SU3Adjoint::generator(b,TRb);
+      TRb=-TRb;
+      cplx = 2.0*trace(ci*tb*Zx); // my convention 1/2 delta ba
+      ZxAd2 = ZxAd2 + cplx * TRb; // is this right? YES - Guido used Anti herm Ta's and with bloody wrong sign.
+      
+    }
+    std::cout << GridLogMessage << " DEBUG: ZxAd " <<smr<<" "<<mu<<" "<<cb<<" "<<" simd "<<AdjMatrix::Nsimd()<<" "<<norm2(ZxAd-ZxAd2)<<" "<<norm2(ZxAd)<<" "<<norm2(ZxAd2)<<std::endl;
+#endif
+    
+    //////////////////////////////////////
+    // J(x) = 1 + Sum_k=1..N (-Zac)^k/(k+1)!
+    //////////////////////////////////////
+    time=-usecond();
+    {GRID_TRACE("JxAd");
+    autoView(JxAd_v,JxAd,AcceleratorWrite);
+    autoView(ZxAd_v,ZxAd,AcceleratorRead);
+    accelerator_for(ss,hgrid->oSites(),Nsimd,{
+        typedef decltype(coalescedRead(JxAd_v[0])) adj_mat;
+        adj_mat X, JxAd_site;
+	RealD kpfac = 1;
+	
+	X=1.0;
+	JxAd_site = X;
+	for(int k=1;k<12;k++){
+	  X=-X*ZxAd_v(ss);
+	  kpfac = kpfac /(k+1);
+	  JxAd_site = JxAd_site + X * kpfac;
+	}
+        coalescedWrite(JxAd_v[ss],JxAd_site);
+      });
+    }
+    time+=usecond();
+    std::cout << GridLogMessage << "Jx took "<<time<< " us"<<std::endl;
+#ifdef DEBUG
+    std::cout << GridLogMessage << " DEBUG: Jx " <<smr<<" "<<mu<<" "<<cb<<" "<<" simd "<<AdjMatrix::Nsimd()<<" "<<Nsimd<<" "<<norm2(JxAd)<<std::endl;
+#endif
+
+    /////////////////////////////////////////////////////////////////
+    // NxxAd
+    /////////////////////////////////////////////////////////////////
+    time=-usecond();
+    {
+      GRID_TRACE("Plaq_for_Nxy");
+    PlaqL = Ident;
+    PlaqR = Ueo*adj(Cmu);
+    }
+    ComputeNxy(PlaqL,PlaqR,NxxAd);
+    time+=usecond();
+    std::cout << GridLogMessage << "ComputeNxy took "<<time<< " us"<<std::endl;
+#ifdef DEBUG
+    AdjMatrixField  NxxAd2(hgrid); NxxAd2.Checkerboard() = cb;
+    ComputeNxy(0,PlaqL,PlaqR,NxxAd2);
+    std::cout << GridLogMessage << " DEBUG: NxxAd " <<smr<<" "<<mu<<" "<<cb<<" "<<" simd "<<AdjMatrix::Nsimd()<<" "<<Nsimd<<" "<<norm2(NxxAd-NxxAd2)<<std::endl;
+#endif
+    ////////////////////////////
+    // Mab
+    ////////////////////////////
+    MpAd = Complex(1.0,0.0);
+    MpAd = MpAd - JxAd * NxxAd; 
+
+    /////////////////////////
+    // invert the 8x8
+    /////////////////////////
+    {GRID_TRACE("MpAdInv");
+    time=-usecond();
+    MpAdInv = Inverse_RealPart(MpAd);
+    MpAdInv.Checkerboard() = cb; //inside, it calls Lattice(GridBase *grid,ViewMode mode=AcceleratorWriteDiscard) in Lattice_base.h & sets checkerboard to 0
+    time+=usecond();
+    }
+    std::cout << GridLogPerformance << "MpAdInv took "<<time<< " us"<<std::endl;
+#ifdef DEBUG
+    std::cout << GridLogMessage << " DEBUG: MpAdInv " <<smr<<" "<<mu<<" "<<cb<<" "<<" simd "<<AdjMatrix::Nsimd()<<" "<<norm2(MpAdInv)<<std::endl;
+#endif
+
+    /////////////////////////////////////////////////////////////////
+    // M'^{-1} J(X)_ad
+    /////////////////////////////////////////////////////////////////
+    AdjMatrixField MpInvJx(hgrid);      MpInvJx.Checkerboard() = cb; 
+    AdjMatrixField MpInvJx_nu(hgrid);   MpInvJx_nu.Checkerboard() = cb;
+    {GRID_TRACE("Mp_inv_Jx");
+    MpInvJx = (-1.0)*MpAdInv * JxAd;// rho is on the plaq factor
+    }
+#ifdef DEBUG
+    std::cout << GridLogMessage << " DEBUG: MpAdInvJx " <<smr<<" "<<mu<<" "<<cb<<" "<<" simd "<<AdjMatrix::Nsimd()<<" "<<norm2(MpInvJx)<<std::endl;
+#endif
+    RealD t3a = usecond();
+    /////////////////////////////////////////////////////////////////
+    // dJ(x)/dxe N M'^{-1}  
+    /////////////////////////////////////////////////////////////////
+    time=-usecond();
+    {GRID_TRACE("dJdX_nMpinv_combined");
+
+    autoView(dJdXe_nMpInv_v,dJdXe_nMpInv,AcceleratorWrite);
+    autoView(ZxAd_v,ZxAd,AcceleratorRead);
+    autoView(NxxAd_v,NxxAd,AcceleratorRead);
+    autoView(MpAdInv_v,MpAdInv,AcceleratorRead);
+    accelerator_for(ss,hgrid->oSites(),Nsimd,{
+	typedef decltype(coalescedRead(ZxAd_v[0]))         adj_mat;
+	typedef decltype(coalescedRead(dJdXe_nMpInv_v[0])) adj_vec;
+      	adj_mat X, t2, dt2, t3, dt3, aunit, nMpInv_site;
+	adj_vec dJdXe_nMpInv_site;
+	iVector<adj_mat,Ngen> iTas;
+	iVector<adj_mat,Ngen> dJdX_b;
+
+	for(int b=0;b<Ngen;b++){
+	  SU3Adjoint::generator(b, iTas(b));
+	  dJdX_b(b) = iTas(b);
+	}
+	aunit = ComplexD(1.0);
+	X  = (-1.0)*ZxAd_v(ss);
+	t2 = X;
+	for (int j = 12; j > 1; --j) {
+	  t3  = t2*(1.0 / (j + 1))  + aunit;
+	  t2  = X * t3;
+	  for(int b=0;b<Ngen;b++){
+	    dJdX_b(b)= iTas(b) * t3 + X * dJdX_b(b)*(1.0 / (j + 1));
+	  }
+	}
+	nMpInv_site= NxxAd_v(ss) * MpAdInv_v(ss);
+	for(int e=0;e<Ngen;e++){
+	  dJdXe_nMpInv_site()()(e) = traceProduct((-1.0)*dJdX_b(e),nMpInv_site)()()();
+	}
+	coalescedWrite(dJdXe_nMpInv_v[ss],dJdXe_nMpInv_site);
+      });
+    }//make sure view object is closed
+    time += usecond();
+    std::cout << GridLogMessage << "dJdX_nMpinv_combined took "<<time<< " us"<<std::endl;
+#ifdef DEBUG
+    std::cout << GridLogMessage << " DEBUG: dJdX_nMpinv_combined " <<smr<<" "<<mu<<" "<<cb<<" "<<" simd "<<AdjMatrix::Nsimd()<<" "<<norm2(dJdXe_nMpInv)<<std::endl;
+#endif
+    RealD t3b = usecond();
+    
+    AdjVectorField  Fdet1_mu(grid);
+    AdjVectorField  Fdet2_mu(grid);
+    AdjVectorField  Fdet1_nu(grid);
+    AdjVectorField  Fdet2_nu(grid);
+
+    AdjVectorField  FdetV(hgrid);       FdetV.Checkerboard() = cb;
+    AdjVectorField  Fdet1_nu_eo(hgrid); Fdet1_nu_eo.Checkerboard() = cb;       
+    AdjVectorField  Fdet1_nu_oe(hgrid); Fdet1_nu_oe.Checkerboard() = (cb+1)%2; 
+    AdjVectorField  Fdet2_nu_eo(hgrid); Fdet2_nu_eo.Checkerboard() = cb;       
+    AdjVectorField  Fdet2_nu_oe(hgrid); Fdet1_nu_oe.Checkerboard() = (cb+1)%2;
+    
+    AdjVectorField  Fdet1_mu_oe(hgrid); Fdet1_mu_oe.Checkerboard() = (cb+1)%2; Fdet1_mu_oe = Zero();
+    AdjVectorField  Fdet2_mu_oe(hgrid); Fdet2_mu_oe.Checkerboard() = (cb+1)%2; Fdet2_mu_oe = Zero();
+
+    // Set even part of Fdet1_mu & Fdet2_mu <- cb is referred to as even here, regardless of actual parity corredponding to cb
+    setCheckerboard(Fdet1_mu, (AdjVectorField) (transpose(NxxAd)*dJdXe_nMpInv)); 
+    Compute_MpInvJx_dNxxdSy(PlaqL,PlaqR,MpInvJx,FdetV);
+    setCheckerboard(Fdet2_mu,FdetV);
+
+    RealD t3c = usecond();
+    
+    //    dJdXe_nMpInv needs to multiply:
+    //       Nxx_mu (site local)                           (1)
+    //       Nxy_mu one site forward  in each nu direction (3)
+    //       Nxy_mu one site backward in each nu direction (3)
+    //       Nxy_nu 0,0  ; +mu,0; 0,-nu; +mu-nu   [ 3x4 = 12]
+    // 19 terms.
+    AdjMatrixField Nxy(hgrid);
+    //// force = Fdet1 + Fdet2
+    GaugeField Fdet1(grid);
+    GaugeField Fdet2(grid);
+
+    RealD t4 = usecond(), tLR = 0, tNxy = 0, tMJx = 0, t_ins=0, t_ck = 0, t_stencil=0, t_cshift=0;
+
+    GaugeLinkField  gPlaqL(ggrid), gPlaqR(ggrid);
+    
+    autoView( gPlaqL_v , gPlaqL, AcceleratorWrite);
+    autoView( gPlaqR_v , gPlaqR, AcceleratorWrite);
+    autoView( gU_mu_v , gUmu[mu], AcceleratorRead);
+    
+    for(int nu=0;nu<Nd;nu++){
+      
+      if (nu!=mu) {
+	GRID_TRACE("MuNuLoopBody");
+	autoView( gU_nu_v , gUmu[nu], AcceleratorRead);
+
+
+	///////////////// +ve nu /////////////////
+	//     __
+	//    |  |
+	//    x==    // nu polarisation -- clockwise
+
+	time=-usecond(); tLR -= usecond();
+
+	PlaqL=Ident;
+	{
+	  GRID_TRACE("Staple");
+	autoView( gStencil_v  , gStencils[mu*(Nd-1)*6+(nu-(mu<=nu))*6], AcceleratorRead);
+	accelerator_for(ss, ggrid->oSites(), ggrid->Nsimd(), {
+          GeneralStencilEntry const* e = gStencil_v.GetEntry(0,ss);
+          auto U_nu_x = coalescedReadGeneralPermute(gU_nu_v[e->_offset], e->_permute, Nd);
+          e = gStencil_v.GetEntry(1,ss);
+          auto U_mu_xpnu = coalescedReadGeneralPermute(gU_mu_v[e->_offset], e->_permute, Nd);
+          e = gStencil_v.GetEntry(2,ss);
+          auto Udag_nu_xpmu = adj(coalescedReadGeneralPermute(gU_nu_v[e->_offset], e->_permute, Nd));
+	  e = gStencil_v.GetEntry(3,ss);
+          auto Udag_mu_x = adj(coalescedReadGeneralPermute(gU_mu_v[e->_offset], e->_permute, Nd));
+
+          auto stencil_ss = (-rho) * U_nu_x * U_mu_xpnu * Udag_nu_xpmu * Udag_mu_x;
+
+          coalescedWrite(gPlaqR_v[ss],stencil_ss);
+	}
+        );
+	t_ck -= usecond();
+	pickCheckerboard(cb,PlaqR,Ghost.Extract(gPlaqR));
+	t_ck += usecond();
+	}
+	time+=usecond(); tLR += usecond();
+	std::cout << GridLogMessage << "PlaqLR took "<<time<< " us "<<" checkerboard_extract "<<t_ck<<" us"<<std::endl;
+#ifdef DEBUG
+	std::cout << GridLogMessage << " DEBUG: munuLoopBody PlqR" <<smr<<" "<<mu<<" "<<cb<<" "<<" simd "<<AdjMatrix::Nsimd()<<" "<< ggrid->Nsimd()<<" "<<norm2(PlaqR)<<std::endl;
+#endif
+	time=-usecond(); tNxy -= usecond();
+	PlaqL.Checkerboard() = cb; Nxy.Checkerboard() = cb; FdetV.Checkerboard() = cb;
+	
+	dJdXe_nMpInv_y = dJdXe_nMpInv;
+	ComputeNxy(PlaqL,PlaqR,Nxy);
+	
+	Fdet1_nu_eo = transpose(Nxy)*dJdXe_nMpInv_y;
+	time+=usecond(); tNxy += usecond();
+	std::cout << GridLogMessage << "ComputeNxy (occurs 6x) took "<<time<< " us"<<std::endl;
+
+	time=-usecond(); tMJx -= usecond();
+	PlaqR=(-1.0)*PlaqR;
+	Compute_MpInvJx_dNxxdSy(PlaqL,PlaqR,MpInvJx,FdetV);
+	Fdet2_nu_eo = FdetV;
+	time+=usecond(); tMJx += usecond();
+	std::cout << GridLogMessage << "Compute_MpInvJx_dNxxSy (occurs 6x) took "<<time<< " us"<<std::endl;
+#ifdef DEBUG
+	std::cout << GridLogMessage << " DEBUG: munuLoopBody FdetV " <<smr<<" "<<mu<<" "<<cb<<" "<<" simd "<<AdjMatrix::Nsimd()<<" "<<norm2(FdetV)<<std::endl;
+#endif
+	//    x==
+	//    |  |
+	//    .__|    // nu polarisation -- anticlockwise
+
+	tLR -= usecond();
+	{
+	  GRID_TRACE("Staple");
+	autoView( gStencil_v  , gStencils[mu*(Nd-1)*6+(nu-(mu<=nu))*6+1], AcceleratorRead);
+        accelerator_for(ss, ggrid->oSites(), ggrid->Nsimd(), {
+          GeneralStencilEntry const* e = gStencil_v.GetEntry(0,ss);
+          auto U_nu_y = coalescedReadGeneralPermute(gU_nu_v[e->_offset], e->_permute, Nd);
+          e = gStencil_v.GetEntry(1,ss);
+          auto Udag_mu_ymmupnu = adj(coalescedReadGeneralPermute(gU_mu_v[e->_offset], e->_permute, Nd));
+          e = gStencil_v.GetEntry(2,ss);
+          auto Udag_nu_ymmu = adj(coalescedReadGeneralPermute(gU_nu_v[e->_offset], e->_permute, Nd));
+	  auto stencil_ss = (rho) * U_nu_y * Udag_mu_ymmupnu * Udag_nu_ymmu;
+	  coalescedWrite(gPlaqR_v[ss],stencil_ss);
+	  
+          e = gStencil_v.GetEntry(3,ss);
+          auto Udag_mu_ymmu = adj(coalescedReadGeneralPermute(gU_mu_v[e->_offset], e->_permute, Nd));
+	  stencil_ss = Udag_mu_ymmu;
+          coalescedWrite(gPlaqL_v[ss],stencil_ss);
+        }
+        );
+	pickCheckerboard((cb+1)%2,PlaqR,Ghost.Extract(gPlaqR));
+	pickCheckerboard((cb+1)%2,PlaqL,Ghost.Extract(gPlaqL));
+	}
+	tLR += usecond();
+
+	tNxy -= usecond();
+	Nxy.Checkerboard() = (cb+1)%2; 	FdetV.Checkerboard() = (cb+1)%2;
+	t_cshift-= usecond();
+	dJdXe_nMpInv_y = Cshift(dJdXe_nMpInv,mu,-1);
+	t_cshift += usecond();
+	ComputeNxy(PlaqL, PlaqR,Nxy);
+	Fdet1_nu_oe = transpose(Nxy)*dJdXe_nMpInv_y;
+	tNxy += usecond();
+
+	tMJx -= usecond();
+	MpInvJx_nu = Cshift(MpInvJx,mu,-1);
+	Compute_MpInvJx_dNxxdSy(PlaqL,PlaqR,MpInvJx_nu,FdetV);
+	Fdet2_nu_oe = FdetV;
+	tMJx += usecond();
+	
+	///////////////// -ve nu /////////////////
+	//  __
+	// |  |
+	// x==          // nu polarisation -- clockwise
+	tLR -= usecond();
+	{ GRID_TRACE("Staple");
+	autoView( gStencil_v  , gStencils[mu*(Nd-1)*6+(nu-(mu<=nu))*6+2], AcceleratorRead);
+        accelerator_for(ss, ggrid->oSites(), ggrid->Nsimd(), {
+          GeneralStencilEntry const* e = gStencil_v.GetEntry(0,ss);
+          auto U_mu_y = coalescedReadGeneralPermute(gU_mu_v[e->_offset], e->_permute, Nd);
+          e = gStencil_v.GetEntry(1,ss);
+          auto U_nu_ypmu = coalescedReadGeneralPermute(gU_nu_v[e->_offset], e->_permute, Nd);
+          e = gStencil_v.GetEntry(2,ss);
+          auto Udag_mu_ypnu = adj(coalescedReadGeneralPermute(gU_mu_v[e->_offset], e->_permute, Nd));
+
+          auto stencil_ss = (rho) * U_mu_y * U_nu_ypmu * Udag_mu_ypnu;
+          coalescedWrite(gPlaqL_v[ss],stencil_ss);
+        }
+        );
+	pickCheckerboard((cb+1)%2,PlaqL,Ghost.Extract(gPlaqL)); 
+        pickCheckerboard((cb+1)%2,PlaqR,Umu[nu]);
+	}
+	tLR += usecond();
+
+	tNxy -= usecond();
+	t_cshift-=usecond();
+	dJdXe_nMpInv_y = Cshift(dJdXe_nMpInv,nu,1);
+	t_cshift+=usecond();
+	ComputeNxy(PlaqL,PlaqR,Nxy);
+	Fdet1_nu_oe = Fdet1_nu_oe + transpose(Nxy)*dJdXe_nMpInv_y;
+	tNxy += usecond();
+
+	tMJx -= usecond();
+	t_cshift-=usecond();
+	MpInvJx_nu = Cshift(MpInvJx,nu,1);
+	t_cshift+=usecond();
+	Compute_MpInvJx_dNxxdSy(PlaqL,PlaqR,MpInvJx_nu,FdetV);
+	Fdet2_nu_oe = Fdet2_nu_oe+FdetV;
+	tMJx += usecond();
+      
+	// x==
+	// |  |
+	// |__|         // nu polarisation
+	tLR -= usecond();
+	{
+	  GRID_TRACE("Staple");
+	autoView( gStencil_v  , gStencils[mu*(Nd-1)*6+(nu-(mu<=nu))*6+3], AcceleratorRead);
+        accelerator_for(ss, ggrid->oSites(), ggrid->Nsimd(), {
+          GeneralStencilEntry const* e = gStencil_v.GetEntry(0,ss);
+          auto U_nu_y = coalescedReadGeneralPermute(gU_nu_v[e->_offset], e->_permute, Nd);
+          e = gStencil_v.GetEntry(1,ss);
+          auto Udag_mu_ymmupnu = adj(coalescedReadGeneralPermute(gU_mu_v[e->_offset], e->_permute, Nd));
+          auto stencil_ss = (-rho) * U_nu_y * Udag_mu_ymmupnu;
+          coalescedWrite(gPlaqL_v[ss],stencil_ss);
+	  
+          e = gStencil_v.GetEntry(2,ss);
+          auto Udag_mu_ymmu = adj(coalescedReadGeneralPermute(gU_mu_v[e->_offset], e->_permute, Nd));
+          e = gStencil_v.GetEntry(3,ss);
+          auto U_nu_ymmu = coalescedReadGeneralPermute(gU_nu_v[e->_offset], e->_permute, Nd);
+          stencil_ss = Udag_mu_ymmu * U_nu_ymmu;
+          coalescedWrite(gPlaqR_v[ss],stencil_ss);
+        }
+        );
+	pickCheckerboard(cb,PlaqL,Ghost.Extract(gPlaqL));
+        pickCheckerboard(cb,PlaqR,Ghost.Extract(gPlaqR));
+	}
+	tLR += usecond();
+
+	tNxy -= usecond();
+	Nxy.Checkerboard() = cb;  FdetV.Checkerboard() = cb;
+	t_cshift-=usecond();
+	dJdXe_nMpInv_y = Cshift(dJdXe_nMpInv,mu,-1);
+	dJdXe_nMpInv_y = Cshift(dJdXe_nMpInv_y,nu,1);
+	t_cshift+=usecond();
+	
+	ComputeNxy(PlaqL,PlaqR,Nxy);
+	Fdet1_nu_eo = Fdet1_nu_eo + transpose(Nxy)*dJdXe_nMpInv_y;
+	tNxy += usecond();
+
+	tMJx -= usecond();
+	t_cshift-=usecond();
+	MpInvJx_nu = Cshift(MpInvJx,mu,-1);
+	MpInvJx_nu = Cshift(MpInvJx_nu,nu,1);
+	t_cshift+=usecond();
+	Compute_MpInvJx_dNxxdSy(PlaqL,PlaqR,MpInvJx_nu,FdetV);
+	Fdet2_nu_eo = Fdet2_nu_eo+FdetV;
+	tMJx += usecond();
+	
+	/////////////////////////////////////////////////////////////////////
+	// Set up the determinant force contribution in 3x3 algebra basis
+	/////////////////////////////////////////////////////////////////////
+	t_ins -= usecond();
+	setCheckerboard(Fdet1_nu, Fdet1_nu_eo);
+	setCheckerboard(Fdet1_nu, Fdet1_nu_oe);
+	InsertForce(Fdet1,Fdet1_nu,nu);
+	setCheckerboard(Fdet2_nu, Fdet2_nu_eo);
+        setCheckerboard(Fdet2_nu, Fdet2_nu_oe);
+	InsertForce(Fdet2,Fdet2_nu,nu);
+	t_ins+= usecond();
+	
+	//////////////////////////////////////////////////
+	// Parallel direction terms
+	//////////////////////////////////////////////////
+
+	Nxy.Checkerboard() = (cb+1)%2; FdetV.Checkerboard() = (cb+1)%2;
+	
+        //     __
+	//    |  "
+	//    |__"x    // mu polarisation
+	tLR -= usecond();
+	{
+	  GRID_TRACE("Staple");
+	autoView( gStencil_v  , gStencils[mu*(Nd-1)*6+(nu-(mu<=nu))*6+4], AcceleratorRead);
+        accelerator_for(ss, ggrid->oSites(), ggrid->Nsimd(), {
+          GeneralStencilEntry const* e = gStencil_v.GetEntry(0,ss);
+          auto U_mu_y = coalescedReadGeneralPermute(gU_mu_v[e->_offset], e->_permute, Nd);
+          e = gStencil_v.GetEntry(1,ss);
+          auto Udag_nu_ypmumnu = adj(coalescedReadGeneralPermute(gU_nu_v[e->_offset], e->_permute, Nd));
+          e = gStencil_v.GetEntry(2,ss);
+          auto Udag_mu_ymnu = adj(coalescedReadGeneralPermute(gU_mu_v[e->_offset], e->_permute, Nd));
+          auto stencil_ss = (-rho) * U_mu_y * Udag_nu_ypmumnu * Udag_mu_ymnu;
+          coalescedWrite(gPlaqL_v[ss],stencil_ss);
+
+          e = gStencil_v.GetEntry(3,ss);
+          auto Udag_nu_ymnu = adj(coalescedReadGeneralPermute(gU_nu_v[e->_offset], e->_permute, Nd));
+          stencil_ss = Udag_nu_ymnu;
+          coalescedWrite(gPlaqR_v[ss],stencil_ss);
+        }
+        );
+	pickCheckerboard((cb+1)%2,PlaqL,Ghost.Extract(gPlaqL));
+        pickCheckerboard((cb+1)%2,PlaqR,Ghost.Extract(gPlaqR));
+	}
+	tLR += usecond();
+
+	tNxy -= usecond();
+	t_cshift-=usecond();
+	dJdXe_nMpInv_y = Cshift(dJdXe_nMpInv,nu,-1);
+	t_cshift+=usecond();
+	
+	ComputeNxy(PlaqL,PlaqR,Nxy);
+	Fdet1_mu_oe = Fdet1_mu_oe + transpose(Nxy)*dJdXe_nMpInv_y;
+	tNxy += usecond();
+
+	tMJx -= usecond();
+	t_cshift-=usecond();
+	MpInvJx_nu = Cshift(MpInvJx,nu,-1);
+	t_cshift+=usecond();
+	Compute_MpInvJx_dNxxdSy(PlaqL,PlaqR,MpInvJx_nu,FdetV);
+	Fdet2_mu_oe = Fdet2_mu_oe+FdetV;
+	tMJx += usecond();
+
+	//  __
+	// "  |
+	// x__|          // mu polarisation
+	tLR -= usecond();
+	{
+	  GRID_TRACE("Staple");
+	autoView( gStencil_v  , gStencils[mu*(Nd-1)*6+(nu-(mu<=nu))*6+5], AcceleratorRead);
+        accelerator_for(ss, ggrid->oSites(), ggrid->Nsimd(), {
+          GeneralStencilEntry const* e = gStencil_v.GetEntry(0,ss);
+          auto U_mu_y = coalescedReadGeneralPermute(gU_mu_v[e->_offset], e->_permute, Nd);
+          e = gStencil_v.GetEntry(1,ss);
+          auto U_nu_ypmu = coalescedReadGeneralPermute(gU_nu_v[e->_offset], e->_permute, Nd);
+          e = gStencil_v.GetEntry(2,ss);
+          auto Udag_mu_ypnu = adj(coalescedReadGeneralPermute(gU_mu_v[e->_offset], e->_permute, Nd));
+
+          auto stencil_ss = (-rho) * U_mu_y * U_nu_ypmu * Udag_mu_ypnu;
+          coalescedWrite(gPlaqL_v[ss],stencil_ss);
+        }
+        );
+        pickCheckerboard((cb+1)%2,PlaqL,Ghost.Extract(gPlaqL));
+	pickCheckerboard((cb+1)%2,PlaqR,Umu[nu]);
+	}
+	tLR += usecond();
+
+	tNxy -= usecond();
+	t_cshift-=usecond();
+	dJdXe_nMpInv_y = Cshift(dJdXe_nMpInv,nu,1);
+	t_cshift+=usecond();
+	ComputeNxy(PlaqL,PlaqR,Nxy);
+	Fdet1_mu_oe = Fdet1_mu_oe + transpose(Nxy)*dJdXe_nMpInv_y;
+	tNxy += usecond();
+
+	tMJx -= usecond();
+	t_cshift-=usecond();
+	MpInvJx_nu = Cshift(MpInvJx,nu,1);
+	t_cshift+=usecond();
+	Compute_MpInvJx_dNxxdSy(PlaqL,PlaqR,MpInvJx_nu,FdetV);
+	Fdet2_mu_oe = Fdet2_mu_oe+FdetV;
+	tMJx += usecond();
+	
+      }
+    }
+    RealD t5 = usecond();
+    setCheckerboard(Fdet1_mu, Fdet1_mu_oe);
+    setCheckerboard(Fdet2_mu, Fdet2_mu_oe);
+    InsertForce(Fdet1,Fdet1_mu,mu);
+    InsertForce(Fdet2,Fdet2_mu,mu);
+
+    force= (-0.5)*( Fdet1 + Fdet2);
+    RealD t1 = usecond();
+    std::cout << GridLogMessage << " logDetJacobianForce t3-t0 "<<t3a-t0<<" us "<<std::endl;
+    std::cout << GridLogMessage << " logDetJacobianForce t4-t3b dJdXe_nMpInv+even part of Fdet1_mu & Fdet2_mu "<<t4-t3b<<" us "<<std::endl;
+    std::cout << GridLogMessage << " logDetJacobianForce t3b-t3a dJdXe_nMpInv  "<<t3b-t3a<<" us "<<std::endl;
+    std::cout << GridLogMessage << " logDetJacobianForce t3c-t3b even part of Fdet1_mu & Fdet2_mu  "<<t3c-t3b<<" us "<<std::endl;
+    std::cout << GridLogMessage << " logDetJacobianForce t5-t4 mu nu loop "<<t5-t4<<" us Plaq "
+	      <<tLR/1e3<<" ms Nxy "<<tNxy/1e3<<" ms MpInvJx_dNxxdSy "<<tMJx/1e3<<" ms "<<" insert_force "<<t_ins/1e3<< " ms Stencil "
+	      <<t_stencil/1e3<<" ms Cshift "<<t_cshift/1e3<<" ms"<<std::endl;
+    std::cout << GridLogMessage << " logDetJacobianForce t1-t5 "<<t1-t5<<" us "<<std::endl; // turn adj vec to SU3 force
+    std::cout << GridLogMessage << " logDetJacobianForce level took "<<t1-t0<<" us "<<std::endl;
+  }
+
+  //Assume: masking is based on red-black checkerboarding
+  RealD logDetJacobianLevel(const GaugeField &U,int smr)
+  {
+    GRID_TRACE("logDetJacobianLevel");
+    GridBase* grid = U.Grid();
+    GridBase* hgrid = UrbGrid; // For now, assume masking is based on red-black checkerboarding
+    assert(grid==UGrid);
+    GaugeField gU(grid);
+    std::vector<GaugeLinkField> gUmu(Nd,grid);
+    GaugeLinkField PlaqL(hgrid);
+    GaugeLinkField Z(hgrid);
+    GaugeLinkField Umu(hgrid), Cmu(hgrid);
+    LatticeComplex ln_det(hgrid); 
+    typedef typename SU3Adjoint::LatticeAdjMatrix  AdjMatrixField;
+    AdjMatrixField  Ncb(hgrid);
+    AdjMatrixField  Zac(hgrid); 
+    ColourMatrix Ident;
+    
+    RealD time=0, tN=0, tZ=0, t_M=0, tJ_lnDet=0, t_det=0, t_ln=0;
+    int mu= (smr/2) %Nd;
+
+    time -= usecond();
+    Ident = ComplexD(1.0);
+
+    int cb = cbs[smr]; //Assume: the applied mask is the cb mask
 
-	if ( mu == mmu )
-	  temp_Sigma = -rho_munu*adj(sh_field)*adj(U_mu)*iLambda_mu*U_nu;
+    Z.Checkerboard() = cb;
+    Cmu.Checkerboard() = cb;
+    PlaqL.Checkerboard() = cb;
+    Ncb.Checkerboard() = cb;
+    ln_det.Checkerboard() = cb; 
 
-	if ( nu == mmu ) {
-	  temp_Sigma += rho_numu*adj(sh_field)*adj(U_mu)*iLambda_nu*U_nu;
+    pickCheckerboard(cb,Umu,peekLorentz(U,mu));
 
-	  u_tmp = adj(U_nu)*iLambda_nu;
-	  sh_field = Cshift(u_tmp, mu, 1);
-	  temp_Sigma += -rho_numu*sh_field*adj(U_mu)*U_nu;
-	}
-	
-	sh_field = Cshift(temp_Sigma, nu, -1);
-	Gimpl::AddLink(SigmaTerm, sh_field, mu);
+    //////////////////////////////////////////////////////////////////
+    // Assemble the N matrix
+    //////////////////////////////////////////////////////////////////
 
-      }
+    tN -= usecond();
+    {GRID_TRACE("ExchangePeriodic");
+    gU = Ghost.ExchangePeriodic(U);
     }
-  }
-  
-  void BaseSmear(GaugeLinkField& Cup, const GaugeField& U,int mu,RealD rho) {
-    GridBase *grid = U.Grid();
-    GaugeLinkField tmp_stpl(grid);
-    WilsonLoops<Gimpl> WL;
-    Cup = Zero();
-    for(int nu=0; nu<Nd; ++nu){
-      if (nu != mu) {
-	// get the staple in direction mu, nu
-	WL.Staple(tmp_stpl, U, mu, nu);  //nb staple conventions of IroIro and Grid differ by a dagger
-	Cup += adj(tmp_stpl*rho);
-      }
+    double rho=this->StoutSmearing->SmearRho[1];
+    PlaqL = Ident;
+    BaseSmear_ghost(Cmu,gU,mu,rho);
+    ComputeNxy(PlaqL,Umu * adj(Cmu),Ncb);//we can expand function body and form a bigger accelerated loop
+    tN += usecond();
+
+    //////////////////////////////////////////////////////////////////
+    // Assemble Luscher exp diff map J matrix 
+    //////////////////////////////////////////////////////////////////
+    tZ -= usecond();
+    {GRID_TRACE("Z");
+    // Ta so Z lives in Lie algabra
+    Z  = Ta(Cmu * adj(Umu));
+    // Move Z to the Adjoint Rep == make_adjoint_representation
+    SU3Adjoint::make_adjoint_rep(Zac, Z);
     }
-  }
-  // Adjoint vector to GaugeField force
-  void InsertForce(GaugeField &Fdet,AdjVectorField &Fdet_nu,int nu)
-  {
-    Complex ci(0,1);
-    GaugeLinkField Fdet_pol(Fdet.Grid());
-    Fdet_pol=Zero();
-    for(int e=0;e<8;e++){
-      ColourMatrix te;
-      SU3::generator(e, te);
-      auto tmp=peekColour(Fdet_nu,e);
-      Fdet_pol=Fdet_pol + ci*tmp*te; // but norm of te is different.. why?
+    tZ += usecond();
+
+    tJ_lnDet -= usecond();
+    t_M -= usecond();
+    {GRID_TRACE("Mab");
+    autoView(ln_det_v,ln_det,AcceleratorWrite);
+    autoView(Zac_v,Zac,AcceleratorRead);
+    autoView(Ncb_v,Ncb,AcceleratorRead);
+    accelerator_for(ss,hgrid->oSites(),hgrid->Nsimd(),{
+	typedef decltype(coalescedRead(Zac_v(0)))    adj_mat;
+
+	adj_mat X, Jac, Mab_ss;
+	RealD kpfac = 1;
+	
+	//////////////////////////////////////
+	// J(x) = 1 + Sum_k=1..N (-Zac)^k/(k+1)!
+	//////////////////////////////////////
+	X=1.0; 
+	Jac = X;
+	for(int k=1;k<12;k++){
+	  X=(-1.0)*X*Zac_v(ss);
+	  kpfac = kpfac /((RealD) (k+1));
+	  Jac = Jac + X * kpfac;
+	}
+	
+	////////////////////////////
+	// Mab
+	////////////////////////////
+	Mab_ss = Complex(1.0,0.0);
+	Mab_ss = Mab_ss - Jac * Ncb_v(ss);
+
+
+	////////////////////////////
+	// ln det
+	////////////////////////////
+	auto detD = Determinant(Mab_ss);
+	coalescedWrite(ln_det_v[ss],log(detD));
+      });
     }
-    pokeLorentz(Fdet, Fdet_pol, nu);
+    t_M+=usecond();
+
+    tJ_lnDet += usecond();
+    Complex result = sum(ln_det);
+    time += usecond();
+    std::cout << GridLogMessage << " logDetJacobianLevel " << time/1e3 <<" ms N "<<tN/1e3<<" ms Z "<<tZ/1e3<<" ms J_lnDetM " <<tJ_lnDet/1e3
+	      <<" ms Det " <<t_det/1e3<<" ms log "<<t_ln/1e3<<" ms M "<<t_M/1e3<<std::endl;
+
+    return result.real();
   }
-  void Compute_MpInvJx_dNxxdSy(const GaugeLinkField &PlaqL,const GaugeLinkField &PlaqR, AdjMatrixField MpInvJx,AdjVectorField &Fdet2 )
+
+  /*------------------    old implementation for testing      #########################   -------------------------------- */
+  void Compute_MpInvJx_dNxxdSy(int old, const GaugeLinkField &PlaqL,const GaugeLinkField &PlaqR, AdjMatrixField MpInvJx,AdjVectorField &Fdet2 )
   {
-    GaugeLinkField UtaU(PlaqL.Grid());
-    GaugeLinkField D(PlaqL.Grid());
-    AdjMatrixField Dbc(PlaqL.Grid());
-    AdjMatrixField Dbc_opt(PlaqL.Grid());
-    LatticeComplex tmp(PlaqL.Grid());
+    GRID_TRACE("Compute_MpInvJx_dNxxdSy_old");
+    int cb = PlaqL.Checkerboard();
+    GaugeLinkField UtaU(PlaqL.Grid());         UtaU.Checkerboard() = cb;
+    GaugeLinkField D(PlaqL.Grid());            D.Checkerboard() = cb;
+    AdjMatrixField Dbc(PlaqL.Grid());          Dbc.Checkerboard() = cb;
+    AdjMatrixField Dbc_opt(PlaqL.Grid());      Dbc_opt.Checkerboard() = cb;
+    LatticeComplex tmp(PlaqL.Grid());          tmp.Checkerboard() = cb;
     const int Ngen = SU3Adjoint::Dimension;
     Complex ci(0,1);
     ColourMatrix   ta,tb,tc;
-    RealD t=0;
-    RealD tp=0;
+    RealD t=0, time;
+    RealD tp=0, tpl=0;
     RealD tta=0;
     RealD tpk=0;
     t-=usecond();
@@ -167,7 +1094,10 @@ class SmearedConfigurationMasked : public SmearedConfiguration<Gimpl>
       SU3::generator(a, ta);
       ta = 2.0 * ci * ta;
       // Qlat Tb = 2i Tb^Grid
-      UtaU= adj(PlaqL)*ta*PlaqR; // 6ms
+      {
+	GRID_TRACE("UtaU");
+	UtaU= adj(PlaqL)*ta*PlaqR; // 6ms
+      }
       tta+=usecond();
       ////////////////////////////////////////////
       // Could add this entire C-loop to a projection routine
@@ -178,9 +1108,19 @@ class SmearedConfigurationMasked : public SmearedConfiguration<Gimpl>
 	SU3::generator(c, tc);
 	tc = 2.0*ci*tc;
 	tp-=usecond(); 
-	D = Ta( tc *UtaU); // 2ms
+	{
+	  GRID_TRACE("D = Ta tc*UtaU");
+	  D = Ta( tc *UtaU); // 2ms
+	}
 #if 1
-	SU3::LieAlgebraProject(Dbc_opt,D,c); // 5.5ms
+	{
+	  GRID_TRACE("LieAlgebraProject");
+	  time=-usecond();
+	  SU3::LieAlgebraProject(Dbc_opt,D,c); // 5.5ms
+	  time+=usecond();
+	  tpl+=time;
+	  //std::cout << GridLogMessage << " LieAlgebraProject_in_Compute_MpInvJx " << a << " " << c << " "<<time/1e3<<" ms"<<std::endl;
+	}
 #else
 	for(int b=0;b<Ngen;b++){
 	  SU3::generator(b, tb);
@@ -190,55 +1130,82 @@ class SmearedConfigurationMasked : public SmearedConfiguration<Gimpl>
 #endif
 	tp+=usecond();
       }
-      //      Dump(Dbc_opt,"Dbc_opt");
-      //      Dump(Dbc,"Dbc");
+      //Dump(Dbc_opt,"Dbc_opt");
+      //Dump(Dbc,"Dbc");
       tpk-=usecond();
-      tmp = trace(MpInvJx * Dbc_opt);
-      PokeIndex<ColourIndex>(Fdet2,tmp,a);
+      {
+	GRID_TRACE("traceMpDbc");
+	tmp = trace(MpInvJx * Dbc_opt);
+      }
+#if 0
+      // Tesing trace_prod
+      if (a==0){
+	GRID_TRACE("traceMpDbc2");
+	LatticeComplex tmp2(PlaqL.Grid()); tmp2.Checkerboard() = cb;
+	SU3::trace_product(tmp2,MpInvJx,Dbc_opt);
+	std::cout << GridLogMessage <<  "DEBUG: Compute_MpInvJx_dNxxdSy trace " << norm2(tmp2-tmp)<<std::endl;
+      }
+	
+#endif	
+      {
+	GRID_TRACE("pokeIndecx");
+	PokeIndex<ColourIndex>(Fdet2,tmp,a);
+      }
       tpk+=usecond();
     }
     t+=usecond();
-    std::cout << GridLogPerformance << " Compute_MpInvJx_dNxxdSy " << t/1e3 << " ms  proj "<<tp/1e3<< " ms"
-	      << " ta "<<tta/1e3<<" ms" << " poke "<<tpk/1e3<< " ms"<<std::endl;
+    std::cout<< GridLogMessage <<"DEBUG: Full Compute_MpInvJx_dNxxdSy "<<norm2(PlaqL)<<" "<<norm2(PlaqR)<<" "<<norm2(MpInvJx)<<" "<<norm2(Fdet2)<<std::endl;
+    std::cout << GridLogPerformance << " Compute_MpInvJx_dNxxdSy_old " << t/1e3 << " ms  proj "<<tp/1e3<< " ms"
+	      << " ta "<<tta/1e3<<" ms" << " poke "<<tpk/1e3<< " ms LieAlgebraProject "<<tpl/1e3<<" ms"<<std::endl;
   }
   
-  void ComputeNxy(const GaugeLinkField &PlaqL,const GaugeLinkField &PlaqR,AdjMatrixField &NxAd)
+  void ComputeNxy(int old, const GaugeLinkField &PlaqL,const GaugeLinkField &PlaqR,AdjMatrixField &NxAd)
   {
+    GRID_TRACE("ComputeNxy_old");
     GaugeLinkField Nx(PlaqL.Grid());
     const int Ngen = SU3Adjoint::Dimension;
     Complex ci(0,1);
     ColourMatrix   tb;
     ColourMatrix   tc;
+    RealD t = 0, tta = 0, tp = 0, tgen = 0;
+
+    t-=usecond();
     for(int b=0;b<Ngen;b++) {
+      tgen-=usecond();
       SU3::generator(b, tb);
+      tgen+=usecond();
       tb = 2.0 * ci * tb;
-      Nx = Ta( adj(PlaqL)*tb * PlaqR );
+      tta-=usecond();
+      {
+        GRID_TRACE("UtaU");
+	Nx = Ta( adj(PlaqL)*tb * PlaqR );
+      }
+      tta+=usecond();
+      tp-=usecond();
 #if 1
-      SU3::LieAlgebraProject(NxAd,Nx,b);
+      {
+	GRID_TRACE("LieAlgebraProject");
+	SU3::LieAlgebraProject(NxAd,Nx,b);
+      }
 #else
       for(int c=0;c<Ngen;c++) {
 	SU3::generator(c, tc);
-	auto tmp =closure( -trace(ci*tc*Nx)); 
-	PokeIndex<ColourIndex>(NxAd,tmp,c,b); 
+	auto tmp =closure( -trace(ci*tc*Nx));
+	{
+	  GRID_TRACE("pokeIndecx");
+	  PokeIndex<ColourIndex>(NxAd,tmp,c,b);
+	}
       }
 #endif
+      tp+=usecond();
     }
+    t+=usecond();
+    std::cout << GridLogPerformance << " ComputeNxy_old " << t/1e3 << " ms  proj "<<tp/1e3<< " ms"
+              << " ta "<<tta/1e3<<" ms tgen "<< tgen/1e3 << std::endl;
   }
-  void ApplyMask(GaugeField &U,int smr)
-  {
-    LatticeComplex tmp(U.Grid());
-    GaugeLinkField Umu(U.Grid());
-    for(int mu=0;mu<Nd;mu++){
-      Umu=PeekIndex<LorentzIndex>(U,mu);
-      tmp=PeekIndex<LorentzIndex>(masks[smr],mu);
-      Umu=Umu*tmp;
-      PokeIndex<LorentzIndex>(U, Umu, mu);
-    }
-  }
-public:
-
-  void logDetJacobianForceLevel(const GaugeField &U, GaugeField &force ,int smr)
+  void logDetJacobianForceLevel(int old, const GaugeField &U, GaugeField &force ,int smr)
   {
+    GRID_TRACE("logDetJacobianForceLevel_old");
     GridBase* grid = U.Grid();
     ColourMatrix   tb;
     ColourMatrix   tc;
@@ -307,10 +1274,12 @@ class SmearedConfigurationMasked : public SmearedConfiguration<Gimpl>
     // Assemble Luscher exp diff map J matrix 
     //////////////////////////////////////////////////////////////////
     // Ta so Z lives in Lie algabra
+    {GRID_TRACE("Zx");
     Zx  = Ta(Cmu * adj(Umu[mu]));
     time+=usecond();
-    std::cout << GridLogMessage << "Z took "<<time<< " us"<<std::endl;
-
+    }
+    std::cout << GridLogMessage << "Full: Z took "<<time<< " us"<<std::endl;
+    {GRID_TRACE("ZxAd");
     time=-usecond();
     // Move Z to the Adjoint Rep == make_adjoint_representation
     ZxAd = Zero();
@@ -323,11 +1292,13 @@ class SmearedConfigurationMasked : public SmearedConfiguration<Gimpl>
       ZxAd = ZxAd + cplx * TRb; // is this right? YES - Guido used Anti herm Ta's and with bloody wrong sign.
     }
     time+=usecond();
-    std::cout << GridLogMessage << "ZxAd took "<<time<< " us"<<std::endl;
+    }
+    std::cout << GridLogMessage << "Full: ZxAd took "<<time<< " us"<<std::endl;
 
     //////////////////////////////////////
     // J(x) = 1 + Sum_k=1..N (-Zac)^k/(k+1)!
     //////////////////////////////////////
+    {GRID_TRACE("Jx");
     time=-usecond();
     X=1.0; 
     JxAd = X;
@@ -339,7 +1310,8 @@ class SmearedConfigurationMasked : public SmearedConfiguration<Gimpl>
       JxAd = JxAd + X * kpfac;
     }
     time+=usecond();
-    std::cout << GridLogMessage << "Jx took "<<time<< " us"<<std::endl;
+    }
+    std::cout << GridLogMessage << "Full: Jx took "<<time<< " us "<< norm2(JxAd)<<std::endl;
 
     //////////////////////////////////////
     // dJ(x)/dxe
@@ -355,7 +1327,8 @@ class SmearedConfigurationMasked : public SmearedConfiguration<Gimpl>
     AdjMatrixField t3(grid);
     AdjMatrixField dt3(grid);
     AdjMatrixField aunit(grid);
-
+    
+    {GRID_TRACE("dJdX");
     for(int b=0;b<8;b++){
       SU3Adjoint::generator(b, TRb_s[b]);
       dJdX[b] = TRb_s[b];
@@ -374,6 +1347,7 @@ class SmearedConfigurationMasked : public SmearedConfiguration<Gimpl>
     for(int b=0;b<8;b++){
       dJdX[b] = -dJdX[b];
     }
+    }
 #else
     std::vector<AdjMatrixField>  dJdX;    dJdX.resize(8,grid);
     AdjMatrixField tbXn(grid);
@@ -400,16 +1374,16 @@ class SmearedConfigurationMasked : public SmearedConfiguration<Gimpl>
     }
 #endif  
     time+=usecond();
-    std::cout << GridLogMessage << "dJx took "<<time<< " us"<<std::endl;
+    std::cout << GridLogMessage << "Full: dJx took "<<time<< " us "<<norm2(dJdX[0])<<std::endl;
     /////////////////////////////////////////////////////////////////
     // Mask Umu for this link
     /////////////////////////////////////////////////////////////////
     time=-usecond();
     PlaqL = Ident;
     PlaqR = Utmp*adj(Cmu);
-    ComputeNxy(PlaqL,PlaqR,NxxAd);
+    ComputeNxy(old,PlaqL,PlaqR,NxxAd);
     time+=usecond();
-    std::cout << GridLogMessage << "ComputeNxy took "<<time<< " us"<<std::endl;
+    std::cout << GridLogMessage << "Full: ComputeNxy took "<<time<< " us"<<std::endl;
     
     ////////////////////////////
     // Mab
@@ -423,7 +1397,7 @@ class SmearedConfigurationMasked : public SmearedConfiguration<Gimpl>
     time=-usecond();
     MpAdInv = Inverse(MpAd);
     time+=usecond();
-    std::cout << GridLogMessage << "MpAdInv took "<<time<< " us"<<std::endl;
+    std::cout << GridLogMessage << "Full: MpAdInv took "<<time<< " us "<<norm2(MpAdInv)<<std::endl;
     
     RealD t3a = usecond();
     /////////////////////////////////////////////////////////////////
@@ -442,7 +1416,7 @@ class SmearedConfigurationMasked : public SmearedConfiguration<Gimpl>
     AdjMatrixField MpInvJx_nu(grid);
     MpInvJx = (-1.0)*MpAdInv * JxAd;// rho is on the plaq factor
 
-    Compute_MpInvJx_dNxxdSy(PlaqL,PlaqR,MpInvJx,FdetV);
+    Compute_MpInvJx_dNxxdSy(old,PlaqL,PlaqR,MpInvJx,FdetV);
     Fdet2_mu=FdetV;
     Fdet1_mu=Zero();
     
@@ -458,7 +1432,7 @@ class SmearedConfigurationMasked : public SmearedConfiguration<Gimpl>
     ///////////////////////////////
     auto tmp=PeekIndex<LorentzIndex>(masks[smr],mu);
     dJdXe_nMpInv = dJdXe_nMpInv*tmp;
-    
+    std::cout << GridLogMessage << " DEBUG: Full dJdX_nMpAdInv " <<smr<<" "<<mu<<" "<<" simd "<<AdjMatrix::Nsimd()<<" "<<norm2(dJdXe_nMpInv)<<std::endl; 
     //    dJdXe_nMpInv needs to multiply:
     //       Nxx_mu (site local)                           (1)
     //       Nxy_mu one site forward  in each nu direction (3)
@@ -471,7 +1445,7 @@ class SmearedConfigurationMasked : public SmearedConfiguration<Gimpl>
     GaugeField Fdet2(grid);
     GaugeLinkField Fdet_pol(grid); // one polarisation
 
-    RealD t4 = usecond();
+    RealD t4 = usecond(), tLR = 0, tNxy = 0, tMJx = 0;
     for(int nu=0;nu<Nd;nu++){
 
       if (nu!=mu) {
@@ -480,89 +1454,107 @@ class SmearedConfigurationMasked : public SmearedConfiguration<Gimpl>
 	//    |  |
 	//    x==    // nu polarisation -- clockwise
 
-	time=-usecond();
+	time=-usecond(); tLR -= usecond();
 	PlaqL=Ident;
 
 	PlaqR=(-rho)*Gimpl::CovShiftForward(Umu[nu], nu,
  	       Gimpl::CovShiftForward(Umu[mu], mu,
 	         Gimpl::CovShiftBackward(Umu[nu], nu,
 		   Gimpl::CovShiftIdentityBackward(Utmp, mu))));
-	time+=usecond();
-	std::cout << GridLogMessage << "PlaqLR took "<<time<< " us"<<std::endl;
-
-	time=-usecond();
+	time+=usecond(); tLR += usecond();
+	std::cout << GridLogMessage << "Full: PlaqLR took "<<time<< " us"<<std::endl;
+	std::cout << GridLogMessage << " DEBUG: munuLoopBody PalqR " <<smr<<" "<<mu<<" "<<" simd "<<AdjMatrix::Nsimd()<<" "<<norm2(PlaqR)<<std::endl;
+	
+	time=-usecond(); tNxy -= usecond();
 	dJdXe_nMpInv_y =   dJdXe_nMpInv;
-	ComputeNxy(PlaqL,PlaqR,Nxy);
+	ComputeNxy(old,PlaqL,PlaqR,Nxy);
 	Fdet1_nu = transpose(Nxy)*dJdXe_nMpInv_y;
-	time+=usecond();
-	std::cout << GridLogMessage << "ComputeNxy (occurs 6x) took "<<time<< " us"<<std::endl;
+	time+=usecond(); tNxy += usecond();
+	std::cout << GridLogMessage << "Full: ComputeNxy (occurs 6x) took "<<time<< " us"<<std::endl;
 
-	time=-usecond();
+	time=-usecond(); tMJx -= usecond();
 	PlaqR=(-1.0)*PlaqR;
-	Compute_MpInvJx_dNxxdSy(PlaqL,PlaqR,MpInvJx,FdetV);
+	Compute_MpInvJx_dNxxdSy(old,PlaqL,PlaqR,MpInvJx,FdetV);
 	Fdet2_nu = FdetV;
-	time+=usecond();
-	std::cout << GridLogMessage << "Compute_MpInvJx_dNxxSy (occurs 6x) took "<<time<< " us"<<std::endl;
-	
+	time+=usecond(); tMJx += usecond();
+	std::cout << GridLogMessage << "Full: Compute_MpInvJx_dNxxSy (occurs 6x) took "<<time<< " us"<<std::endl;
+	std::cout << GridLogMessage << " DEBUG: Full munuLoopBody FdetV " <<smr<<" "<<mu<<" simd "<<AdjMatrix::Nsimd()<<" "<<norm2(FdetV)<<std::endl;
 	//    x==
 	//    |  |
 	//    .__|    // nu polarisation -- anticlockwise
 
+	tLR -= usecond();
 	PlaqR=(rho)*Gimpl::CovShiftForward(Umu[nu], nu,
 		      Gimpl::CovShiftBackward(Umu[mu], mu,
     	 	        Gimpl::CovShiftIdentityBackward(Umu[nu], nu)));
 
 	PlaqL=Gimpl::CovShiftIdentityBackward(Utmp, mu);
+	tLR += usecond();
 
+	tNxy -= usecond();
 	dJdXe_nMpInv_y = Cshift(dJdXe_nMpInv,mu,-1);
-	ComputeNxy(PlaqL, PlaqR,Nxy);
+	ComputeNxy(old,PlaqL, PlaqR,Nxy);
 	Fdet1_nu = Fdet1_nu+transpose(Nxy)*dJdXe_nMpInv_y;
-	
+	tNxy += usecond();
 
+	tMJx -= usecond();
 	MpInvJx_nu = Cshift(MpInvJx,mu,-1);
-	Compute_MpInvJx_dNxxdSy(PlaqL,PlaqR,MpInvJx_nu,FdetV);
+	Compute_MpInvJx_dNxxdSy(old,PlaqL,PlaqR,MpInvJx_nu,FdetV);
 	Fdet2_nu = Fdet2_nu+FdetV;
+	tMJx += usecond();
 	
 	///////////////// -ve nu /////////////////
 	//  __
 	// |  |
 	// x==          // nu polarisation -- clockwise
 
+	tLR -= usecond();
 	PlaqL=(rho)* Gimpl::CovShiftForward(Umu[mu], mu,
 		       Gimpl::CovShiftForward(Umu[nu], nu,
 			 Gimpl::CovShiftIdentityBackward(Utmp, mu)));
 
         PlaqR = Gimpl::CovShiftIdentityForward(Umu[nu], nu);
+	tLR += usecond();
 
+	tNxy -= usecond();
 	dJdXe_nMpInv_y = Cshift(dJdXe_nMpInv,nu,1);
-	ComputeNxy(PlaqL,PlaqR,Nxy);
+	ComputeNxy(old,PlaqL,PlaqR,Nxy);
 	Fdet1_nu = Fdet1_nu + transpose(Nxy)*dJdXe_nMpInv_y;
+	tNxy += usecond();
 
+	tMJx -= usecond();
 	MpInvJx_nu = Cshift(MpInvJx,nu,1);
-	Compute_MpInvJx_dNxxdSy(PlaqL,PlaqR,MpInvJx_nu,FdetV);
+	Compute_MpInvJx_dNxxdSy(old,PlaqL,PlaqR,MpInvJx_nu,FdetV);
 	Fdet2_nu = Fdet2_nu+FdetV;
+	tMJx += usecond();
 	
 	// x==
 	// |  |
 	// |__|         // nu polarisation
 
+	tLR -= usecond();
 	PlaqL=(-rho)*Gimpl::CovShiftForward(Umu[nu], nu,
  	        Gimpl::CovShiftIdentityBackward(Utmp, mu));
 
 	PlaqR=Gimpl::CovShiftBackward(Umu[mu], mu,
 	        Gimpl::CovShiftIdentityForward(Umu[nu], nu));
+	tLR += usecond();
 
+	tNxy -= usecond();
 	dJdXe_nMpInv_y = Cshift(dJdXe_nMpInv,mu,-1);
 	dJdXe_nMpInv_y = Cshift(dJdXe_nMpInv_y,nu,1);
 
-	ComputeNxy(PlaqL,PlaqR,Nxy);
+	ComputeNxy(old,PlaqL,PlaqR,Nxy);
 	Fdet1_nu = Fdet1_nu + transpose(Nxy)*dJdXe_nMpInv_y;
+	tNxy += usecond();
 
+	tMJx -= usecond();
 	MpInvJx_nu = Cshift(MpInvJx,mu,-1);
 	MpInvJx_nu = Cshift(MpInvJx_nu,nu,1);
-	Compute_MpInvJx_dNxxdSy(PlaqL,PlaqR,MpInvJx_nu,FdetV);
+	Compute_MpInvJx_dNxxdSy(old,PlaqL,PlaqR,MpInvJx_nu,FdetV);
 	Fdet2_nu = Fdet2_nu+FdetV;
-
+	tMJx += usecond();
+	
 	/////////////////////////////////////////////////////////////////////
 	// Set up the determinant force contribution in 3x3 algebra basis
 	/////////////////////////////////////////////////////////////////////
@@ -576,42 +1568,52 @@ class SmearedConfigurationMasked : public SmearedConfiguration<Gimpl>
         //     __
 	//    |  "
 	//    |__"x    // mu polarisation
+	tLR -= usecond();
 	PlaqL=(-rho)*Gimpl::CovShiftForward(Umu[mu], mu,
 		      Gimpl::CovShiftBackward(Umu[nu], nu,
    		        Gimpl::CovShiftIdentityBackward(Utmp, mu)));
 
 	PlaqR=Gimpl::CovShiftIdentityBackward(Umu[nu], nu);
-	
+	tLR += usecond();
+
+	tNxy -= usecond();
 	dJdXe_nMpInv_y = Cshift(dJdXe_nMpInv,nu,-1);
 
-	ComputeNxy(PlaqL,PlaqR,Nxy);
+	ComputeNxy(old,PlaqL,PlaqR,Nxy);
 	Fdet1_mu = Fdet1_mu + transpose(Nxy)*dJdXe_nMpInv_y;
+	tNxy += usecond();
 
+	tMJx -= usecond();
 	MpInvJx_nu = Cshift(MpInvJx,nu,-1);
 
-	Compute_MpInvJx_dNxxdSy(PlaqL,PlaqR,MpInvJx_nu,FdetV);
+	Compute_MpInvJx_dNxxdSy(old,PlaqL,PlaqR,MpInvJx_nu,FdetV);
 	Fdet2_mu = Fdet2_mu+FdetV;
+	tMJx += usecond();
 
 	//  __
 	// "  |
 	// x__|          // mu polarisation
-
+	tLR -= usecond();
 	PlaqL=(-rho)*Gimpl::CovShiftForward(Umu[mu], mu,
 		       Gimpl::CovShiftForward(Umu[nu], nu,
 		 	 Gimpl::CovShiftIdentityBackward(Utmp, mu)));
 
         PlaqR=Gimpl::CovShiftIdentityForward(Umu[nu], nu);
+	tLR += usecond();
 
+	tNxy -= usecond();
 	dJdXe_nMpInv_y = Cshift(dJdXe_nMpInv,nu,1);
 
-	ComputeNxy(PlaqL,PlaqR,Nxy);
+	ComputeNxy(old,PlaqL,PlaqR,Nxy);
 	Fdet1_mu = Fdet1_mu + transpose(Nxy)*dJdXe_nMpInv_y;
+	tNxy += usecond();
 
+	tMJx -= usecond();
 	MpInvJx_nu = Cshift(MpInvJx,nu,1);
 
-	Compute_MpInvJx_dNxxdSy(PlaqL,PlaqR,MpInvJx_nu,FdetV);
+	Compute_MpInvJx_dNxxdSy(old,PlaqL,PlaqR,MpInvJx_nu,FdetV);
 	Fdet2_mu = Fdet2_mu+FdetV;
-	
+	tMJx += usecond();
       }
     }
     RealD t5 = usecond();
@@ -623,14 +1625,14 @@ class SmearedConfigurationMasked : public SmearedConfiguration<Gimpl>
 
     force= (-0.5)*( Fdet1 + Fdet2);
     RealD t1 = usecond();
-    std::cout << GridLogMessage << " logDetJacobianForce level took "<<t1-t0<<" us "<<std::endl;
-    std::cout << GridLogMessage << " logDetJacobianForce t3-t0 "<<t3a-t0<<" us "<<std::endl;
-    std::cout << GridLogMessage << " logDetJacobianForce t4-t3 dJdXe_nMpInv "<<t4-t3a<<" us "<<std::endl;
-    std::cout << GridLogMessage << " logDetJacobianForce t5-t4 mu nu loop "<<t5-t4<<" us "<<std::endl;
-    std::cout << GridLogMessage << " logDetJacobianForce t1-t5 "<<t1-t5<<" us "<<std::endl;
-    std::cout << GridLogMessage << " logDetJacobianForce level took "<<t1-t0<<" us "<<std::endl;
+    std::cout << GridLogMessage << " Full: logDetJacobianForce t3-t0 "<<t3a-t0<<" us "<<std::endl;
+    std::cout << GridLogMessage << " Full: logDetJacobianForce t4-t3 dJdXe_nMpInv "<<t4-t3a<<" us "<<std::endl;
+    std::cout << GridLogMessage << " Full: logDetJacobianForce t5-t4 mu nu loop "<<t5-t4<<" us Plaq "
+	      <<tLR/1e3<<" ms Nxy "<<tNxy/1e3<<" ms MpInvJx_dNxxdSy "<<tMJx/1e3<<" ms"<<std::endl;
+    std::cout << GridLogMessage << " Full: logDetJacobianForce t1-t5 "<<t1-t5<<" us "<<std::endl; // turn adj vec to SU3 force
+    std::cout << GridLogMessage << " Full: logDetJacobianForce level took "<<t1-t0<<" us "<<std::endl;
   }
-  RealD logDetJacobianLevel(const GaugeField &U,int smr)
+  RealD logDetJacobianLevel(int old, const GaugeField &U,int smr)
   {
     GridBase* grid = U.Grid();
     GaugeField C(grid);
@@ -647,43 +1649,78 @@ class SmearedConfigurationMasked : public SmearedConfiguration<Gimpl>
     LatticeComplex       cplx(grid); 
     AdjVectorField  AlgV(grid); 
     AdjMatrixField  Mab(grid);
-    AdjMatrixField  Ncb(grid);
+    AdjMatrixField  Ncb(grid), Ncb_opt(grid);
     AdjMatrixField  Jac(grid);
     AdjMatrixField  Zac(grid);
     AdjMatrixField  mZac(grid);
     AdjMatrixField  X(grid);
 
+    RealD time=0, tta=0, tpk=0, tN=0, tZ=0, tJ=0, tlnDetM=0;
     int mu= (smr/2) %Nd;
 
+    time -= usecond();
     auto mask=PeekIndex<LorentzIndex>(masks[smr],mu); // the cb mask
 
     //////////////////////////////////////////////////////////////////
     // Assemble the N matrix
     //////////////////////////////////////////////////////////////////
+
+    tN -= usecond();
     double rho=this->StoutSmearing->SmearRho[1];
     BaseSmear(Cmu, U,mu,rho);
 
     Umu = peekLorentz(U, mu);
     Complex ci(0,1);
+#if 1
     for(int b=0;b<Ngen;b++) {
       SU3::generator(b, Tb);
       // Qlat Tb = 2i Tb^Grid
+      tta -= usecond();
       Nb = (2.0)*Ta( ci*Tb * Umu * adj(Cmu));
+      tta += usecond();
       // FIXME -- replace this with LieAlgebraProject
 #if 0
-      SU3::LieAlgebraProject(Ncb,tmp,b);
+      // Fixed it
+      SU3::LieAlgebraProject(Ncb_opt,Nb,b);
 #else
       for(int c=0;c<Ngen;c++) {
 	SU3::generator(c, Tc);
 	auto tmp = -trace(ci*Tc*Nb); // Luchang's norm: (2Tc) (2Td) N^db = -2 delta cd N^db // - was important
-	PokeIndex<ColourIndex>(Ncb,tmp,c,b); 
+	tpk -= usecond();
+	PokeIndex<ColourIndex>(Ncb,tmp,c,b);
+	tpk += usecond();
       }
-#endif
-    }      
 
+#endif
+    }
+#else
+      autoView(NxAd_v,NxAd,AcceleratorWrite);
+    autoView(PlaqL_v,PlaqL,AcceleratorRead);
+    autoView(PlaqR_v,PlaqR,AcceleratorRead);
+    const int nsimd = vAlgebraMatrix::Nsimd();
+    accelerator_for(ss,grid->oSites(),nsimd,{
+        typedef decltype(coalescedRead(PlaqL_v[0])) SU3_mat;
+        typedef decltype(coalescedRead(NxAd_v[0]))  adj_mat;
+        SU3_mat Nx;
+        adj_mat NxAd_site;
+        for(int b=0;b<Ngen;b++) {
+          SU3::generator(b, tb);
+          tb = 2.0 * ci * tb;
+          //auto Nx =Ta( adj(PlaqL_v(ss)) * tb * PlaqR_v(ss) );
+          auto Nx = 0.5*( adj(PlaqL_v(ss))*tb*PlaqR_v(ss) + adj(PlaqR_v(ss))*tb*PlaqL_v(ss) );//-tr(*) part does not contribute
+          SU3::LieAlgebraProject(NxAd_site,Nx,b);
+        }
+        coalescedWrite(NxAd_v[ss],NxAd_site);
+      });
+#endif
+    //Dump(Ncb_opt,"Ncb_opt");
+    //Dump(Ncb,"Ncb");
+    tN += usecond();
+    
     //////////////////////////////////////////////////////////////////
     // Assemble Luscher exp diff map J matrix 
     //////////////////////////////////////////////////////////////////
+    tZ -= usecond();
     // Ta so Z lives in Lie algabra
     Z  = Ta(Cmu * adj(Umu));
 
@@ -699,10 +1736,12 @@ class SmearedConfigurationMasked : public SmearedConfiguration<Gimpl>
       cplx = 2.0*trace(ci*Tb*Z); // my convention 1/2 delta ba
       Zac = Zac + cplx * TRb; // is this right? YES - Guido used Anti herm Ta's and with bloody wrong sign.
     }
-
+    tZ += usecond();
+    
     //////////////////////////////////////
     // J(x) = 1 + Sum_k=1..N (-Zac)^k/(k+1)!
     //////////////////////////////////////
+    tJ -= usecond();
     X=1.0; 
     Jac = X;
     mZac = (-1.0)*Zac; 
@@ -712,7 +1751,9 @@ class SmearedConfigurationMasked : public SmearedConfiguration<Gimpl>
       kpfac = kpfac /(k+1);
       Jac = Jac + X * kpfac;
     }
+    tJ += usecond();
 
+    tlnDetM -= usecond();
     ////////////////////////////
     // Mab
     ////////////////////////////
@@ -735,10 +1776,102 @@ class SmearedConfigurationMasked : public SmearedConfiguration<Gimpl>
     // Masked sum
     ////////////////////////////
     ln_det = ln_det * mask;
+    tlnDetM += usecond();
+    time += usecond();
+    std::cout << GridLogMessage << " Full: logDetJacobianLevel " << time/1e3 << " ms ta "<<tta/1e3<<" ms" << " poke "<<tpk/1e3<< " ms"
+	      <<" N "<<tN/1e3<<" ms Z "<<tZ/1e3<<" ms J " <<tJ/1e3<<" ms lnDetM "<<tlnDetM/1e3<<" ms" <<std::endl;
+    
     Complex result = sum(ln_det);
     return result.real();
   }
+
+  
 public:
+  RealD logDetJacobian(int old)
+  {
+    RealD ln_det = 0;
+    if (this->smearingLevels > 0)
+    {
+      double start = usecond();
+      for (int ismr = this->smearingLevels - 1; ismr > 0; --ismr) {
+	ln_det+= logDetJacobianLevel(old,this->get_smeared_conf(ismr-1),ismr);
+      }
+      ln_det +=logDetJacobianLevel(old,*(this->ThinLinks),0);
+
+      double end = usecond();
+      double time = (end - start)/ 1e3;
+      std::cout << GridLogMessage << "Full: GaugeConfigurationMasked: logDetJacobian took " << time << " ms" << std::endl;  
+    }
+    std::cout << GridLogMessage << " DEBUG: logDetJacobian Full " << std::endl;
+    return ln_det;
+  }
+  void logDetJacobianForce(int old, GaugeField &force)
+  {
+    GRID_TRACE("logDetJacobianForce");
+    force =Zero();
+    GaugeField force_det(force.Grid());
+
+    if (this->smearingLevels > 0)
+    {
+      double start = usecond(), tas = 0;
+
+      GaugeLinkField tmp_mu(force.Grid());
+
+      for (int ismr = this->smearingLevels - 1; ismr > 0; --ismr) {
+
+	// remove U in UdSdU...
+	for (int mu = 0; mu < Nd; mu++) {
+	  tmp_mu = adj(peekLorentz(this->get_smeared_conf(ismr), mu)) * peekLorentz(force, mu);
+	  pokeLorentz(force, tmp_mu, mu);
+	}
+
+	tas -= usecond();
+      	// Propagate existing force
+        force = this->AnalyticSmearedForce(force, this->get_smeared_conf(ismr - 1), ismr);
+	tas += usecond();
+	
+	// Add back U in UdSdU...
+	for (int mu = 0; mu < Nd; mu++) {
+	  tmp_mu = peekLorentz(this->get_smeared_conf(ismr - 1), mu) * peekLorentz(force, mu);
+	  pokeLorentz(force, tmp_mu, mu);
+	}
+    	
+	// Get this levels determinant force
+	force_det = Zero();
+	logDetJacobianForceLevel(old, this->get_smeared_conf(ismr-1),force_det,ismr);
+
+	// Sum the contributions
+	force = force + force_det;
+      }
+    
+      // remove U in UdSdU...
+      for (int mu = 0; mu < Nd; mu++) {
+	tmp_mu = adj(peekLorentz(this->get_smeared_conf(0), mu)) * peekLorentz(force, mu);
+	pokeLorentz(force, tmp_mu, mu);
+      }
+
+      tas -= usecond();
+      force = this->AnalyticSmearedForce(force, *this->ThinLinks,0);
+      tas += usecond();
+      
+      for (int mu = 0; mu < Nd; mu++) {
+	tmp_mu = peekLorentz(*this->ThinLinks, mu) * peekLorentz(force, mu);
+	pokeLorentz(force, tmp_mu, mu);
+      }
+
+      force_det = Zero();
+
+      logDetJacobianForceLevel(old, *this->ThinLinks,force_det,0);
+
+      force = force + force_det;
+
+      force=Ta(force); // Ta
+      
+    }  // if smearingLevels = 0 do nothing
+    std::cout << GridLogMessage << " DEBUG: logDetJacobianForce Full " << std::endl;
+  }
+  /*------------------------------- OLD IMPLEMENTATION ----------------------------------------------------------*/
+  
   RealD logDetJacobian(void)
   {
     RealD ln_det = 0;
@@ -752,18 +1885,23 @@ class SmearedConfigurationMasked : public SmearedConfiguration<Gimpl>
 
       double end = usecond();
       double time = (end - start)/ 1e3;
-      std::cout << GridLogMessage << "GaugeConfigurationMasked: logDetJacobian took " << time << " ms" << std::endl;  
+      std::cout << GridLogMessage << "GaugeConfigurationMasked: logDetJacobian took " << time << " ms" << std::endl;
+#ifdef DEBUG
+      RealD ln_det2 = logDetJacobian(1);
+      std::cout << GridLogMessage << " DEBUG: logDetJacobian_diff " << abs(ln_det2-ln_det) <<":"<<ln_det<<" "<<ln_det2<<std::endl;
+#endif
     }
     return ln_det;
   }
   void logDetJacobianForce(GaugeField &force)
   {
-    force =Zero();
+    GRID_TRACE("logDetJacobianForce");
+    force = Zero();
     GaugeField force_det(force.Grid());
 
     if (this->smearingLevels > 0)
     {
-      double start = usecond();
+      double start = usecond(), tas = 0;
 
       GaugeLinkField tmp_mu(force.Grid());
 
@@ -774,10 +1912,12 @@ class SmearedConfigurationMasked : public SmearedConfiguration<Gimpl>
 	  tmp_mu = adj(peekLorentz(this->get_smeared_conf(ismr), mu)) * peekLorentz(force, mu);
 	  pokeLorentz(force, tmp_mu, mu);
 	}
-	
+
+	tas -= usecond();
       	// Propagate existing force
         force = this->AnalyticSmearedForce(force, this->get_smeared_conf(ismr - 1), ismr);
-
+	tas += usecond();
+	
 	// Add back U in UdSdU...
 	for (int mu = 0; mu < Nd; mu++) {
 	  tmp_mu = peekLorentz(this->get_smeared_conf(ismr - 1), mu) * peekLorentz(force, mu);
@@ -798,8 +1938,10 @@ class SmearedConfigurationMasked : public SmearedConfiguration<Gimpl>
 	pokeLorentz(force, tmp_mu, mu);
       }
 
+      tas -= usecond();
       force = this->AnalyticSmearedForce(force, *this->ThinLinks,0);
-
+      tas += usecond();
+      
       for (int mu = 0; mu < Nd; mu++) {
 	tmp_mu = peekLorentz(*this->ThinLinks, mu) * peekLorentz(force, mu);
 	pokeLorentz(force, tmp_mu, mu);
@@ -813,8 +1955,14 @@ class SmearedConfigurationMasked : public SmearedConfiguration<Gimpl>
 
       force=Ta(force); // Ta
       
+#ifdef DEBUG // debug
+      GaugeField force_debug(force.Grid()); 
+      logDetJacobianForce(1,force_debug);
+      std::cout << GridLogMessage << " DEBUG: logDetJacobianForce_diff " << norm2(force-force_debug) <<":"<<norm2(force)<<" "<<norm2(force_debug)<< std::endl;
+#endif
       double end = usecond();
       double time = (end - start)/ 1e3;
+      std::cout << GridLogMessage << "GaugeConfigurationMasked: AnalyticSmearedForce in lnDetJacobianForce took " << tas/1e3 << " ms" << std::endl;
       std::cout << GridLogMessage << "GaugeConfigurationMasked: lnDetJacobianForce took " << time << " ms" << std::endl;  
     }  // if smearingLevels = 0 do nothing
   }
@@ -856,10 +2004,20 @@ class SmearedConfigurationMasked : public SmearedConfiguration<Gimpl>
         // For debug purposes
         RealD impl_plaq = WilsonLoops<Gimpl>::avgPlaquette(previous_u);
         std::cout << GridLogMessage << "[SmearedConfigurationMasked] smeared Plaq: " << impl_plaq << std::endl;
+#ifndef PRINT_SNAPSHOTS
+      }
+#else
+      // overhead not small
+      //RealD impl_5Li = WilsonLoops<PeriodicGimplR>::TopologicalCharge5Li(previous_u);
+      //std::cout << GridLogMessage << "[SmearedConfigurationMasked] smeared 5Li at " << smearLvl << " : " << impl_5Li << std::endl;
       }
+      RealD intmdt_lnDetJ = logDetJacobian();
+      std::cout << GridLogMessage << "GaugeConfigurationMasked: Intermediate Jacobian " << intmdt_lnDetJ << std::endl;
+#endif
+
       double end = usecond();
       double time = (end - start)/ 1e3;
-      std::cout << GridLogMessage << "GaugeConfigurationMasked: Link smearing took " << time << " ms" << std::endl;  
+      std::cout << GridLogMessage << "GaugeConfigurationMasked: Link smearing took " << time << " ms" << std::endl;
     }
   }
   //====================================================================
@@ -867,6 +2025,7 @@ class SmearedConfigurationMasked : public SmearedConfiguration<Gimpl>
   virtual GaugeField AnalyticSmearedForce(const GaugeField& SigmaKPrime,
 					  const GaugeField& GaugeK,int level) 
   {
+    GRID_TRACE("AnalyticSmearedForce");
     GridBase* grid = GaugeK.Grid();
     GaugeField SigmaK(grid), iLambda(grid);
     GaugeField SigmaKPrimeA(grid);
@@ -880,6 +2039,9 @@ class SmearedConfigurationMasked : public SmearedConfiguration<Gimpl>
     int cb= (level%2);
     double rho=this->StoutSmearing->SmearRho[1];
 
+    RealD time = 0;
+
+    time -= usecond();
     // Can override this to do one direction only.
     SigmaK = Zero();
     iLambda = Zero();
@@ -914,7 +2076,7 @@ class SmearedConfigurationMasked : public SmearedConfiguration<Gimpl>
       this->set_iLambda(iLambda_mu, e_iQ, iQ, SigmaKPrime_mu, GaugeKmu);
       pokeLorentz(SigmaK, SigmaKPrime_mu * e_iQ + adj(Cmu) * iLambda_mu, mu);
       pokeLorentz(iLambda, iLambda_mu, mu);
-      std::cout << " mu "<<mu<<" SigmaKPrime_mu"<<norm2(SigmaKPrime_mu)<< " iLambda_mu " <<norm2(iLambda_mu)<<std::endl;
+      std::cout << " mu "<<mu<<" SigmaKPrime_mu "<<norm2(SigmaKPrime_mu)<< " iLambda_mu " <<norm2(iLambda_mu)<<std::endl;
     }
     //    GaugeField SigmaKcopy(grid);
     //    SigmaKcopy = SigmaK;
@@ -927,26 +2089,35 @@ class SmearedConfigurationMasked : public SmearedConfiguration<Gimpl>
     // propagate the rest of the force as identity map, just add back
     ////////////////////////////////////////////////////////////////////////////////////
     SigmaK = SigmaK+SigmaKPrimeB;
+    time += usecond();
 
+    std::cout << GridLogMessage << "GaugeConfigurationMasked: Analytic smearing force  took " << time/1e3 << " ms" << std::endl;
     return SigmaK;
   }
 
 public:
 
   /* Standard constructor */
+  virtual ~SmearedConfigurationMasked()
+  {
+    delete UrbGrid;
+    gStencils.clear();
+    gStencils_smear.clear();
+  }
   SmearedConfigurationMasked(GridCartesian* _UGrid, unsigned int Nsmear, Smear_Stout<Gimpl>& Stout)
-    : SmearedConfiguration<Gimpl>(_UGrid, Nsmear,Stout)
+    : SmearedConfiguration<Gimpl>(_UGrid, Nsmear,Stout),
+      UGrid(_UGrid),
+      Ghost(1,_UGrid)
   {
     assert(Nsmear%(2*Nd)==0); // Or multiply by 8??
 
     // was resized in base class
     assert(this->SmearedSet.size()==Nsmear);
     
-    GridRedBlackCartesian * UrbGrid;
     UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(_UGrid);
+    
     LatticeComplex one(_UGrid); one = ComplexD(1.0,0.0);
     LatticeComplex tmp(_UGrid);
-
     for (unsigned int i = 0; i < this->smearingLevels; ++i) {
 
       masks.push_back(*(new LatticeLorentzComplex(_UGrid)));
@@ -954,7 +2125,9 @@ class SmearedConfigurationMasked : public SmearedConfiguration<Gimpl>
       int mu= (i/2) %Nd;
       int cb= (i%2);
       LatticeComplex tmpcb(UrbGrid);
-	
+
+      cbs.push_back(cb);
+      
       masks[i]=Zero();
       ////////////////////
       // Setup the mask
@@ -965,11 +2138,82 @@ class SmearedConfigurationMasked : public SmearedConfiguration<Gimpl>
       PokeIndex<LorentzIndex>(masks[i],tmp, mu);
 	
     }
-    delete UrbGrid;
+    ///////////////////////////////
+    // Setup stencils for staples
+    ///////////////////////////////
+    int depth            = 1;
+    Coordinate local     =UGrid->LocalDimensions();
+    Coordinate simd      =UGrid->_simd_layout;
+    Coordinate processors=UGrid->_processors;
+    Coordinate plocal    =UGrid->LocalDimensions();
+    Coordinate global(Nd);
+
+    // pre-stencil calculation for force calculation
+    GridBase *ggrid = Ghost.grids[Nd-1];
+    std::vector<Coordinate> shifts;
+
+    gStencils.clear();
+    Coordinate shift_0(Nd,0);
+    for(int mu=0;mu<Nd;mu++){
+      Coordinate shift_mu(Nd,0);  shift_mu[mu]=1;
+      Coordinate shift_mmu(Nd,0); shift_mmu[mu]=-1;
+      for(int nu=0;nu<Nd;nu++){
+	if (nu!=mu) {
+	  Coordinate shift_nu(Nd,0);  shift_nu[nu]=1;
+	  Coordinate shift_mnu(Nd,0); shift_mnu[nu]=-1;
+	  Coordinate shift_pmu_pnu(Nd,0); shift_pmu_pnu[mu]= 1;  shift_pmu_pnu[nu]= 1;
+	  Coordinate shift_pmu_mnu(Nd,0); shift_pmu_mnu[mu]= 1;  shift_pmu_mnu[nu]=-1;
+	  Coordinate shift_mmu_pnu(Nd,0); shift_mmu_pnu[mu]=-1;  shift_mmu_pnu[nu]= 1;
+	  shifts.clear();
+	  shifts.push_back(shift_0);shifts.push_back(shift_nu);shifts.push_back(shift_mu);shifts.push_back(shift_0);
+	  gStencils.push_back(GeneralLocalStencil(ggrid,shifts));
+
+	  shifts.clear();
+	  shifts.push_back(shift_0);shifts.push_back(shift_mmu_pnu);shifts.push_back(shift_mmu);shifts.push_back(shift_mmu);
+	  gStencils.push_back(GeneralLocalStencil(ggrid,shifts));
+
+	  shifts.clear();
+	  shifts.push_back(shift_0);shifts.push_back(shift_mu);shifts.push_back(shift_nu);shifts.push_back(shift_0);
+	  gStencils.push_back(GeneralLocalStencil(ggrid,shifts));
+	  
+	  shifts.clear();
+	  shifts.push_back(shift_0);shifts.push_back(shift_mmu_pnu);shifts.push_back(shift_mmu);shifts.push_back(shift_mmu);
+	  gStencils.push_back(GeneralLocalStencil(ggrid,shifts));
+	  
+	  shifts.clear();
+	  shifts.push_back(shift_0);shifts.push_back(shift_pmu_mnu);shifts.push_back(shift_mnu);shifts.push_back(shift_mnu);
+	  gStencils.push_back(GeneralLocalStencil(ggrid,shifts));
+	  
+	  shifts.clear();
+	  shifts.push_back(shift_0);shifts.push_back(shift_mu);shifts.push_back(shift_nu);shifts.push_back(shift_0);
+	  gStencils.push_back(GeneralLocalStencil(ggrid,shifts));
+	}
+      }
+    }
+    // pre-stencil calculation for BaseSmear
+    gStencils_smear.clear();
+    for(int mu=0;mu<Nd;mu++){
+      Coordinate shift_mu(Nd,0);  shift_mu[mu]=1;
+      shifts.clear();
+      for(int nu=0;nu<Nd;nu++){
+        if (nu!=mu) {
+	  Coordinate shift_nu(Nd,0);  shift_nu[nu]=1;
+          Coordinate shift_mnu(Nd,0); shift_mnu[nu]=-1;
+	  Coordinate shift_pmu_mnu(Nd,0); shift_pmu_mnu[mu]=1; shift_pmu_mnu[nu]=-1;
+
+	  //U_nu(x+mu)U^dag_mu(x+nu) U^dag_nu(x)
+	  shifts.push_back(shift_0); shifts.push_back(shift_nu); shifts.push_back(shift_mu);
+	  //U_nu^dag(x-nu+mu) U_mu^dag(x-nu) U_nu(x-nu)
+	  shifts.push_back(shift_mnu); shifts.push_back(shift_mnu); shifts.push_back(shift_pmu_mnu);
+	}
+      }
+      gStencils_smear.push_back(GeneralLocalStencil(ggrid,shifts));
+    }
   }
   
   virtual void smeared_force(GaugeField &SigmaTilde) 
   {
+    GRID_TRACE("smeared_force");
     if (this->smearingLevels > 0)
     {
       double start = usecond();
@@ -1009,4 +2253,4 @@ class SmearedConfigurationMasked : public SmearedConfiguration<Gimpl>
 };
 
 NAMESPACE_END(Grid);
-
+#endif
diff --git a/Grid/qcd/smearing/JacobianAction.h b/Grid/qcd/smearing/JacobianAction.h
index bc691581ac..a16a87ec61 100644
--- a/Grid/qcd/smearing/JacobianAction.h
+++ b/Grid/qcd/smearing/JacobianAction.h
@@ -78,10 +78,12 @@ class JacobianAction : public Action<typename Gimpl::GaugeField> {
   {
     GRID_ASSERT( &U == smearer );
     smearer->logDetJacobianForce(dSdU);
+#ifdef PRINT_SNAPSHOTS
+    writeField(dSdU, "F_"+action_name()+"_smr."+std::to_string(this->deriv_num));
+#endif
   }
 
 private:
- };
+};
 
 NAMESPACE_END(Grid);
-
diff --git a/Grid/qcd/utils/GaugeFix.h b/Grid/qcd/utils/GaugeFix.h
index 95c4dd481e..170c5c56b1 100644
--- a/Grid/qcd/utils/GaugeFix.h
+++ b/Grid/qcd/utils/GaugeFix.h
@@ -24,7 +24,6 @@ Author: Peter Boyle <paboyle@ph.ed.ac.uk>
     See the full license in the file "LICENSE" in the top level distribution directory
 *************************************************************************************/
 /*  END LEGAL */
-//#include <Grid/Grid.h>
 
 #ifndef GRID_QCD_GAUGE_FIX_H
 #define GRID_QCD_GAUGE_FIX_H
diff --git a/Grid/qcd/utils/GaugeGroup.h b/Grid/qcd/utils/GaugeGroup.h
index 6811d247e9..f0f9ef1f06 100644
--- a/Grid/qcd/utils/GaugeGroup.h
+++ b/Grid/qcd/utils/GaugeGroup.h
@@ -46,7 +46,19 @@ directory
                           std::is_same<dummy_name, group_name>::value && \
                       is_sp<dummy_name>::value >
 
+
 NAMESPACE_BEGIN(Grid);
+
+namespace CplxMultExt {
+  template< typename T, typename SCALAR > inline
+  typename std::enable_if< std::is_integral<SCALAR>::value, std::complex<T> >::type
+  operator* ( const std::complex<T>& c, SCALAR n ) { return c * static_cast<T>(n) ; }
+
+  template< typename T, typename SCALAR > inline
+  typename std::enable_if< std::is_integral<SCALAR>::value, std::complex<T> >::type
+  operator* ( SCALAR n, const std::complex<T>& c ) { return static_cast<T>(n) * c ; }
+}
+
 namespace GroupName {
 class SU {};
 class Sp {};
@@ -135,7 +147,11 @@ class GaugeGroup {
   typedef Lattice<vAlgebraVector> LatticeAlgebraVector;
   typedef Lattice<vAlgebraVectorF> LatticeAlgebraVectorF;
   typedef Lattice<vAlgebraVectorD> LatticeAlgebraVectorD;
-   
+
+  typedef iSUnAlgebraMatrix<Complex>  AlgebraMatrix;
+  typedef iSUnAlgebraMatrix<ComplexF> AlgebraMatrixF;
+  typedef iSUnAlgebraMatrix<ComplexD> AlgebraMatrixD;
+  
   typedef iSUnAlgebraMatrix<vComplex>  vAlgebraMatrix;
   typedef iSUnAlgebraMatrix<vComplexF> vAlgebraMatrixF;
   typedef iSUnAlgebraMatrix<vComplexD> vAlgebraMatrixD;
@@ -168,7 +184,7 @@ class GaugeGroup {
 
  public:
   template <class cplx>
-  static void generator(int lieIndex, iGroupMatrix<cplx> &ta) {
+  static accelerator_inline void generator(int lieIndex, iGroupMatrix<cplx> &ta) {
     return generator(lieIndex, ta, group_name());
   }
 
@@ -401,52 +417,355 @@ class GaugeGroup {
     }
   }
 
-// Ta are hermitian (?)
-// Anti herm is i Ta basis
-static void LieAlgebraProject(LatticeAlgebraMatrix &out,const LatticeMatrix &in, int b)
-{
-  conformable(in, out);
-  GridBase *grid = out.Grid();
-  LatticeComplex tmp(grid);
-  Matrix ta;
-  // Using Luchang's projection convention
-  //  2 Tr{Ta Tb} A_b= 2/2 delta ab A_b = A_a
-  autoView(out_v,out,AcceleratorWrite);
-  autoView(in_v,in,AcceleratorRead);
-  int N = ncolour;
-  int NNm1 = N * (N - 1);
-  int hNNm1= NNm1/2;
-  RealD sqrt_2 = sqrt(2.0);
-  Complex ci(0.0,1.0);
-  for(int su2Index=0;su2Index<hNNm1;su2Index++){
-    int i1, i2;
-    su2SubGroupIndex(i1, i2, su2Index);
-    int ax = su2Index*2;
-    int ay = su2Index*2+1;
-    accelerator_for(ss,grid->oSites(),1,{
-	// in is traceless ANTI-hermitian whereas Grid generators are Hermitian.
-	// trace( Ta x Ci in)
-	// Bet I need to move to real part with mult by -i
-	out_v[ss]()()(ax,b) = 0.5*(real(in_v[ss]()()(i2,i1)) - real(in_v[ss]()()(i1,i2)));
-	out_v[ss]()()(ay,b) = 0.5*(imag(in_v[ss]()()(i1,i2)) + imag(in_v[ss]()()(i2,i1)));
+  // Site-local operation
+  // Work with traceless anti-hermitian matrices for Lie algebra elements with Luscher's normalization convention
+  // Explicitly, T^a = -i t^a where t^a is hermitian generators in Grid convention, c.f., Compute_MpInvJx_dNxxdSy
+  //  in GaugeConfigurationMasked.h
+  // Used only in old implementation left for debugging purposes; to be deleted
+  template<typename T_out, typename T_in>
+  static accelerator_inline void LieAlgebraProject(T_out &out_v,const T_in &in_v, int b)
+  {
+    int N = ncolour;
+    int NNm1 = N * (N - 1);
+    int hNNm1= NNm1/2;
+    
+    for(int su2Index=0;su2Index<hNNm1;su2Index++){
+      int i1, i2;
+      su2SubGroupIndex(i1, i2, su2Index);
+      int ax = su2Index*2;
+      int ay = su2Index*2+1;
+      // in is traceless ANTI-hermitian
+      out_v()()(ax,b)=real(in_v()()(i2,i1));
+      out_v()()(ay,b)=imag(in_v()()(i1,i2));
+    }
+    for(int diagIndex=0;diagIndex<N-1;diagIndex++){
+      int k = diagIndex + 1; // diagIndex starts from 0
+      int a = NNm1+diagIndex;
+      RealD scale = 1.0/sqrt(2.0*k*(k+1));
+      auto tmp = in_v()()(0,0);
+      for(int i=1;i<k;i++){
+	tmp=tmp+in_v()()(i,i);
+      }
+      tmp = tmp - in_v()()(k,k)*((RealD) k);
+      out_v()()(a,b)=imag(tmp) * scale;
+    }
+  }
+  
+  // Lattice-wide operator
+  // Used only in old implementation left for debugging purposes; to be deleted
+  static void LieAlgebraProject(LatticeAlgebraMatrix &out,const LatticeMatrix &in, int b)
+  {
+    conformable(in, out);
+    GridBase *grid = out.Grid();
+    LatticeComplex tmp(grid);
+    Matrix ta;
+    
+    // for convention for algebra basis, see the comments in the LieAlgebraProject below.
+    autoView(out_v,out,AcceleratorWrite);
+    autoView(in_v,in,AcceleratorRead);
+    int N = ncolour;
+    int NNm1 = N * (N - 1);
+    int hNNm1= NNm1/2;
+    RealD sqrt_2 = sqrt(2.0);
+    Complex ci(0.0,1.0);
+    
+    const int nsimd = vMatrix::Nsimd();
+    accelerator_for(ss,grid->oSites(),nsimd,{
+	for(int su2Index=0;su2Index<hNNm1;su2Index++){
+	  int i1, i2;
+	  su2SubGroupIndex(i1, i2, su2Index);
+	  int ax = su2Index*2;
+	  int ay = su2Index*2+1;
+	  // in is traceless ANTI-hermitian
+	  coalescedWrite(out_v[ss]()()(ax,b),real(in_v(ss)()()(i2,i1)));
+	  coalescedWrite(out_v[ss]()()(ay,b),imag(in_v(ss)()()(i1,i2)));  
+	}
+	for(int diagIndex=0;diagIndex<N-1;diagIndex++){
+	  int k = diagIndex + 1; // diagIndex starts from 0
+	  int a = NNm1+diagIndex;
+	  RealD scale = 1.0/sqrt(2.0*k*(k+1));
+	  auto tmp = in_v(ss)()()(0,0);
+	  for(int i=1;i<k;i++){
+	    tmp=tmp+in_v(ss)()()(i,i);
+	  }
+	  tmp = tmp - in_v(ss)()()(k,k)*((RealD) k);
+	  coalescedWrite(out_v[ss]()()(a,b),imag(tmp) * scale);
+	}
       });
   }
-  for(int diagIndex=0;diagIndex<N-1;diagIndex++){
-    int k = diagIndex + 1; // diagIndex starts from 0
-    int a = NNm1+diagIndex;
-    RealD scale = 1.0/sqrt(2.0*k*(k+1));
-    accelerator_for(ss,grid->oSites(),vComplex::Nsimd(),{
-	auto tmp = in_v[ss]()()(0,0);
-	for(int i=1;i<k;i++){
-	  tmp=tmp+in_v[ss]()()(i,i);
+
+  // Site-local operation
+  // Work with traceless anti-hermitian matrices for Lie algebra elements with Luscher's normalization convention
+  // Explicitly, T^a = -i t^a where t^a is hermitian generators in Grid convention, c.f., Compute_MpInvJx_dNxxdSy
+  //  in GaugeConfigurationMasked.h
+  template<typename T_out, typename T_in>
+  static accelerator_inline void LieAlgebraProject(T_out &out_v,const T_in &in_v){
+    int N = ncolour;
+    int NNm1 = N * (N - 1);
+    int hNNm1= NNm1/2;
+
+#if 1
+    using namespace CplxMultExt;
+
+    // Compute N_{ab}
+    for(int su2Index_b=0;su2Index_b<hNNm1;su2Index_b++){
+      int ib1, ib2;
+      su2SubGroupIndex(ib1, ib2, su2Index_b);
+      int bx = su2Index_b*2;
+      int by = su2Index_b*2+1;
+      for(int su2Index=0;su2Index<hNNm1;su2Index++){
+        int i1, i2;
+        su2SubGroupIndex(i1, i2, su2Index); //i1<i2
+        int ax = su2Index*2;
+        int ay = su2Index*2+1;
+        // Compute: tr(T^a P) where P is traceless anti-hermitian
+        // Here, P = Ta(M) where M is defined below
+        // Note: T^b appearing in M (c.f. below) is 2it^b where t^b in Grid's convention, c.f. original def UtaU in the earlier version
+	
+        // real( 0.5*[M - M^\dag]_(i2, i1) )=real( 0.5*[T^b*in + in^\dag*T^b]_(i2, i1) ) <- M = T^b*in <- T^b: suN matrix in Luchang's convention; in: input matrix
+        //    sign flip of the complex part of the 2nd term is neglected, as real is taken
+        out_v()()(ax,bx) = 0.5*real( !(ib2^i2)*in_v()()(ib1,i1)-!(ib1^i2)*in_v()()(ib2,i1)
+                                    +!(ib1^i1)*in_v()()(ib2,i2)-!(ib2^i1)*in_v()()(ib1,i2) );
+        // imag( 0.5*[M - M^\dag]_(i1, i2) )
+        //    sign flip of the complex part of the 2nd term due to dagger
+        out_v()()(ay,bx) = 0.5*imag( !(ib2^i1)*in_v()()(ib1,i2)-!(ib1^i1)*in_v()()(ib2,i2)
+                                    -!(ib1^i2)*in_v()()(ib2,i1)+!(ib2^i2)*in_v()()(ib1,i1) );
+	
+	//    sign flip of the complex part of the 2nd term due to dagger <- imag turns into real due to i from T^b
+        out_v()()(ax,by) = -0.5*imag(  !(ib2^i2)*in_v()()(ib1,i1)+!(ib1^i2)*in_v()()(ib2,i1)
+                                      -!(ib1^i1)*in_v()()(ib2,i2)-!(ib2^i1)*in_v()()(ib1,i2) );
+
+        out_v()()(ay,by) = 0.5*real(  !(ib2^i1)*in_v()()(ib1,i2)+!(ib1^i1)*in_v()()(ib2,i2)
+                                     +!(ib1^i2)*in_v()()(ib2,i1)+!(ib2^i2)*in_v()()(ib1,i1) );
+      }
+      for(int diagIndex=0;diagIndex<N-1;diagIndex++){
+        int k = diagIndex + 1; // diagIndex starts from 0
+        int a = NNm1+diagIndex;
+        RealD scale = 1.0/sqrt(2.0*k*(k+1));
+        out_v()()(a,bx) = scale*imag(-((ib1<k)-k*(ib1==k))*in_v()()(ib2,ib1) + ((ib2<k)-k*(ib2==k))*in_v()()(ib1,ib2) );
+        out_v()()(a,by) = scale*real( ((ib1<k)-k*(ib1==k))*in_v()()(ib2,ib1) + ((ib2<k)-k*(ib2==k))*in_v()()(ib1,ib2) );
+      }
+    }
+    for(int diagIndex_b=0;diagIndex_b<N-1;diagIndex_b++){
+      int k_b = diagIndex_b + 1; // diagIndex starts from 0
+      int b = NNm1+diagIndex_b;
+      RealD scale_b = 2.0/sqrt(2.0*k_b*(k_b+1));
+      for(int su2Index=0;su2Index<hNNm1;su2Index++){
+        int i1, i2;
+        su2SubGroupIndex(i1, i2, su2Index);
+        int ax = su2Index*2;
+        int ay = su2Index*2+1;
+        out_v()()(ax,b)=-0.5*scale_b*imag( ((i2<k_b)-k_b*(i2==k_b))*in_v()()(i2,i1) - ((i1<k_b)-k_b*(i1==k_b))*in_v()()(i1,i2) );
+        out_v()()(ay,b)= 0.5*scale_b*real( ((i1<k_b)-k_b*(i1==k_b))*in_v()()(i1,i2) + ((i2<k_b)-k_b*(i2==k_b))*in_v()()(i2,i1) );
+      }
+      for(int diagIndex=0;diagIndex<N-1;diagIndex++){
+        int k = diagIndex + 1; // diagIndex starts from 0
+        int a = NNm1+diagIndex;
+        RealD scale = 1.0/sqrt(2.0*k*(k+1));
+        auto tmp = in_v()()(0,0);
+        int k_min = (k < k_b)? k:k_b;
+        for(int i=1;i<k_min;i++){
+          tmp=tmp+in_v()()(i,i);
+        }
+        int c = (k==k_b)? -k:1;
+        tmp = tmp - in_v()()(k_min,k_min)*k_min*c;
+        out_v()()(a,b)=real(tmp) * scale * scale_b;
+      }
+    }
+#else
+    // Compute N_{ab}
+    for(int su2Index_b=0;su2Index_b<hNNm1;su2Index_b++){
+      int ib1, ib2;
+      su2SubGroupIndex(ib1, ib2, su2Index_b);
+      int bx = su2Index_b*2;
+      int by = su2Index_b*2+1;
+      for(int su2Index=0;su2Index<hNNm1;su2Index++){
+	int i1, i2;
+	su2SubGroupIndex(i1, i2, su2Index); //i1<i2
+	int ax = su2Index*2;
+	int ay = su2Index*2+1;
+	// Compute: tr(T^a P) where P is traceless anti-hermitian
+	// Here, P = Ta(M) where M is defined below
+	// Note: T^b appearing in M (c.f. below) is 2it^b where t^b in Grid's convention, c.f. original def UtaU in the earlier version
+	
+	// real( 0.5*[M - M^\dag]_(i2, i1) )=real( 0.5*[T^b*in + in^\dag*T^b]_(i2, i1) ) <- M = T^b*in <- T^b: suN matrix in Luchang's convention; in: input matrix
+	//    sign flip of the complex part of the 2nd term is neglected, as real is taken
+	out_v()()(ax,bx) = 0.5*( !(ib2^i2)*real(in_v()()(ib1,i1))-!(ib1^i2)*real(in_v()()(ib2,i1))
+				 +!(ib1^i1)*real(in_v()()(ib2,i2))-!(ib2^i1)*real(in_v()()(ib1,i2)) );
+	// imag( 0.5*[M - M^\dag]_(i1, i2) )
+	//    sign flip of the complex part of the 2nd term due to dagger
+	out_v()()(ay,bx) = 0.5*( !(ib2^i1)*imag(in_v()()(ib1,i2))-!(ib1^i1)*imag(in_v()()(ib2,i2))
+				 -!(ib1^i2)*imag(in_v()()(ib2,i1))+!(ib2^i2)*imag(in_v()()(ib1,i1)) );
+	
+	//    sign flip of the complex part of the 2nd term due to dagger <- imag turns into real due to i from T^b
+	out_v()()(ax,by) = -0.5*(  !(ib2^i2)*imag(in_v()()(ib1,i1))+!(ib1^i2)*imag(in_v()()(ib2,i1))
+				   -!(ib1^i1)*imag(in_v()()(ib2,i2))-!(ib2^i1)*imag(in_v()()(ib1,i2)) );
+	
+	out_v()()(ay,by) = 0.5*(  !(ib2^i1)*real(in_v()()(ib1,i2))+!(ib1^i1)*real(in_v()()(ib2,i2))
+				  +!(ib1^i2)*real(in_v()()(ib2,i1))+!(ib2^i2)*real(in_v()()(ib1,i1)) );  
+      }
+      for(int diagIndex=0;diagIndex<N-1;diagIndex++){
+	int k = diagIndex + 1; // diagIndex starts from 0
+	int a = NNm1+diagIndex;
+	RealD scale = 1.0/sqrt(2.0*k*(k+1));
+	out_v()()(a,bx) = scale*(-((ib1<k)-k*(ib1==k))*imag(in_v()()(ib2,ib1)) + ((ib2<k)-k*(ib2==k))*imag(in_v()()(ib1,ib2)) );
+	out_v()()(a,by) = scale*( ((ib1<k)-k*(ib1==k))*real(in_v()()(ib2,ib1)) + ((ib2<k)-k*(ib2==k))*real(in_v()()(ib1,ib2)) );
+      }
+    }
+    for(int diagIndex_b=0;diagIndex_b<N-1;diagIndex_b++){
+      int k_b = diagIndex_b + 1; // diagIndex starts from 0
+      int b = NNm1+diagIndex_b;
+      RealD scale_b = 2.0/sqrt(2.0*k_b*(k_b+1));
+      for(int su2Index=0;su2Index<hNNm1;su2Index++){
+	int i1, i2;
+	su2SubGroupIndex(i1, i2, su2Index);
+	int ax = su2Index*2;
+	int ay = su2Index*2+1;
+	out_v()()(ax,b)=-0.5*scale_b*( ((i2<k_b)-k_b*(i2==k_b))*imag(in_v()()(i2,i1)) - ((i1<k_b)-k_b*(i1==k_b))*imag(in_v()()(i1,i2)) );
+	out_v()()(ay,b)= 0.5*scale_b*( ((i1<k_b)-k_b*(i1==k_b))*real(in_v()()(i1,i2)) + ((i2<k_b)-k_b*(i2==k_b))*real(in_v()()(i2,i1)) );
+      }
+      for(int diagIndex=0;diagIndex<N-1;diagIndex++){
+	int k = diagIndex + 1; // diagIndex starts from 0
+	int a = NNm1+diagIndex;
+	RealD scale = 1.0/sqrt(2.0*k*(k+1));
+	auto tmp = in_v()()(0,0);
+	int k_min = (k < k_b)? k:k_b;
+	for(int i=1;i<k_min;i++){
+	  tmp=tmp+in_v()()(i,i);
 	}
-	tmp = tmp - in_v[ss]()()(k,k)*k;
-	out_v[ss]()()(a,b) =imag(tmp) * scale;
+	int c = (k==k_b)? -k:1;
+	tmp = tmp - in_v()()(k_min,k_min)*((RealD) k_min*c);
+	out_v()()(a,b)=real(tmp) * scale * scale_b;
+      }
+    }
+#endif
+  }
+
+#if 0
+  // Lattice-wide operator
+  // Work with traceless anti-hermitian matrices for Lie algebra elements with Luscher's normalization convention
+  // Explicitly, T^a = -i t^a where t^a is hermitian generators in Grid convention, c.f., Compute_MpInvJx_dNxxdSy
+  //  in GaugeConfigurationMasked.h
+  // not used any more; to be removed
+  static void LieAlgebraProject(LatticeAlgebraMatrix &out,const LatticeMatrix &in)
+  {
+    conformable(in, out);
+    GridBase *grid = out.Grid();
+    autoView(out_v,out,AcceleratorWrite);
+    autoView(in_v,in,AcceleratorRead);
+    const int nsimd = vMatrix::Nsimd();
+    accelerator_for(ss,grid->oSites(),nsimd,{
+	typedef decltype(coalescedRead(out_v[0])) adj_mat;
+	adj_mat out;
+	LieAlgebraProject(out,in_v(ss));
+	coalescedWrite(out_v[ss],out);
       });
+  }
+  #endif
+
+  // Site-local operation
+  // Work with traceless anti-hermitian matrices for Lie algebra elements with Luscher's normalization convention
+  // Explicitly, T^a = -i t^a where t^a is hermitian generators in Grid convention, c.f., Compute_MpInvJx_dNxxdSy
+  //  in GaugeConfigurationMasked.h
+  // To be used in ComputeNxy
+  template<typename T_out, typename T_in>
+  static accelerator_inline void LieAlgebraProject(T_out &out_v,const T_in &inL_v,const T_in &inR_v){
+    int N = ncolour;
+    int NNm1 = N * (N - 1);
+    int hNNm1= NNm1/2;
+
+    //    using namespace CplxMultExt;
+    
+    // Compute N_{ab} = tr(T^a P^b) where P^b is traceless anti-hermitian
+    // Here, P^b = Ta(M^b) where M^b = adj(inL_v) * T'^b * inR_v => P^b = 0.5*( adj(inL_v(ss))*T'^b*inR_v(ss) + adj(inR_v(ss))*T'^b*inL_v(ss) )
+    // Note: 
+    //     * -tr(*)/2N part does not contribute and thus neglected
+    //     * T'^b appearing in M (c.f. below) is 2i t^b where t^b in Grid's convention, c.f. original def UtaU in the earlier version, diff from normalization of T^a
+    for(int su2Index_b=0;su2Index_b<hNNm1;su2Index_b++){
+      int ib1, ib2;
+      su2SubGroupIndex(ib1, ib2, su2Index_b);
+      int bx = su2Index_b*2;
+      int by = su2Index_b*2+1;
+      for(int su2Index=0;su2Index<hNNm1;su2Index++){
+	int i1, i2;
+	su2SubGroupIndex(i1, i2, su2Index); //i1<i2
+	int ax = su2Index*2;
+	int ay = su2Index*2+1;
+	// real( P^b[i2,i1])
+	out_v()()(ax,bx)=0.5*real( adj(inL_v()()(ib2,i2))*inR_v()()(ib1,i1) - adj(inL_v()()(ib1,i2))*inR_v()()(ib2,i1)
+				  +adj(inR_v()()(ib2,i2))*inL_v()()(ib1,i1) - adj(inR_v()()(ib1,i2))*inL_v()()(ib2,i1));
+	// imag( P^b[i1,i2])
+	out_v()()(ay,bx)=0.5*imag( adj(inL_v()()(ib2,i1))*inR_v()()(ib1,i2) - adj(inL_v()()(ib1,i1))*inR_v()()(ib2,i2)
+				  +adj(inR_v()()(ib2,i1))*inL_v()()(ib1,i2) - adj(inR_v()()(ib1,i1))*inL_v()()(ib2,i2));
+	// real( P^b[i2,i1])
+	out_v()()(ax,by)=-0.5*imag( adj(inL_v()()(ib2,i2))*inR_v()()(ib1,i1) + adj(inL_v()()(ib1,i2))*inR_v()()(ib2,i1)
+				   +adj(inR_v()()(ib2,i2))*inL_v()()(ib1,i1) + adj(inR_v()()(ib1,i2))*inL_v()()(ib2,i1));
+	// imag( P^b[i1,i2])) 	
+	out_v()()(ay,by)= 0.5*real( adj(inL_v()()(ib2,i1))*inR_v()()(ib1,i2) + adj(inL_v()()(ib1,i1))*inR_v()()(ib2,i2)
+				   +adj(inR_v()()(ib2,i1))*inL_v()()(ib1,i2) + adj(inR_v()()(ib1,i1))*inL_v()()(ib2,i2));
+      }
+      for(int diagIndex=0;diagIndex<N-1;diagIndex++){
+	int k = diagIndex + 1; // diagIndex starts from 0
+	int a = NNm1+diagIndex;
+	RealD scale = 1.0/sqrt(2.0*k*(k+1));
+	auto tmpx = imag( adj(inL_v()()(ib2,0))*inR_v()()(ib1,0) - adj(inL_v()()(ib1,0))*inR_v()()(ib2,0));
+	auto tmpy = real( adj(inL_v()()(ib2,0))*inR_v()()(ib1,0) + adj(inL_v()()(ib1,0))*inR_v()()(ib2,0));
+	for(int i=1;i<k;i++){
+	  tmpx = tmpx + imag( adj(inL_v()()(ib2,i))*inR_v()()(ib1,i) - adj(inL_v()()(ib1,i))*inR_v()()(ib2,i));
+	  tmpy = tmpy + real( adj(inL_v()()(ib2,i))*inR_v()()(ib1,i) + adj(inL_v()()(ib1,i))*inR_v()()(ib2,i));
+        }
+	tmpx = tmpx - k*imag( adj(inL_v()()(ib2,k))*inR_v()()(ib1,k) - adj(inL_v()()(ib1,k))*inR_v()()(ib2,k));
+	tmpy = tmpy - k*real( adj(inL_v()()(ib2,k))*inR_v()()(ib1,k) + adj(inL_v()()(ib1,k))*inR_v()()(ib2,k));
+	out_v()()(a,bx)= scale*tmpx;
+	out_v()()(a,by)= scale*tmpy;
+      }
     }
-}
+    for(int diagIndex_b=0;diagIndex_b<N-1;diagIndex_b++){
+      int k_b = diagIndex_b + 1; // diagIndex starts from 0
+      int b = NNm1+diagIndex_b;
+      RealD scale_b = 1.0/sqrt(2.0*k_b*(k_b+1));
+      
+      for(int su2Index=0;su2Index<hNNm1;su2Index++){
+	int i1, i2;
+	su2SubGroupIndex(i1, i2, su2Index);
+	int ax = su2Index*2;
+	int ay = su2Index*2+1;
+	
+	auto tmpx = adj(inL_v()()(0,i2))*inR_v()()(0,i1) + adj(inR_v()()(0,i2))*inL_v()()(0,i1);
+	auto tmpy = adj(inL_v()()(0,i1))*inR_v()()(0,i2) + adj(inR_v()()(0,i1))*inL_v()()(0,i2);
+	for(int i=1;i<k_b;i++){
+	  tmpx = tmpx + adj(inL_v()()(i,i2))*inR_v()()(i,i1) + adj(inR_v()()(i,i2))*inL_v()()(i,i1);
+	  tmpy = tmpy + adj(inL_v()()(i,i1))*inR_v()()(i,i2) + adj(inR_v()()(i,i1))*inL_v()()(i,i2);
+	}
+	//	tmpx = tmpx - ((RealD) k_b)*(adj(inL_v()()(k_b,i2))*inR_v()()(k_b,i1) + adj(inR_v()()(k_b,i2))*inL_v()()(k_b,i1));
+	//	tmpy = tmpy - ((RealD) k_b)*(adj(inL_v()()(k_b,i1))*inR_v()()(k_b,i2) + adj(inR_v()()(k_b,i1))*inL_v()()(k_b,i2));
+	tmpx = tmpx - ((RealD) k_b)*(adj(inL_v()()(k_b,i2))*inR_v()()(k_b,i1) + adj(inR_v()()(k_b,i2))*inL_v()()(k_b,i1));
+	tmpy = tmpy - ((RealD) k_b)*(adj(inL_v()()(k_b,i1))*inR_v()()(k_b,i2) + adj(inR_v()()(k_b,i1))*inL_v()()(k_b,i2));
+	
+	out_v()()(ax,b)=-scale_b*imag(tmpx);
+	out_v()()(ay,b)= scale_b*real(tmpy);
+      }
+      for(int diagIndex=0;diagIndex<N-1;diagIndex++){
+	int k = diagIndex + 1; // diagIndex starts from 0
+	int a = NNm1+diagIndex;
+	RealD scale = 1.0/sqrt(2.0*k*(k+1));
+	
+	auto tmp = -k_b*real(adj(inL_v()()(k_b,k))*inR_v()()(k_b,k));
+	for(int i=0;i<k_b;i++){
+	  tmp = tmp + real( adj(inL_v()()(i,k))*inR_v()()(i,k));
+	}
+	tmp = -k*tmp;
 
-  
+	for(int i=0;i<k;i++){
+	  for(int j=0;j<k_b;j++){
+	    tmp=tmp+real(adj(inL_v()()(j,i))*inR_v()()(j,i));
+	  }
+	  tmp=tmp-k_b*real(adj(inL_v()()(k_b,i))*inR_v()()(k_b,i));
+	}
+	out_v()()(a,b)=2.0 * tmp * scale * scale_b;
+      }
+    }
+  }
 };
     
 template <int ncolour>
diff --git a/Grid/qcd/utils/SUn.impl.h b/Grid/qcd/utils/SUn.impl.h
index bcccc7c48e..a2156bd1ff 100644
--- a/Grid/qcd/utils/SUn.impl.h
+++ b/Grid/qcd/utils/SUn.impl.h
@@ -62,7 +62,7 @@ static int su2subgroups(GroupName::SU) { return (ncolour * (ncolour - 1)) / 2; }
 //
 ////////////////////////////////////////////////////////////////////////
 template <class cplx, ONLY_IF_SU>
-static void generator(int lieIndex, iGroupMatrix<cplx> &ta, GroupName::SU) {
+static accelerator_inline void generator(int lieIndex, iGroupMatrix<cplx> &ta, GroupName::SU) {
   // map lie index to which type of generator
   int diagIndex;
   int su2Index;
@@ -82,7 +82,7 @@ static void generator(int lieIndex, iGroupMatrix<cplx> &ta, GroupName::SU) {
 }
 
 template <class cplx, ONLY_IF_SU>
-static void generatorSigmaY(int su2Index, iGroupMatrix<cplx> &ta) {
+static accelerator_inline void generatorSigmaY(int su2Index, iGroupMatrix<cplx> &ta) {
   ta = Zero();
   int i1, i2;
   su2SubGroupIndex(i1, i2, su2Index);
@@ -92,9 +92,9 @@ static void generatorSigmaY(int su2Index, iGroupMatrix<cplx> &ta) {
 }
 
 template <class cplx, ONLY_IF_SU>
-static void generatorSigmaX(int su2Index, iGroupMatrix<cplx> &ta) {
+static accelerator_inline void generatorSigmaX(int su2Index, iGroupMatrix<cplx> &ta) {
   ta = Zero();
-  cplx i(0.0, 1.0);
+  cplx i(Complex(0.0, 1.0));
   int i1, i2;
   su2SubGroupIndex(i1, i2, su2Index);
   ta()()(i1, i2) = i;
@@ -103,7 +103,7 @@ static void generatorSigmaX(int su2Index, iGroupMatrix<cplx> &ta) {
 }
 
 template <class cplx, ONLY_IF_SU>
-static void generatorDiagonal(int diagIndex, iGroupMatrix<cplx> &ta) {
+static accelerator_inline void generatorDiagonal(int diagIndex, iGroupMatrix<cplx> &ta) {
   // diag ({1, 1, ..., 1}(k-times), -k, 0, 0, ...)
   ta = Zero();
   int k = diagIndex + 1;                  // diagIndex starts from 0
diff --git a/Grid/qcd/utils/SUnAdjoint.h b/Grid/qcd/utils/SUnAdjoint.h
index 5924760c4a..374d4fcfa3 100644
--- a/Grid/qcd/utils/SUnAdjoint.h
+++ b/Grid/qcd/utils/SUnAdjoint.h
@@ -22,6 +22,7 @@
 //
 ////////////////////////////////////////////////////////////////////////
 
+
 NAMESPACE_BEGIN(Grid);
 
 template <int ncolour>
@@ -58,23 +59,21 @@ class SU_Adjoint : public SU<ncolour> {
   typedef Lattice<iScalar<iScalar<iVector<vComplex, Dimension> > > >  LatticeAdjVector;
 
   template <class cplx>
-  static void generator(int Index, iSUnAdjointMatrix<cplx> &iAdjTa) {
+  static accelerator_inline void generator(int Index, iSUnAdjointMatrix<cplx> &iAdjTa) {
     // returns i(T_Adj)^index necessary for the projectors
     // see definitions above
     iAdjTa = Zero();
-    iSUnMatrix<cplx> ta[ncolour * ncolour - 1];
+    iVector<iSUnMatrix<cplx>,Dimension> ta;
     iSUnMatrix<cplx> tmp;
 
     // FIXME not very efficient to get all the generators everytime
-    for (int a = 0; a < Dimension; a++) SU<ncolour>::generator(a, ta[a]);
+    for (int a = 0; a < Dimension; a++) SU<ncolour>::generator(a, ta(a));
 
     for (int a = 0; a < Dimension; a++) {
-      tmp = ta[a] * ta[Index] - ta[Index] * ta[a];
-      for (int b = 0; b < (ncolour * ncolour - 1); b++) {
-        iSUnMatrix<cplx> tmp1 = 2.0 * tmp * ta[b];  // 2.0 from the normalization
-        Complex iTr = TensorRemove(timesI(trace(tmp1)));
-        //iAdjTa()()(b, a) = iTr;
-        iAdjTa()()(a, b) = iTr;
+      tmp = ta(a) * ta(Index) - ta(Index) * ta(a);
+      for (int b = 0; b < Dimension; b++) {
+        iSUnMatrix<cplx> tmp1 = 2.0 * tmp * ta(b);  // 2.0 from the normalization
+        iAdjTa()()(a, b) = TensorRemove(timesI(trace(tmp1)));
       }
     }
   }
@@ -161,7 +160,65 @@ class SU_Adjoint : public SU<ncolour> {
     }
   }
 
+  // Turn complex 3x3 Lie algebra elem into its real 8x8 adj rep
+  static void make_adjoint_rep(LatticeAdjMatrix &out, const typename SU<ncolour>::LatticeMatrix &in) {
+    // Use real and totally anti-symmetric matrices as adjoint rep
+    //  => basis matrices for su(N) are anti-hermitian
+    // Take: T^a = i t^a where t^a: basis elem in Grid's convention; T^a satisfies Luscher's normalization convention -2tr[T^a,T^b] = \delta_{ab}
+    // Recall: adj rep generated above, iAdjTa, is multiplied by i
+    //   This makes adj rep here is consistent with the one written in terms of T^a
+    // So for adj rep, we follow Luscher's convention
+    
+    // in: traceless-anti-hermitian
 
+    GridBase *grid = out.Grid();
+      
+    autoView(out_v,out,AcceleratorWrite);
+    autoView(in_v,in,AcceleratorRead);
+    int N = ncolour;
+    int NNm1 = N * (N - 1);
+    int hNNm1= NNm1/2;
+    const int nsimd = vAMatrix::Nsimd();
+    accelerator_for(ss,grid->oSites(),nsimd,{
+	typedef decltype(coalescedRead(out_v[0])) adj_mat;
+	typedef decltype(coalescedRead(vTComplex())) computeComplex;
+	
+	adj_mat iTa, Tadj=Zero();
+	computeComplex c;
+	
+        for(int su2Index=0;su2Index<hNNm1;su2Index++){
+          int i1, i2;
+          SU<ncolour>::su2SubGroupIndex(i1, i2, su2Index);
+          int ax = su2Index*2;
+          int ay = su2Index*2+1;
+
+	  // multiply t^a by -i
+	  // cplx = -2.0*trace(ci*tb*Zx);
+	  // ZxAd = ZxAd + cplx * TRb;
+	  generator(ax,iTa);
+	  c()()() = 2.0*real(in_v(ss)()()(i2,i1));
+	  Tadj = Tadj + c*iTa;
+	  generator(ay,iTa);
+          c()()() = 2.0*imag(in_v(ss)()()(i1,i2));
+	  Tadj = Tadj + c*iTa;
+        }
+        for(int diagIndex=0;diagIndex<N-1;diagIndex++){
+          int k = diagIndex + 1; // diagIndex starts from 0
+          int a = NNm1+diagIndex;
+          RealD scale = 1.0/sqrt(2.0*k*(k+1));
+
+	  generator(a,iTa);
+          auto tmp = in_v(ss)()()(0,0);
+          for(int i=1;i<k;i++){
+            tmp=tmp+in_v(ss)()()(i,i);
+          }
+          tmp = tmp - in_v(ss)()()(k,k)*((RealD) k);
+	  c()()() = 2.0 * scale * imag(tmp);
+          Tadj = Tadj + c*iTa;
+        }
+	coalescedWrite(out_v[ss],Tadj);
+      });
+  }
 };
 
 
diff --git a/Grid/qcd/utils/WilsonLoops.h b/Grid/qcd/utils/WilsonLoops.h
index 23ad237c78..2af6e799ee 100644
--- a/Grid/qcd/utils/WilsonLoops.h
+++ b/Grid/qcd/utils/WilsonLoops.h
@@ -642,7 +642,28 @@ template <class Gimpl> class WilsonLoops : public Gimpl {
     return TensorRemove(Tq).real();
   }
 
+  // cp'ed from the above; added for convinience; used in Grid/HMC/analyze_snapshots.cc
+  static ComplexField TopologicalChargeDensity(const GaugeLorentz &U){
+    // 4d topological charge
+    assert(Nd==4);
+    // Bx = -iF(y,z), By = -iF(z,y), Bz = -iF(x,y)
+    GaugeMat Bx(U.Grid()), By(U.Grid()), Bz(U.Grid());
+    FieldStrength(Bx, U, Ydir, Zdir);
+    FieldStrength(By, U, Zdir, Xdir);
+    FieldStrength(Bz, U, Xdir, Ydir);
 
+    // Ex = -iF(t,x), Ey = -iF(t,y), Ez = -iF(t,z)
+    GaugeMat Ex(U.Grid()), Ey(U.Grid()), Ez(U.Grid());
+    FieldStrength(Ex, U, Tdir, Xdir);
+    FieldStrength(Ey, U, Tdir, Ydir);
+    FieldStrength(Ez, U, Tdir, Zdir);
+
+    double coeff = 8.0/(32.0*M_PI*M_PI);
+
+    ComplexField qfield = coeff*trace(Bx*Ex + By*Ey + Bz*Ez);
+    return qfield;
+  }
+  
   //Clover-leaf Wilson loop combination for arbitrary mu-extent M and nu extent N,  mu >= nu
   //cf  https://arxiv.org/pdf/hep-lat/9701012.pdf Eq 7  for 1x2 Wilson loop    
   //Clockwise ordering
diff --git a/Grid/tensors/Tensor_determinant.h b/Grid/tensors/Tensor_determinant.h
index a3862c94d9..8f958f7e13 100644
--- a/Grid/tensors/Tensor_determinant.h
+++ b/Grid/tensors/Tensor_determinant.h
@@ -49,7 +49,8 @@ template<class vtype,int N, typename std::enable_if< GridTypeMapper<vtype>::Tens
 accelerator_inline iScalar<vtype> Determinant(const iMatrix<vtype,N> &arg)
 {
   iMatrix<vtype,N> ret(arg);
-  iScalar<vtype> det = vtype(1.0);
+  iScalar<vtype> det;
+  det._internal = vtype(1.0);
   /* Conversion of matrix to upper triangular */
   for(int i = 0; i < N; i++){
     for(int j = 0; j < N; j++){
@@ -63,7 +64,7 @@ accelerator_inline iScalar<vtype> Determinant(const iMatrix<vtype,N> &arg)
   }      
 
   for(int i = 0; i < N; i++)
-    det *= ret._internal[i][i];   
+    det._internal *= ret._internal[i][i];   
 
   return det;
 }
diff --git a/Grid/threads/Accelerator.h b/Grid/threads/Accelerator.h
index 2c18796d5a..697c8b9f46 100644
--- a/Grid/threads/Accelerator.h
+++ b/Grid/threads/Accelerator.h
@@ -324,6 +324,8 @@ accelerator_inline int acceleratorSIMTlane(int Nsimd) {
  return 0;
 #endif
 } // SYCL specific
+//  printf("accelerator_for\n");					\
+  BACKTRACEFP(stdout); fflush(stdout);					\
 
 #define accelerator_for2dNB( iter1, num1, iter2, num2, nsimd, ... )	\
   theGridAccelerator->submit([&](sycl::handler &cgh) {		\
@@ -525,10 +527,11 @@ inline void *acceleratorAllocDevice(size_t bytes)
 inline void acceleratorFreeHost(void *ptr){ auto discard=hipFree(ptr);};
 inline void acceleratorFreeShared(void *ptr){ auto discard=hipFree(ptr);};
 inline void acceleratorFreeDevice(void *ptr){ auto discard=hipFree(ptr);};
-inline void acceleratorCopyToDevice(const void *from,void *to,size_t bytes)  { auto discard=hipMemcpy(to,from,bytes, hipMemcpyHostToDevice);}
-inline void acceleratorCopyFromDevice(const void *from,void *to,size_t bytes){ auto discard=hipMemcpy(to,from,bytes, hipMemcpyDeviceToHost);}
+inline void acceleratorCopyToDevice(const void *from,void *to,size_t bytes)  { GRID_TRACE("acceleratorCopyToDevice"); auto discard=hipMemcpy(to,from,bytes, hipMemcpyHostToDevice);}
+inline void acceleratorCopyFromDevice(const void *from,void *to,size_t bytes){ GRID_TRACE("acceleratorCopyFromDevice"); auto discard=hipMemcpy(to,from,bytes, hipMemcpyDeviceToHost);}
+inline void acceleratorMemSet(void *base,int value,size_t bytes) { GRID_TRACE("acceleratorMemset"); auto discard=hipMemset(base,value,bytes);}
 
-inline void acceleratorMemSet(void *base,int value,size_t bytes) { auto discard=hipMemset(base,value,bytes);}
+inline void acceleratorCopyDeviceToDeviceAsynch(const void *from,void *to,size_t bytes) // Asynch
 
 typedef int acceleratorEvent_t;
 
diff --git a/Grid/util/Init.cc b/Grid/util/Init.cc
index 23061372a1..65efd83ede 100644
--- a/Grid/util/Init.cc
+++ b/Grid/util/Init.cc
@@ -863,6 +863,15 @@ void Grid_heartbeat(void)
   it_val.it_interval = it_val.it_value;
   setitimer(ITIMER_REAL, &it_val, NULL);
 }
+void Grid_heartbeat_off(void)
+{
+  struct itimerval it_val;
+  it_val.it_value.tv_sec = 0;
+  it_val.it_value.tv_usec = 0;
+  it_val.it_interval = it_val.it_value;
+  setitimer(ITIMER_REAL, &it_val, NULL);
+}
+
 void Grid_exit_handler(void)
 {
   BACKTRACEFP(stdout);
diff --git a/Grid/util/Init.h b/Grid/util/Init.h
index b370077a83..1cdcff400d 100644
--- a/Grid/util/Init.h
+++ b/Grid/util/Init.h
@@ -43,6 +43,7 @@ void Grid_empty_signal_handler(int sig,siginfo_t *si,void * ptr);
 void Grid_debug_handler_init(void);
 void Grid_debug_heartbeat(void);
 void Grid_heartbeat(void);
+void Grid_heartbeat_off(void);
 void Grid_quiesce_nodes(void);
 void Grid_unquiesce_nodes(void);
 
diff --git a/HMC/ACC.hpp b/HMC/ACC.hpp
new file mode 100644
index 0000000000..13b83803f4
--- /dev/null
+++ b/HMC/ACC.hpp
@@ -0,0 +1,214 @@
+namespace Grid{
+  struct WFParameters: Serializable {
+    GRID_SERIALIZABLE_CLASS_MEMBERS(WFParameters,
+	    int, tau,
+	    std::string, data_name,
+	    std::string, path);
+       
+
+    template <class ReaderClass >
+    WFParameters(Reader<ReaderClass>& Reader){
+      read(Reader, "WilsonFlow", *this);
+    }
+
+  };
+
+  struct ConfParameters: Serializable {
+    GRID_SERIALIZABLE_CLASS_MEMBERS(ConfParameters,
+           std::string, conf_prefix,
+	   int, StartConfiguration,
+	   int, EndConfiguration);
+  
+    template <class ReaderClass >
+    ConfParameters(Reader<ReaderClass>& Reader){
+      read(Reader, "Configurations", *this);
+    }
+
+  };
+
+  struct ACFParameters: Serializable {
+    GRID_SERIALIZABLE_CLASS_MEMBERS(ACFParameters,
+	   int, MDtime_div_fac,				    
+	   int, MScut,	          // Cutoff time for Madras-Sokal approx. to error in ACC
+	   std::vector<int>, space_block_sizes,
+	   int, R,                // Summation radius of Master field tecnnique
+	   int, isFullTimeAvg);
+
+    template <class ReaderClass >
+    ACFParameters(Reader<ReaderClass>& Reader){
+      read(Reader, "Autocorrelations", *this);
+    }
+
+  };
+}
+
+template <class T> void readFile(T& out, std::string const fname){
+  Grid::emptyUserRecord record;
+  Grid::ScidacReader RD;
+  RD.open(fname);
+  RD.readScidacFieldRecord(out,record);
+  RD.close();
+}
+
+template <class T> void writeFile(T& in, std::string const fname){  
+#ifdef HAVE_LIME
+  // Ref: https://github.com/paboyle/Grid/blob/feature/scidac-wp1/tests/debug/Test_general_coarse_hdcg_phys48.cc#L111
+  // the same function also in ActionBase.h
+  std::cout << Grid::GridLogMessage << "Writes to: " << fname << std::endl;
+  Grid::emptyUserRecord record;
+  Grid::ScidacWriter WR(in.Grid()->IsBoss());
+  WR.open(fname);
+  WR.writeScidacFieldRecord(in,record,0);
+  WR.close();
+#endif
+}
+
+template <class L, typename A> void MS_approx(Grid::GridBase *Coarse, std::vector<L,A>  const& in, int T, int W, int block_size, int tau,
+					      std::string const e_name, std::string ACC_Type ){
+  using namespace Grid;
+  // applicable only if blocking the src fields
+  // Otherwise, for block averaged G_B(t), its error is hard to estimate as we need to consider correlation of G_x(t) at diff. sites wihting a block
+
+  RealD sumG[in.size()], ACC;
+  for (int t=0; t<T; t++) {
+    sumG[t] = TensorRemove(sum(in[t])).real();
+    std::cout << GridLogMessage << ACC_Type + " MFACC " + e_name + " (Madras-Sokal Approx): " << tau << " " << block_size << " " << 1 << " " << t << " "
+              << sumG[t]/sumG[0] << std::endl;
+    std::cout << GridLogMessage << ACC_Type + " MFCOV "  + e_name + " (Madras-Sokal Approx): " << tau << " " << block_size << " " << 1 << " " << t << " "
+	      << sumG[t] << std::endl;
+  }
+  
+  L var(Coarse), tmp(Coarse), Gt(Coarse), G0(Coarse);
+  for (int t=0; 2*t<T-W; t++){ // \sigma_{\rho}(t) uses \rho(t') up to 2t+W
+    ACC = sumG[t]/sumG[0];
+    var = (1/sumG[t])*in[t] - (1/sumG[0])*in[0];
+    var = (2.0/RealD(T))*var*var;
+    for (int k=1; k<W+t; k++) {
+      tmp = (1.0/sumG[t]/sumG[t]+2.0/sumG[0]/sumG[0])*in[k]*in[k] +
+	in[std::abs(k-t)]*( (1.0/sumG[t]/sumG[t])*in[k+t] - (2.0/sumG[0]/sumG[t])*in[k]) -(2.0/sumG[0]/sumG[t])*in[k]*in[k+t];
+      var = var + (2.0/RealD(T))*tmp;
+    }
+    std::cout << GridLogMessage << ACC_Type + " Variance " + e_name + " (Madras-Sokal Approx): " << tau << " " << block_size << " " << 1 << " " << t << " "
+	      << TensorRemove(sum(var)).real()*ACC*ACC << std::endl;
+  }
+}
+
+template <class L, typename A> void binning(Grid::GridBase *Coarse, std::vector<L,A>  const& in, int T, int n_bin, int block_size, int tau,
+					    std::string const e_name, std::string ACC_Type){
+  using namespace Grid;
+  // error of G_b(t) is first computed via binning for each block b
+  // G_b(t) with diff. b is considered as indep. measurement of G(t) => sigma_{G_(t)}^2 is reduced by V_b for the avg over blocks
+  // Then, proceed to compute error of G(t)/G(0)
+  // Computing ACC_b(t) = G_b(t)/G_b(0), 
+  
+  L avg(Coarse), avg0(Coarse), var(Coarse), var0(Coarse), cov(Coarse), tmp(Coarse);
+  for (int t=0; t<T; t++){
+    
+    avg = Zero(); var = Zero(); cov = Zero();
+    
+    // Find avg field over binns
+    for (int b=0; b<n_bin; b++)
+      avg = avg + (1/RealD(n_bin))*in[b*T + t];
+    if ( t== 0 ) avg0 = avg;
+    RealD ACC = TensorRemove(sum(avg)).real()/TensorRemove(sum(avg0)).real();
+    std::cout << GridLogMessage << ACC_Type + " MFACC " + e_name + " (Binning): " << tau << " " << block_size << " " << n_bin << " " << t << " "
+              << ACC << std::endl;
+
+    // Find variance of Gt/G0
+    for	(int b=0; b<n_bin; b++){
+      tmp = in[b*T + t] - avg;
+      var = var + (1/RealD((n_bin-1)*n_bin))*(tmp*tmp);
+      cov = cov + (1/RealD((n_bin-1)*n_bin))*tmp*(in[b*T] - avg0);
+    }
+    RealD G0 = real(sum(avg0)), Gt = real(sum(avg)); // the factor of 1/N_B is cancel out by 1/V^2 from Cov[G0,Gt] in the expression of var of Gt/G0
+    if ( t == 0 ) var0 = var;
+    std::cout << GridLogMessage << ACC_Type + " Variance " + e_name + " (Binning): " << tau << " " << block_size << " " << n_bin << " " << t << " "
+              << (real(sum(var))/Gt/Gt + real(sum(var0))/G0/G0 - 2.0*real(sum(cov))/G0/Gt)*ACC*ACC 
+	      << " " << (real(sum(var))/Gt/Gt + real(sum(var0))/G0/G0)*ACC*ACC << " " << - 2.0*real(sum(cov))/G0/Gt*ACC*ACC << std::endl;
+
+  }
+}
+
+template <class L, typename A> void binning2(Grid::GridBase *Coarse, std::vector<L,A>  const& in, int T, int n_bin, int block_size, int tau,
+					     std::string const e_name, std::string ACC_Type){
+  using namespace Grid;
+  // if bin_size > 30, the mean over binns approx dist. like Gaussian; the validity of this  estimation rests on Central Limit Theorem
+  // we avg G_x(t) over blocks first and then take the ratio G(t)/G(0) before avg over bins
+  // block avg might give a better sig as it incorpolates more stats; both should be unbiased
+  
+  RealD avg, avg0, var, tmp;
+  std::vector<RealD> in_sum(in.size());
+
+  for (int i=0; i<in.size(); i++) {
+    in_sum[i] = TensorRemove(sum(in[i])).real();
+    std::cout << GridLogMessage << ACC_Type + " MFCOV " + e_name + " (Binning2): " << tau << " " << block_size << " " << n_bin << " " << i << " "
+	      << in_sum[i] << std::endl;
+  }
+  
+  // Find avg field over binns
+  for (int t=0; t<T; t++){
+    avg = 0.0; var = 0.0;
+    for (int b=0; b<n_bin; b++) 
+      avg += (1/RealD(n_bin))*in_sum[b*T + t]/in_sum[b*T];
+    std::cout << GridLogMessage << ACC_Type + " MFACC " + e_name + " (Binning2): " << tau << " " << block_size << " " << n_bin << " " << t << " "
+              << avg << std::endl;
+
+    // Find variance of Gt/G0 
+    for (int b=0; b<n_bin; b++){
+      tmp = in_sum[b*T + t]/in_sum[b*T] - avg;
+      var = var + (1/RealD((n_bin-1)*n_bin))*(tmp*tmp);
+    }
+    std::cout << GridLogMessage << ACC_Type + " Variance " + e_name + " (Binning2): " << tau << " " << block_size << " " << n_bin << " " << t << " "
+              << var << " " << (1-avg*avg)*(1-avg*avg)/RealD(Coarse->gSites()-3)*RealD(n_bin) << std::endl;
+  }
+}
+
+template <class L, typename A> void MF_approx(Grid::GridBase *Coarse, std::vector<L,A>  const& in, int T, int R, int block_size, int tau,
+                                              std::string const e_name, std::string ACC_Type ){
+  using namespace Grid;
+
+  // Assume: block averaging
+  
+  int R_b = R/block_size;
+  RealD sumG[in.size()], ACC;
+  for (int t=0; t<T; t++) {
+    sumG[t] = TensorRemove(sum(in[t])).real();
+    std::cout << GridLogMessage << ACC_Type + " MFACC " + e_name + " (Master-Field Approx): " << tau << " " << block_size << " " << R << " " << t << " "
+              << sumG[t]/sumG[0] << std::endl;
+  }
+
+  RealD var, cov;
+  L tmp(Coarse), Gt(Coarse), G0(Coarse), one(Coarse); one = typename L::scalar_type(1.0,0.0);
+  G0 = in[0];
+  G0 = G0 - (sumG[0]/RealD(Coarse->gSites())) * one;
+  //std::cout << GridLogMessage << "in MF" << T <<  " " << block_size << " "<< R_b << " " << R << std::endl;
+  for(int t=0; t<T; t++) {
+    Gt = in[t] - (sumG[t]/RealD(Coarse->gSites())) * one;
+    tmp = Gt;
+    var = TensorRemove(sum(Gt*Gt)).real(); // C_tt
+    cov = TensorRemove(sum(Gt*G0)).real(); // C_t0
+    
+    for (int r=1; r<=R_b; r++){
+      // Shift Gt by y s.t. |y| <= r
+      for(int x=-r; x<=r; x++)
+	for(int y=-r; y<=r; y++)
+	  for(int z=-r; z<=r; z++)
+	    for(int s=-r; s<=r; s++) {
+	      int r2_c = x*x+y*y+z*z+s*s;
+	      if ( (r-1)*(r-1) < r2_c && r2_c <= r*r) {
+		
+		int d[4] = {x,y,z,s};
+		Gt = tmp;
+		for(int mu=0; mu<Nd; mu++) 
+		  Gt = Cshift(Gt, mu, d[mu]);
+		var += TensorRemove(sum(Gt*tmp)).real();
+		cov += TensorRemove(sum(Gt*G0)).real();
+	      }
+	    }
+      
+      // Note: No division by the volume of the coarse lattice
+      std::cout << GridLogMessage << ACC_Type + " Variance " + e_name + " (Master-Field Approx): " << tau << " " << block_size << " " << r << " " << t << " "
+		<< sumG[t] << " " << var << " " << cov << std::endl;
+    }
+  }
+}
diff --git a/HMC/ComputeWilsonFlow.cc b/HMC/ComputeWilsonFlow.cc
index 9bc935361d..043c2e619a 100644
--- a/HMC/ComputeWilsonFlow.cc
+++ b/HMC/ComputeWilsonFlow.cc
@@ -71,6 +71,7 @@ namespace Grid{
 template <class T> void writeFile(T& in, std::string const fname){  
 #ifdef HAVE_LIME
   // Ref: https://github.com/paboyle/Grid/blob/feature/scidac-wp1/tests/debug/Test_general_coarse_hdcg_phys48.cc#L111
+  // the same function also in ActionBase.h
   std::cout << Grid::GridLogMessage << "Writes to: " << fname << std::endl;
   Grid::emptyUserRecord record;
   Grid::ScidacWriter WR(in.Grid()->IsBoss());
@@ -91,7 +92,7 @@ int main(int argc, char **argv) {
   auto latt_size   = GridDefaultLatt();
   auto simd_layout = GridDefaultSimd(Nd, vComplex::Nsimd());
   auto mpi_layout  = GridDefaultMpi();
-  GridCartesian               Grid(latt_size, simd_layout, mpi_layout);
+  GridCartesian Grid(latt_size, simd_layout, mpi_layout);
   
   std::vector<int> seeds({1, 2, 3, 4, 5});
   GridSerialRNG sRNG;
@@ -153,7 +154,6 @@ int main(int argc, char **argv) {
     //However as  F_numu = -F_munu, only need to sum the trace of the squares of the following 6 field strengths:
     //F_01 F_02 F_03   F_12 F_13  F_23
     GaugeMat F(U.Grid());
-    //LatticeComplexD R(U.Grid());
     ComplexField R(U.Grid());
     R = Zero();
   
@@ -186,17 +186,13 @@ int main(int argc, char **argv) {
     //double coeff = 2.0 / (1.0 * Nd * (Nd - 1)) / 3.0;
     //Plq = coeff * Plq;
 
-
-    RealD WFlow_TC5Li   = WilsonLoops<PeriodicGimplR>::TopologicalCharge5Li(U);
-
+    RealD WFlow_TC5Li = WilsonLoops<PeriodicGimplR>::TopologicalCharge5Li(U);
+    
     int tau = std::round(t);
-
     std::string efile = file_pre + "E_dnsty_" + std::to_string(tau) + "_" + file_post;
-    //    writeFile(R,efile);
-
+    writeFile(R,efile);
     std::string tfile = file_pre + "Top_dnsty_" + std::to_string(tau) + "_" + file_post;
-    //    writeFile(qfield,tfile);
-
+    writeFile(qfield,tfile);
     std::string ufile = file_pre + "U_" + std::to_string(tau) + "_" + file_post;
     {
       //      PeriodicGimplR::GaugeField Ucopy = U;
@@ -209,7 +205,7 @@ int main(int argc, char **argv) {
     RealD E0 = real(peekSite(R,scoor));
     RealD T0 = real(peekSite(qfield,scoor));
     std::cout << GridLogMessage << "[WilsonFlow] Saved energy density (clover) & topo. charge density: "  << conf << " " << step << "  " << tau << "  "
-	      << "(E_avg,T_sum) " << E << " " << T << " (E, T at origin) " << E0 << " " << T0 << " Q5Li "<< WFlow_TC5Li << std::endl;
+	      << "(E_avg,T_sum) " << E << " " << T << " (E, T at origin) " << E0 << " " << T0 << " TC5Li " <<  WFlow_TC5Li << std::endl;
     
   });
   
diff --git a/HMC/autocovariance_NVS.cc b/HMC/autocovariance_NVS.cc
new file mode 100644
index 0000000000..d9af5ded5a
--- /dev/null
+++ b/HMC/autocovariance_NVS.cc
@@ -0,0 +1,220 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid
+
+Source file: ./tests/hmc/Test_WilsonFlow.cc
+
+Copyright (C) 2017
+
+Author: Guido Cossu <guido.cossu@ed.ac.uk>
+Author: Shuhei Yamamoto <syamamoto@bnl.gov>
+
+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 */
+#include <Grid/Grid.h>
+#include <string>
+#include "ACC.hpp"
+
+
+int main(int argc, char **argv) {
+  using namespace Grid;
+  
+  Grid_init(&argc, &argv);
+  GridLogLayout();
+
+  auto latt_size   = GridDefaultLatt();
+  auto simd_layout = GridDefaultSimd(Nd, vComplex::Nsimd());
+  auto mpi_layout  = GridDefaultMpi();
+  GridCartesian Grid(latt_size, simd_layout, mpi_layout);
+
+  typedef typename PeriodicGimplR::ComplexField ComplexField;
+
+  typedef Grid::XmlReader Serialiser;
+  Serialiser              Reader("input_ACF.xml", false, "root");
+  WFParameters            WFPar(Reader);
+  ConfParameters          CPar(Reader);
+  ACFParameters           APar(Reader);
+
+  std::string fname;
+  std::string file_path  = WFPar.path;
+  int W = APar.MScut, tau = WFPar.tau; //W: data[t>W] is not used in Madras-Sokal approximation of the error  
+  ComplexD coeff0, coeff1;
+  int total_configs = CPar.EndConfiguration - CPar.StartConfiguration + 1;
+  int T             = total_configs/APar.MDtime_div_fac; // #processed configs per bin
+  int arr_size      = T*APar.MDtime_div_fac;
+  if ( APar.MDtime_div_fac == 1 ) assert( APar.R>0 || W >= 100);
+
+  ComplexField A0(&Grid), A1(&Grid), avg(&Grid), one(&Grid); one=ComplexField::scalar_type(1.0,0.0);
+  std::vector<ComplexField> F(total_configs,&Grid), G(total_configs,&Grid), G2(total_configs,&Grid);
+
+  
+  std::cout << std::setprecision(15);
+
+  ////////////////////////
+  // Retrieve the fields
+  ////////////////////////
+  for (int i=0; i<total_configs; i++) {
+    int conf = CPar.StartConfiguration + i;
+    fname = file_path + WFPar.data_name + "_" + std::to_string(tau) + "_" + CPar.conf_prefix + ".";
+    readFile(F[i], fname + std::to_string(conf));
+#if 1 // DEBUG
+    RealD out = real(sum(F[i]));
+    std::cout << GridLogMessage << "NVS " + WFPar.data_name + " (conf, tau, val):   " << " " << conf << " " << tau << " "
+              << out/Real(Grid.gSites()) << std::endl;
+#endif
+  }
+
+  /////////////////////////////////////////
+  // Compute NVS autocovariance function
+  /////////////////////////////////////////
+  for (int i_bin=0, it=0; it<arr_size; i_bin++){  
+
+    // Compute estimation of ensemble average
+    avg = Zero();
+    for (int t=0 ; t<T; t++) {
+      avg = avg + (1/RealD(T))*F[i_bin*T + t];
+    }
+
+    // Compute NVS Auto-COV in two ways (<- corresponds to 2 diff estimation methods)
+    for (int t=0; t<T; t++, it++){
+
+      G[it] = Zero(); G2[it] = Zero();
+      
+      int n_src_configs = APar.isFullTimeAvg ? T-t : 1;
+      for (int i=0; i<n_src_configs; i++){ // if we store G_x(t) in a vector, we can loop over i first and then over t
+	  
+	A0 = F[i_bin*T + i];
+	A1 = F[i_bin*T + i + t];
+	
+	coeff0 = TensorRemove(sum(avg))/RealD(Grid.gSites());
+	// We replace \bar a(x) by << \bar a >>
+	G[it] = G[it] + ( A0 - coeff0*one )*( A1 - coeff0*one );
+	// Use \bar a(x) which is inaccurate est. of <a>
+	G2[it] = G2[it] + ( A0 - avg )*( A1 - avg );
+	}
+
+      G[it] = (1.0/((RealD) n_src_configs)) * G[it]; G2[it] = (1.0/((RealD) n_src_configs)) * G2[it];
+
+#if 0 // DEBUG: Monitor variation of COV
+      Coordinate scoor; for (int mu=0; mu < Nd; mu++) scoor[mu] = 0;
+      RealD Gt0 = real(peekSite(G[it],scoor));
+      ComplexD coeff0 = real(sum(G[it]))/RealD(Grid.gSites());
+      A0 = G[it] - coeff0*one;
+      A0 = A0*A0;
+      std::cout << GridLogMessage << "NVS COV at origin for " + WFPar.data_name + " " << t << " " << Gt0 << " " << real(sum(A0))/RealD(Grid.gSites()) << std::endl;
+#endif
+      
+      }// END: loop within a bin for MD time
+    }// END: loop over bins
+
+  ///////////////////////////
+  // Data for Error Estimate
+  ///////////////////////////
+  for(const int& bs : APar.space_block_sizes){
+    ///////////////////////////////////////////////////////////////////////////////////////
+    // Apply blocking on autocovariance function fields (block average or sparse sampling) 
+    ///////////////////////////////////////////////////////////////////////////////////////
+    Coordinate clatt_size(Nd);
+    for(int i=0;i<Nd;i++) clatt_size[i] = Grid.FullDimensions()[i]/bs;
+    GridCartesian Coarse(clatt_size,simd_layout,mpi_layout);
+    
+    // block averaging
+    std::vector<ComplexField> G_B(arr_size,&Coarse), G2_B(arr_size,&Coarse);
+    for (int i=0; i<arr_size; i++){
+      blockSum(G_B[i],G[i]);
+      G_B[i] = (1.0/RealD(std::pow(bs,Nd)))*G_B[i];
+      blockSum(G2_B[i],G2[i]);
+      G2_B[i] = (1.0/RealD(std::pow(bs,Nd)))*G2_B[i];
+    }
+    
+    // sparse sampling
+    std::vector<ComplexField> G_s(arr_size,&Coarse), G2_s(arr_size,&Coarse);
+    for (int i=0; i<arr_size; i++){
+      
+    LatticeInteger coor(&Grid);
+    ComplexField tmp(&Grid), filter(&Grid), zero(&Grid); filter = one; zero = Zero();
+    for (int d=0; d<Nd; d++) {
+      LatticeCoordinate(coor,d);
+      filter = where(mod(coor,bs)==Integer(0),filter,zero);
+    }
+    tmp = filter*G[i]; blockSum(G_s[i],tmp);
+    tmp = filter*G2[i]; blockSum(G2_s[i],tmp);
+    }
+    
+    /***********   Generate data  ************************************************
+      No Binning => Use Madras-Sokal approximation for error estimation
+      Otherwise  => error estimate via sample variance by binning over MD time
+    ******************************************************************************/
+    // TODO: add MF estimate (c.f. VS file)
+    int n_bin = total_configs/T;
+    if ( total_configs == T ){
+      // Madras-Sokal Approximation
+      //   Valid: when t << T
+      assert( APar.R>0 || T >= W );
+      MS_approx(&Coarse, G_B, T, W, bs, tau, "Blocked " + WFPar.data_name + " ACC", "NVS");
+      MS_approx(&Coarse, G_s, T, W, bs, tau, "Sparsed " + WFPar.data_name + " ACC", "NVS");
+      
+      MS_approx(&Coarse, G2_B, T, W, bs, tau, "Blocked " + WFPar.data_name + " ACC2", "NVS");
+      MS_approx(&Coarse, G2_s, T, W, bs, tau, "Sparsed " + WFPar.data_name + " ACC2", "NVS");
+    }
+    else {
+      // Binning
+      binning(&Coarse, G_B, T, n_bin, bs, tau, "Blocked " + WFPar.data_name + " ACC", "NVS");
+      binning(&Coarse, G_s, T, n_bin, bs, tau, "Sparsed " + WFPar.data_name + " ACC", "NVS");
+      
+      binning2(&Coarse, G_B, T, n_bin, bs, tau, "Blocked " + WFPar.data_name + " ACC", "NVS");
+      binning2(&Coarse, G_s, T, n_bin, bs, tau, "Sparsed " + WFPar.data_name + " ACC", "NVS");
+
+      
+      binning(&Coarse, G2_B, T, n_bin, bs, tau, "Blocked " + WFPar.data_name + " ACC2", "NVS");
+      binning(&Coarse, G2_s, T, n_bin, bs, tau, "Sparsed " + WFPar.data_name + " ACC2", "NVS");
+      
+      binning2(&Coarse, G2_B, T, n_bin, bs, tau, "Blocked " + WFPar.data_name + " ACC2", "NVS");
+      binning2(&Coarse, G2_s, T, n_bin, bs, tau, "Sparsed " + WFPar.data_name + " ACC2", "NVS");
+    }
+  }
+  
+  Grid_finalize();
+}  // main
+
+
+/*
+Input file example
+
+
+JSON
+
+{
+    "WilsonFlow":{
+	"steps": 200,
+	"step_size": 0.01,
+	"meas_interval": 50,
+  "maxTau": 2.0
+    },
+    "Configurations":{
+	"conf_prefix": "ckpoint_lat",
+	"rng_prefix": "ckpoint_rng",
+	"StartConfiguration": 3000,
+	"EndConfiguration": 3000,
+	"Skip": 5
+    }
+}
+
+
+*/
diff --git a/HMC/autocovariance_VS.cc b/HMC/autocovariance_VS.cc
new file mode 100644
index 0000000000..ba33bf7045
--- /dev/null
+++ b/HMC/autocovariance_VS.cc
@@ -0,0 +1,237 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid
+
+Source file: ./tests/hmc/Test_WilsonFlow.cc
+
+Copyright (C) 2017
+
+Author: Guido Cossu <guido.cossu@ed.ac.uk>
+Author: Shuhei Yamamoto <syamamoto@bnl.gov>
+
+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 */
+#include <Grid/Grid.h>
+#include <string>
+#include "ACC.hpp"
+
+
+int main(int argc, char **argv) {
+  using namespace Grid;
+  
+  Grid_init(&argc, &argv);
+  GridLogLayout();
+
+  auto latt_size   = GridDefaultLatt();
+  auto simd_layout = GridDefaultSimd(Nd, vComplex::Nsimd());
+  auto mpi_layout  = GridDefaultMpi();
+  GridCartesian Grid(latt_size, simd_layout, mpi_layout);
+
+  typedef typename PeriodicGimplR::ComplexField ComplexField;
+  
+  typedef Grid::XmlReader Serialiser;
+  Serialiser              Reader("input_ACF.xml", false, "root");
+  WFParameters            WFPar(Reader);
+  ConfParameters          CPar(Reader);
+  ACFParameters           APar(Reader);
+
+  std::string fname;
+  std::string file_path  = WFPar.path;
+  int W = APar.MScut, tau = WFPar.tau; // W: data[t>W] is not used in Madras-Sokal approximation of the error
+  ComplexD coeff0, coeff1;
+  int total_configs = CPar.EndConfiguration - CPar.StartConfiguration + 1;
+  int T             = total_configs/APar.MDtime_div_fac; // #processed configs per bin
+  int arr_size      = T*APar.MDtime_div_fac;  
+  if ( APar.MDtime_div_fac == 1 ) assert(APar.R < 0 || W >= 100 );
+
+  ComplexField A0(&Grid), A1(&Grid), one(&Grid); one=ComplexField::scalar_type(1.0,0.0);
+  std::vector<ComplexField> F(total_configs,&Grid), G(total_configs,&Grid), G2(total_configs,&Grid);
+
+  
+#if 0 // DEBUG: Test to see Cshift shifts the lattice cyclically
+  Coordinate latt_size0({4,4,4,4});
+  GridCartesian Grid0(latt_size0, simd_layout, mpi_layout);
+  ComplexField one0(&Grid0); one0=ComplexField::scalar_type(1.0,0.0);
+  LatticeInteger coor(&Grid0);
+  ComplexField tmp(&Grid0), filter(&Grid0), zero(&Grid0); filter = one0; zero=Zero();
+  int d[4] = {1,0,0,0};
+  for (int d=0; d<Nd; d++) {
+    LatticeCoordinate(coor,d);
+    filter= where(mod(coor,2)==Integer(0),filter,zero);
+  }
+  tmp = filter*one0;
+  for(int mu=0; mu<Nd; mu++)
+    tmp = Cshift(tmp, mu, d[mu]);
+  std::cout<< GridLogMessage << "Sfhited Sparsed Field " << tmp << std::endl;
+#endif
+  
+  std::cout << std::setprecision(15);
+
+  ////////////////////////////////
+  // Retrieve the fields
+  ///////////////////////////////
+  for (int i=0; i<total_configs; i++) {
+    int conf = CPar.StartConfiguration + i;
+    fname = file_path + WFPar.data_name + "_" + std::to_string(tau) + "_" + CPar.conf_prefix + ".";
+    readFile(A0, fname + std::to_string(conf));
+    F[i] = A0;
+#if 1 // DEBUG
+    RealD out = real(sum(F[i]));
+    std::cout << GridLogMessage << "LVS " + WFPar.data_name + " (conf, tau, val):   " << " " << conf << " " << tau << " "
+	      <<  out/Real(Grid.gSites()) << std::endl;
+#endif
+  }
+
+  ////////////////////////////////////////////////////////////////////////////////////////////////
+  // Compute VS autocovariance function in two ways (<- corresponds to 2 diff estimation methods)     
+  ////////////////////////////////////////////////////////////////////////////////////////////////
+  for (int i_bin=0, it=0; it<arr_size; i_bin++){
+    for (int t=0; t<T; t++, it++){
+	
+      G[it] = Zero(); G2[it] = Zero();
+
+      int n_src_configs = APar.isFullTimeAvg ? T-t : 1;
+      for (int i=0; i<n_src_configs; i++){
+
+	A0 = F[i_bin*T + i];
+	A1 = F[i_bin*T + i + t];
+
+	coeff0 = TensorRemove(sum(A0))/RealD(Grid.gSites());
+	coeff1 = TensorRemove(sum(A1))/RealD(Grid.gSites());
+	G2[it] = G2[it] + A0*A1 - (coeff0*coeff1)*one;
+	  
+	A0 = A0 - coeff0*one;
+	A1 = A1 - coeff1*one;
+	G[it] = G[it] + A0*A1;
+
+      }
+
+      G[it] = (1.0/((RealD) n_src_configs )) * G[it]; G2[it] = (1.0/((RealD) n_src_configs )) * G2[it];
+
+#if 0 // DEBUG: Monitor variation of COV
+      Coordinate scoor; for (int mu=0; mu < Nd; mu++) scoor[mu] = 0;
+      RealD Gt0 = real(peekSite(G[it],scoor));
+      coeff0 = real(sum(G[it]))/RealD(Grid.gSites());
+      A0 = G[it] - coeff0*one;
+      A0 = A0*A0;
+      std::cout << GridLogMessage << "LVS COV at origin for " + WFPar.data_name + " " << t << " " << Gt0 << " " << real(sum(A0))/RealD(Grid.gSites()) << std::endl;
+#endif
+      
+    }// END(t): loop within a bin for MD time
+  }// END(cong_s): loop over bins
+	
+  /////////////////////////////////
+  // Data for Error Estimate
+  ////////////////////////////////
+  for(const int& bs : APar.space_block_sizes){
+    ///////////////////////////////////////////////////////////////////////////////////////
+    // Apply blocking on autocovariance function fields (block average or sparse sampling)
+    ///////////////////////////////////////////////////////////////////////////////////////
+    Coordinate clatt_size(Nd);
+    for(int i=0;i<Nd;i++) clatt_size[i] = Grid.FullDimensions()[i]/bs;
+    GridCartesian Coarse(clatt_size,simd_layout,mpi_layout);
+    
+    // blocking
+    std::vector<ComplexField> G_B(arr_size,&Coarse), G2_B(arr_size,&Coarse);
+    for (int i=0; i<arr_size; i++){
+      blockSum(G_B[i],G[i]); blockSum(G2_B[i],G2[i]);
+      G_B[i] = (1/RealD(std::pow(bs,Nd)))*G_B[i]; G2_B[i] = (1/RealD(std::pow(bs,Nd)))*G2_B[i];
+    }
+    
+    // sparse sampling
+    std::vector<ComplexField> G_s(arr_size,&Coarse), G2_s(arr_size,&Coarse);
+    for (int i=0; i<arr_size; i++){
+      
+      LatticeInteger coor(&Grid);
+      ComplexField tmp(&Grid), filter(&Grid), zero(&Grid); filter = one; zero = Zero();
+      for (int d=0; d<Nd; d++) {
+	LatticeCoordinate(coor,d);
+	filter = where(mod(coor,bs)==Integer(0),filter,zero);
+      }
+      tmp = filter*G[i]; blockSum(G_s[i],tmp);
+      tmp = filter*G2[i]; blockSum(G2_s[i],tmp);
+    }
+
+    /***********   Generate data  ****************************************************************
+      No Binning => Use Madras-Sokal approximation or Master-field technique for error estimation
+      Otherwise  => error estimate via sample variance by binning over MD time
+    **********************************************************************************************/
+    int n_bin = total_configs/T;
+    if ( total_configs == T ){
+
+      // Madras-Sokal Approximation
+      if ( APar.R < 0 ){
+	//   Valid: when t << T
+	assert( T > W );
+	MS_approx(&Coarse, G_B, T, W, bs, tau, "Blocked " + WFPar.data_name + " ACC", "LVS");
+	MS_approx(&Coarse, G2_B, T, W, bs, tau, "Blocked " + WFPar.data_name + " ACC2", "LVS");
+	
+	MS_approx(&Coarse, G_s, T, W, bs, tau, "Sparsed " + WFPar.data_name + " ACC", "LVS");
+	MS_approx(&Coarse, G2_s, T, W, bs, tau, "Sparsed " + WFPar.data_name + " ACC2", "LVS");
+      }
+
+      // Master-Field Approximation
+      MF_approx(&Coarse, G_B, T, APar.R, bs, tau, "Blocked " + WFPar.data_name + " ACC", "LVS");
+      //MF_approx(&Coarse, G2_B, T, APar.R, bs, tau, "Blocked " + WFPar.data_name + " ACC2", "LVS");
+    }
+    else {
+      // Binning
+      binning(&Coarse, G_B, T, n_bin, bs, tau, "Blocked " + WFPar.data_name + " ACC", "LVS");
+      binning(&Coarse, G2_B, T, n_bin, bs, tau, "Blocked " + WFPar.data_name + " ACC2", "LVS");
+      
+      binning(&Coarse, G_s, T, n_bin, bs, tau, "Sparsed " + WFPar.data_name + " ACC", "LVS");
+      binning(&Coarse, G2_s, T, n_bin, bs, tau, "Sparsed " + WFPar.data_name + " ACC2", "LVS");
+
+      
+      binning2(&Coarse, G_B, T, n_bin, bs, tau, "Blocked " + WFPar.data_name + " ACC", "LVS");
+      binning2(&Coarse, G2_B, T, n_bin, bs, tau, "Blocked " + WFPar.data_name + " ACC2", "LVS");
+      
+      binning2(&Coarse, G_s, T, n_bin, bs, tau, "Sparsed " + WFPar.data_name + " ACC", "LVS");
+      binning2(&Coarse, G2_s, T, n_bin, bs, tau, "Sparsed " + WFPar.data_name + " ACC2", "LVS");
+    }
+  }
+    
+  Grid_finalize();
+}  // main
+
+
+/*
+Input file example
+
+
+JSON
+
+{
+    "WilsonFlow":{
+	"steps": 200,
+	"step_size": 0.01,
+	"meas_interval": 50,
+  "maxTau": 2.0
+    },
+    "Configurations":{
+	"conf_prefix": "ckpoint_lat",
+	"rng_prefix": "ckpoint_rng",
+	"StartConfiguration": 3000,
+	"EndConfiguration": 3000,
+	"Skip": 5
+    }
+}
+
+
+*/
diff --git a/benchmarks/Benchmark_comms.cc b/benchmarks/Benchmark_comms.cc
index 7755635a21..65f63bfa12 100644
--- a/benchmarks/Benchmark_comms.cc
+++ b/benchmarks/Benchmark_comms.cc
@@ -166,7 +166,7 @@ int main (int argc, char ** argv)
   }  
 
 
-
+#ifdef GRID_OMP
   std::cout<<GridLogMessage << "===================================================================================================="<<std::endl;
   std::cout<<GridLogMessage << "= Benchmarking sequential STENCIL halo exchange in "<<nmu<<" dimensions"<<std::endl;
   std::cout<<GridLogMessage << "===================================================================================================="<<std::endl;
@@ -249,7 +249,7 @@ int main (int argc, char ** argv)
  
     }
   }    
-
+#endif
 
 #ifdef GRID_OMP
   std::cout<<GridLogMessage << "===================================================================================================="<<std::endl;
diff --git a/systems/Aurora/config-command b/systems/Aurora/config-command
index 02d57c0dd3..1516996912 100644
--- a/systems/Aurora/config-command
+++ b/systems/Aurora/config-command
@@ -1,8 +1,23 @@
-export MPFR=`spack find --paths mpfr    | grep ^mpfr  | awk '{print $2}' `
-export GMP=`spack find --paths gmp   | grep ^gmp  | awk '{print $2}' `
-export CLIME=`spack find --paths c-lime  | grep ^c-lime  | awk '{print $2}' `
-export UNWIND=`spack find --paths libunwind  | grep ^libunwind  | awk '{print $2}' `
-export HDF5=/opt/cray/pe/hdf5/1.12.2.3/gnu/9.1
+#Ahead of time compile for PVC
+
+export LDFLAGS="-fiopenmp -fsycl -fsycl-device-code-split=per_kernel -fsycl-targets=spir64 -Xs -device -Xs pvc -fsycl-device-lib=all -lze_loader -L${MKLROOT}/lib -qmkl=parallel  -fsycl  -lsycl -lnuma -L/opt/aurora/24.180.3/spack/unified/0.8.0/install/linux-sles15-x86_64/oneapi-2024.07.30.002/numactl-2.0.14-7v6edad/lib -fPIC -fsycl-max-parallel-link-jobs=16 -fno-sycl-rdc" 
+export CXXFLAGS="-O3 -fiopenmp -fsycl-unnamed-lambda -fsycl -Wno-narrowing -Wno-tautological-compare -qmkl=parallel  -fsycl -fno-exceptions -I/opt/aurora/24.180.3/spack/unified/0.8.0/install/linux-sles15-x86_64/oneapi-2024.07.30.002/numactl-2.0.14-7v6edad/include/ -fPIC "
+
+### Original
+#export LDFLAGS="-fiopenmp -fsycl -fsycl-device-code-split=per_kernel -fsycl-targets=spir64_gen -Xs -device -Xs pvc -fsycl-device-lib=all -lze_loader -L${MKLROOT}/lib -qmkl=parallel  -fsycl  -lsycl -lnuma -L/opt/aurora/24.180.3/spack/unified/0.8.0/install/linux-sles15-x86_64/oneapi-2024.07.30.002/numactl-2.0.14-7v6edad/lib -fPIC -fsycl-max-parallel-link-jobs=16 -fno-sycl-rdc" 
+#export CXXFLAGS="-O3 -fiopenmp -fsycl-unnamed-lambda -fsycl -Wno-narrowing -Wno-tautological-compare -qmkl=parallel  -fsycl -fno-exceptions -I/opt/aurora/24.180.3/spack/unified/0.8.0/install/linux-sles15-x86_64/oneapi-2024.07.30.002/numactl-2.0.14-7v6edad/include/ -fPIC"
+
+### JIT compile 
+#export LDFLAGS="-fiopenmp -fsycl -fsycl-device-code-split=per_kernel  -fsycl-device-lib=all -lze_loader -L${MKLROOT}/lib -qmkl=parallel  -fsycl  -lsycl " 
+#export CXXFLAGS="-O3 -fiopenmp -fsycl-unnamed-lambda -fsycl -Wno-tautological-compare -qmkl=parallel  -fsycl -fno-exceptions "
+
+### upstream
+#export MPFR=`spack find --paths mpfr    | grep ^mpfr  | awk '{print $2}' `
+#export GMP=`spack find --paths gmp   | grep ^gmp  | awk '{print $2}' `
+#export CLIME=`spack find --paths c-lime  | grep ^c-lime  | awk '{print $2}' `
+#export UNWIND=`spack find --paths libunwind  | grep ^libunwind  | awk '{print $2}' `
+#export HDF5=/opt/cray/pe/hdf5/1.12.2.3/gnu/9.1
+
 ../../configure \
 	--enable-simd=GPU \
 	--enable-gen-simd-width=64 \
@@ -16,14 +31,14 @@ export HDF5=/opt/cray/pe/hdf5/1.12.2.3/gnu/9.1
 	--enable-accelerator-aware-mpi=no \
 	--enable-unified=no \
 	--with-lime=$CLIME \
-	--with-hdf5=${HDF5} \
 	--with-gmp=$GMP \
 	--with-mpfr=$MPFR \
-	--with-unwind=$UNWIND \
 	MPICXX=mpicxx \
 	CXX=icpx \
 	LDFLAGS="-fiopenmp -fsycl -fsycl-device-code-split=per_kernel -fsycl-device-lib=all -lze_loader -L${MKLROOT}/lib -qmkl=parallel -lsycl -fsycl-max-parallel-link-jobs=16 -fno-sycl-rdc -lnuma" \
 	CXXFLAGS="-fPIC -fiopenmp -fsycl-unnamed-lambda -fsycl -Wno-tautological-compare -qmkl=parallel"
 
 
-
+#	--with-unwind=$UNWIND \
+#	--prefix $HOME/gpt-install \
+#	--enable-debug \
diff --git a/systems/Aurora/sourceme.sh b/systems/Aurora/sourceme.sh
index 4009609c07..f21256f5f9 100644
--- a/systems/Aurora/sourceme.sh
+++ b/systems/Aurora/sourceme.sh
@@ -1,13 +1,43 @@
-export HTTP_PROXY=http://proxy.alcf.anl.gov:3128
-export HTTPS_PROXY=http://proxy.alcf.anl.gov:3128
-export http_proxy=http://proxy.alcf.anl.gov:3128
-export https_proxy=http://proxy.alcf.anl.gov:3128
-git config --global http.proxy http://proxy.alcf.anl.gov:3128
+
+### After Update (Hopefully Working) Setup
+module reset
+module use /opt/aurora/24.347.0/spack/unified/0.9.2/install/modulefiles/Core
+module use /opt/aurora/24.347.0/spack/unified/0.9.2/install/modulefiles/oneapi/2025.0.5
+#module load py-numpy
+module unload mpich
+module unload oneapi
+module use /soft/compilers/oneapi/2025.1.0/modulefiles
+module load oneapi/public/2025.1.0
+module use /home/bertoni/mpich_module/
+module load aurora_test_2025.1
+
+### Before Update Setup
+#module restore
+#module load oneapi/eng-compiler/2024.07.30.002
+#module load mpich/opt/4.3.0rc3
+#module load cmake
+#module unload gcc/12.2.0
+#module use /opt/aurora/24.347.0/spack/unified/0.9.2/install/modulefiles/Core
+#module load gcc/13.3.0
+
+##########
+##module load oneapi/release/2023.12.15.001
+##module load mpich/icc-all-debug-pmix-gpu/52.2
+##module load mpich-config/mode/deterministic
+##module load intel_compute_runtime/release/821.35
 
 source ~/spack/share/spack/setup-env.sh 
 spack load c-lime
-spack load openssl@3.3.1%gcc@12.2.0
-spack load unwind
-export UNWIND=`spack find --paths libunwind  | grep ^libunwind  | awk '{print $2}' `
+spack load openssl
+spack load gmp
+spack load mpfr
 export CLIME=`spack find --paths c-lime | grep ^c-lime | awk '{print $2}' `
-export SYCL_PROGRAM_COMPILE_OPTIONS="-ze-opt-large-register-file"
+export GMP=`spack find --paths gmp | grep ^gmp | awk '{print $2}' `
+export MPFR=`spack find --paths mpfr | grep ^mpfr | awk '{print $2}' `
+
+git config --global http.proxy http://proxy.alcf.anl.gov:3128
+
+#spack load openssl@3.3.1%gcc@12.2.0
+#spack load unwind
+#export UNWIND=`spack find --paths libunwind  | grep ^libunwind  | awk '{print $2}' `
+#export SYCL_PROGRAM_COMPILE_OPTIONS="-ze-opt-large-register-file"
diff --git a/systems/Frontier/config-command b/systems/Frontier/config-command
index 7561fb155c..049283c168 100644
--- a/systems/Frontier/config-command
+++ b/systems/Frontier/config-command
@@ -1,4 +1,5 @@
 CLIME=`spack find --paths c-lime@2-3-9 | grep c-lime| cut -c 15-`
+GMP=`spack find --paths gmp | grep gmp| cut -c 12-`
 ../../configure --enable-comms=mpi-auto \
 --with-lime=$CLIME \
 --enable-unified=no \
@@ -9,13 +10,13 @@ CLIME=`spack find --paths c-lime@2-3-9 | grep c-lime| cut -c 15-`
 --disable-gparity \
 --disable-fermion-reps \
 --enable-simd=GPU \
---with-gmp=$OLCF_GMP_ROOT \
+--enable-accelerator-cshift \
+--with-gmp=$GMP \
+--with-fftw=$FFTW_DIR/.. \
 --with-mpfr=/opt/cray/pe/gcc/mpfr/3.1.4/ \
 --disable-fermion-reps \
 CXX=hipcc MPICXX=mpicxx \
-CXXFLAGS="-fPIC -I${ROCM_PATH}/include/ -I${MPICH_DIR}/include -L/lib64 " \
- LDFLAGS="-L/lib64 -L${ROCM_PATH}/lib -L${MPICH_DIR}/lib -lmpi -L${CRAY_MPICH_ROOTDIR}/gtl/lib -lmpi_gtl_hsa -lhipblas -lrocblas -lhipfft"
-
-
-
-
+CXXFLAGS="-fPIC -I{$ROCM_PATH}/include/ -I${MPICH_DIR}/include -L/lib64 " \
+ LDFLAGS="-L/lib64 -L${MPICH_DIR}/lib -lmpi -L${CRAY_MPICH_ROOTDIR}/gtl/lib -lmpi_gtl_hsa -lamdhip64 -lhipblas -lrocblas"
+#hipcc
+#--enable-tracing=roctx or timer
diff --git a/systems/Frontier/sourceme.sh b/systems/Frontier/sourceme.sh
index 71028ff785..d7f1afe2b3 100644
--- a/systems/Frontier/sourceme.sh
+++ b/systems/Frontier/sourceme.sh
@@ -1,9 +1,15 @@
 
 echo spack
 . /autofs/nccs-svm1_home1/paboyle/Crusher/Grid/spack/share/spack/setup-env.sh
+spack load c-lime gmp@6.2.1
+module load emacs 
+module load PrgEnv-cray-amd/8.5.0
+module load rocm #/6.2.0
+#module load PrgEnv-gnu #present in the original
+#module load rocm/6.0.0 #present in the original
 
-module load cce/15.0.1
-module load rocm/5.3.0
+module load cray-mpich
+#module load gmp-6.2.1-gcc-7.5.0-h5wpx5n  #gmp
 module load cray-fftw
 module load craype-accel-amd-gfx90a
 
diff --git a/tests/Test_meson_field.cc b/tests/Test_meson_field.cc
index fa428d6a2c..4ca1d9c64e 100644
--- a/tests/Test_meson_field.cc
+++ b/tests/Test_meson_field.cc
@@ -146,11 +146,13 @@ int main(int argc, char *argv[])
   using Default_Writer = Grid::BinaryWriter;
   FileName.append(".bin");
 #endif
+
   {
     Default_Writer w(FileName);
     write(w,"MesonField",Mpp);
     write(w,"AslashField",App);
   }
+
   // epilogue
   std::cout << GridLogMessage << "Grid is finalizing now" << std::endl;
   Grid_finalize();
diff --git a/tests/forces/Test_fthmc.cc b/tests/forces/Test_fthmc.cc
index 2b0d20d5af..704b334e2c 100644
--- a/tests/forces/Test_fthmc.cc
+++ b/tests/forces/Test_fthmc.cc
@@ -61,11 +61,12 @@ void ForceTest(Action<LatticeGaugeField> &action,ConfigurationBase<LatticeGaugeF
   std::cout << GridLogMessage << "+++++++++++++++++++++++++++++++++++++++++++++++++++++++++"<<std::endl;
   
   Gimpl::generate_momenta(P,sRNG,RNG4);
+  std::cout << GridLogMessage << "P mag: " << Gimpl::FieldSquareNorm(P) << std::endl;
   //  Filter.applyFilter(P);
   std::cout << GridLogMessage << "Initial momenta " << norm2(P) << std::endl;
 
   action.refresh(smU,sRNG,RNG4);
-
+  
   std::cout << GridLogMessage << "+++++++++++++++++++++++++++++++++++++++++++++++++++++++++"<<std::endl;
   std::cout << GridLogMessage << " Action "<<action.action_name()<<std::endl;
   std::cout << GridLogMessage << "+++++++++++++++++++++++++++++++++++++++++++++++++++++++++"<<std::endl;
@@ -88,6 +89,7 @@ void ForceTest(Action<LatticeGaugeField> &action,ConfigurationBase<LatticeGaugeF
   DumpSliceNorm("Force",UdSdU,Nd-1);
   
   Gimpl::update_field(P,U,eps);
+  std::cout << GridLogMessage << "U mag after update: " << Gimpl::FieldSquareNorm(smU.get_U(false)) << " " << Gimpl::FieldSquareNorm(smU.get_U(true))<<std::endl;
   smU.set_Field(U);
 
   std::cout << GridLogMessage << "+++++++++++++++++++++++++++++++++++++++++++++++++++++++++"<<std::endl;
@@ -154,10 +156,10 @@ int main (int argc, char ** argv)
 
 
   
-  WilsonGaugeActionR  PlaqAction(6.0);
-  IwasakiGaugeActionR RectAction(2.13);
-  PlaqAction.is_smeared = true;  
-  RectAction.is_smeared = true;  
+  WilsonGaugeActionR  PlaqAction(6.0); PlaqAction.reset_timer();
+  IwasakiGaugeActionR RectAction(2.13); RectAction.reset_timer();
+  PlaqAction.is_smeared = true; //true;  
+  RectAction.is_smeared = true; //true;  
 
   ////////////////////////////////////
   // Fermion Action
@@ -178,7 +180,7 @@ int main (int argc, char ** argv)
   double MaxCGIterations = 50000;
   ConjugateGradient<LatticeFermion>  CG(StoppingCondition,MaxCGIterations);
 
-  TwoFlavourRatioPseudoFermionAction<FimplD> Nf2(PVPeriodic, DdwfPeriodic,CG,CG);
+  TwoFlavourRatioPseudoFermionAction<FimplD> Nf2(PVPeriodic, DdwfPeriodic,CG,CG); Nf2.reset_timer();
   Nf2.is_smeared = true;  
   
   ////////////////////////////////////////////////
@@ -189,7 +191,7 @@ int main (int argc, char ** argv)
   SmearedConfigurationMasked<PeriodicGimplR> SmartConfig(UGrid,2*Nd,Smearer);
   SmearedConfiguration<PeriodicGimplR> StoutConfig(UGrid,1,Smearer);
 
-  JacobianAction<PeriodicGimplR> Jacobian(&SmartConfig);
+  JacobianAction<PeriodicGimplR> Jacobian(&SmartConfig); Jacobian.reset_timer();
   
   ////////////////////////////////////////////////
   // Run some tests
@@ -210,6 +212,15 @@ int main (int argc, char ** argv)
   SmartConfig.set_Field(U);
   ForceTest<GimplTypesR>(Nf2,SmartConfig,FilterNone);
 
+  std::cout << " *********  FIELD TRANSFORM SMEARING: no smearing ***** "<<std::endl;
+
+  PlaqAction.is_smeared = false;
+  RectAction.is_smeared = false;
+  SmartConfig.set_Field(U);
+  ForceTest<GimplTypesR>(PlaqAction,SmartConfig,FilterNone);
+  SmartConfig.set_Field(U);
+  ForceTest<GimplTypesR>(RectAction,SmartConfig,FilterNone);
+  
   std::cout << " *********    STOUT SMEARING ***** "<<std::endl;
 
   StoutConfig.set_Field(U);