diff --git a/lib/GridConfig.h.in b/lib/GridConfig.h.in
index 1ec0ecf5..bcb8d811 100644
--- a/lib/GridConfig.h.in
+++ b/lib/GridConfig.h.in
@@ -33,9 +33,6 @@
 /* Support AVX2 (Advanced Vector Extensions 2) instructions */
 #undef HAVE_AVX2
 
-/* define if the compiler supports basic C++11 syntax */
-#undef HAVE_CXX11
-
 /* Define to 1 if you have the declaration of `be64toh', and to 0 if you
    don't. */
 #undef HAVE_DECL_BE64TOH
diff --git a/lib/lattice/Lattice_base.h b/lib/lattice/Lattice_base.h
index 2a365cdd..9c9ef044 100644
--- a/lib/lattice/Lattice_base.h
+++ b/lib/lattice/Lattice_base.h
@@ -60,6 +60,11 @@ public:
     GridBase *_grid;
     int checkerboard;
     std::vector<vobj,alignedAllocator<vobj> > _odata;
+    
+    // to pthread need a computable loop where loop induction is not required
+    int begin(void) { return 0;};
+    int end(void)   { return _odata.size(); }
+    vobj & operator[](int i) { return _odata[i]; };
 
 public:
     typedef typename vobj::scalar_type scalar_type;
diff --git a/lib/qcd/action/fermion/FermionOperator.h b/lib/qcd/action/fermion/FermionOperator.h
index 1aa2f8f0..a8f10380 100644
--- a/lib/qcd/action/fermion/FermionOperator.h
+++ b/lib/qcd/action/fermion/FermionOperator.h
@@ -39,10 +39,23 @@ namespace Grid {
       virtual void Dhop  (const FermionField &in, FermionField &out,int dag)=0;
       virtual void DhopOE(const FermionField &in, FermionField &out,int dag)=0;
       virtual void DhopEO(const FermionField &in, FermionField &out,int dag)=0;
+      virtual void DhopDir(const FermionField &in, FermionField &out,int dir,int disp)=0; // implemented by WilsonFermion and WilsonFermion5D
+
+      // force terms; five routines; default to Dhop on diagonal
+      virtual void MDeriv  (LatticeGaugeField &mat,const FermionField &U,const FermionField &V,int dag){DhopDeriv(mat,U,V,dag);};
+      virtual void MoeDeriv(LatticeGaugeField &mat,const FermionField &U,const FermionField &V,int dag){DhopDerivOE(mat,U,V,dag);};
+      virtual void MeoDeriv(LatticeGaugeField &mat,const FermionField &U,const FermionField &V,int dag){DhopDerivEO(mat,U,V,dag);};
+      virtual void MooDeriv(LatticeGaugeField &mat,const FermionField &U,const FermionField &V,int dag){mat=zero;};
+      virtual void MeeDeriv(LatticeGaugeField &mat,const FermionField &U,const FermionField &V,int dag){mat=zero;};
+
+      virtual void DhopDeriv  (LatticeGaugeField &mat,const FermionField &U,const FermionField &V,int dag)=0;
+      virtual void DhopDerivEO(LatticeGaugeField &mat,const FermionField &U,const FermionField &V,int dag)=0;
+      virtual void DhopDerivOE(LatticeGaugeField &mat,const FermionField &U,const FermionField &V,int dag)=0;
+
+
+      virtual void  Mdiag  (const FermionField &in, FermionField &out) { Mooee(in,out);};   // Same as Mooee applied to both CB's
+      virtual void  Mdir   (const FermionField &in, FermionField &out,int dir,int disp)=0;   // case by case Wilson, Clover, Cayley, ContFrac, PartFrac
 
-      virtual void  Mdiag(const FermionField &in, FermionField &out) { Mooee(in,out);};   // Same as Mooee applied to both CB's
-      virtual void  Mdir (const FermionField &in, FermionField &out,int dir,int disp)=0;   // case by case Wilson, Clover, Cayley, ContFrac, PartFrac
-      virtual void  DhopDir(const FermionField &in, FermionField &out,int dir,int disp)=0; // implemented by WilsonFermion and WilsonFermion5D
 
     };
 
diff --git a/lib/qcd/action/fermion/WilsonFermion.cc b/lib/qcd/action/fermion/WilsonFermion.cc
index f5ff7cb0..afd7edd7 100644
--- a/lib/qcd/action/fermion/WilsonFermion.cc
+++ b/lib/qcd/action/fermion/WilsonFermion.cc
@@ -98,7 +98,9 @@ void WilsonFermion::Mdir (const LatticeFermion &in, LatticeFermion &out,int dir,
   DhopDir(in,out,dir,disp);
 }
 void WilsonFermion::DhopDir(const LatticeFermion &in, LatticeFermion &out,int dir,int disp){
+
   WilsonCompressor compressor(DaggerNo);
+
   Stencil.HaloExchange<vSpinColourVector,vHalfSpinColourVector,WilsonCompressor>(in,comm_buf,compressor);
   
   assert( (disp==1)||(disp==-1) );
@@ -109,9 +111,22 @@ void WilsonFermion::DhopDir(const LatticeFermion &in, LatticeFermion &out,int di
 
 PARALLEL_FOR_LOOP
   for(int sss=0;sss<in._grid->oSites();sss++){
-    DiracOptDhopDir(Stencil,Umu,comm_buf,sss,sss,in,out,dirdisp);
+    DiracOptDhopDir(Stencil,Umu,comm_buf,sss,sss,in,out,dirdisp,dirdisp);
   }
+ 
+};
+void WilsonFermion::DhopDirDisp(const LatticeFermion &in, LatticeFermion &out,int dirdisp,int gamma,int dag)
+{
 
+  WilsonCompressor compressor(dag);
+
+  Stencil.HaloExchange<vSpinColourVector,vHalfSpinColourVector,WilsonCompressor>(in,comm_buf,compressor);
+
+PARALLEL_FOR_LOOP
+  for(int sss=0;sss<in._grid->oSites();sss++){
+    DiracOptDhopDir(Stencil,Umu,comm_buf,sss,sss,in,out,dirdisp,gamma);
+  }
+ 
 };
 
 void WilsonFermion::DhopInternal(CartesianStencil & st,LatticeDoubledGaugeField & U,
@@ -177,6 +192,77 @@ void WilsonFermion::Dhop(const LatticeFermion &in, LatticeFermion &out,int dag)
   DhopInternal(Stencil,Umu,in,out,dag);
 }
 
+void WilsonFermion::DerivInternal(CartesianStencil & st,LatticeDoubledGaugeField & U,
+				  LatticeGaugeField &mat,const LatticeFermion &A,const LatticeFermion &B,int dag)
+{
+  assert((dag==DaggerNo) ||(dag==DaggerYes));
+
+  WilsonCompressor compressor(dag);
+
+  LatticeColourMatrix tmp(B._grid);
+  LatticeFermion Btilde(B._grid);
+
+  st.HaloExchange<vSpinColourVector,vHalfSpinColourVector,WilsonCompressor>(B,comm_buf,compressor);
+  
+  for(int mu=0;mu<Nd;mu++){
+      
+    ////////////////////////////////////////////////////////////////////////
+    // Flip gamma (1+g)<->(1-g) if dag
+    ////////////////////////////////////////////////////////////////////////
+    int gamma = mu;
+    if ( dag ) gamma+= Nd;
+
+    ////////////////////////
+    // Call the single hop
+    ////////////////////////
+PARALLEL_FOR_LOOP
+    for(int sss=0;sss<B._grid->oSites();sss++){
+      DiracOptDhopDir(st,U,comm_buf,sss,sss,B,Btilde,mu,gamma);
+    }
+
+    //////////////////////////////////////////////////
+    // spin trace outer product
+    //////////////////////////////////////////////////
+    tmp = TraceIndex<SpinIndex>(outerProduct(Btilde,A)); 
+    PokeIndex<LorentzIndex>(mat,tmp,mu);
+
+  }
+}
+void WilsonFermion::DhopDeriv(LatticeGaugeField &mat,const LatticeFermion &U,const LatticeFermion &V,int dag)
+{
+  conformable(U._grid,_grid);  
+  conformable(U._grid,V._grid);
+  conformable(U._grid,mat._grid);
+
+  mat.checkerboard = U.checkerboard;
+
+  DerivInternal(Stencil,Umu,mat,U,V,dag);
+}
+void WilsonFermion::DhopDerivOE(LatticeGaugeField &mat,const LatticeFermion &U,const LatticeFermion &V,int dag)
+{
+  conformable(U._grid,_cbgrid);  
+  conformable(U._grid,V._grid);
+  conformable(U._grid,mat._grid);
+
+  assert(V.checkerboard==Even);
+  assert(U.checkerboard==Odd);
+  mat.checkerboard = Odd;
+
+  DerivInternal(StencilEven,UmuOdd,mat,U,V,dag);
+}
+void WilsonFermion::DhopDerivEO(LatticeGaugeField &mat,const LatticeFermion &U,const LatticeFermion &V,int dag)
+{
+  conformable(U._grid,_cbgrid);  
+  conformable(U._grid,V._grid);
+  conformable(U._grid,mat._grid);
+
+  assert(V.checkerboard==Odd);
+  assert(U.checkerboard==Even);
+  mat.checkerboard = Even;
+
+  DerivInternal(StencilOdd,UmuEven,mat,U,V,dag);
+}
+
 }}
 
 
diff --git a/lib/qcd/action/fermion/WilsonFermion.h b/lib/qcd/action/fermion/WilsonFermion.h
index 60726fdd..5422ca52 100644
--- a/lib/qcd/action/fermion/WilsonFermion.h
+++ b/lib/qcd/action/fermion/WilsonFermion.h
@@ -24,11 +24,95 @@ namespace Grid {
       // half checkerboard operaions
       void   Meooe       (const LatticeFermion &in, LatticeFermion &out);
       void   MeooeDag    (const LatticeFermion &in, LatticeFermion &out);
+
       virtual void   Mooee       (const LatticeFermion &in, LatticeFermion &out); // remain virtual so we 
       virtual void   MooeeDag    (const LatticeFermion &in, LatticeFermion &out); // can derive Clover
       virtual void   MooeeInv    (const LatticeFermion &in, LatticeFermion &out); // from Wilson base
       virtual void   MooeeInvDag (const LatticeFermion &in, LatticeFermion &out);
 
+      ////////////////////////
+      // 
+      // Force term: d/dtau S = 0
+      //
+      // It is simplest to consider the two flavour force term 
+      // 
+      //       S[U,phi] = phidag (MdagM)^-1 phi
+      //
+      // But simplify even this to
+      //
+      //       S[U,phi] = phidag MdagM phi
+      //
+      // (other options exist depending on nature of action fragment.)
+      // 
+      // Require momentum be traceless anti-hermitian to move within group manifold [ P = i P^a T^a ]
+      //
+      // Define the HMC hamiltonian
+      //
+      //       H = 1/2 Tr P^2 + S(U,phi)
+      // 
+      // .
+      // U =  P U    (lorentz & color indices multiplied)
+      //
+      // Hence
+      //
+      // .c    c  c       c
+      // U =  U  P   = - U  P        (c == dagger)
+      //
+      // So, taking some liberty with implicit indices
+      //                  .      .          .c      c
+      // dH/dt = 0 = Tr P P +Tr[ U  dS/dU + U  dS/dU  ]
+      //
+      //                .                          c      c
+      //           = Tr P P + i Tr[  P U dS/dU  - U   P dS/dU  ]
+      //
+      //                   .                        c  c
+      //           = Tr P (P + i ( U dS/dU - P dS/dU  U ]
+      //
+      //              .                        c  c
+      //           => P  = -i [ U dS/dU - dS/dU  U ]      generates HMC EoM
+      //
+      // Simple case work this out using S = phi^dag MdagM phi for wilson:
+      //                               c       c
+      // dSdt     = dU_xdt  dSdUx  + dUxdt dSdUx
+      //          
+      //         = Tr i P U_x [ (\phi^\dag)_x (1+g) (M \phi)_x+\mu   +(\phi^\dag M^\dag)_x (1-g) \phi_{x+\mu} ]
+      //                 c
+      //            - i U_x P [ (\phi^\dag)_x+mu (1-g) (M \phi)_x    +(\phi^\dag M^\dag)_(x+\mu) (1+g) \phi_{x} ]
+      //
+      //         = i [(\phi^\dag)_x      ]_j P_jk  [U_x(1+g) (M \phi)_x+\mu]_k        (1)
+      //         + i [(\phi^\dagM^\dag)_x]_j P_jk  [U_x(1-g) (\phi)_x+\mu]_k          (2)
+      //         - i [(\phi^\dag)_x+mu (1-g) U^dag_x]_j P_jk  [(M \phi)_xk            (3)
+      //         - i [(\phi^\dagM^\dag)_x+mu (1+g) U^dag_x]_j P_jk  [ \phi]_xk        (4)
+      //
+      // Observe that (1)* = (4)
+      //              (2)* = (3)
+      //
+      // Write as    .
+      //             P_{kj}  = - i (  [U_x(1+g) (M \phi)_x+\mu] (x) [(\phi^\dag)_x] + [U_x(1-g) (\phi)_x+\mu] (x) [(\phi^\dagM^\dag)_x] - h.c )
+      //
+      // where (x) denotes outer product in colour and spins are traced.
+      //
+      // Need only evaluate (1) and (2)  [Chroma] or (2) and (4) [IroIro] and take the 
+      // traceless anti hermitian part (of term in brackets w/o the "i")
+      // 
+      // Generalisation to S=phi^dag (MdagM)^{-1} phi is simple:
+      //
+      // For more complicated DWF etc... apply product rule in differentiation
+      //
+      ////////////////////////
+      void DhopDeriv  (LatticeGaugeField &mat,const LatticeFermion &U,const LatticeFermion &V,int dag);
+      void DhopDerivEO(LatticeGaugeField &mat,const LatticeFermion &U,const LatticeFermion &V,int dag);
+      void DhopDerivOE(LatticeGaugeField &mat,const LatticeFermion &U,const LatticeFermion &V,int dag);
+
+      // Extra support internal
+      void DerivInternal(CartesianStencil & st,
+			 LatticeDoubledGaugeField & U,
+			 LatticeGaugeField &mat,
+			 const LatticeFermion &A,
+			 const LatticeFermion &B,
+			 int dag);
+
+
       // non-hermitian hopping term; half cb or both
       void Dhop  (const LatticeFermion &in, LatticeFermion &out,int dag);
       void DhopOE(const LatticeFermion &in, LatticeFermion &out,int dag);
@@ -37,6 +121,7 @@ namespace Grid {
       // Multigrid assistance
       void   Mdir (const LatticeFermion &in, LatticeFermion &out,int dir,int disp);
       void DhopDir(const LatticeFermion &in, LatticeFermion &out,int dir,int disp);
+      void DhopDirDisp(const LatticeFermion &in, LatticeFermion &out,int dirdisp,int gamma,int dag);
 
       ///////////////////////////////////////////////////////////////
       // Extra methods added by derived
@@ -59,7 +144,8 @@ namespace Grid {
       static int HandOptDslash; // these are a temporary hack
       static int MortonOrder;
 
-    protected:
+      //    protected:
+    public:
 
       RealD                        mass;
 
diff --git a/lib/qcd/action/fermion/WilsonFermion5D.cc b/lib/qcd/action/fermion/WilsonFermion5D.cc
index 43829ea7..523934c4 100644
--- a/lib/qcd/action/fermion/WilsonFermion5D.cc
+++ b/lib/qcd/action/fermion/WilsonFermion5D.cc
@@ -105,14 +105,105 @@ PARALLEL_FOR_LOOP
     for(int s=0;s<Ls;s++){
       int sU=ss;
       int sF = s+Ls*sU; 
-      DiracOptDhopDir(Stencil,Umu,comm_buf,sF,sU,in,out,dirdisp);
+      DiracOptDhopDir(Stencil,Umu,comm_buf,sF,sU,in,out,dirdisp,dirdisp);
     }
   }
 };
 
+void WilsonFermion5D::DerivInternal(CartesianStencil & st,
+				    LatticeDoubledGaugeField & U,
+				    LatticeGaugeField &mat,
+				    const LatticeFermion &A,
+				    const LatticeFermion &B,
+				    int dag)
+{
+  assert((dag==DaggerNo) ||(dag==DaggerYes));
+
+  WilsonCompressor compressor(dag);
+  
+  LatticeColourMatrix tmp(B._grid);
+  LatticeFermion Btilde(B._grid);
+
+  st.HaloExchange<vSpinColourVector,vHalfSpinColourVector,WilsonCompressor>(B,comm_buf,compressor);
+  
+  for(int mu=0;mu<Nd;mu++){
+      
+    ////////////////////////////////////////////////////////////////////////
+    // Flip gamma if dag
+    ////////////////////////////////////////////////////////////////////////
+    int gamma = mu;
+    if ( dag ) gamma+= Nd;
+
+    ////////////////////////
+    // Call the single hop
+    ////////////////////////
+PARALLEL_FOR_LOOP
+    for(int sss=0;sss<B._grid->oSites();sss++){
+      for(int s=0;s<Ls;s++){
+	int sU=sss;
+	int sF = s+Ls*sU;
+	DiracOptDhopDir(st,U,comm_buf,sF,sU,B,Btilde,mu,gamma);
+      }
+    }
+
+    ////////////////////////////
+    // spin trace outer product
+    ////////////////////////////
+    tmp = TraceIndex<SpinIndex>(outerProduct(Btilde,A)); // ordering here
+    PokeIndex<LorentzIndex>(mat,tmp,mu);
+
+  }
+}
+
+void WilsonFermion5D::DhopDeriv(      LatticeGaugeField &mat,
+				const LatticeFermion &A,
+				const LatticeFermion &B,
+				int dag)
+{
+  conformable(A._grid,FermionGrid());  
+  conformable(A._grid,B._grid);
+  conformable(GaugeGrid(),mat._grid);
+
+  mat.checkerboard = A.checkerboard;
+
+  DerivInternal(Stencil,Umu,mat,A,B,dag);
+}
+
+void WilsonFermion5D::DhopDerivEO(LatticeGaugeField &mat,
+				  const LatticeFermion &A,
+				  const LatticeFermion &B,
+				  int dag)
+{
+  conformable(A._grid,FermionRedBlackGrid());
+  conformable(GaugeRedBlackGrid(),mat._grid);
+  conformable(A._grid,B._grid);
+
+  assert(B.checkerboard==Odd);
+  assert(A.checkerboard==Even);
+  mat.checkerboard = Even;
+
+  DerivInternal(StencilOdd,UmuEven,mat,A,B,dag);
+}
+
+void WilsonFermion5D::DhopDerivOE(LatticeGaugeField &mat,
+				  const LatticeFermion &A,
+				  const LatticeFermion &B,
+				  int dag)
+{
+  conformable(A._grid,FermionRedBlackGrid());
+  conformable(GaugeRedBlackGrid(),mat._grid);
+  conformable(A._grid,B._grid);
+
+  assert(B.checkerboard==Even);
+  assert(A.checkerboard==Odd);
+  mat.checkerboard = Odd;
+
+  DerivInternal(StencilEven,UmuOdd,mat,A,B,dag);
+}
+
 void WilsonFermion5D::DhopInternal(CartesianStencil & st, LebesgueOrder &lo,
 				   LatticeDoubledGaugeField & U,
-			   const LatticeFermion &in, LatticeFermion &out,int dag)
+				   const LatticeFermion &in, LatticeFermion &out,int dag)
 {
   //  assert((dag==DaggerNo) ||(dag==DaggerYes));
 
diff --git a/lib/qcd/action/fermion/WilsonFermion5D.h b/lib/qcd/action/fermion/WilsonFermion5D.h
index e345d2a4..7ff5d1c3 100644
--- a/lib/qcd/action/fermion/WilsonFermion5D.h
+++ b/lib/qcd/action/fermion/WilsonFermion5D.h
@@ -44,12 +44,18 @@ namespace Grid {
 
       // half checkerboard operations; leave unimplemented as abstract for now
       virtual void   Meooe       (const LatticeFermion &in, LatticeFermion &out){assert(0);};
-      virtual void   MeooeDag    (const LatticeFermion &in, LatticeFermion &out){assert(0);};
       virtual void   Mooee       (const LatticeFermion &in, LatticeFermion &out){assert(0);};
-      virtual void   MooeeDag    (const LatticeFermion &in, LatticeFermion &out){assert(0);};
       virtual void   MooeeInv    (const LatticeFermion &in, LatticeFermion &out){assert(0);};
+
+      virtual void   MeooeDag    (const LatticeFermion &in, LatticeFermion &out){assert(0);};
+      virtual void   MooeeDag    (const LatticeFermion &in, LatticeFermion &out){assert(0);};
       virtual void   MooeeInvDag (const LatticeFermion &in, LatticeFermion &out){assert(0);};
 
+      // These can be overridden by fancy 5d chiral actions
+      virtual void DhopDeriv  (LatticeGaugeField &mat,const LatticeFermion &U,const LatticeFermion &V,int dag);
+      virtual void DhopDerivEO(LatticeGaugeField &mat,const LatticeFermion &U,const LatticeFermion &V,int dag);
+      virtual void DhopDerivOE(LatticeGaugeField &mat,const LatticeFermion &U,const LatticeFermion &V,int dag);
+
       // Implement hopping term non-hermitian hopping term; half cb or both
       // Implement s-diagonal DW
       void DW    (const LatticeFermion &in, LatticeFermion &out,int dag);
@@ -64,6 +70,14 @@ namespace Grid {
       ///////////////////////////////////////////////////////////////
       // New methods added 
       ///////////////////////////////////////////////////////////////
+
+      void DerivInternal(CartesianStencil & st,
+			 LatticeDoubledGaugeField & U,
+			 LatticeGaugeField &mat,
+			 const LatticeFermion &A,
+			 const LatticeFermion &B,
+			 int dag);
+
       void DhopInternal(CartesianStencil & st,
 			LebesgueOrder &lo,
 			LatticeDoubledGaugeField &U,
diff --git a/lib/qcd/action/fermion/WilsonKernels.cc b/lib/qcd/action/fermion/WilsonKernels.cc
index bdc18eee..43ccff63 100644
--- a/lib/qcd/action/fermion/WilsonKernels.cc
+++ b/lib/qcd/action/fermion/WilsonKernels.cc
@@ -294,8 +294,8 @@ void DiracOptDhopSiteDag(CartesianStencil &st,LatticeDoubledGaugeField &U,
 }
 
 void DiracOptDhopDir(CartesianStencil &st,LatticeDoubledGaugeField &U,
-			std::vector<vHalfSpinColourVector,alignedAllocator<vHalfSpinColourVector> >  &buf,
-		       int sF,int sU,const LatticeFermion &in, LatticeFermion &out,int dirdisp)
+		     std::vector<vHalfSpinColourVector,alignedAllocator<vHalfSpinColourVector> >  &buf,
+		     int sF,int sU,const LatticeFermion &in, LatticeFermion &out,int dir,int gamma)
 {
     vHalfSpinColourVector  tmp;    
     vHalfSpinColourVector  chi;    
@@ -304,13 +304,13 @@ void DiracOptDhopDir(CartesianStencil &st,LatticeDoubledGaugeField &U,
     int offset,local,perm, ptype;
     int ss=sF;
     
-    offset = st._offsets [dirdisp][ss];
-    local  = st._is_local[dirdisp][ss];
-    perm   = st._permute[dirdisp][ss];
-    ptype  = st._permute_type[dirdisp];
+    offset = st._offsets [dir][ss];
+    local  = st._is_local[dir][ss];
+    perm   = st._permute[dir][ss];
+    ptype  = st._permute_type[dir];
 
     // Xp
-    if(dirdisp==Xp){
+    if(gamma==Xp){
       if ( local && perm ) {
 	spProjXp(tmp,in._odata[offset]);
 	permute(chi,tmp,ptype);
@@ -319,12 +319,12 @@ void DiracOptDhopDir(CartesianStencil &st,LatticeDoubledGaugeField &U,
       } else { 
 	chi=buf[offset];
       }
-      mult(&Uchi(),&U._odata[sU](Xp),&chi());
+      mult(&Uchi(),&U._odata[sU](dir),&chi());
       spReconXp(result,Uchi);
     }
 
     // Yp
-    if ( dirdisp==Yp ){
+    if ( gamma==Yp ){
       if ( local && perm ) {
 	spProjYp(tmp,in._odata[offset]);
 	permute(chi,tmp,ptype);
@@ -333,12 +333,12 @@ void DiracOptDhopDir(CartesianStencil &st,LatticeDoubledGaugeField &U,
       } else { 
 	chi=buf[offset];
       }
-      mult(&Uchi(),&U._odata[sU](Yp),&chi());
+      mult(&Uchi(),&U._odata[sU](dir),&chi());
       spReconYp(result,Uchi);
     }
 
     // Zp
-    if ( dirdisp ==Zp ){
+    if ( gamma ==Zp ){
       if ( local && perm ) {
 	spProjZp(tmp,in._odata[offset]);
 	permute(chi,tmp,ptype);
@@ -347,12 +347,12 @@ void DiracOptDhopDir(CartesianStencil &st,LatticeDoubledGaugeField &U,
       } else { 
 	chi=buf[offset];
       }
-      mult(&Uchi(),&U._odata[sU](Zp),&chi());
+      mult(&Uchi(),&U._odata[sU](dir),&chi());
       spReconZp(result,Uchi);
     }
 
     // Tp
-    if ( dirdisp ==Tp ){
+    if ( gamma ==Tp ){
       if ( local && perm ) {
 	spProjTp(tmp,in._odata[offset]);
 	permute(chi,tmp,ptype);
@@ -361,12 +361,12 @@ void DiracOptDhopDir(CartesianStencil &st,LatticeDoubledGaugeField &U,
       } else { 
 	chi=buf[offset];
       }
-      mult(&Uchi(),&U._odata[sU](Tp),&chi());
+      mult(&Uchi(),&U._odata[sU](dir),&chi());
       spReconTp(result,Uchi);
     }
 
     // Xm
-    if ( dirdisp==Xm ){
+    if ( gamma==Xm ){
       if ( local && perm ) {
 	spProjXm(tmp,in._odata[offset]);
 	permute(chi,tmp,ptype);
@@ -375,12 +375,12 @@ void DiracOptDhopDir(CartesianStencil &st,LatticeDoubledGaugeField &U,
       } else { 
 	chi=buf[offset];
       }
-      mult(&Uchi(),&U._odata[sU](Xm),&chi());
+      mult(&Uchi(),&U._odata[sU](dir),&chi());
       spReconXm(result,Uchi);
     }
 
     // Ym
-    if ( dirdisp == Ym ){
+    if ( gamma == Ym ){
       if ( local && perm ) {
 	spProjYm(tmp,in._odata[offset]);
 	permute(chi,tmp,ptype);
@@ -389,12 +389,12 @@ void DiracOptDhopDir(CartesianStencil &st,LatticeDoubledGaugeField &U,
       } else { 
 	chi=buf[offset];
       }
-      mult(&Uchi(),&U._odata[sU](Ym),&chi());
+      mult(&Uchi(),&U._odata[sU](dir),&chi());
       spReconYm(result,Uchi);
     }
 
     // Zm
-    if ( dirdisp == Zm ){
+    if ( gamma == Zm ){
       if ( local && perm ) {
 	spProjZm(tmp,in._odata[offset]);
 	permute(chi,tmp,ptype);
@@ -403,12 +403,12 @@ void DiracOptDhopDir(CartesianStencil &st,LatticeDoubledGaugeField &U,
       } else { 
 	chi=buf[offset];
       }
-      mult(&Uchi(),&U._odata[sU](Zm),&chi());
+      mult(&Uchi(),&U._odata[sU](dir),&chi());
       spReconZm(result,Uchi);
     }
 
     // Tm
-    if ( dirdisp==Tm ) {
+    if ( gamma==Tm ) {
       if ( local && perm ) {
 	spProjTm(tmp,in._odata[offset]);
 	permute(chi,tmp,ptype);
@@ -417,7 +417,7 @@ void DiracOptDhopDir(CartesianStencil &st,LatticeDoubledGaugeField &U,
       } else { 
 	chi=buf[offset];
       }
-      mult(&Uchi(),&U._odata[sU](Tm),&chi());
+      mult(&Uchi(),&U._odata[sU](dir),&chi());
       spReconTm(result,Uchi);
     }
 
diff --git a/lib/qcd/action/fermion/WilsonKernels.h b/lib/qcd/action/fermion/WilsonKernels.h
index ec807af5..cb0e07bb 100644
--- a/lib/qcd/action/fermion/WilsonKernels.h
+++ b/lib/qcd/action/fermion/WilsonKernels.h
@@ -22,7 +22,7 @@ namespace Grid {
 			       int sF,int sU,const LatticeFermion &in, LatticeFermion &out);
       void DiracOptDhopDir(CartesianStencil &st,LatticeDoubledGaugeField &U,
 			   std::vector<vHalfSpinColourVector,alignedAllocator<vHalfSpinColourVector> >  &buf,
-			   int sF,int sU,const LatticeFermion &in, LatticeFermion &out,int dirdisp);
+			   int sF,int sU,const LatticeFermion &in, LatticeFermion &out,int dirdisp,int gamma);
       
       //  };
 
diff --git a/lib/qcd/hmc/HMC.h b/lib/qcd/hmc/HMC.h
index e3fe57ed..21c6bb0d 100644
--- a/lib/qcd/hmc/HMC.h
+++ b/lib/qcd/hmc/HMC.h
@@ -68,14 +68,12 @@ namespace Grid{
 		     Integrator<Algorithm>& MolDyn):
       Params(Pms),MD(MolDyn){
 	//FIXME
-
 	// initialize RNGs also with seed
 	sRNG.SeedRandomDevice();
       }
       ~HybridMonteCarlo(){};
 
 
-
       void evolve(LatticeLorentzColourMatrix& Uin){
 	Real DeltaH;
 	
diff --git a/lib/simd/Grid_vector_types.h b/lib/simd/Grid_vector_types.h
index 034d5314..71472407 100644
--- a/lib/simd/Grid_vector_types.h
+++ b/lib/simd/Grid_vector_types.h
@@ -414,7 +414,7 @@ namespace Grid {
   template<class S, class V >
   inline Grid_simd< S, V> outerProduct(const Grid_simd< S, V> &l, const Grid_simd< S, V> & r)
   {
-    return l*r;
+    return l*conjugate(r);
   }
 
   template<class S, class V >
diff --git a/lib/tensors/Tensor_outer.h b/lib/tensors/Tensor_outer.h
index 816f6ad8..d0150709 100644
--- a/lib/tensors/Tensor_outer.h
+++ b/lib/tensors/Tensor_outer.h
@@ -28,11 +28,13 @@ auto outerProduct (const iScalar<l>& lhs,const iScalar<r>& rhs) -> iScalar<declt
 
 inline ComplexF outerProduct(const ComplexF &l, const ComplexF& r)
 {
-  return l*r;
+  std::cout << "outer product taking conj "<<r<<" "<<conj(r)<<std::endl;
+  return l*conj(r);
 }
 inline ComplexD outerProduct(const ComplexD &l, const ComplexD& r)
 {
-  return l*r;
+  std::cout << "outer product taking conj "<<r<<" "<<conj(r)<<std::endl;
+  return l*conj(r);
 }
 inline RealF outerProduct(const RealF &l, const RealF& r)
 {
diff --git a/lib/tensors/Tensor_trace.h b/lib/tensors/Tensor_trace.h
index 4cf5edc2..bd75bb91 100644
--- a/lib/tensors/Tensor_trace.h
+++ b/lib/tensors/Tensor_trace.h
@@ -31,6 +31,16 @@ inline auto trace(const iScalar<vtype> &arg) -> iScalar<decltype(trace(arg._inte
     return ret;
 }
 
+template<class vtype,int N>
+  inline auto trace(const iVector<vtype,N> &arg) -> iVector<decltype(trace(arg._internal[0])),N>
+{
+  iVector<decltype(trace(arg._internal[0])),N> ret;
+  for(int i=0;i<N;i++){
+    ret._internal[i]=trace(arg._internal[i]);
+  }
+  return ret;
+}
+
 
 }
 #endif
diff --git a/tests/Make.inc b/tests/Make.inc
index f04da4eb..b78b4d16 100644
--- a/tests/Make.inc
+++ b/tests/Make.inc
@@ -1,5 +1,5 @@
 
-bin_PROGRAMS = Test_GaugeAction Test_cayley_cg Test_cayley_coarsen_support Test_cayley_even_odd Test_cayley_ldop_cr Test_cf_coarsen_support Test_cf_cr_unprec Test_contfrac_cg Test_contfrac_even_odd Test_cshift Test_cshift_red_black Test_dwf_cg_prec Test_dwf_cg_schur Test_dwf_cg_unprec Test_dwf_cr_unprec Test_dwf_even_odd Test_dwf_fpgcr Test_dwf_hdcr Test_gamma Test_hmc_WilsonGauge Test_lie_generators Test_main Test_multishift_sqrt Test_nersc_io Test_quenched_update Test_remez Test_rng Test_rng_fixed Test_simd Test_stencil Test_wilson_cg_prec Test_wilson_cg_schur Test_wilson_cg_unprec Test_wilson_cr_unprec Test_wilson_even_odd
+bin_PROGRAMS = Test_GaugeAction Test_cayley_cg Test_cayley_coarsen_support Test_cayley_even_odd Test_cayley_ldop_cr Test_cf_coarsen_support Test_cf_cr_unprec Test_contfrac_cg Test_contfrac_even_odd Test_cshift Test_cshift_red_black Test_dwf_cg_prec Test_dwf_cg_schur Test_dwf_cg_unprec Test_dwf_cr_unprec Test_dwf_even_odd Test_dwf_fpgcr Test_dwf_hdcr Test_gamma Test_hmc_WilsonGauge Test_lie_generators Test_main Test_multishift_sqrt Test_nersc_io Test_quenched_update Test_remez Test_rng Test_rng_fixed Test_simd Test_stencil Test_wilson_cg_prec Test_wilson_cg_schur Test_wilson_cg_unprec Test_wilson_cr_unprec Test_wilson_even_odd Test_wilson_force Test_wilson_force_phiMphi
 
 
 Test_GaugeAction_SOURCES=Test_GaugeAction.cc
@@ -141,3 +141,11 @@ Test_wilson_cr_unprec_LDADD=-lGrid
 Test_wilson_even_odd_SOURCES=Test_wilson_even_odd.cc
 Test_wilson_even_odd_LDADD=-lGrid
 
+
+Test_wilson_force_SOURCES=Test_wilson_force.cc
+Test_wilson_force_LDADD=-lGrid
+
+
+Test_wilson_force_phiMphi_SOURCES=Test_wilson_force_phiMphi.cc
+Test_wilson_force_phiMphi_LDADD=-lGrid
+
diff --git a/tests/Test_dwf_cg_prec.cc b/tests/Test_dwf_cg_prec.cc
index d429ee92..633d280f 100644
--- a/tests/Test_dwf_cg_prec.cc
+++ b/tests/Test_dwf_cg_prec.cc
@@ -34,7 +34,9 @@ int main (int argc, char ** argv)
 
   LatticeFermion    src(FGrid); random(RNG5,src);
   LatticeFermion result(FGrid); result=zero;
-  LatticeGaugeField Umu(UGrid); random(RNG4,Umu);
+  LatticeGaugeField Umu(UGrid); 
+
+  SU3::HotConfiguration(RNG4,Umu);
 
   std::vector<LatticeColourMatrix> U(4,UGrid);
   for(int mu=0;mu<Nd;mu++){
diff --git a/tests/Test_wilson_force.cc b/tests/Test_wilson_force.cc
new file mode 100644
index 00000000..6a8d7611
--- /dev/null
+++ b/tests/Test_wilson_force.cc
@@ -0,0 +1,102 @@
+#include <Grid.h>
+
+using namespace std;
+using namespace Grid;
+using namespace Grid::QCD;
+
+#define parallel_for PARALLEL_FOR_LOOP for
+
+int main (int argc, char ** argv)
+{
+  Grid_init(&argc,&argv);
+
+  std::vector<int> latt_size   = GridDefaultLatt();
+  std::vector<int> simd_layout = GridDefaultSimd(Nd,vComplex::Nsimd());
+  std::vector<int> mpi_layout  = GridDefaultMpi();
+
+  GridCartesian               Grid(latt_size,simd_layout,mpi_layout);
+  GridRedBlackCartesian     RBGrid(latt_size,simd_layout,mpi_layout);
+
+  int threads = GridThread::GetThreads();
+  std::cout<<GridLogMessage << "Grid is setup to use "<<threads<<" threads"<<std::endl;
+
+  std::vector<int> seeds({1,2,3,4});
+
+  GridParallelRNG          pRNG(&Grid);
+  pRNG.SeedRandomDevice();
+
+  LatticeFermion phi        (&Grid); gaussian(pRNG,phi);
+  LatticeFermion Mphi       (&Grid); 
+  LatticeFermion MphiPrime  (&Grid); 
+
+  LatticeGaugeField U(&Grid);
+
+  SU3::HotConfiguration(pRNG,U);
+  //  SU3::ColdConfiguration(pRNG,U);
+  
+  ////////////////////////////////////
+  // Unmodified matrix element
+  ////////////////////////////////////
+  RealD mass=-4.0; //kills the diagonal term
+  WilsonFermion Dw     (U,     Grid,RBGrid,mass);
+  Dw.M   (phi,Mphi);
+
+  ComplexD S    = innerProduct(Mphi,Mphi); // pdag MdagM p
+
+  // get the deriv of phidag MdagM phi with respect to "U"
+  LatticeGaugeField UdSdU(&Grid);
+  LatticeGaugeField tmp(&Grid);
+
+  Dw.MDeriv(tmp , Mphi,  phi,DaggerNo );  UdSdU=tmp;
+  Dw.MDeriv(tmp , phi,  Mphi,DaggerYes ); UdSdU=UdSdU+tmp;
+
+
+  LatticeFermion Ftmp      (&Grid);
+
+  ////////////////////////////////////
+  // Modify the gauge field a little 
+  ////////////////////////////////////
+  RealD dt = 1.0e-6;
+
+  LatticeColourMatrix mommu(&Grid); 
+  LatticeGaugeField mom(&Grid); 
+  LatticeGaugeField Uprime(&Grid); 
+
+  for(int mu=0;mu<Nd;mu++){
+
+    SU3::GaussianLieAlgebraMatrix(pRNG, mommu); // Traceless antihermitian momentum; gaussian in lie alg
+
+    PokeIndex<LorentzIndex>(mom,mommu,mu);
+    parallel_for(auto i=mom.begin();i<mom.end();i++){
+      Uprime[i](mu) =U[i](mu)+ mom[i](mu)*U[i](mu)*dt;
+    }
+
+  }
+
+  Dw.DoubleStore(Dw.Umu,Uprime);
+  Dw.M          (phi,MphiPrime);
+
+  ComplexD Sprime    = innerProduct(MphiPrime   ,MphiPrime);
+
+  //////////////////////////////////////////////
+  // Use derivative to estimate dS
+  //////////////////////////////////////////////
+  LatticeComplex dS(&Grid); dS = zero;
+
+  parallel_for(auto i=mom.begin();i<mom.end();i++){
+    for(int mu=0;mu<Nd;mu++){
+      //      dS[i]() = dS[i]()+trace(mom[i](mu) * UdSdU[i](mu) - mom[i](mu)* adj( UdSdU[i](mu)) )*dt;
+      dS[i]() =    dS[i]()+trace(mom[i](mu) * (UdSdU[i](mu)))*dt;
+      dS[i]() =    dS[i]()-trace(mom[i](mu) * adj(UdSdU[i](mu)))*dt;
+    }
+  }
+  Complex dSpred    = sum(dS);
+
+  std::cout << GridLogMessage << " S      "<<S<<std::endl;
+  std::cout << GridLogMessage << " Sprime "<<Sprime<<std::endl;
+  std::cout << GridLogMessage << "dS      "<<Sprime-S<<std::endl;
+  std::cout << GridLogMessage << "predict dS    "<< dSpred <<std::endl;
+
+  std::cout<< GridLogMessage << "Done" <<std::endl;
+  Grid_finalize();
+}
diff --git a/tests/Test_wilson_force_phiMdagMphi.cc b/tests/Test_wilson_force_phiMdagMphi.cc
new file mode 100644
index 00000000..33c5f6e1
--- /dev/null
+++ b/tests/Test_wilson_force_phiMdagMphi.cc
@@ -0,0 +1,140 @@
+#include <Grid.h>
+
+using namespace std;
+using namespace Grid;
+using namespace Grid::QCD;
+
+#define parallel_for PARALLEL_FOR_LOOP for
+
+int main (int argc, char ** argv)
+{
+  Grid_init(&argc,&argv);
+
+  std::vector<int> latt_size   = GridDefaultLatt();
+  std::vector<int> simd_layout = GridDefaultSimd(Nd,vComplex::Nsimd());
+  std::vector<int> mpi_layout  = GridDefaultMpi();
+
+  GridCartesian               Grid(latt_size,simd_layout,mpi_layout);
+  GridRedBlackCartesian     RBGrid(latt_size,simd_layout,mpi_layout);
+
+  int threads = GridThread::GetThreads();
+  std::cout<<GridLogMessage << "Grid is setup to use "<<threads<<" threads"<<std::endl;
+
+  std::vector<int> seeds({1,2,3,4});
+
+  GridParallelRNG          pRNG(&Grid);
+  pRNG.SeedRandomDevice();
+
+  LatticeFermion phi        (&Grid); gaussian(pRNG,phi);
+  LatticeFermion Mphi       (&Grid); 
+  LatticeFermion Mdagphi       (&Grid); 
+  LatticeFermion MphiPrime  (&Grid); 
+  LatticeFermion MdagphiPrime  (&Grid); 
+  LatticeFermion dMphi      (&Grid); 
+
+  LatticeGaugeField U(&Grid);
+
+
+  SU3::HotConfiguration(pRNG,U);
+  //  SU3::ColdConfiguration(pRNG,U);
+  
+  ////////////////////////////////////
+  // Unmodified matrix element
+  ////////////////////////////////////
+  RealD mass=-4.0; //kills the diagonal term
+  WilsonFermion Dw     (U,     Grid,RBGrid,mass);
+  Dw.M   (phi,Mphi);
+  Dw.Mdag(phi,Mdagphi);
+
+  ComplexD S    = innerProduct(Mphi,Mphi); // pdag MdagM p
+  ComplexD Sdag = innerProduct(Mdagphi,Mdagphi); // pdag MMdag p
+
+  // get the deriv of phidag MdagM phi with respect to "U"
+  LatticeGaugeField UdSdU(&Grid);
+  LatticeGaugeField UdSdUdag(&Grid);
+  LatticeGaugeField tmp(&Grid);
+
+  Dw.MDeriv(tmp , Mphi,  phi,DaggerNo );  UdSdU=tmp;
+
+  Dw.MDeriv(tmp , Mdagphi,  phi,DaggerYes );  UdSdUdag=tmp;
+
+
+  LatticeFermion dMdagphi      (&Grid);  dMdagphi=zero;
+  LatticeFermion Ftmp      (&Grid);
+
+
+  //  Dw.MDeriv(UdSdU,Mdagphi,  phi,DaggerYes );// UdSdU =UdSdU +tmp;
+
+  ////////////////////////////////////
+  // Modify the gauge field a little in one dir
+  ////////////////////////////////////
+  RealD dt = 1.0e-3;
+
+  LatticeColourMatrix mommu(&Grid); 
+  LatticeGaugeField mom(&Grid); 
+  LatticeGaugeField Uprime(&Grid); 
+
+  for(int mu=0;mu<Nd;mu++){
+
+    SU3::GaussianLieAlgebraMatrix(pRNG, mommu); // Traceless antihermitian momentum; gaussian in lie alg
+
+    //    Dw.DoubleStore(Dw.Umu,Uprime); // update U _and_ Udag
+    Dw.DhopDirDisp(phi,Ftmp,mu,mu+4,DaggerYes); 
+    dMdagphi=dMdagphi+mommu*Ftmp*dt;
+    
+    PokeIndex<LorentzIndex>(mom,mommu,mu);
+
+    parallel_for(auto i=mom.begin();i<mom.end();i++){
+      Uprime[i](mu) =U[i](mu)+ mom[i](mu)*U[i](mu)*dt;
+      Dw.Umu[i](mu) =Uprime[i](mu); // update U but _not_ Udag
+    }
+
+  }
+
+  Dw.Mdag(phi,MdagphiPrime);
+  Dw.M   (phi,MphiPrime);
+
+  std::cout << GridLogMessage << "deltaMdag phi    "<< norm2(dMdagphi) <<std::endl;
+  Ftmp=MdagphiPrime - Mdagphi;
+  std::cout << GridLogMessage << "diff Mdag phi    "<< norm2(Ftmp) <<std::endl;
+  Ftmp = Ftmp - dMdagphi;
+  std::cout << GridLogMessage << "err  Mdag phi    "<< norm2(Ftmp) <<std::endl;
+  std::cout << dMdagphi<<std::endl;
+  Ftmp=MdagphiPrime - Mdagphi;
+  std::cout << Ftmp<<std::endl;
+
+
+  ComplexD Sprime    = innerProduct(Mphi   ,MphiPrime);
+  ComplexD Sprimedag = innerProduct(Mdagphi,MdagphiPrime);
+
+  ComplexD deltaSdag = innerProduct(Mdagphi,dMdagphi);
+  std::cout << GridLogMessage << "deltaSdag from inner prod of mom* M[u]     "<<deltaSdag<<std::endl;
+
+  //////////////////////////////////////////////
+  // Use derivative to estimate dS
+  //////////////////////////////////////////////
+  LatticeComplex dS(&Grid); dS = zero;
+  LatticeComplex dSdag(&Grid); dSdag = zero;
+  parallel_for(auto i=mom.begin();i<mom.end();i++){
+    for(int mu=0;mu<Nd;mu++){
+      //      dS[i]() = dS[i]()+trace(mom[i](mu) * UdSdU[i](mu) - mom[i](mu)* adj( UdSdU[i](mu)) )*dt;
+      dS[i]()    =    dS[i]()+trace(mom[i](mu) * UdSdU[i](mu) )*dt;
+      dSdag[i]() = dSdag[i]()+trace(mom[i](mu) * UdSdUdag[i](mu) )*dt;
+    }
+  }
+  Complex dSpred    = sum(dS);
+  Complex dSdagpred = sum(dSdag);
+
+  std::cout << GridLogMessage << " S      "<<S<<std::endl;
+  std::cout << GridLogMessage << " Sprime "<<Sprime<<std::endl;
+  std::cout << GridLogMessage << "dS      "<<Sprime-S<<std::endl;
+  std::cout << GridLogMessage << "predict dS    "<< dSpred <<std::endl;
+  std::cout << "\n\n"<<std::endl;
+  std::cout << GridLogMessage << " Sdag      "<<Sdag<<std::endl;
+  std::cout << GridLogMessage << " Sprimedag "<<Sprimedag<<std::endl;
+  std::cout << GridLogMessage << "dSdag      "<<Sprimedag-Sdag<<std::endl;
+  std::cout << GridLogMessage << "predict dSdag    "<< dSdagpred <<std::endl;
+
+  std::cout<< GridLogMessage << "Done" <<std::endl;
+  Grid_finalize();
+}
diff --git a/tests/Test_wilson_force_phiMphi.cc b/tests/Test_wilson_force_phiMphi.cc
new file mode 100644
index 00000000..f1e42195
--- /dev/null
+++ b/tests/Test_wilson_force_phiMphi.cc
@@ -0,0 +1,162 @@
+#include <Grid.h>
+
+using namespace std;
+using namespace Grid;
+using namespace Grid::QCD;
+
+#define parallel_for PARALLEL_FOR_LOOP for
+
+int main (int argc, char ** argv)
+{
+  Grid_init(&argc,&argv);
+
+  std::vector<int> latt_size   = GridDefaultLatt();
+  std::vector<int> simd_layout = GridDefaultSimd(Nd,vComplex::Nsimd());
+  std::vector<int> mpi_layout  = GridDefaultMpi();
+
+  GridCartesian               Grid(latt_size,simd_layout,mpi_layout);
+  GridRedBlackCartesian     RBGrid(latt_size,simd_layout,mpi_layout);
+
+  int threads = GridThread::GetThreads();
+  std::cout<<GridLogMessage << "Grid is setup to use "<<threads<<" threads"<<std::endl;
+
+  std::vector<int> seeds({1,2,3,4});
+
+  GridParallelRNG          pRNG(&Grid);
+  pRNG.SeedRandomDevice();
+
+  LatticeFermion phi        (&Grid); gaussian(pRNG,phi);
+  LatticeFermion Mphi       (&Grid); 
+  LatticeFermion MphiPrime  (&Grid); 
+  LatticeFermion dMphi      (&Grid); 
+
+  LatticeGaugeField U(&Grid);
+
+
+  SU3::HotConfiguration(pRNG,U);
+
+  ////////////////////////////////////
+  // Unmodified matrix element
+  ////////////////////////////////////
+  RealD mass=-4.0; //kills the diagonal term
+  WilsonFermion Dw     (U,     Grid,RBGrid,mass);
+  Dw.M(phi,Mphi);
+
+  ComplexD S = innerProduct(phi,Mphi);
+
+  // get the deriv
+  LatticeGaugeField UdSdU(&Grid);
+  Dw.MDeriv(UdSdU,phi, phi,DaggerNo ); 
+
+
+  ////////////////////////////////////
+  // Modify the gauge field a little in one dir
+  ////////////////////////////////////
+  RealD dt = 1.0e-3;
+  Complex Complex_i(0,1);
+
+  LatticeColourMatrix Umu(&Grid);
+  LatticeColourMatrix Umu_save(&Grid);
+  LatticeColourMatrix dU (&Grid);
+  LatticeColourMatrix mom(&Grid); 
+  SU3::GaussianLieAlgebraMatrix(pRNG, mom); // Traceless antihermitian momentum; gaussian in lie alg
+
+
+  // check mom is as i expect
+  LatticeColourMatrix tmpmom(&Grid); 
+  tmpmom = mom+adj(mom);
+  std::cout << GridLogMessage << "mom anti-herm check "<< norm2(tmpmom)<<std::endl;
+  std::cout << GridLogMessage << "mom tr check "<< norm2(trace(mom))<<std::endl;
+  
+  const int mu=0;
+  Umu = PeekIndex<LorentzIndex>(U,mu);
+  Umu_save=Umu;
+  dU = mom * Umu * dt;
+  Umu= Umu+dU;
+  PokeIndex<LorentzIndex>(Dw.Umu,Umu,mu);
+
+  Dw.M(phi,MphiPrime);
+
+  ComplexD Sprime = innerProduct(phi,MphiPrime);
+
+  std::cout << GridLogMessage << " S      "<<S<<std::endl;
+  std::cout << GridLogMessage << " Sprime "<<Sprime<<std::endl;
+  std::cout << GridLogMessage << "dS      "<<Sprime-S<<std::endl;
+
+  Dw.Umu=zero;
+  PokeIndex<LorentzIndex>(Dw.Umu,dU,mu);
+  Dw.M(phi,dMphi);
+
+
+  ComplexD deltaS = innerProduct(phi,dMphi);
+  std::cout << GridLogMessage << "deltaS      "<<deltaS<<std::endl;
+
+  Dw.Umu=zero;
+  PokeIndex<LorentzIndex>(Dw.Umu,Umu_save,mu);
+  Dw.Mdir(phi,dMphi,mu,1);
+  dMphi = dt*mom*dMphi;
+
+  deltaS = innerProduct(phi,dMphi);
+  std::cout << GridLogMessage << "deltaS from inner prod of mom* M[u]     "<<deltaS<<std::endl;
+
+  deltaS = sum(trace(outerProduct(dMphi,phi)));
+  std::cout << GridLogMessage << "deltaS from trace outer prod of deltaM      "<<deltaS<<std::endl;
+
+/*
+  LatticeComplex lip(&Grid);
+  lip  = localInnerProduct(phi,dMphi);
+  
+  LatticeComplex trop(&Grid);
+  trop = trace(outerProduct(dMphi,phi));
+
+  LatticeSpinColourMatrix op(&Grid);
+  op = outerProduct(dMphi,phi);
+
+  LatticeSpinColourMatrix hop(&Grid);
+  LatticeComplex op_cpt(&Grid);
+  for(int s1=0;s1<Ns;s1++){
+  for(int s2=0;s2<Ns;s2++){
+  for(int c1=0;c1<Nc;c1++){
+  for(int c2=0;c2<Nc;c2++){
+
+    op_cpt = peekColour(peekSpin(dMphi,s1),c1) * adj(peekColour(peekSpin(phi,s2),c2));
+
+    parallel_for(auto i=hop.begin();i<hop.end();i++){
+      hop[i]()(s1,s2)(c1,c2) = op_cpt[i]()()();
+    }
+    
+  }}}}
+
+  LatticeSpinColourMatrix diffop(&Grid);
+
+  diffop = hop - op;
+  std::cout << GridLogMessage << "hand outer prod diff   "<<norm2(diffop)<<std::endl;
+
+  deltaS = sum(trace(hop));
+  std::cout << GridLogMessage << "deltaS hop   "<<deltaS<<std::endl;
+
+  std::cout << GridLogMessage<< "  phi[0] : "<<  phi._odata[0]<<std::endl;
+  std::cout << GridLogMessage<< "dMphi[0] : "<<dMphi._odata[0]<<std::endl;
+  std::cout << GridLogMessage<< "hop[0]   : "<<  hop._odata[0]<<std::endl;
+  std::cout << GridLogMessage<< " op[0]   : "<<   op._odata[0]<<std::endl;
+
+
+  std::cout << GridLogMessage << "lip      "<<lip<<std::endl;
+  std::cout << GridLogMessage << "trop     "<<trop<<std::endl;
+
+*/  
+  
+  //  std::cout << GridLogMessage << " UdSdU " << UdSdU << std::endl;
+
+  LatticeComplex dS(&Grid); dS = zero;
+  parallel_for(auto i=mom.begin();i<mom.end();i++){
+    dS[i]() = trace(mom[i]() * UdSdU[i](mu) )*dt;
+  }
+  Complex dSpred = sum(dS);
+
+  std::cout << GridLogMessage << "predict dS    "<< dSpred <<std::endl;
+
+
+  cout<< GridLogMessage << "Done" <<std::endl;
+  Grid_finalize();
+}