diff --git a/Grid/GridQCDcore.h b/Grid/GridQCDcore.h
index cae6f43f..065b62cd 100644
--- a/Grid/GridQCDcore.h
+++ b/Grid/GridQCDcore.h
@@ -36,6 +36,7 @@ Author: paboyle <paboyle@ph.ed.ac.uk>
 #include <Grid/GridCore.h>
 #include <Grid/qcd/QCD.h>
 #include <Grid/qcd/spin/Spin.h>
+#include <Grid/qcd/gparity/Gparity.h>
 #include <Grid/qcd/utils/Utils.h>
 #include <Grid/qcd/representations/Representations.h>
 NAMESPACE_CHECK(GridQCDCore);
diff --git a/Grid/algorithms/Algorithms.h b/Grid/algorithms/Algorithms.h
index 7f27784b..47a0a92b 100644
--- a/Grid/algorithms/Algorithms.h
+++ b/Grid/algorithms/Algorithms.h
@@ -54,6 +54,7 @@ NAMESPACE_CHECK(BiCGSTAB);
 #include <Grid/algorithms/iterative/SchurRedBlack.h>
 #include <Grid/algorithms/iterative/ConjugateGradientMultiShift.h>
 #include <Grid/algorithms/iterative/ConjugateGradientMixedPrec.h>
+#include <Grid/algorithms/iterative/ConjugateGradientMultiShiftMixedPrec.h>
 #include <Grid/algorithms/iterative/BiCGSTABMixedPrec.h>
 #include <Grid/algorithms/iterative/BlockConjugateGradient.h>
 #include <Grid/algorithms/iterative/ConjugateGradientReliableUpdate.h>
diff --git a/Grid/algorithms/iterative/ConjugateGradientMixedPrec.h b/Grid/algorithms/iterative/ConjugateGradientMixedPrec.h
index cd2e4374..31ac55e0 100644
--- a/Grid/algorithms/iterative/ConjugateGradientMixedPrec.h
+++ b/Grid/algorithms/iterative/ConjugateGradientMixedPrec.h
@@ -49,6 +49,7 @@ NAMESPACE_BEGIN(Grid);
     Integer TotalInnerIterations; //Number of inner CG iterations
     Integer TotalOuterIterations; //Number of restarts
     Integer TotalFinalStepIterations; //Number of CG iterations in final patch-up step
+    RealD TrueResidual;
 
     //Option to speed up *inner single precision* solves using a LinearFunction that produces a guess
     LinearFunction<FieldF> *guesser;
@@ -68,6 +69,7 @@ NAMESPACE_BEGIN(Grid);
     }
   
   void operator() (const FieldD &src_d_in, FieldD &sol_d){
+    std::cout << GridLogMessage << "MixedPrecisionConjugateGradient: Starting mixed precision CG with outer tolerance " << Tolerance << " and inner tolerance " << InnerTolerance << std::endl;
     TotalInnerIterations = 0;
 	
     GridStopWatch TotalTimer;
@@ -97,6 +99,7 @@ NAMESPACE_BEGIN(Grid);
     FieldF sol_f(SinglePrecGrid);
     sol_f.Checkerboard() = cb;
     
+    std::cout<<GridLogMessage<<"MixedPrecisionConjugateGradient: Starting initial inner CG with tolerance " << inner_tol << std::endl;
     ConjugateGradient<FieldF> CG_f(inner_tol, MaxInnerIterations);
     CG_f.ErrorOnNoConverge = false;
 
@@ -130,6 +133,7 @@ NAMESPACE_BEGIN(Grid);
 	(*guesser)(src_f, sol_f);
 
       //Inner CG
+      std::cout<<GridLogMessage<<"MixedPrecisionConjugateGradient: Outer iteration " << outer_iter << " starting inner CG with tolerance " << inner_tol << std::endl;
       CG_f.Tolerance = inner_tol;
       InnerCGtimer.Start();
       CG_f(Linop_f, src_f, sol_f);
@@ -150,6 +154,7 @@ NAMESPACE_BEGIN(Grid);
     ConjugateGradient<FieldD> CG_d(Tolerance, MaxInnerIterations);
     CG_d(Linop_d, src_d_in, sol_d);
     TotalFinalStepIterations = CG_d.IterationsToComplete;
+    TrueResidual = CG_d.TrueResidual;
 
     TotalTimer.Stop();
     std::cout<<GridLogMessage<<"MixedPrecisionConjugateGradient: Inner CG iterations " << TotalInnerIterations << " Restarts " << TotalOuterIterations << " Final CG iterations " << TotalFinalStepIterations << std::endl;
diff --git a/Grid/algorithms/iterative/ConjugateGradientMultiShift.h b/Grid/algorithms/iterative/ConjugateGradientMultiShift.h
index 63a6c55f..3b079e99 100644
--- a/Grid/algorithms/iterative/ConjugateGradientMultiShift.h
+++ b/Grid/algorithms/iterative/ConjugateGradientMultiShift.h
@@ -52,7 +52,7 @@ public:
   MultiShiftFunction shifts;
   std::vector<RealD> TrueResidualShift;
 
-  ConjugateGradientMultiShift(Integer maxit,MultiShiftFunction &_shifts) : 
+  ConjugateGradientMultiShift(Integer maxit, const MultiShiftFunction &_shifts) : 
     MaxIterations(maxit),
     shifts(_shifts)
   { 
@@ -182,6 +182,9 @@ public:
     for(int s=0;s<nshift;s++) {
       axpby(psi[s],0.,-bs[s]*alpha[s],src,src);
     }
+
+    std::cout << GridLogIterative << "ConjugateGradientMultiShift: initial rn (|src|^2) =" << rn << " qq (|MdagM src|^2) =" << qq << " d ( dot(src, [MdagM + m_0]src) ) =" << d << " c=" << c << std::endl;
+    
   
   ///////////////////////////////////////
   // Timers
diff --git a/Grid/algorithms/iterative/ConjugateGradientMultiShiftMixedPrec.h b/Grid/algorithms/iterative/ConjugateGradientMultiShiftMixedPrec.h
new file mode 100644
index 00000000..1470a229
--- /dev/null
+++ b/Grid/algorithms/iterative/ConjugateGradientMultiShiftMixedPrec.h
@@ -0,0 +1,409 @@
+/*************************************************************************************
+
+    Grid physics library, www.github.com/paboyle/Grid 
+
+    Source file: ./lib/algorithms/iterative/ConjugateGradientMultiShift.h
+
+    Copyright (C) 2015
+
+Author: Azusa Yamaguchi <ayamaguc@staffmail.ed.ac.uk>
+Author: Peter Boyle <paboyle@ph.ed.ac.uk>
+Author: Christopher Kelly <ckelly@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 */
+#ifndef GRID_CONJUGATE_GRADIENT_MULTI_SHIFT_MIXEDPREC_H
+#define GRID_CONJUGATE_GRADIENT_MULTI_SHIFT_MIXEDPREC_H
+
+NAMESPACE_BEGIN(Grid);
+
+//CK 2020: A variant of the multi-shift conjugate gradient with the matrix multiplication in single precision. 
+//The residual is stored in single precision, but the search directions and solution are stored in double precision. 
+//Every update_freq iterations the residual is corrected in double precision. 
+    
+//For safety the a final regular CG is applied to clean up if necessary
+
+//Linop to add shift to input linop, used in cleanup CG
+namespace ConjugateGradientMultiShiftMixedPrecSupport{
+template<typename Field>
+class ShiftedLinop: public LinearOperatorBase<Field>{
+public:
+  LinearOperatorBase<Field> &linop_base;
+  RealD shift;
+
+  ShiftedLinop(LinearOperatorBase<Field> &_linop_base, RealD _shift): linop_base(_linop_base), shift(_shift){}
+
+  void OpDiag (const Field &in, Field &out){ assert(0); }
+  void OpDir  (const Field &in, Field &out,int dir,int disp){ assert(0); }
+  void OpDirAll  (const Field &in, std::vector<Field> &out){ assert(0); }
+  
+  void Op     (const Field &in, Field &out){ assert(0); }
+  void AdjOp  (const Field &in, Field &out){ assert(0); }
+
+  void HermOp(const Field &in, Field &out){
+    linop_base.HermOp(in, out);
+    axpy(out, shift, in, out);
+  }    
+
+  void HermOpAndNorm(const Field &in, Field &out,RealD &n1,RealD &n2){
+    HermOp(in,out);
+    ComplexD dot = innerProduct(in,out);
+    n1=real(dot);
+    n2=norm2(out);
+  }
+};
+};
+
+
+template<class FieldD, class FieldF,
+	 typename std::enable_if< getPrecision<FieldD>::value == 2, int>::type = 0,
+	 typename std::enable_if< getPrecision<FieldF>::value == 1, int>::type = 0> 
+class ConjugateGradientMultiShiftMixedPrec : public OperatorMultiFunction<FieldD>,
+					     public OperatorFunction<FieldD>
+{
+public:                                                
+
+  using OperatorFunction<FieldD>::operator();
+
+  RealD   Tolerance;
+  Integer MaxIterations;
+  Integer IterationsToComplete; //Number of iterations the CG took to finish. Filled in upon completion
+  std::vector<int> IterationsToCompleteShift;  // Iterations for this shift
+  int verbose;
+  MultiShiftFunction shifts;
+  std::vector<RealD> TrueResidualShift;
+
+  int ReliableUpdateFreq; //number of iterations between reliable updates
+
+  GridBase* SinglePrecGrid; //Grid for single-precision fields
+  LinearOperatorBase<FieldF> &Linop_f; //single precision
+
+  ConjugateGradientMultiShiftMixedPrec(Integer maxit, const MultiShiftFunction &_shifts,
+				       GridBase* _SinglePrecGrid, LinearOperatorBase<FieldF> &_Linop_f,
+				       int _ReliableUpdateFreq
+				       ) : 
+    MaxIterations(maxit),  shifts(_shifts), SinglePrecGrid(_SinglePrecGrid), Linop_f(_Linop_f), ReliableUpdateFreq(_ReliableUpdateFreq)
+  { 
+    verbose=1;
+    IterationsToCompleteShift.resize(_shifts.order);
+    TrueResidualShift.resize(_shifts.order);
+  }
+
+  void operator() (LinearOperatorBase<FieldD> &Linop, const FieldD &src, FieldD &psi)
+  {
+    GridBase *grid = src.Grid();
+    int nshift = shifts.order;
+    std::vector<FieldD> results(nshift,grid);
+    (*this)(Linop,src,results,psi);
+  }
+  void operator() (LinearOperatorBase<FieldD> &Linop, const FieldD &src, std::vector<FieldD> &results, FieldD &psi)
+  {
+    int nshift = shifts.order;
+
+    (*this)(Linop,src,results);
+  
+    psi = shifts.norm*src;
+    for(int i=0;i<nshift;i++){
+      psi = psi + shifts.residues[i]*results[i];
+    }
+
+    return;
+  }
+
+  void operator() (LinearOperatorBase<FieldD> &Linop_d, const FieldD &src_d, std::vector<FieldD> &psi_d)
+  { 
+    GridBase *DoublePrecGrid = src_d.Grid();
+
+    ////////////////////////////////////////////////////////////////////////
+    // Convenience references to the info stored in "MultiShiftFunction"
+    ////////////////////////////////////////////////////////////////////////
+    int nshift = shifts.order;
+
+    std::vector<RealD> &mass(shifts.poles); // Make references to array in "shifts"
+    std::vector<RealD> &mresidual(shifts.tolerances);
+    std::vector<RealD> alpha(nshift,1.0);
+
+    //Double precision search directions
+    FieldD p_d(DoublePrecGrid);
+    std::vector<FieldD> ps_d(nshift, DoublePrecGrid);// Search directions (double precision)
+
+    FieldD tmp_d(DoublePrecGrid);
+    FieldD r_d(DoublePrecGrid);
+    FieldD mmp_d(DoublePrecGrid);
+
+    assert(psi_d.size()==nshift);
+    assert(mass.size()==nshift);
+    assert(mresidual.size()==nshift);
+  
+    // dynamic sized arrays on stack; 2d is a pain with vector
+    RealD  bs[nshift];
+    RealD  rsq[nshift];
+    RealD  z[nshift][2];
+    int     converged[nshift];
+  
+    const int       primary =0;
+  
+    //Primary shift fields CG iteration
+    RealD a,b,c,d;
+    RealD cp,bp,qq; //prev
+  
+    // Matrix mult fields
+    FieldF r_f(SinglePrecGrid);
+    FieldF p_f(SinglePrecGrid);
+    FieldF tmp_f(SinglePrecGrid);
+    FieldF mmp_f(SinglePrecGrid);
+    FieldF src_f(SinglePrecGrid);
+    precisionChange(src_f, src_d);
+
+    // Check lightest mass
+    for(int s=0;s<nshift;s++){
+      assert( mass[s]>= mass[primary] );
+      converged[s]=0;
+    }
+  
+    // Wire guess to zero
+    // Residuals "r" are src
+    // First search direction "p" is also src
+    cp = norm2(src_d);
+
+    // Handle trivial case of zero src.
+    if( cp == 0. ){
+      for(int s=0;s<nshift;s++){
+	psi_d[s] = Zero();
+	IterationsToCompleteShift[s] = 1;
+	TrueResidualShift[s] = 0.;
+      }
+      return;
+    }
+
+    for(int s=0;s<nshift;s++){
+      rsq[s] = cp * mresidual[s] * mresidual[s];
+      std::cout<<GridLogMessage<<"ConjugateGradientMultiShiftMixedPrec: shift "<< s <<" target resid "<<rsq[s]<<std::endl;
+      ps_d[s] = src_d;
+    }
+    // r and p for primary
+    r_f=src_f; //residual maintained in single
+    p_f=src_f;
+    p_d = src_d; //primary copy --- make this a reference to ps_d to save axpys
+  
+    //MdagM+m[0]
+    Linop_f.HermOpAndNorm(p_f,mmp_f,d,qq); // mmp = MdagM p        d=real(dot(p, mmp)),  qq=norm2(mmp)
+    axpy(mmp_f,mass[0],p_f,mmp_f);
+    RealD rn = norm2(p_f);
+    d += rn*mass[0];
+
+    b = -cp /d;
+  
+    // Set up the various shift variables
+    int       iz=0;
+    z[0][1-iz] = 1.0;
+    z[0][iz]   = 1.0;
+    bs[0]      = b;
+    for(int s=1;s<nshift;s++){
+      z[s][1-iz] = 1.0;
+      z[s][iz]   = 1.0/( 1.0 - b*(mass[s]-mass[0]));
+      bs[s]      = b*z[s][iz]; 
+    }
+  
+    // r += b[0] A.p[0]
+    // c= norm(r)
+    c=axpy_norm(r_f,b,mmp_f,r_f);
+  
+    for(int s=0;s<nshift;s++) {
+      axpby(psi_d[s],0.,-bs[s]*alpha[s],src_d,src_d);
+    }
+  
+    ///////////////////////////////////////
+    // Timers
+    ///////////////////////////////////////
+    GridStopWatch AXPYTimer, ShiftTimer, QRTimer, MatrixTimer, SolverTimer, PrecChangeTimer, CleanupTimer;
+
+    SolverTimer.Start();
+  
+    // Iteration loop
+    int k;
+  
+    for (k=1;k<=MaxIterations;k++){    
+      a = c /cp;
+
+      //Update double precision search direction by residual
+      PrecChangeTimer.Start();
+      precisionChange(r_d, r_f);
+      PrecChangeTimer.Stop();
+
+      AXPYTimer.Start();
+      axpy(p_d,a,p_d,r_d); 
+
+      for(int s=0;s<nshift;s++){
+	if ( ! converged[s] ) { 
+	  if (s==0){
+	    axpy(ps_d[s],a,ps_d[s],r_d);
+	  } else{
+	    RealD as =a *z[s][iz]*bs[s] /(z[s][1-iz]*b);
+	    axpby(ps_d[s],z[s][iz],as,r_d,ps_d[s]);
+	  }
+	}
+      }
+      AXPYTimer.Stop();
+
+      PrecChangeTimer.Start();
+      precisionChange(p_f, p_d); //get back single prec search direction for linop
+      PrecChangeTimer.Stop();
+
+      cp=c;
+      MatrixTimer.Start();  
+      Linop_f.HermOp(p_f,mmp_f); 
+      d=real(innerProduct(p_f,mmp_f));    
+      MatrixTimer.Stop();  
+
+      AXPYTimer.Start();
+      axpy(mmp_f,mass[0],p_f,mmp_f);
+      AXPYTimer.Stop();
+      RealD rn = norm2(p_f);
+      d += rn*mass[0];
+    
+      bp=b;
+      b=-cp/d;
+    
+      // Toggle the recurrence history
+      bs[0] = b;
+      iz = 1-iz;
+      ShiftTimer.Start();
+      for(int s=1;s<nshift;s++){
+	if((!converged[s])){
+	  RealD z0 = z[s][1-iz];
+	  RealD z1 = z[s][iz];
+	  z[s][iz] = z0*z1*bp
+	    / (b*a*(z1-z0) + z1*bp*(1- (mass[s]-mass[0])*b)); 
+	  bs[s] = b*z[s][iz]/z0; // NB sign  rel to Mike
+	}
+      }
+      ShiftTimer.Stop();
+
+      //Update double precision solutions
+      AXPYTimer.Start();
+      for(int s=0;s<nshift;s++){
+	int ss = s;
+	if( (!converged[s]) ) { 
+	  axpy(psi_d[ss],-bs[s]*alpha[s],ps_d[s],psi_d[ss]);
+	}
+      }
+
+      //Perform reliable update if necessary; otherwise update residual from single-prec mmp
+      RealD c_f = axpy_norm(r_f,b,mmp_f,r_f);
+      AXPYTimer.Stop();
+
+      c = c_f;
+
+      if(k % ReliableUpdateFreq == 0){
+	//Replace r with true residual
+	MatrixTimer.Start();  
+	Linop_d.HermOp(psi_d[0],mmp_d); 
+	MatrixTimer.Stop();  
+
+	AXPYTimer.Start();
+	axpy(mmp_d,mass[0],psi_d[0],mmp_d);
+
+	RealD c_d = axpy_norm(r_d, -1.0, mmp_d, src_d);
+	AXPYTimer.Stop();
+
+	std::cout<<GridLogMessage<<"ConjugateGradientMultiShiftMixedPrec k="<<k<< ", replaced |r|^2 = "<<c_f <<" with |r|^2 = "<<c_d<<std::endl;
+	
+	PrecChangeTimer.Start();
+	precisionChange(r_f, r_d);
+	PrecChangeTimer.Stop();
+	c = c_d;
+      }
+    
+      // Convergence checks
+      int all_converged = 1;
+      for(int s=0;s<nshift;s++){
+      
+	if ( (!converged[s]) ){
+	  IterationsToCompleteShift[s] = k;
+	
+	  RealD css  = c * z[s][iz]* z[s][iz];
+	
+	  if(css<rsq[s]){
+	    if ( ! converged[s] )
+	      std::cout<<GridLogMessage<<"ConjugateGradientMultiShiftMixedPrec k="<<k<<" Shift "<<s<<" has converged"<<std::endl;
+	    converged[s]=1;
+	  } else {
+	    all_converged=0;
+	  }
+
+	}
+      }
+
+      if ( all_converged ){
+
+	SolverTimer.Stop();
+	std::cout<<GridLogMessage<< "ConjugateGradientMultiShiftMixedPrec: All shifts have converged iteration "<<k<<std::endl;
+	std::cout<<GridLogMessage<< "ConjugateGradientMultiShiftMixedPrec: Checking solutions"<<std::endl;
+      
+	// Check answers 
+	for(int s=0; s < nshift; s++) { 
+	  Linop_d.HermOpAndNorm(psi_d[s],mmp_d,d,qq);
+	  axpy(tmp_d,mass[s],psi_d[s],mmp_d);
+	  axpy(r_d,-alpha[s],src_d,tmp_d);
+	  RealD rn = norm2(r_d);
+	  RealD cn = norm2(src_d);
+	  TrueResidualShift[s] = std::sqrt(rn/cn);
+	  std::cout<<GridLogMessage<<"ConjugateGradientMultiShiftMixedPrec: shift["<<s<<"] true residual "<< TrueResidualShift[s] << " target " << mresidual[s] << std::endl;
+
+	  //If we have not reached the desired tolerance, do a (mixed precision) CG cleanup
+	  if(rn >= rsq[s]){
+	    CleanupTimer.Start();
+	    std::cout<<GridLogMessage<<"ConjugateGradientMultiShiftMixedPrec: performing cleanup step for shift " << s << std::endl;
+
+	    //Setup linear operators for final cleanup
+	    ConjugateGradientMultiShiftMixedPrecSupport::ShiftedLinop<FieldD> Linop_shift_d(Linop_d, mass[s]);
+	    ConjugateGradientMultiShiftMixedPrecSupport::ShiftedLinop<FieldF> Linop_shift_f(Linop_f, mass[s]);
+					       
+	    MixedPrecisionConjugateGradient<FieldD,FieldF> cg(mresidual[s], MaxIterations, MaxIterations, SinglePrecGrid, Linop_shift_f, Linop_shift_d); 
+	    cg(src_d, psi_d[s]);
+	    
+	    TrueResidualShift[s] = cg.TrueResidual;
+	    CleanupTimer.Stop();
+	  }
+	}
+
+	std::cout << GridLogMessage << "ConjugateGradientMultiShiftMixedPrec: Time Breakdown for body"<<std::endl;
+	std::cout << GridLogMessage << "\tSolver    " << SolverTimer.Elapsed()     <<std::endl;
+	std::cout << GridLogMessage << "\t\tAXPY    " << AXPYTimer.Elapsed()     <<std::endl;
+	std::cout << GridLogMessage << "\t\tMatrix    " << MatrixTimer.Elapsed()     <<std::endl;
+	std::cout << GridLogMessage << "\t\tShift    " << ShiftTimer.Elapsed()     <<std::endl;
+	std::cout << GridLogMessage << "\t\tPrecision Change " << PrecChangeTimer.Elapsed()     <<std::endl;
+	std::cout << GridLogMessage << "\tFinal Cleanup " << CleanupTimer.Elapsed()     <<std::endl;
+	std::cout << GridLogMessage << "\tSolver+Cleanup " << SolverTimer.Elapsed() + CleanupTimer.Elapsed() << std::endl;
+
+	IterationsToComplete = k;	
+
+	return;
+      }
+
+   
+    }
+    // ugly hack
+    std::cout<<GridLogMessage<<"CG multi shift did not converge"<<std::endl;
+    //  assert(0);
+  }
+
+};
+NAMESPACE_END(Grid);
+#endif
diff --git a/Grid/algorithms/iterative/LocalCoherenceLanczos.h b/Grid/algorithms/iterative/LocalCoherenceLanczos.h
index dc82134a..0a2fe55c 100644
--- a/Grid/algorithms/iterative/LocalCoherenceLanczos.h
+++ b/Grid/algorithms/iterative/LocalCoherenceLanczos.h
@@ -44,6 +44,7 @@ public:
 				  int, MinRes);    // Must restart
 };
 
+//This class is the input parameter class for some testing programs
 struct LocalCoherenceLanczosParams : Serializable {
 public:
   GRID_SERIALIZABLE_CLASS_MEMBERS(LocalCoherenceLanczosParams,
@@ -155,6 +156,7 @@ public:
       _coarse_relax_tol(coarse_relax_tol)  
   {    };
 
+  //evalMaxApprox: approximation of largest eval of the fine Chebyshev operator (suitably wrapped by block projection)
   int TestConvergence(int j,RealD eresid,CoarseField &B, RealD &eval,RealD evalMaxApprox)
   {
     CoarseField v(B);
@@ -181,8 +183,16 @@ public:
     if( (vv<eresid*eresid) ) conv = 1;
     return conv;
   }
-  int ReconstructEval(int j,RealD eresid,CoarseField &B, RealD &eval,RealD evalMaxApprox)
+
+  //This function is called at the end of the coarse grid Lanczos. It promotes the coarse eigenvector 'B' to the fine grid,
+  //applies a smoother to the result then computes the computes the *fine grid* eigenvalue (output as 'eval').
+
+  //evalMaxApprox should be the approximation of the largest eval of the fine Hermop. However when this function is called by IRL it actually passes the largest eval of the *Chebyshev* operator (as this is the max approx used for the TestConvergence above)
+  //As the largest eval of the Chebyshev is typically several orders of magnitude larger this makes the convergence test pass even when it should not.
+  //We therefore ignore evalMaxApprox here and use a value of 1.0 (note this value is already used by TestCoarse)
+  int ReconstructEval(int j,RealD eresid,CoarseField &B, RealD &eval,RealD evalMaxApprox)  
   {
+    evalMaxApprox = 1.0; //cf above
     GridBase *FineGrid = _subspace[0].Grid();    
     int checkerboard   = _subspace[0].Checkerboard();
     FineField fB(FineGrid);fB.Checkerboard() =checkerboard;
@@ -201,13 +211,13 @@ public:
     eval   = vnum/vden;
     fv -= eval*fB;
     RealD vv = norm2(fv) / ::pow(evalMaxApprox,2.0);
-
+    if ( j > nbasis ) eresid = eresid*_coarse_relax_tol;
+    
     std::cout.precision(13);
     std::cout<<GridLogIRL  << "[" << std::setw(3)<<j<<"] "
 	     <<"eval = "<<std::setw(25)<< eval << " (" << eval_poly << ")"
-	     <<" |H B[i] - eval[i]B[i]|^2 / evalMaxApprox^2 " << std::setw(25) << vv
+	     <<" |H B[i] - eval[i]B[i]|^2 / evalMaxApprox^2 " << std::setw(25) << vv << " target " << eresid*eresid
 	     <<std::endl;
-    if ( j > nbasis ) eresid = eresid*_coarse_relax_tol;
     if( (vv<eresid*eresid) ) return 1;
     return 0;
   }
@@ -285,6 +295,10 @@ public:
     evals_coarse.resize(0);
   };
 
+  //The block inner product is the inner product on the fine grid locally summed over the blocks
+  //to give a Lattice<Scalar> on the coarse grid. This function orthnormalizes the fine-grid subspace
+  //vectors under the block inner product. This step must be performed after computing the fine grid
+  //eigenvectors and before computing the coarse grid eigenvectors.    
   void Orthogonalise(void ) {
     CoarseScalar InnerProd(_CoarseGrid);
     std::cout << GridLogMessage <<" Gramm-Schmidt pass 1"<<std::endl;
@@ -328,6 +342,8 @@ public:
     }
   }
 
+  //While this method serves to check the coarse eigenvectors, it also recomputes the eigenvalues from the smoothed reconstructed eigenvectors
+  //hence the smoother can be tuned after running the coarse Lanczos by using a different smoother here
   void testCoarse(RealD resid,ChebyParams cheby_smooth,RealD relax) 
   {
     assert(evals_fine.size() == nbasis);
@@ -376,25 +392,31 @@ public:
     evals_fine.resize(nbasis);
     subspace.resize(nbasis,_FineGrid);
   }
+
+
+  //cheby_op: Parameters of the fine grid Chebyshev polynomial used for the Lanczos acceleration
+  //cheby_smooth: Parameters of a separate Chebyshev polynomial used after the Lanczos has completed to smooth out high frequency noise in the reconstructed fine grid eigenvectors prior to computing the eigenvalue
+  //relax: Reconstructed eigenvectors (post smoothing) are naturally not as precise as true eigenvectors. This factor acts as a multiplier on the stopping condition when determining whether the results satisfy the user provided stopping condition
   void calcCoarse(ChebyParams cheby_op,ChebyParams cheby_smooth,RealD relax,
 		  int Nstop, int Nk, int Nm,RealD resid, 
 		  RealD MaxIt, RealD betastp, int MinRes)
   {
-    Chebyshev<FineField>                          Cheby(cheby_op);
-    ProjectedHermOp<Fobj,CComplex,nbasis>         Op(_FineOp,subspace);
-    ProjectedFunctionHermOp<Fobj,CComplex,nbasis> ChebyOp (Cheby,_FineOp,subspace);
+    Chebyshev<FineField>                          Cheby(cheby_op); //Chebyshev of fine operator on fine grid
+    ProjectedHermOp<Fobj,CComplex,nbasis>         Op(_FineOp,subspace); //Fine operator on coarse grid with intermediate fine grid conversion
+    ProjectedFunctionHermOp<Fobj,CComplex,nbasis> ChebyOp (Cheby,_FineOp,subspace); //Chebyshev of fine operator on coarse grid with intermediate fine grid conversion
     //////////////////////////////////////////////////////////////////////////////////////////////////
     // create a smoother and see if we can get a cheap convergence test and smooth inside the IRL
     //////////////////////////////////////////////////////////////////////////////////////////////////
 
-    Chebyshev<FineField>                                           ChebySmooth(cheby_smooth);
-    ImplicitlyRestartedLanczosSmoothedTester<Fobj,CComplex,nbasis> ChebySmoothTester(ChebyOp,ChebySmooth,_FineOp,subspace,relax);
+    Chebyshev<FineField>                                           ChebySmooth(cheby_smooth); //lower order Chebyshev of fine operator on fine grid used to smooth regenerated eigenvectors
+    ImplicitlyRestartedLanczosSmoothedTester<Fobj,CComplex,nbasis> ChebySmoothTester(ChebyOp,ChebySmooth,_FineOp,subspace,relax); 
 
     evals_coarse.resize(Nm);
     evec_coarse.resize(Nm,_CoarseGrid);
 
     CoarseField src(_CoarseGrid);     src=1.0; 
 
+    //Note the "tester" here is also responsible for generating the fine grid eigenvalues which are output into the "evals_coarse" array
     ImplicitlyRestartedLanczos<CoarseField> IRL(ChebyOp,ChebyOp,ChebySmoothTester,Nstop,Nk,Nm,resid,MaxIt,betastp,MinRes);
     int Nconv=0;
     IRL.calc(evals_coarse,evec_coarse,src,Nconv,false);
@@ -405,6 +427,14 @@ public:
       std::cout << i << " Coarse eval = " << evals_coarse[i]  << std::endl;
     }
   }
+
+  //Get the fine eigenvector 'i' by reconstruction
+  void getFineEvecEval(FineField &evec, RealD &eval, const int i) const{
+    blockPromote(evec_coarse[i],evec,subspace);  
+    eval = evals_coarse[i];
+  }
+    
+    
 };
 
 NAMESPACE_END(Grid);
diff --git a/Grid/algorithms/iterative/PowerMethod.h b/Grid/algorithms/iterative/PowerMethod.h
index 6aa8e923..027ea68c 100644
--- a/Grid/algorithms/iterative/PowerMethod.h
+++ b/Grid/algorithms/iterative/PowerMethod.h
@@ -29,6 +29,8 @@ template<class Field> class PowerMethod
       RealD vnum = real(innerProduct(src_n,tmp)); // HermOp. 
       RealD vden = norm2(src_n); 
       RealD na = vnum/vden; 
+
+      std::cout << GridLogIterative << "PowerMethod: Current approximation of largest eigenvalue " << na << std::endl;
       
       if ( (fabs(evalMaxApprox/na - 1.0) < 0.001) || (i==_MAX_ITER_EST_-1) ) { 
  	evalMaxApprox = na; 
diff --git a/Grid/lattice/Lattice_transfer.h b/Grid/lattice/Lattice_transfer.h
index ef489ea6..aee55e93 100644
--- a/Grid/lattice/Lattice_transfer.h
+++ b/Grid/lattice/Lattice_transfer.h
@@ -855,7 +855,7 @@ void ExtractSliceLocal(Lattice<vobj> &lowDim,const Lattice<vobj> & higherDim,int
 
 
 template<class vobj>
-void Replicate(Lattice<vobj> &coarse,Lattice<vobj> & fine)
+void Replicate(const Lattice<vobj> &coarse,Lattice<vobj> & fine)
 {
   typedef typename vobj::scalar_object sobj;
 
diff --git a/Grid/parallelIO/NerscIO.h b/Grid/parallelIO/NerscIO.h
index 99011e25..88278131 100644
--- a/Grid/parallelIO/NerscIO.h
+++ b/Grid/parallelIO/NerscIO.h
@@ -39,9 +39,11 @@ using namespace Grid;
 ////////////////////////////////////////////////////////////////////////////////
 class NerscIO : public BinaryIO { 
 public:
-
   typedef Lattice<vLorentzColourMatrixD> GaugeField;
 
+  // Enable/disable exiting if the plaquette in the header does not match the value computed (default true)
+  static bool & exitOnReadPlaquetteMismatch(){ static bool v=true; return v; }
+
   static inline void truncate(std::string file){
     std::ofstream fout(file,std::ios::out);
   }
@@ -198,7 +200,7 @@ public:
       std::cerr << " nersc_csum  " <<std::hex<< nersc_csum << " " << header.checksum<< std::dec<< std::endl;
       exit(0);
     }
-    assert(fabs(clone.plaquette -header.plaquette ) < 1.0e-5 );
+    if(exitOnReadPlaquetteMismatch()) assert(fabs(clone.plaquette -header.plaquette ) < 1.0e-5 );
     assert(fabs(clone.link_trace-header.link_trace) < 1.0e-6 );
     assert(nersc_csum == header.checksum );
       
diff --git a/Grid/qcd/QCD.h b/Grid/qcd/QCD.h
index 858aead7..81356a66 100644
--- a/Grid/qcd/QCD.h
+++ b/Grid/qcd/QCD.h
@@ -63,6 +63,7 @@ static constexpr int Ngp=2; // gparity index range
 #define ColourIndex  (2)
 #define SpinIndex    (1)
 #define LorentzIndex (0)
+#define GparityFlavourIndex (0)
 
 // Also should make these a named enum type
 static constexpr int DaggerNo=0;
@@ -87,6 +88,8 @@ template<typename T> struct isCoarsened {
 template <typename T> using IfCoarsened    = Invoke<std::enable_if< isCoarsened<T>::value,int> > ;
 template <typename T> using IfNotCoarsened = Invoke<std::enable_if<!isCoarsened<T>::value,int> > ;
 
+const int GparityFlavourTensorIndex = 3; //TensorLevel counts from the bottom!
+
 // ChrisK very keen to add extra space for Gparity doubling.
 //
 // Also add domain wall index, in a way where Wilson operator 
@@ -110,8 +113,10 @@ template<typename vtype> using iHalfSpinColourVector      = iScalar<iVector<iVec
     template<typename vtype> using iSpinColourSpinColourMatrix  = iScalar<iMatrix<iMatrix<iMatrix<iMatrix<vtype, Nc>, Ns>, Nc>, Ns> >;
 
 
+template<typename vtype> using iGparityFlavourVector                = iVector<iScalar<iScalar<vtype> >, Ngp>;
 template<typename vtype> using iGparitySpinColourVector       = iVector<iVector<iVector<vtype, Nc>, Ns>, Ngp >;
 template<typename vtype> using iGparityHalfSpinColourVector   = iVector<iVector<iVector<vtype, Nc>, Nhs>, Ngp >;
+template<typename vtype> using iGparityFlavourMatrix = iMatrix<iScalar<iScalar<vtype> >, Ngp>;
 
 // Spin matrix
 typedef iSpinMatrix<Complex  >          SpinMatrix;
@@ -176,6 +181,16 @@ typedef iDoubleStoredColourMatrix<vComplex > vDoubleStoredColourMatrix;
 typedef iDoubleStoredColourMatrix<vComplexF> vDoubleStoredColourMatrixF;
 typedef iDoubleStoredColourMatrix<vComplexD> vDoubleStoredColourMatrixD;
 
+//G-parity flavour matrix
+typedef iGparityFlavourMatrix<Complex> GparityFlavourMatrix;
+typedef iGparityFlavourMatrix<ComplexF> GparityFlavourMatrixF;
+typedef iGparityFlavourMatrix<ComplexD> GparityFlavourMatrixD;
+
+typedef iGparityFlavourMatrix<vComplex> vGparityFlavourMatrix;
+typedef iGparityFlavourMatrix<vComplexF> vGparityFlavourMatrixF;
+typedef iGparityFlavourMatrix<vComplexD> vGparityFlavourMatrixD;
+
+
 // Spin vector
 typedef iSpinVector<Complex >           SpinVector;
 typedef iSpinVector<ComplexF>           SpinVectorF;
@@ -220,6 +235,16 @@ typedef iHalfSpinColourVector<ComplexD> HalfSpinColourVectorD;
 typedef iHalfSpinColourVector<vComplex > vHalfSpinColourVector;
 typedef iHalfSpinColourVector<vComplexF> vHalfSpinColourVectorF;
 typedef iHalfSpinColourVector<vComplexD> vHalfSpinColourVectorD;
+
+//G-parity flavour vector
+typedef iGparityFlavourVector<Complex >         GparityFlavourVector;
+typedef iGparityFlavourVector<ComplexF>         GparityFlavourVectorF;
+typedef iGparityFlavourVector<ComplexD>         GparityFlavourVectorD;
+
+typedef iGparityFlavourVector<vComplex >         vGparityFlavourVector;
+typedef iGparityFlavourVector<vComplexF>         vGparityFlavourVectorF;
+typedef iGparityFlavourVector<vComplexD>         vGparityFlavourVectorD;
+
     
 // singlets
 typedef iSinglet<Complex >         TComplex;     // FIXME This is painful. Tensor singlet complex type.
diff --git a/Grid/qcd/action/fermion/GparityWilsonImpl.h b/Grid/qcd/action/fermion/GparityWilsonImpl.h
index fd627aed..8017bc76 100644
--- a/Grid/qcd/action/fermion/GparityWilsonImpl.h
+++ b/Grid/qcd/action/fermion/GparityWilsonImpl.h
@@ -30,6 +30,18 @@ directory
 
 NAMESPACE_BEGIN(Grid);
 
+/*
+  Policy implementation for G-parity boundary conditions
+
+  Rather than treating the gauge field as a flavored field, the Grid implementation of G-parity treats the gauge field as a regular
+  field with complex conjugate boundary conditions. In order to ensure the second flavor interacts with the conjugate links and the first
+  with the regular links we overload the functionality of doubleStore, whose purpose is to store the gauge field and the barrel-shifted gauge field
+  to avoid communicating links when applying the Dirac operator, such that the double-stored field contains also a flavor index which maps to
+  either the link or the conjugate link. This flavored field is then used by multLink to apply the correct link to a spinor.
+
+  Here the first Nd-1 directions are treated as "spatial", and a twist value of 1 indicates G-parity BCs in that direction. 
+  mu=Nd-1 is assumed to be the time direction and a twist value of 1 indicates antiperiodic BCs
+ */
 template <class S, class Representation = FundamentalRepresentation, class Options=CoeffReal>
 class GparityWilsonImpl : public ConjugateGaugeImpl<GaugeImplTypes<S, Representation::Dimension> > {
 public:
@@ -113,7 +125,7 @@ public:
     || ((distance== 1)&&(icoor[direction]==1))
     || ((distance==-1)&&(icoor[direction]==0));
 
-    permute_lane = permute_lane && SE->_around_the_world && St.parameters.twists[mmu]; //only if we are going around the world
+    permute_lane = permute_lane && SE->_around_the_world && St.parameters.twists[mmu] && mmu < Nd-1; //only if we are going around the world in a spatial direction
 
     //Apply the links
     int f_upper = permute_lane ? 1 : 0;
@@ -139,10 +151,10 @@ public:
     assert((distance == 1) || (distance == -1));  // nearest neighbour stencil hard code
     assert((sl == 1) || (sl == 2));
 
-    if ( SE->_around_the_world && St.parameters.twists[mmu] ) {
-
+    //If this site is an global boundary site, perform the G-parity flavor twist
+    if ( mmu < Nd-1 && SE->_around_the_world && St.parameters.twists[mmu] ) {
       if ( sl == 2 ) {
-       
+	//Only do the twist for lanes on the edge of the physical node
 	ExtractBuffer<sobj> vals(Nsimd);
 
 	extract(chi,vals);
@@ -197,6 +209,19 @@ public:
     reg = memory;
   }
 
+
+  //Poke 'poke_f0' onto flavor 0 and 'poke_f1' onto flavor 1 in direction mu of the doubled gauge field Uds
+  inline void pokeGparityDoubledGaugeField(DoubledGaugeField &Uds, const GaugeLinkField &poke_f0, const GaugeLinkField &poke_f1, const int mu){
+    autoView(poke_f0_v, poke_f0, CpuRead);
+    autoView(poke_f1_v, poke_f1, CpuRead);
+    autoView(Uds_v, Uds, CpuWrite);
+    thread_foreach(ss,poke_f0_v,{
+	Uds_v[ss](0)(mu) = poke_f0_v[ss]();
+	Uds_v[ss](1)(mu) = poke_f1_v[ss]();
+      });
+  }
+    
+
   inline void DoubleStore(GridBase *GaugeGrid,DoubledGaugeField &Uds,const GaugeField &Umu)
   {
     conformable(Uds.Grid(),GaugeGrid);
@@ -207,14 +232,19 @@ public:
     GaugeLinkField Uconj(GaugeGrid);
    
     Lattice<iScalar<vInteger> > coor(GaugeGrid);
-        
-    for(int mu=0;mu<Nd;mu++){
-          
-      LatticeCoordinate(coor,mu);
+
+    //Here the first Nd-1 directions are treated as "spatial", and a twist value of 1 indicates G-parity BCs in that direction. 
+    //mu=Nd-1 is assumed to be the time direction and a twist value of 1 indicates antiperiodic BCs        
+    for(int mu=0;mu<Nd-1;mu++){
+
+      if( Params.twists[mu] ){
+	LatticeCoordinate(coor,mu);
+      }
           
       U     = PeekIndex<LorentzIndex>(Umu,mu);
       Uconj = conjugate(U);
      
+      // Implement the isospin rotation sign on the boundary between f=1 and f=0
       // This phase could come from a simple bc 1,1,-1,1 ..
       int neglink = GaugeGrid->GlobalDimensions()[mu]-1;
       if ( Params.twists[mu] ) { 
@@ -229,7 +259,7 @@ public:
 	thread_foreach(ss,U_v,{
 	    Uds_v[ss](0)(mu) = U_v[ss]();
 	    Uds_v[ss](1)(mu) = Uconj_v[ss]();
-	  });
+	});
       }
           
       U     = adj(Cshift(U    ,mu,-1));      // correct except for spanning the boundary
@@ -260,6 +290,38 @@ public:
         });
       }
     }
+
+    { //periodic / antiperiodic temporal BCs
+      int mu = Nd-1;
+      int L   = GaugeGrid->GlobalDimensions()[mu];
+      int Lmu = L - 1;
+
+      LatticeCoordinate(coor, mu);
+
+      U = PeekIndex<LorentzIndex>(Umu, mu); //Get t-directed links
+      
+      GaugeLinkField *Upoke = &U;
+
+      if(Params.twists[mu]){ //antiperiodic
+	Utmp =  where(coor == Lmu, -U, U);
+	Upoke = &Utmp;
+      }
+    
+      Uconj = conjugate(*Upoke); //second flavor interacts with conjugate links      
+      pokeGparityDoubledGaugeField(Uds, *Upoke, Uconj, mu);
+
+      //Get the barrel-shifted field
+      Utmp = adj(Cshift(U, mu, -1)); //is a forward shift!
+      Upoke = &Utmp;
+
+      if(Params.twists[mu]){
+	U = where(coor == 0, -Utmp, Utmp);  //boundary phase
+	Upoke = &U;
+      }
+      
+      Uconj = conjugate(*Upoke);
+      pokeGparityDoubledGaugeField(Uds, *Upoke, Uconj, mu + 4);
+    }
   }
       
   inline void InsertForce4D(GaugeField &mat, FermionField &Btilde, FermionField &A, int mu) {
@@ -298,28 +360,48 @@ public:
   inline void extractLinkField(std::vector<GaugeLinkField> &mat, DoubledGaugeField &Uds){
     assert(0);
   }
-  
+ 
   inline void InsertForce5D(GaugeField &mat, FermionField &Btilde, FermionField &Atilde, int mu) {
-
-    int Ls = Btilde.Grid()->_fdimensions[0];
-        
-    GaugeLinkField tmp(mat.Grid());
-    tmp = Zero();
+    int Ls=Btilde.Grid()->_fdimensions[0];
+    
     {
-      autoView( tmp_v , tmp, CpuWrite);
-      autoView( Atilde_v , Atilde, CpuRead);
-      autoView( Btilde_v , Btilde, CpuRead);
-      thread_for(ss,tmp.Grid()->oSites(),{
-	  for (int s = 0; s < Ls; s++) {
-	    int sF = s + Ls * ss;
-	    auto ttmp = traceIndex<SpinIndex>(outerProduct(Btilde_v[sF], Atilde_v[sF]));
-	    tmp_v[ss]() = tmp_v[ss]() + ttmp(0, 0) + conjugate(ttmp(1, 1));
-	  }
-	});
+      GridBase *GaugeGrid = mat.Grid();
+      Lattice<iScalar<vInteger> > coor(GaugeGrid);
+
+      if( Params.twists[mu] ){
+	LatticeCoordinate(coor,mu);
+      }
+
+      autoView( mat_v , mat, AcceleratorWrite);
+      autoView( Btilde_v , Btilde, AcceleratorRead);
+      autoView( Atilde_v , Atilde, AcceleratorRead);
+      accelerator_for(sss,mat.Grid()->oSites(), FermionField::vector_type::Nsimd(),{	  
+  	  int sU=sss;
+  	  typedef decltype(coalescedRead(mat_v[sU](mu)() )) ColorMatrixType;
+  	  ColorMatrixType sum;
+  	  zeroit(sum);
+  	  for(int s=0;s<Ls;s++){
+  	    int sF = s+Ls*sU;
+  	    for(int spn=0;spn<Ns;spn++){ //sum over spin
+	      //Flavor 0
+  	      auto bb = coalescedRead(Btilde_v[sF](0)(spn) ); //color vector
+  	      auto aa = coalescedRead(Atilde_v[sF](0)(spn) );
+  	      sum = sum + outerProduct(bb,aa);
+
+  	      //Flavor 1
+  	      bb = coalescedRead(Btilde_v[sF](1)(spn) );
+  	      aa = coalescedRead(Atilde_v[sF](1)(spn) );
+  	      sum = sum + conjugate(outerProduct(bb,aa));
+  	    }
+  	  }	    
+  	  coalescedWrite(mat_v[sU](mu)(), sum);
+  	});
     }
-    PokeIndex<LorentzIndex>(mat, tmp, mu);
-    return;
   }
+
+
+  
+
   
 };
 
diff --git a/Grid/qcd/action/filters/MomentumFilter.h b/Grid/qcd/action/filters/MomentumFilter.h
index 2a15d80c..864166f5 100644
--- a/Grid/qcd/action/filters/MomentumFilter.h
+++ b/Grid/qcd/action/filters/MomentumFilter.h
@@ -37,7 +37,7 @@ NAMESPACE_BEGIN(Grid);
 
 template<typename MomentaField>
 struct MomentumFilterBase{
-  virtual void applyFilter(MomentaField &P) const;
+  virtual void applyFilter(MomentaField &P) const = 0;
 };
 
 //Do nothing
diff --git a/Grid/qcd/action/gauge/GaugeImplementations.h b/Grid/qcd/action/gauge/GaugeImplementations.h
index 16147c77..f518b236 100644
--- a/Grid/qcd/action/gauge/GaugeImplementations.h
+++ b/Grid/qcd/action/gauge/GaugeImplementations.h
@@ -69,6 +69,11 @@ public:
     return PeriodicBC::ShiftStaple(Link,mu);
   }
 
+  //Same as Cshift for periodic BCs
+  static inline GaugeLinkField CshiftLink(const GaugeLinkField &Link, int mu, int shift){
+    return PeriodicBC::CshiftLink(Link,mu,shift);
+  }
+
   static inline bool isPeriodicGaugeField(void) { return true; }
 };
 
@@ -110,6 +115,11 @@ public:
       return PeriodicBC::CovShiftBackward(Link, mu, field);
   }
 
+  //If mu is a conjugate BC direction
+  //Out(x) = U^dag_\mu(x-mu)  | x_\mu != 0
+  //       = U^T_\mu(L-1)  | x_\mu == 0
+  //else
+  //Out(x) = U^dag_\mu(x-mu mod L)
   static inline GaugeLinkField
   CovShiftIdentityBackward(const GaugeLinkField &Link, int mu)
   {
@@ -129,6 +139,13 @@ public:
       return PeriodicBC::CovShiftIdentityForward(Link,mu);
   }
 
+
+  //If mu is a conjugate BC direction
+  //Out(x) = S_\mu(x+mu)  | x_\mu != L-1
+  //       = S*_\mu(x+mu)  | x_\mu == L-1
+  //else
+  //Out(x) = S_\mu(x+mu mod L)
+  //Note: While this is used for Staples it is also applicable for shifting gauge links or gauge transformation matrices
   static inline GaugeLinkField ShiftStaple(const GaugeLinkField &Link, int mu)
   {
     assert(_conjDirs.size() == Nd);
@@ -138,6 +155,27 @@ public:
       return PeriodicBC::ShiftStaple(Link,mu);
   }
 
+  //Boundary-aware C-shift of gauge links / gauge transformation matrices
+  //For conjugate BC direction
+  //shift = 1
+  //Out(x) = U_\mu(x+\hat\mu)  | x_\mu != L-1
+  //       = U*_\mu(0)  | x_\mu == L-1
+  //shift = -1
+  //Out(x) = U_\mu(x-mu)  | x_\mu != 0
+  //       = U*_\mu(L-1)  | x_\mu == 0
+  //else
+  //shift = 1
+  //Out(x) = U_\mu(x+\hat\mu mod L)
+  //shift = -1
+  //Out(x) = U_\mu(x-\hat\mu mod L)
+  static inline GaugeLinkField CshiftLink(const GaugeLinkField &Link, int mu, int shift){
+    assert(_conjDirs.size() == Nd);
+    if(_conjDirs[mu]) 
+      return ConjugateBC::CshiftLink(Link,mu,shift);
+    else     
+      return PeriodicBC::CshiftLink(Link,mu,shift);
+  }
+
   static inline void       setDirections(std::vector<int> &conjDirs) { _conjDirs=conjDirs; }
   static inline std::vector<int> getDirections(void) { return _conjDirs; }
   static inline bool isPeriodicGaugeField(void) { return false; }
diff --git a/Grid/qcd/gparity/Gparity.h b/Grid/qcd/gparity/Gparity.h
new file mode 100644
index 00000000..ce1c70eb
--- /dev/null
+++ b/Grid/qcd/gparity/Gparity.h
@@ -0,0 +1,6 @@
+#ifndef GRID_GPARITY_H_
+#define GRID_GPARITY_H_
+
+#include<Grid/qcd/gparity/GparityFlavour.h>
+
+#endif
diff --git a/Grid/qcd/gparity/GparityFlavour.cc b/Grid/qcd/gparity/GparityFlavour.cc
new file mode 100644
index 00000000..4596f96b
--- /dev/null
+++ b/Grid/qcd/gparity/GparityFlavour.cc
@@ -0,0 +1,34 @@
+#include <Grid/Grid.h>
+
+NAMESPACE_BEGIN(Grid);
+
+const std::array<const GparityFlavour, 3> GparityFlavour::sigma_mu = {{
+    GparityFlavour(GparityFlavour::Algebra::SigmaX),
+    GparityFlavour(GparityFlavour::Algebra::SigmaY),
+    GparityFlavour(GparityFlavour::Algebra::SigmaZ)
+    }};
+
+const std::array<const GparityFlavour, 6> GparityFlavour::sigma_all = {{
+  GparityFlavour(GparityFlavour::Algebra::Identity),
+  GparityFlavour(GparityFlavour::Algebra::SigmaX),
+  GparityFlavour(GparityFlavour::Algebra::SigmaY),
+  GparityFlavour(GparityFlavour::Algebra::SigmaZ),
+  GparityFlavour(GparityFlavour::Algebra::ProjPlus),
+  GparityFlavour(GparityFlavour::Algebra::ProjMinus)
+}};
+
+const std::array<const char *, GparityFlavour::nSigma> GparityFlavour::name = {{
+    "SigmaX",
+    "MinusSigmaX",
+    "SigmaY",
+    "MinusSigmaY",
+    "SigmaZ",
+    "MinusSigmaZ",
+    "Identity",
+    "MinusIdentity",
+    "ProjPlus",
+    "MinusProjPlus",
+    "ProjMinus",
+    "MinusProjMinus"}};
+
+NAMESPACE_END(Grid);
diff --git a/Grid/qcd/gparity/GparityFlavour.h b/Grid/qcd/gparity/GparityFlavour.h
new file mode 100644
index 00000000..b2009235
--- /dev/null
+++ b/Grid/qcd/gparity/GparityFlavour.h
@@ -0,0 +1,475 @@
+#ifndef GRID_QCD_GPARITY_FLAVOUR_H
+#define GRID_QCD_GPARITY_FLAVOUR_H
+
+//Support for flavour-matrix operations acting on the G-parity flavour index
+
+#include <array>
+
+NAMESPACE_BEGIN(Grid);
+
+class GparityFlavour {
+  public:
+    GRID_SERIALIZABLE_ENUM(Algebra, undef,
+                           SigmaX, 0,
+			   MinusSigmaX, 1,
+                           SigmaY, 2,
+			   MinusSigmaY, 3,
+                           SigmaZ, 4,
+			   MinusSigmaZ, 5,
+			   Identity, 6,
+			   MinusIdentity, 7,
+			   ProjPlus, 8,
+			   MinusProjPlus, 9,
+			   ProjMinus, 10,
+			   MinusProjMinus, 11
+			   );
+    static constexpr unsigned int nSigma = 12;
+    static const std::array<const char *, nSigma>                name;
+    static const std::array<const GparityFlavour, 3>             sigma_mu;
+    static const std::array<const GparityFlavour, 6>            sigma_all;
+    Algebra                                                      g;
+  public:
+  accelerator GparityFlavour(Algebra initg): g(initg) {}  
+};
+
+
+
+// 0 1  x   vector
+// 1 0
+template<class vtype>
+accelerator_inline void multFlavourSigmaX(iVector<vtype, Ngp> &ret, const iVector<vtype, Ngp> &rhs)
+{
+  ret(0) = rhs(1);
+  ret(1) = rhs(0);
+};
+template<class vtype>
+accelerator_inline void lmultFlavourSigmaX(iMatrix<vtype, Ngp> &ret, const iMatrix<vtype, Ngp> &rhs)
+{
+  ret(0,0) = rhs(1,0);
+  ret(0,1) = rhs(1,1);
+  ret(1,0) = rhs(0,0);
+  ret(1,1) = rhs(0,1);
+};
+template<class vtype>
+accelerator_inline void rmultFlavourSigmaX(iMatrix<vtype, Ngp> &ret, const iMatrix<vtype, Ngp> &rhs)
+{
+  ret(0,0) = rhs(0,1);
+  ret(0,1) = rhs(0,0);
+  ret(1,0) = rhs(1,1);
+  ret(1,1) = rhs(1,0);
+};
+
+
+template<class vtype>
+accelerator_inline void multFlavourMinusSigmaX(iVector<vtype, Ngp> &ret, const iVector<vtype, Ngp> &rhs)
+{
+  ret(0) = -rhs(1);
+  ret(1) = -rhs(0);
+};
+template<class vtype>
+accelerator_inline void lmultFlavourMinusSigmaX(iMatrix<vtype, Ngp> &ret, const iMatrix<vtype, Ngp> &rhs)
+{
+  ret(0,0) = -rhs(1,0);
+  ret(0,1) = -rhs(1,1);
+  ret(1,0) = -rhs(0,0);
+  ret(1,1) = -rhs(0,1);
+};
+template<class vtype>
+accelerator_inline void rmultFlavourMinusSigmaX(iMatrix<vtype, Ngp> &ret, const iMatrix<vtype, Ngp> &rhs)
+{
+  ret(0,0) = -rhs(0,1);
+  ret(0,1) = -rhs(0,0);
+  ret(1,0) = -rhs(1,1);
+  ret(1,1) = -rhs(1,0);
+};
+
+
+
+
+
+// 0 -i  x   vector
+// i 0
+template<class vtype>
+accelerator_inline void multFlavourSigmaY(iVector<vtype, Ngp> &ret, const iVector<vtype, Ngp> &rhs)
+{
+  ret(0) = timesMinusI(rhs(1));
+  ret(1) = timesI(rhs(0));
+};
+template<class vtype>
+accelerator_inline void lmultFlavourSigmaY(iMatrix<vtype, Ngp> &ret, const iMatrix<vtype, Ngp> &rhs)
+{
+  ret(0,0) = timesMinusI(rhs(1,0));
+  ret(0,1) = timesMinusI(rhs(1,1));
+  ret(1,0) = timesI(rhs(0,0));
+  ret(1,1) = timesI(rhs(0,1));
+};
+template<class vtype>
+accelerator_inline void rmultFlavourSigmaY(iMatrix<vtype, Ngp> &ret, const iMatrix<vtype, Ngp> &rhs)
+{
+  ret(0,0) = timesI(rhs(0,1));
+  ret(0,1) = timesMinusI(rhs(0,0));
+  ret(1,0) = timesI(rhs(1,1));
+  ret(1,1) = timesMinusI(rhs(1,0));
+};
+
+template<class vtype>
+accelerator_inline void multFlavourMinusSigmaY(iVector<vtype, Ngp> &ret, const iVector<vtype, Ngp> &rhs)
+{
+  ret(0) = timesI(rhs(1));
+  ret(1) = timesMinusI(rhs(0));
+};
+template<class vtype>
+accelerator_inline void lmultFlavourMinusSigmaY(iMatrix<vtype, Ngp> &ret, const iMatrix<vtype, Ngp> &rhs)
+{
+  ret(0,0) = timesI(rhs(1,0));
+  ret(0,1) = timesI(rhs(1,1));
+  ret(1,0) = timesMinusI(rhs(0,0));
+  ret(1,1) = timesMinusI(rhs(0,1));
+};
+template<class vtype>
+accelerator_inline void rmultFlavourMinusSigmaY(iMatrix<vtype, Ngp> &ret, const iMatrix<vtype, Ngp> &rhs)
+{
+  ret(0,0) = timesMinusI(rhs(0,1));
+  ret(0,1) = timesI(rhs(0,0));
+  ret(1,0) = timesMinusI(rhs(1,1));
+  ret(1,1) = timesI(rhs(1,0));
+};
+
+
+
+
+
+// 1 0  x   vector
+// 0 -1
+template<class vtype>
+accelerator_inline void multFlavourSigmaZ(iVector<vtype, Ngp> &ret, const iVector<vtype, Ngp> &rhs)
+{
+  ret(0) = rhs(0);
+  ret(1) = -rhs(1);
+};
+template<class vtype>
+accelerator_inline void lmultFlavourSigmaZ(iMatrix<vtype, Ngp> &ret, const iMatrix<vtype, Ngp> &rhs)
+{
+  ret(0,0) = rhs(0,0);
+  ret(0,1) = rhs(0,1);
+  ret(1,0) = -rhs(1,0);
+  ret(1,1) = -rhs(1,1);
+};
+template<class vtype>
+accelerator_inline void rmultFlavourSigmaZ(iMatrix<vtype, Ngp> &ret, const iMatrix<vtype, Ngp> &rhs)
+{
+  ret(0,0) = rhs(0,0);
+  ret(0,1) = -rhs(0,1);
+  ret(1,0) = rhs(1,0);
+  ret(1,1) = -rhs(1,1);
+};
+
+
+template<class vtype>
+accelerator_inline void multFlavourMinusSigmaZ(iVector<vtype, Ngp> &ret, const iVector<vtype, Ngp> &rhs)
+{
+  ret(0) = -rhs(0);
+  ret(1) = rhs(1);
+};
+template<class vtype>
+accelerator_inline void lmultFlavourMinusSigmaZ(iMatrix<vtype, Ngp> &ret, const iMatrix<vtype, Ngp> &rhs)
+{
+  ret(0,0) = -rhs(0,0);
+  ret(0,1) = -rhs(0,1);
+  ret(1,0) = rhs(1,0);
+  ret(1,1) = rhs(1,1);
+};
+template<class vtype>
+accelerator_inline void rmultFlavourMinusSigmaZ(iMatrix<vtype, Ngp> &ret, const iMatrix<vtype, Ngp> &rhs)
+{
+  ret(0,0) = -rhs(0,0);
+  ret(0,1) = rhs(0,1);
+  ret(1,0) = -rhs(1,0);
+  ret(1,1) = rhs(1,1);
+};
+
+
+
+
+
+
+template<class vtype>
+accelerator_inline void multFlavourIdentity(iVector<vtype, Ngp> &ret, const iVector<vtype, Ngp> &rhs)
+{
+  ret(0) = rhs(0);
+  ret(1) = rhs(1);
+};
+template<class vtype>
+accelerator_inline void lmultFlavourIdentity(iMatrix<vtype, Ngp> &ret, const iMatrix<vtype, Ngp> &rhs)
+{
+  ret(0,0) = rhs(0,0);
+  ret(0,1) = rhs(0,1);
+  ret(1,0) = rhs(1,0);
+  ret(1,1) = rhs(1,1);
+};
+template<class vtype>
+accelerator_inline void rmultFlavourIdentity(iMatrix<vtype, Ngp> &ret, const iMatrix<vtype, Ngp> &rhs)
+{
+  ret(0,0) = rhs(0,0);
+  ret(0,1) = rhs(0,1);
+  ret(1,0) = rhs(1,0);
+  ret(1,1) = rhs(1,1);
+};
+
+template<class vtype>
+accelerator_inline void multFlavourMinusIdentity(iVector<vtype, Ngp> &ret, const iVector<vtype, Ngp> &rhs)
+{
+  ret(0) = -rhs(0);
+  ret(1) = -rhs(1);
+};
+template<class vtype>
+accelerator_inline void lmultFlavourMinusIdentity(iMatrix<vtype, Ngp> &ret, const iMatrix<vtype, Ngp> &rhs)
+{
+  ret(0,0) = -rhs(0,0);
+  ret(0,1) = -rhs(0,1);
+  ret(1,0) = -rhs(1,0);
+  ret(1,1) = -rhs(1,1);
+};
+template<class vtype>
+accelerator_inline void rmultFlavourMinusIdentity(iMatrix<vtype, Ngp> &ret, const iMatrix<vtype, Ngp> &rhs)
+{
+  ret(0,0) = -rhs(0,0);
+  ret(0,1) = -rhs(0,1);
+  ret(1,0) = -rhs(1,0);
+  ret(1,1) = -rhs(1,1);
+};
+
+
+
+
+
+//G-parity flavour projection 1/2(1+\sigma_2)
+//1 -i
+//i  1
+template<class vtype>
+accelerator_inline void multFlavourProjPlus(iVector<vtype, Ngp> &ret, const iVector<vtype, Ngp> &rhs)
+{
+  ret(0) = 0.5*rhs(0) + 0.5*timesMinusI(rhs(1));
+  ret(1) = 0.5*timesI(rhs(0)) + 0.5*rhs(1);
+};
+template<class vtype>
+accelerator_inline void lmultFlavourProjPlus(iMatrix<vtype, Ngp> &ret, const iMatrix<vtype, Ngp> &rhs)
+{
+  ret(0,0) = 0.5*rhs(0,0) + 0.5*timesMinusI(rhs(1,0));
+  ret(0,1) = 0.5*rhs(0,1) + 0.5*timesMinusI(rhs(1,1));
+  ret(1,0) = 0.5*timesI(rhs(0,0)) + 0.5*rhs(1,0);
+  ret(1,1) = 0.5*timesI(rhs(0,1)) + 0.5*rhs(1,1);
+};
+template<class vtype>
+accelerator_inline void rmultFlavourProjPlus(iMatrix<vtype, Ngp> &ret, const iMatrix<vtype, Ngp> &rhs)
+{
+  ret(0,0) = 0.5*rhs(0,0) + 0.5*timesI(rhs(0,1));
+  ret(0,1) = 0.5*timesMinusI(rhs(0,0)) + 0.5*rhs(0,1);
+  ret(1,0) = 0.5*rhs(1,0) + 0.5*timesI(rhs(1,1));
+  ret(1,1) = 0.5*timesMinusI(rhs(1,0)) + 0.5*rhs(1,1);
+};
+
+
+template<class vtype>
+accelerator_inline void multFlavourMinusProjPlus(iVector<vtype, Ngp> &ret, const iVector<vtype, Ngp> &rhs)
+{
+  ret(0) = -0.5*rhs(0) + 0.5*timesI(rhs(1));
+  ret(1) = 0.5*timesMinusI(rhs(0)) - 0.5*rhs(1);
+};
+template<class vtype>
+accelerator_inline void lmultFlavourMinusProjPlus(iMatrix<vtype, Ngp> &ret, const iMatrix<vtype, Ngp> &rhs)
+{
+  ret(0,0) = -0.5*rhs(0,0) + 0.5*timesI(rhs(1,0));
+  ret(0,1) = -0.5*rhs(0,1) + 0.5*timesI(rhs(1,1));
+  ret(1,0) = 0.5*timesMinusI(rhs(0,0)) - 0.5*rhs(1,0);
+  ret(1,1) = 0.5*timesMinusI(rhs(0,1)) - 0.5*rhs(1,1);
+};
+template<class vtype>
+accelerator_inline void rmultFlavourMinusProjPlus(iMatrix<vtype, Ngp> &ret, const iMatrix<vtype, Ngp> &rhs)
+{
+  ret(0,0) = -0.5*rhs(0,0) + 0.5*timesMinusI(rhs(0,1));
+  ret(0,1) = 0.5*timesI(rhs(0,0)) - 0.5*rhs(0,1);
+  ret(1,0) = -0.5*rhs(1,0) + 0.5*timesMinusI(rhs(1,1));
+  ret(1,1) = 0.5*timesI(rhs(1,0)) - 0.5*rhs(1,1);
+};
+
+
+
+
+
+//G-parity flavour projection 1/2(1-\sigma_2)
+//1 i
+//-i  1
+template<class vtype>
+accelerator_inline void multFlavourProjMinus(iVector<vtype, Ngp> &ret, const iVector<vtype, Ngp> &rhs)
+{
+  ret(0) = 0.5*rhs(0) + 0.5*timesI(rhs(1));
+  ret(1) = 0.5*timesMinusI(rhs(0)) + 0.5*rhs(1);
+};
+template<class vtype>
+accelerator_inline void lmultFlavourProjMinus(iMatrix<vtype, Ngp> &ret, const iMatrix<vtype, Ngp> &rhs)
+{
+  ret(0,0) = 0.5*rhs(0,0) + 0.5*timesI(rhs(1,0));
+  ret(0,1) = 0.5*rhs(0,1) + 0.5*timesI(rhs(1,1));
+  ret(1,0) = 0.5*timesMinusI(rhs(0,0)) + 0.5*rhs(1,0);
+  ret(1,1) = 0.5*timesMinusI(rhs(0,1)) + 0.5*rhs(1,1);
+};
+template<class vtype>
+accelerator_inline void rmultFlavourProjMinus(iMatrix<vtype, Ngp> &ret, const iMatrix<vtype, Ngp> &rhs)
+{
+  ret(0,0) = 0.5*rhs(0,0) + 0.5*timesMinusI(rhs(0,1));
+  ret(0,1) = 0.5*timesI(rhs(0,0)) + 0.5*rhs(0,1);
+  ret(1,0) = 0.5*rhs(1,0) + 0.5*timesMinusI(rhs(1,1));
+  ret(1,1) = 0.5*timesI(rhs(1,0)) + 0.5*rhs(1,1);
+};
+
+
+template<class vtype>
+accelerator_inline void multFlavourMinusProjMinus(iVector<vtype, Ngp> &ret, const iVector<vtype, Ngp> &rhs)
+{
+  ret(0) = -0.5*rhs(0) + 0.5*timesMinusI(rhs(1));
+  ret(1) = 0.5*timesI(rhs(0)) - 0.5*rhs(1);
+};
+template<class vtype>
+accelerator_inline void lmultFlavourMinusProjMinus(iMatrix<vtype, Ngp> &ret, const iMatrix<vtype, Ngp> &rhs)
+{
+  ret(0,0) = -0.5*rhs(0,0) + 0.5*timesMinusI(rhs(1,0));
+  ret(0,1) = -0.5*rhs(0,1) + 0.5*timesMinusI(rhs(1,1));
+  ret(1,0) = 0.5*timesI(rhs(0,0)) - 0.5*rhs(1,0);
+  ret(1,1) = 0.5*timesI(rhs(0,1)) - 0.5*rhs(1,1);
+};
+template<class vtype>
+accelerator_inline void rmultFlavourMinusProjMinus(iMatrix<vtype, Ngp> &ret, const iMatrix<vtype, Ngp> &rhs)
+{
+  ret(0,0) = -0.5*rhs(0,0) + 0.5*timesI(rhs(0,1));
+  ret(0,1) = 0.5*timesMinusI(rhs(0,0)) - 0.5*rhs(0,1);
+  ret(1,0) = -0.5*rhs(1,0) + 0.5*timesI(rhs(1,1));
+  ret(1,1) = 0.5*timesMinusI(rhs(1,0)) - 0.5*rhs(1,1);
+};
+
+
+
+
+
+
+
+
+
+
+template<class vtype> 
+accelerator_inline auto operator*(const GparityFlavour &G, const iVector<vtype, Ngp> &arg)
+->typename std::enable_if<matchGridTensorIndex<iVector<vtype, Ngp>, GparityFlavourTensorIndex>::value, iVector<vtype, Ngp>>::type
+{
+  iVector<vtype, Ngp> ret;
+
+  switch (G.g) 
+  {
+  case GparityFlavour::Algebra::SigmaX:
+    multFlavourSigmaX(ret, arg); break;
+  case GparityFlavour::Algebra::MinusSigmaX:
+    multFlavourMinusSigmaX(ret, arg); break;
+  case GparityFlavour::Algebra::SigmaY:
+    multFlavourSigmaY(ret, arg); break;
+  case GparityFlavour::Algebra::MinusSigmaY:
+    multFlavourMinusSigmaY(ret, arg); break;
+  case GparityFlavour::Algebra::SigmaZ:
+    multFlavourSigmaZ(ret, arg); break;
+  case GparityFlavour::Algebra::MinusSigmaZ:
+    multFlavourMinusSigmaZ(ret, arg); break;
+  case GparityFlavour::Algebra::Identity:
+    multFlavourIdentity(ret, arg); break;
+  case GparityFlavour::Algebra::MinusIdentity:
+    multFlavourMinusIdentity(ret, arg); break;
+  case GparityFlavour::Algebra::ProjPlus:
+    multFlavourProjPlus(ret, arg); break;
+  case GparityFlavour::Algebra::MinusProjPlus:
+    multFlavourMinusProjPlus(ret, arg); break;
+  case GparityFlavour::Algebra::ProjMinus:
+    multFlavourProjMinus(ret, arg); break;
+  case GparityFlavour::Algebra::MinusProjMinus:
+    multFlavourMinusProjMinus(ret, arg); break;
+  default: assert(0);
+  }
+ 
+  return ret;
+}
+
+template<class vtype> 
+accelerator_inline auto operator*(const GparityFlavour &G, const iMatrix<vtype, Ngp> &arg)
+->typename std::enable_if<matchGridTensorIndex<iMatrix<vtype, Ngp>, GparityFlavourTensorIndex>::value, iMatrix<vtype, Ngp>>::type
+{
+  iMatrix<vtype, Ngp> ret;
+
+  switch (G.g) 
+  {
+  case GparityFlavour::Algebra::SigmaX:
+    lmultFlavourSigmaX(ret, arg); break;
+  case GparityFlavour::Algebra::MinusSigmaX:
+    lmultFlavourMinusSigmaX(ret, arg); break;
+  case GparityFlavour::Algebra::SigmaY:
+    lmultFlavourSigmaY(ret, arg); break;
+  case GparityFlavour::Algebra::MinusSigmaY:
+    lmultFlavourMinusSigmaY(ret, arg); break;
+  case GparityFlavour::Algebra::SigmaZ:
+    lmultFlavourSigmaZ(ret, arg); break;
+  case GparityFlavour::Algebra::MinusSigmaZ:
+    lmultFlavourMinusSigmaZ(ret, arg); break;
+  case GparityFlavour::Algebra::Identity:
+    lmultFlavourIdentity(ret, arg); break;
+  case GparityFlavour::Algebra::MinusIdentity:
+    lmultFlavourMinusIdentity(ret, arg); break;
+  case GparityFlavour::Algebra::ProjPlus:
+    lmultFlavourProjPlus(ret, arg); break;
+  case GparityFlavour::Algebra::MinusProjPlus:
+    lmultFlavourMinusProjPlus(ret, arg); break;
+  case GparityFlavour::Algebra::ProjMinus:
+    lmultFlavourProjMinus(ret, arg); break;
+  case GparityFlavour::Algebra::MinusProjMinus:
+    lmultFlavourMinusProjMinus(ret, arg); break;  
+  default: assert(0);
+  }
+  
+  return ret;
+}
+
+template<class vtype> 
+accelerator_inline auto operator*(const iMatrix<vtype, Ngp> &arg, const GparityFlavour &G)
+->typename std::enable_if<matchGridTensorIndex<iMatrix<vtype, Ngp>, GparityFlavourTensorIndex>::value, iMatrix<vtype, Ngp>>::type
+{
+  iMatrix<vtype, Ngp> ret;
+
+  switch (G.g) 
+  {
+  case GparityFlavour::Algebra::SigmaX:
+    rmultFlavourSigmaX(ret, arg); break;
+  case GparityFlavour::Algebra::MinusSigmaX:
+    rmultFlavourMinusSigmaX(ret, arg); break;
+  case GparityFlavour::Algebra::SigmaY:
+    rmultFlavourSigmaY(ret, arg); break;
+  case GparityFlavour::Algebra::MinusSigmaY:
+    rmultFlavourMinusSigmaY(ret, arg); break;
+  case GparityFlavour::Algebra::SigmaZ:
+    rmultFlavourSigmaZ(ret, arg); break;
+  case GparityFlavour::Algebra::MinusSigmaZ:
+    rmultFlavourMinusSigmaZ(ret, arg); break;
+  case GparityFlavour::Algebra::Identity:
+    rmultFlavourIdentity(ret, arg); break;
+  case GparityFlavour::Algebra::MinusIdentity:
+    rmultFlavourMinusIdentity(ret, arg); break;
+  case GparityFlavour::Algebra::ProjPlus:
+    rmultFlavourProjPlus(ret, arg); break;
+  case GparityFlavour::Algebra::MinusProjPlus:
+    rmultFlavourMinusProjPlus(ret, arg); break;
+  case GparityFlavour::Algebra::ProjMinus:
+    rmultFlavourProjMinus(ret, arg); break;
+  case GparityFlavour::Algebra::MinusProjMinus:
+    rmultFlavourMinusProjMinus(ret, arg); break;
+  default: assert(0);
+  }
+
+  return ret;
+}
+
+NAMESPACE_END(Grid);
+
+#endif // include guard
diff --git a/Grid/qcd/utils/CovariantCshift.h b/Grid/qcd/utils/CovariantCshift.h
index 6c70706f..79cf8e0f 100644
--- a/Grid/qcd/utils/CovariantCshift.h
+++ b/Grid/qcd/utils/CovariantCshift.h
@@ -88,6 +88,12 @@ namespace PeriodicBC {
     return CovShiftBackward(Link,mu,arg);
   }
 
+  //Boundary-aware C-shift of gauge links / gauge transformation matrices
+  template<class gauge> Lattice<gauge>
+  CshiftLink(const Lattice<gauge> &Link, int mu, int shift)
+  {
+    return Cshift(Link, mu, shift);
+  }
 
 }
 
@@ -158,6 +164,9 @@ namespace ConjugateBC {
     //    std::cout<<"Gparity::CovCshiftBackward mu="<<mu<<std::endl;
     return Cshift(tmp,mu,-1);// moves towards positive mu
   }
+
+  //Out(x) = U^dag_\mu(x-mu)  | x_\mu != 0
+  //       = U^T_\mu(L-1)  | x_\mu == 0
   template<class gauge> Lattice<gauge>
   CovShiftIdentityBackward(const Lattice<gauge> &Link, int mu) {
     GridBase *grid = Link.Grid();
@@ -176,6 +185,9 @@ namespace ConjugateBC {
     return Link;
   }
 
+  //Out(x) = S_\mu(x+\hat\mu)  | x_\mu != L-1
+  //       = S*_\mu(0)  | x_\mu == L-1
+  //Note: While this is used for Staples it is also applicable for shifting gauge links or gauge transformation matrices
   template<class gauge> Lattice<gauge>
   ShiftStaple(const Lattice<gauge> &Link, int mu)
   {
@@ -208,6 +220,35 @@ namespace ConjugateBC {
     return CovShiftBackward(Link,mu,arg);
   }
 
+  //Boundary-aware C-shift of gauge links / gauge transformation matrices
+  //shift = 1
+  //Out(x) = U_\mu(x+\hat\mu)  | x_\mu != L-1
+  //       = U*_\mu(0)  | x_\mu == L-1
+  //shift = -1
+  //Out(x) = U_\mu(x-mu)  | x_\mu != 0
+  //       = U*_\mu(L-1)  | x_\mu == 0
+  template<class gauge> Lattice<gauge>
+  CshiftLink(const Lattice<gauge> &Link, int mu, int shift)
+  {
+    GridBase *grid = Link.Grid();
+    int Lmu = grid->GlobalDimensions()[mu] - 1;
+
+    Lattice<iScalar<vInteger>> coor(grid);
+    LatticeCoordinate(coor, mu);
+
+    Lattice<gauge> tmp(grid);
+    if(shift == 1){
+      tmp = Cshift(Link, mu, 1);
+      tmp = where(coor == Lmu, conjugate(tmp), tmp);
+      return tmp;
+    }else if(shift == -1){
+      tmp = Link;
+      tmp = where(coor == Lmu, conjugate(tmp), tmp);
+      return Cshift(tmp, mu, -1);
+    }else assert(0 && "Invalid shift value");
+    return tmp; //shuts up the compiler fussing about the return type
+  }
+
 }
 
 
diff --git a/Grid/qcd/utils/GaugeFix.h b/Grid/qcd/utils/GaugeFix.h
index 2b3384da..c0bc2c83 100644
--- a/Grid/qcd/utils/GaugeFix.h
+++ b/Grid/qcd/utils/GaugeFix.h
@@ -40,27 +40,46 @@ public:
   typedef typename Gimpl::GaugeLinkField GaugeMat;
   typedef typename Gimpl::GaugeField GaugeLorentz;
 
-  static void GaugeLinkToLieAlgebraField(const std::vector<GaugeMat> &U,std::vector<GaugeMat> &A) {
-    for(int mu=0;mu<Nd;mu++){
-      Complex cmi(0.0,-1.0);
-      A[mu] = Ta(U[mu]) * cmi;
-    }
+  //A_\mu(x) = -i Ta(U_\mu(x) )   where Ta(U) = 1/2( U - U^dag ) - 1/2N tr(U - U^dag)  is the traceless antihermitian part. This is an O(A^3) approximation to the logarithm of U
+  static void GaugeLinkToLieAlgebraField(const GaugeMat &U, GaugeMat &A) {
+    Complex cmi(0.0,-1.0);
+    A = Ta(U) * cmi;
   }
-  static void DmuAmu(const std::vector<GaugeMat> &A,GaugeMat &dmuAmu,int orthog) {
+  
+  //The derivative of the Lie algebra field
+  static void DmuAmu(const std::vector<GaugeMat> &U, GaugeMat &dmuAmu,int orthog) {
+    GridBase* grid = U[0].Grid();
+    GaugeMat Ax(grid);
+    GaugeMat Axm1(grid);
+    GaugeMat Utmp(grid);
+
     dmuAmu=Zero();
     for(int mu=0;mu<Nd;mu++){
       if ( mu != orthog ) {
-	dmuAmu = dmuAmu + A[mu] - Cshift(A[mu],mu,-1);
+	//Rather than define functionality to work out how the BCs apply to A_\mu we simply use the BC-aware Cshift to the gauge links and compute A_\mu(x) and A_\mu(x-1) separately
+	//Ax = A_\mu(x)
+	GaugeLinkToLieAlgebraField(U[mu], Ax);
+	
+	//Axm1 = A_\mu(x_\mu-1)
+	Utmp = Gimpl::CshiftLink(U[mu], mu, -1);
+	GaugeLinkToLieAlgebraField(Utmp, Axm1);
+	
+	//Derivative
+	dmuAmu = dmuAmu + Ax - Axm1;
       }
     }
   }  
 
-  static void SteepestDescentGaugeFix(GaugeLorentz &Umu,Real & alpha,int maxiter,Real Omega_tol, Real Phi_tol,bool Fourier=false,int orthog=-1) {
+  //Fix the gauge field Umu
+  //0 < alpha < 1 is related to the step size, cf https://arxiv.org/pdf/1405.5812.pdf
+  static void SteepestDescentGaugeFix(GaugeLorentz &Umu, Real alpha,int maxiter,Real Omega_tol, Real Phi_tol,bool Fourier=false,int orthog=-1) {
     GridBase *grid = Umu.Grid();
     GaugeMat xform(grid);
     SteepestDescentGaugeFix(Umu,xform,alpha,maxiter,Omega_tol,Phi_tol,Fourier,orthog);
   }
-  static void SteepestDescentGaugeFix(GaugeLorentz &Umu,GaugeMat &xform,Real & alpha,int maxiter,Real Omega_tol, Real Phi_tol,bool Fourier=false,int orthog=-1) {
+
+  //Fix the gauge field Umu and also return the gauge transformation from the original gauge field, xform
+  static void SteepestDescentGaugeFix(GaugeLorentz &Umu,GaugeMat &xform, Real alpha,int maxiter,Real Omega_tol, Real Phi_tol,bool Fourier=false,int orthog=-1) {
 
     GridBase *grid = Umu.Grid();
 
@@ -122,27 +141,24 @@ public:
 
       }
     }
+    assert(0 && "Gauge fixing did not converge within the specified number of iterations");
   };
-  static Real SteepestDescentStep(std::vector<GaugeMat> &U,GaugeMat &xform,Real & alpha, GaugeMat & dmuAmu,int orthog) {
+  static Real SteepestDescentStep(std::vector<GaugeMat> &U,GaugeMat &xform, Real alpha, GaugeMat & dmuAmu,int orthog) {
     GridBase *grid = U[0].Grid();
 
-    std::vector<GaugeMat> A(Nd,grid);
     GaugeMat g(grid);
-
-    GaugeLinkToLieAlgebraField(U,A);
-    ExpiAlphaDmuAmu(A,g,alpha,dmuAmu,orthog);
-
+    ExpiAlphaDmuAmu(U,g,alpha,dmuAmu,orthog);
 
     Real vol = grid->gSites();
     Real trG = TensorRemove(sum(trace(g))).real()/vol/Nc;
 
     xform = g*xform ;
-    SU<Nc>::GaugeTransform(U,g);
+    SU<Nc>::GaugeTransform<Gimpl>(U,g);
 
     return trG;
   }
 
-  static Real FourierAccelSteepestDescentStep(std::vector<GaugeMat> &U,GaugeMat &xform,Real & alpha, GaugeMat & dmuAmu,int orthog) {
+  static Real FourierAccelSteepestDescentStep(std::vector<GaugeMat> &U,GaugeMat &xform, Real alpha, GaugeMat & dmuAmu,int orthog) {
 
     GridBase *grid = U[0].Grid();
 
@@ -157,11 +173,7 @@ public:
 
     GaugeMat g(grid);
     GaugeMat dmuAmu_p(grid);
-    std::vector<GaugeMat> A(Nd,grid);
-
-    GaugeLinkToLieAlgebraField(U,A);
-
-    DmuAmu(A,dmuAmu,orthog);
+    DmuAmu(U,dmuAmu,orthog);
 
     std::vector<int> mask(Nd,1);
     for(int mu=0;mu<Nd;mu++) if (mu==orthog) mask[mu]=0;
@@ -205,16 +217,16 @@ public:
     Real trG = TensorRemove(sum(trace(g))).real()/vol/Nc;
 
     xform = g*xform ;
-    SU<Nc>::GaugeTransform(U,g);
+    SU<Nc>::GaugeTransform<Gimpl>(U,g);
 
     return trG;
   }
 
-  static void ExpiAlphaDmuAmu(const std::vector<GaugeMat> &A,GaugeMat &g,Real & alpha, GaugeMat &dmuAmu,int orthog) {
+  static void ExpiAlphaDmuAmu(const std::vector<GaugeMat> &U,GaugeMat &g, Real alpha, GaugeMat &dmuAmu,int orthog) {
     GridBase *grid = g.Grid();
     Complex cialpha(0.0,-alpha);
     GaugeMat ciadmam(grid);
-    DmuAmu(A,dmuAmu,orthog);
+    DmuAmu(U,dmuAmu,orthog);
     ciadmam = dmuAmu*cialpha;
     SU<Nc>::taExp(ciadmam,g);
   }  
diff --git a/Grid/qcd/utils/SUn.h b/Grid/qcd/utils/SUn.h
index 675493b3..b9660c65 100644
--- a/Grid/qcd/utils/SUn.h
+++ b/Grid/qcd/utils/SUn.h
@@ -694,32 +694,32 @@ public:
  * Adjoint rep gauge xform
  */
 
-  template<typename GaugeField,typename GaugeMat>
-  static void GaugeTransform( GaugeField &Umu, GaugeMat &g){
+  template<typename Gimpl>
+  static void GaugeTransform(typename Gimpl::GaugeField &Umu, typename Gimpl::GaugeLinkField &g){
     GridBase *grid = Umu.Grid();
     conformable(grid,g.Grid());
 
-    GaugeMat U(grid);
-    GaugeMat ag(grid); ag = adj(g);
+    typename Gimpl::GaugeLinkField U(grid);
+    typename Gimpl::GaugeLinkField ag(grid); ag = adj(g);
 
     for(int mu=0;mu<Nd;mu++){
       U= PeekIndex<LorentzIndex>(Umu,mu);
-      U = g*U*Cshift(ag, mu, 1);
+      U = g*U*Gimpl::CshiftLink(ag, mu, 1); //BC-aware
       PokeIndex<LorentzIndex>(Umu,U,mu);
     }
   }
-  template<typename GaugeMat>
-  static void GaugeTransform( std::vector<GaugeMat> &U, GaugeMat &g){
+  template<typename Gimpl>
+  static void GaugeTransform( std::vector<typename Gimpl::GaugeLinkField> &U, typename Gimpl::GaugeLinkField &g){
     GridBase *grid = g.Grid();
-    GaugeMat ag(grid); ag = adj(g);
+    typename Gimpl::GaugeLinkField ag(grid); ag = adj(g);
     for(int mu=0;mu<Nd;mu++){
-      U[mu] = g*U[mu]*Cshift(ag, mu, 1);
+      U[mu] = g*U[mu]*Gimpl::CshiftLink(ag, mu, 1); //BC-aware
     }
   }
-  template<typename GaugeField,typename GaugeMat>
-  static void RandomGaugeTransform(GridParallelRNG &pRNG, GaugeField &Umu, GaugeMat &g){
+  template<typename Gimpl>
+  static void RandomGaugeTransform(GridParallelRNG &pRNG, typename Gimpl::GaugeField &Umu, typename Gimpl::GaugeLinkField &g){
     LieRandomize(pRNG,g,1.0);
-    GaugeTransform(Umu,g);
+    GaugeTransform<Gimpl>(Umu,g);
   }
 
   // Projects the algebra components a lattice matrix (of dimension ncol*ncol -1 )
diff --git a/Grid/qcd/utils/WilsonLoops.h b/Grid/qcd/utils/WilsonLoops.h
index 0367c9fa..da1f5ac8 100644
--- a/Grid/qcd/utils/WilsonLoops.h
+++ b/Grid/qcd/utils/WilsonLoops.h
@@ -125,6 +125,56 @@ public:
     return sumplaq / vol / faces / Nc; // Nd , Nc dependent... FIXME
   }
 
+  //////////////////////////////////////////////////
+  // sum over all spatial planes of plaquette
+  //////////////////////////////////////////////////
+  static void siteSpatialPlaquette(ComplexField &Plaq,
+                            const std::vector<GaugeMat> &U) {
+    ComplexField sitePlaq(U[0].Grid());
+    Plaq = Zero();
+    for (int mu = 1; mu < Nd-1; mu++) {
+      for (int nu = 0; nu < mu; nu++) {
+        traceDirPlaquette(sitePlaq, U, mu, nu);
+        Plaq = Plaq + sitePlaq;
+      }
+    }
+  }
+
+  ////////////////////////////////////
+  // sum over all x,y,z and over all spatial planes of plaquette
+  //////////////////////////////////////////////////
+  static std::vector<RealD> timesliceSumSpatialPlaquette(const GaugeLorentz &Umu) {
+    std::vector<GaugeMat> U(Nd, Umu.Grid());
+    // inefficient here
+    for (int mu = 0; mu < Nd; mu++) {
+      U[mu] = PeekIndex<LorentzIndex>(Umu, mu);
+    }
+
+    ComplexField Plaq(Umu.Grid());
+
+    siteSpatialPlaquette(Plaq, U);
+    typedef typename ComplexField::scalar_object sobj;
+    std::vector<sobj> Tq;
+    sliceSum(Plaq, Tq, Nd-1);
+
+    std::vector<Real> out(Tq.size());
+    for(int t=0;t<Tq.size();t++) out[t] = TensorRemove(Tq[t]).real();
+    return out;
+  }
+  
+  //////////////////////////////////////////////////
+  // average over all x,y,z and over all spatial planes of plaquette
+  //////////////////////////////////////////////////
+  static std::vector<RealD> timesliceAvgSpatialPlaquette(const GaugeLorentz &Umu) {
+    std::vector<RealD> sumplaq = timesliceSumSpatialPlaquette(Umu);
+    int Lt = Umu.Grid()->FullDimensions()[Nd-1];
+    assert(sumplaq.size() == Lt);
+    double vol = Umu.Grid()->gSites() / Lt;
+    double faces = (1.0 * (Nd - 1)* (Nd - 2)) / 2.0;
+    for(int t=0;t<Lt;t++)
+      sumplaq[t] = sumplaq[t] / vol / faces / Nc; // Nd , Nc dependent... FIXME
+    return sumplaq;
+  }
 
   //////////////////////////////////////////////////
   // average over all x,y,z the temporal loop
@@ -363,11 +413,11 @@ public:
     GaugeMat u = PeekIndex<LorentzIndex>(Umu, mu);  // some redundant copies
     GaugeMat vu = v*u;
       //FS = 0.25*Ta(u*v + Cshift(vu, mu, -1));
-      FS = (u*v + Cshift(vu, mu, -1));
+      FS = (u*v + Gimpl::CshiftLink(vu, mu, -1));
       FS = 0.125*(FS - adj(FS));
   }
 
-  static Real TopologicalCharge(GaugeLorentz &U){
+  static Real TopologicalCharge(const GaugeLorentz &U){
     // 4d topological charge
     assert(Nd==4);
     // Bx = -iF(y,z), By = -iF(z,y), Bz = -iF(x,y)
@@ -390,6 +440,203 @@ public:
   }
 
 
+  //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
+  static void CloverleafMxN(GaugeMat &FS, const GaugeMat &Umu, const GaugeMat &Unu, int mu, int nu, int M, int N){  
+#define Fmu(A) Gimpl::CovShiftForward(Umu, mu, A)
+#define Bmu(A) Gimpl::CovShiftBackward(Umu, mu, A)
+#define Fnu(A) Gimpl::CovShiftForward(Unu, nu, A)
+#define Bnu(A) Gimpl::CovShiftBackward(Unu, nu, A)
+#define FmuI Gimpl::CovShiftIdentityForward(Umu, mu)
+#define BmuI Gimpl::CovShiftIdentityBackward(Umu, mu)
+#define FnuI Gimpl::CovShiftIdentityForward(Unu, nu)
+#define BnuI Gimpl::CovShiftIdentityBackward(Unu, nu)
+
+    //Upper right loop
+    GaugeMat tmp = BmuI;
+    for(int i=1;i<M;i++)
+      tmp = Bmu(tmp);
+    for(int j=0;j<N;j++)
+      tmp = Bnu(tmp);
+    for(int i=0;i<M;i++)
+      tmp = Fmu(tmp);
+    for(int j=0;j<N;j++)
+      tmp = Fnu(tmp);
+      
+    FS = tmp;
+
+    //Upper left loop
+    tmp = BnuI;
+    for(int j=1;j<N;j++)
+      tmp = Bnu(tmp);
+    for(int i=0;i<M;i++)
+      tmp = Fmu(tmp);
+    for(int j=0;j<N;j++)
+      tmp = Fnu(tmp);
+    for(int i=0;i<M;i++)
+      tmp = Bmu(tmp);
+      
+    FS = FS + tmp;
+
+    //Lower right loop
+    tmp = FnuI;
+    for(int j=1;j<N;j++)
+      tmp = Fnu(tmp);
+    for(int i=0;i<M;i++)
+      tmp = Bmu(tmp);
+    for(int j=0;j<N;j++)
+      tmp = Bnu(tmp);
+    for(int i=0;i<M;i++)
+      tmp = Fmu(tmp);
+      
+    FS = FS + tmp;
+
+    //Lower left loop
+    tmp = FmuI;
+    for(int i=1;i<M;i++)
+      tmp = Fmu(tmp);
+    for(int j=0;j<N;j++)
+      tmp = Fnu(tmp);
+    for(int i=0;i<M;i++)
+      tmp = Bmu(tmp);
+    for(int j=0;j<N;j++)
+      tmp = Bnu(tmp);
+
+    FS = FS + tmp;
+
+#undef Fmu
+#undef Bmu
+#undef Fnu
+#undef Bnu
+#undef FmuI
+#undef BmuI
+#undef FnuI
+#undef BnuI
+  }
+
+  //Field strength from MxN Wilson loop
+  //Note F_numu = - F_munu
+  static void FieldStrengthMxN(GaugeMat &FS, const GaugeLorentz &U, int mu, int nu, int M, int N){  
+    GaugeMat Umu = PeekIndex<LorentzIndex>(U, mu);
+    GaugeMat Unu = PeekIndex<LorentzIndex>(U, nu);
+    if(M == N){
+      GaugeMat F(Umu.Grid());
+      CloverleafMxN(F, Umu, Unu, mu, nu, M, N);
+      FS = 0.125 * ( F - adj(F) );
+    }else{
+      //Average over both orientations
+      GaugeMat horizontal(Umu.Grid()), vertical(Umu.Grid());
+      CloverleafMxN(horizontal, Umu, Unu, mu, nu, M, N);
+      CloverleafMxN(vertical, Umu, Unu, mu, nu, N, M);
+      FS = 0.0625 * ( horizontal - adj(horizontal) + vertical - adj(vertical) );
+    }
+  }
+
+  //Topological charge contribution from MxN Wilson loops
+  //cf  https://arxiv.org/pdf/hep-lat/9701012.pdf  Eq 6
+  //output is the charge by timeslice: sum over timeslices to obtain the total
+  static std::vector<Real> TimesliceTopologicalChargeMxN(const GaugeLorentz &U, int M, int N){
+    assert(Nd == 4);
+    std::vector<std::vector<GaugeMat*> > F(Nd,std::vector<GaugeMat*>(Nd,nullptr));
+    //Note F_numu = - F_munu
+    //hence we only need to loop over mu,nu,rho,sigma that aren't related by permuting mu,nu  or rho,sigma
+    //Use nu > mu
+    for(int mu=0;mu<Nd-1;mu++){
+      for(int nu=mu+1; nu<Nd; nu++){
+	F[mu][nu] = new GaugeMat(U.Grid());
+	FieldStrengthMxN(*F[mu][nu], U, mu, nu, M, N);
+      }
+    }
+    Real coeff = -1./(32 * M_PI*M_PI * M*M * N*N); //overall sign to match CPS and Grid conventions, possibly related to time direction = 3 vs 0
+
+    static const int combs[3][4] = { {0,1,2,3}, {0,2,1,3}, {0,3,1,2} };
+    static const int signs[3] = { 1, -1, 1 }; //epsilon_{mu nu rho sigma}
+
+    ComplexField fsum(U.Grid());
+    fsum = Zero();
+    for(int c=0;c<3;c++){
+      int mu = combs[c][0], nu = combs[c][1], rho = combs[c][2], sigma = combs[c][3];
+      int eps = signs[c];
+      fsum = fsum + (8. * coeff * eps) * trace( (*F[mu][nu]) * (*F[rho][sigma]) ); 
+    }
+
+    for(int mu=0;mu<Nd-1;mu++)
+      for(int nu=mu+1; nu<Nd; nu++)
+	delete F[mu][nu];
+    
+    typedef typename ComplexField::scalar_object sobj;
+    std::vector<sobj> Tq;
+    sliceSum(fsum, Tq, Nd-1);
+
+    std::vector<Real> out(Tq.size());
+    for(int t=0;t<Tq.size();t++) out[t] = TensorRemove(Tq[t]).real();
+    return out;
+  }
+  static Real TopologicalChargeMxN(const GaugeLorentz &U, int M, int N){
+    std::vector<Real> Tq = TimesliceTopologicalChargeMxN(U,M,N);
+    Real out(0);
+    for(int t=0;t<Tq.size();t++) out += Tq[t];
+    return out;
+  }
+
+  //Generate the contributions to the 5Li topological charge from Wilson loops of the following sizes
+  //Use coefficients from hep-lat/9701012
+  //1x1 : c1=(19.-55.*c5)/9.
+  //2x2 : c2=(1-64.*c5)/9.
+  //1x2 : c3=(-64.+640.*c5)/45.
+  //1x3 : c4=1./5.-2.*c5
+  //3x3 : c5=1./20.
+  //Output array outer index contains the loops in the above order
+  //Inner index is the time coordinate
+  static std::vector<std::vector<Real> > TimesliceTopologicalCharge5LiContributions(const GaugeLorentz &U){
+    static const int exts[5][2] = { {1,1}, {2,2}, {1,2}, {1,3}, {3,3} };       
+    std::vector<std::vector<Real> > out(5);
+    for(int i=0;i<5;i++){	
+      out[i] = TimesliceTopologicalChargeMxN(U,exts[i][0],exts[i][1]);
+    }
+    return out;
+  }   
+
+  static std::vector<Real> TopologicalCharge5LiContributions(const GaugeLorentz &U){   
+    static const int exts[5][2] = { {1,1}, {2,2}, {1,2}, {1,3}, {3,3} };
+    std::vector<Real> out(5);
+    std::cout << GridLogMessage << "Computing topological charge" << std::endl;
+    for(int i=0;i<5;i++){
+      out[i] = TopologicalChargeMxN(U,exts[i][0],exts[i][1]);
+      std::cout << GridLogMessage << exts[i][0] << "x" << exts[i][1] << " Wilson loop contribution " << out[i] << std::endl;
+    }
+    return out;
+  }
+
+  //Compute the 5Li topological charge
+  static std::vector<Real> TimesliceTopologicalCharge5Li(const GaugeLorentz &U){
+    std::vector<std::vector<Real> > loops = TimesliceTopologicalCharge5LiContributions(U);
+
+    double c5=1./20.;
+    double c4=1./5.-2.*c5;
+    double c3=(-64.+640.*c5)/45.;
+    double c2=(1-64.*c5)/9.;
+    double c1=(19.-55.*c5)/9.;
+
+    int Lt = loops[0].size();
+    std::vector<Real> out(Lt,0.);
+    for(int t=0;t<Lt;t++)
+      out[t] += c1*loops[0][t] + c2*loops[1][t] + c3*loops[2][t] + c4*loops[3][t] + c5*loops[4][t];
+    return out;
+  }
+
+  static Real TopologicalCharge5Li(const GaugeLorentz &U){
+    std::vector<Real> Qt = TimesliceTopologicalCharge5Li(U);
+    Real Q = 0.;
+    for(int t=0;t<Qt.size();t++) Q += Qt[t];
+    std::cout << GridLogMessage << "5Li Topological charge: " << Q << std::endl;
+    return Q;
+  }
+
+
+
+
   //////////////////////////////////////////////////////
   // Similar to above for rectangle is required
   //////////////////////////////////////////////////////
diff --git a/examples/Example_Laplacian.cc b/examples/Example_Laplacian.cc
index fa8466cf..402f7235 100644
--- a/examples/Example_Laplacian.cc
+++ b/examples/Example_Laplacian.cc
@@ -324,7 +324,7 @@ int main(int argc, char ** argv)
 
   U_GT = U;
   // Make a random xform to teh gauge field
-  SU<Nc>::RandomGaugeTransform(RNG,U_GT,g); // Unit gauge
+  SU<Nc>::RandomGaugeTransform<PeriodicGimplR>(RNG,U_GT,g); // Unit gauge
 
   Field in_GT(&Grid); 
   Field out_GT(&Grid);
diff --git a/tests/core/Test_fft_gfix.cc b/tests/core/Test_fft_gfix.cc
index 87dbc242..6d617e25 100644
--- a/tests/core/Test_fft_gfix.cc
+++ b/tests/core/Test_fft_gfix.cc
@@ -29,14 +29,10 @@ Author: Peter Boyle <paboyle@ph.ed.ac.uk>
 #include <Grid/Grid.h>
 
 using namespace Grid;
- ;
 
-int main (int argc, char ** argv)
-{
+template<typename Gimpl>
+void run(double alpha, bool do_fft_gfix){
   std::vector<int> seeds({1,2,3,4});
-
-  Grid_init(&argc,&argv);
-
   int threads = GridThread::GetThreads();
 
   Coordinate latt_size   = GridDefaultLatt();
@@ -55,10 +51,7 @@ int main (int argc, char ** argv)
   FFT theFFT(&GRID);
 
   std::cout<<GridLogMessage << "Grid is setup to use "<<threads<<" threads"<<std::endl;
-
-  std::cout<< "*****************************************************************" <<std::endl;
-  std::cout<< "* Testing we can gauge fix steep descent a RGT of Unit gauge    *" <<std::endl;
-  std::cout<< "*****************************************************************" <<std::endl;
+  std::cout<<GridLogMessage << "Using alpha=" << alpha << std::endl;
 
   //  int coulomb_dir = -1;
   int coulomb_dir = Nd-1;
@@ -72,81 +65,167 @@ int main (int argc, char ** argv)
   LatticeColourMatrix   xform1(&GRID); // Gauge xform
   LatticeColourMatrix   xform2(&GRID); // Gauge xform
   LatticeColourMatrix   xform3(&GRID); // Gauge xform
+
+  //#########################################################################################
+
+  std::cout<< "*********************************************************************************************************" <<std::endl;
+  std::cout<< "* Testing steepest descent fixing to Landau gauge with randomly transformed unit gauge configuration    *" <<std::endl;
+  std::cout<< "*********************************************************************************************************" <<std::endl;
   
   SU<Nc>::ColdConfiguration(pRNG,Umu); // Unit gauge
   Uorg=Umu;
+
+  Real init_plaq=WilsonLoops<Gimpl>::avgPlaquette(Umu);
+  std::cout << " Initial plaquette "<< init_plaq << std::endl;
+
+  //Apply a random gauge transformation to the unit gauge config
   Urnd=Umu;
+  SU<Nc>::RandomGaugeTransform<Gimpl>(pRNG,Urnd,g);
 
-  SU<Nc>::RandomGaugeTransform(pRNG,Urnd,g); // Unit gauge
-
-  Real plaq=WilsonLoops<PeriodicGimplR>::avgPlaquette(Umu);
-  std::cout << " Initial plaquette "<<plaq << std::endl;
-
-  Real alpha=0.1;
-
+  //Gauge fix the randomly transformed field 
   Umu = Urnd;
-  FourierAcceleratedGaugeFixer<PeriodicGimplR>::SteepestDescentGaugeFix(Umu,xform1,alpha,10000,1.0e-12, 1.0e-12,false);
+  FourierAcceleratedGaugeFixer<Gimpl>::SteepestDescentGaugeFix(Umu,xform1,alpha,10000,1.0e-12, 1.0e-12,false);
 
   // Check the gauge xform matrices
   Utmp=Urnd;
-  SU<Nc>::GaugeTransform(Utmp,xform1);
+  SU<Nc>::GaugeTransform<Gimpl>(Utmp,xform1);
   Utmp = Utmp - Umu;
-  std::cout << " Norm Difference of xformed gauge "<< norm2(Utmp) << std::endl;
+  std::cout << " Check the output gauge transformation matrices applied to the original field produce the xformed field "<< norm2(Utmp) << " (expect 0)" << std::endl;
   
 
-  plaq=WilsonLoops<PeriodicGimplR>::avgPlaquette(Umu);
-  std::cout << " Final plaquette "<<plaq << std::endl;
+  Real plaq=WilsonLoops<Gimpl>::avgPlaquette(Umu);
+  std::cout << " Final plaquette "<<plaq << " diff " << plaq - init_plaq << " (expect 0)" << std::endl;
 
   Uorg = Uorg - Umu;
-  std::cout << " Norm Difference "<< norm2(Uorg) << std::endl;
-  std::cout << " Norm "<< norm2(Umu) << std::endl;
+  std::cout << " Norm difference between a unit gauge configuration and the gauge fixed configuration "<< norm2(Uorg) << " (expect 0)" << std::endl;
+  std::cout << " Norm of gauge fixed configuration "<< norm2(Umu) << std::endl;
+
+  //#########################################################################################
+  if(do_fft_gfix){
+    std::cout<< "*************************************************************************************" <<std::endl;
+    std::cout<< "* Testing Fourier accelerated fixing to Landau gauge with unit gauge configuration  *" <<std::endl;
+    std::cout<< "*************************************************************************************" <<std::endl;
+    Umu=Urnd;
+    FourierAcceleratedGaugeFixer<Gimpl>::SteepestDescentGaugeFix(Umu,xform2,alpha,10000,1.0e-12, 1.0e-12,true);
+
+    Utmp=Urnd;
+    SU<Nc>::GaugeTransform<Gimpl>(Utmp,xform2);
+    Utmp = Utmp - Umu;
+    std::cout << " Check the output gauge transformation matrices applied to the original field produce the xformed field "<< norm2(Utmp) << " (expect 0)" << std::endl;
 
 
-  std::cout<< "*****************************************************************" <<std::endl;
-  std::cout<< "* Testing Fourier accelerated fixing                            *" <<std::endl;
-  std::cout<< "*****************************************************************" <<std::endl;
-  Umu=Urnd;
-  FourierAcceleratedGaugeFixer<PeriodicGimplR>::SteepestDescentGaugeFix(Umu,xform2,alpha,10000,1.0e-12, 1.0e-12,true);
+    plaq=WilsonLoops<Gimpl>::avgPlaquette(Umu);
+    std::cout << " Final plaquette "<<plaq << " diff " << plaq - init_plaq << " (expect 0)" << std::endl;
+  }
+  //#########################################################################################
 
-  Utmp=Urnd;
-  SU<Nc>::GaugeTransform(Utmp,xform2);
-  Utmp = Utmp - Umu;
-  std::cout << " Norm Difference of xformed gauge "<< norm2(Utmp) << std::endl;
+  std::cout<< "******************************************************************************************" <<std::endl;
+  std::cout<< "* Testing steepest descent fixing to Landau gauge with random configuration             **" <<std::endl;
+  std::cout<< "******************************************************************************************" <<std::endl;
 
+  SU<Nc>::HotConfiguration(pRNG,Umu);
 
-  plaq=WilsonLoops<PeriodicGimplR>::avgPlaquette(Umu);
-  std::cout << " Final plaquette "<<plaq << std::endl;
+  init_plaq=WilsonLoops<Gimpl>::avgPlaquette(Umu);
+  std::cout << " Initial plaquette "<< init_plaq << std::endl;
 
-  std::cout<< "*****************************************************************" <<std::endl;
-  std::cout<< "* Testing non-unit configuration                                *" <<std::endl;
-  std::cout<< "*****************************************************************" <<std::endl;
+  FourierAcceleratedGaugeFixer<Gimpl>::SteepestDescentGaugeFix(Umu,alpha,10000,1.0e-12, 1.0e-12,false);
 
-  SU<Nc>::HotConfiguration(pRNG,Umu); // Unit gauge
+  plaq=WilsonLoops<Gimpl>::avgPlaquette(Umu);
+  std::cout << " Final plaquette "<<plaq << " diff " << plaq - init_plaq << " (expect 0)" << std::endl;
 
-  plaq=WilsonLoops<PeriodicGimplR>::avgPlaquette(Umu);
-  std::cout << " Initial plaquette "<<plaq << std::endl;
+  //#########################################################################################
+  if(do_fft_gfix){
+    std::cout<< "******************************************************************************************" <<std::endl;
+    std::cout<< "* Testing Fourier accelerated fixing to Landau gauge with random configuration          **" <<std::endl;
+    std::cout<< "******************************************************************************************" <<std::endl;
 
-  FourierAcceleratedGaugeFixer<PeriodicGimplR>::SteepestDescentGaugeFix(Umu,alpha,10000,1.0e-12, 1.0e-12,true);
+    SU<Nc>::HotConfiguration(pRNG,Umu);
 
-  plaq=WilsonLoops<PeriodicGimplR>::avgPlaquette(Umu);
-  std::cout << " Final plaquette "<<plaq << std::endl;
+    init_plaq=WilsonLoops<Gimpl>::avgPlaquette(Umu);
+    std::cout << " Initial plaquette "<< init_plaq << std::endl;
 
-  std::cout<< "*****************************************************************" <<std::endl;
-  std::cout<< "* Testing Fourier accelerated fixing to coulomb gauge           *" <<std::endl;
-  std::cout<< "*****************************************************************" <<std::endl;
+    FourierAcceleratedGaugeFixer<Gimpl>::SteepestDescentGaugeFix(Umu,alpha,10000,1.0e-12, 1.0e-12,true);
+
+    plaq=WilsonLoops<Gimpl>::avgPlaquette(Umu);
+    std::cout << " Final plaquette "<<plaq << " diff " << plaq - init_plaq << " (expect 0)" << std::endl;
+  }
+  //#########################################################################################
+  
+  std::cout<< "*******************************************************************************************" <<std::endl;
+  std::cout<< "* Testing steepest descent fixing to coulomb gauge with random configuration           *" <<std::endl;
+  std::cout<< "*******************************************************************************************" <<std::endl;
 
   Umu=Urnd;
-  SU<Nc>::HotConfiguration(pRNG,Umu); // Unit gauge
+  SU<Nc>::HotConfiguration(pRNG,Umu);
 
-  plaq=WilsonLoops<PeriodicGimplR>::avgPlaquette(Umu);
-  std::cout << " Initial plaquette "<<plaq << std::endl;
+  init_plaq=WilsonLoops<Gimpl>::avgPlaquette(Umu);
+  std::cout << " Initial plaquette "<< init_plaq << std::endl;
 
-  FourierAcceleratedGaugeFixer<PeriodicGimplR>::SteepestDescentGaugeFix(Umu,xform3,alpha,10000,1.0e-12, 1.0e-12,true,coulomb_dir);
+  FourierAcceleratedGaugeFixer<Gimpl>::SteepestDescentGaugeFix(Umu,xform3,alpha,10000,1.0e-12, 1.0e-12,false,coulomb_dir);
 
-  std::cout << Umu<<std::endl;
+  plaq=WilsonLoops<Gimpl>::avgPlaquette(Umu);
+  std::cout << " Final plaquette "<<plaq << " diff " << plaq - init_plaq << " (expect 0)" << std::endl;
 
-  plaq=WilsonLoops<PeriodicGimplR>::avgPlaquette(Umu);
-  std::cout << " Final plaquette "<<plaq << std::endl;
 
+  //#########################################################################################
+  if(do_fft_gfix){
+    std::cout<< "*******************************************************************************************" <<std::endl;
+    std::cout<< "* Testing Fourier accelerated fixing to coulomb gauge with random configuration           *" <<std::endl;
+    std::cout<< "*******************************************************************************************" <<std::endl;
+
+    Umu=Urnd;
+    SU<Nc>::HotConfiguration(pRNG,Umu);
+
+    init_plaq=WilsonLoops<Gimpl>::avgPlaquette(Umu);
+    std::cout << " Initial plaquette "<< init_plaq << std::endl;
+
+    FourierAcceleratedGaugeFixer<Gimpl>::SteepestDescentGaugeFix(Umu,xform3,alpha,10000,1.0e-12, 1.0e-12,true,coulomb_dir);
+
+    plaq=WilsonLoops<Gimpl>::avgPlaquette(Umu);
+    std::cout << " Final plaquette "<<plaq << " diff " << plaq - init_plaq << " (expect 0)" << std::endl;
+  }
+}
+
+int main (int argc, char ** argv)
+{
+  Grid_init(&argc,&argv);
+
+  double alpha=0.1; //step size
+  std::string gimpl = "periodic";
+  bool do_fft_gfix = true; //test fourier transformed gfix as well as steepest descent
+  for(int i=1;i<argc;i++){
+    std::string sarg(argv[i]);
+    if(sarg == "--gimpl"){
+      assert(i<argc-1 && "--gimpl option requires an argument");
+      gimpl = argv[i+1];
+      if(gimpl != "periodic" && gimpl != "conjugate")
+	assert(0 && "Invalid gimpl");
+      if(gimpl == "conjugate")
+	alpha = 0.025; //default alpha too large for CCBC
+    }else if(sarg == "--no-fft-gfix"){
+      std::cout << "Not doing the Fourier accelerated gauge fixing tests" << std::endl;
+      do_fft_gfix = false;
+    }else if(sarg == "--alpha"){
+      assert(i<argc-1 && "--alpha option requires an argument");
+      std::istringstream ss(argv[i+1]); ss >> alpha;
+    }
+  }
+
+
+  if(gimpl == "periodic"){
+    std::cout << GridLogMessage << "Using periodic boundary condition" << std::endl;
+    run<PeriodicGimplR>(alpha, do_fft_gfix);
+  }else{
+    std::vector<int> conjdirs = {1,1,0,0}; //test with 2 conjugate dirs and 2 not
+    std::cout << GridLogMessage << "Using complex conjugate boundary conditions in dimensions ";
+    for(int i=0;i<Nd;i++)
+      if(conjdirs[i])
+	std::cout << i << " ";   
+    std::cout << std::endl;
+
+    ConjugateGimplR::setDirections(conjdirs);
+    run<ConjugateGimplR>(alpha, do_fft_gfix);
+  }
+  
   Grid_finalize();
 }
diff --git a/tests/core/Test_gamma.cc b/tests/core/Test_gamma.cc
index e52049fe..05f8c505 100644
--- a/tests/core/Test_gamma.cc
+++ b/tests/core/Test_gamma.cc
@@ -228,6 +228,59 @@ void checkGammaL(const Gamma::Algebra a, GridSerialRNG &rng)
   std::cout << std::endl;
 }
 
+void checkChargeConjMatrix(){
+  //Check the properties of the charge conjugation matrix
+  //In the Grid basis C = -\gamma^2 \gamma^4
+  SpinMatrix C = testAlgebra[Gamma::Algebra::MinusGammaY] * testAlgebra[Gamma::Algebra::GammaT];
+  SpinMatrix mC = -C;
+  SpinMatrix one = testAlgebra[Gamma::Algebra::Identity];
+
+  std::cout << "Testing properties of charge conjugation matrix C = -\\gamma^2 \\gamma^4 (in Grid's basis)" << std::endl;
+
+  //C^T = -C
+  SpinMatrix Ct = transpose(C);
+  std::cout << GridLogMessage << "C^T=-C ";
+  test(Ct, mC);
+  std::cout << std::endl;
+
+  //C^\dagger = -C
+  SpinMatrix Cdag = adj(C);
+  std::cout << GridLogMessage << "C^dag=-C ";
+  test(Cdag, mC);
+  std::cout << std::endl;
+
+  //C^* = C
+  SpinMatrix Cstar = conjugate(C);
+  std::cout << GridLogMessage << "C^*=C ";
+  test(Cstar, C);
+  std::cout << std::endl;
+
+  //C^{-1} = -C
+  SpinMatrix CinvC = mC * C;
+  std::cout << GridLogMessage << "C^{-1}=-C ";
+  test(CinvC, one);
+  std::cout << std::endl;
+
+  // C^{-1} \gamma^\mu C = -[\gamma^\mu]^T
+  Gamma::Algebra gmu_a[4] = { Gamma::Algebra::GammaX, Gamma::Algebra::GammaY, Gamma::Algebra::GammaZ, Gamma::Algebra::GammaT };
+  for(int mu=0;mu<4;mu++){
+    SpinMatrix gmu = testAlgebra[gmu_a[mu]];
+    SpinMatrix Cinv_gmu_C = mC * gmu * C;
+    SpinMatrix mgmu_T = -transpose(gmu);
+    std::cout << GridLogMessage << "C^{-1} \\gamma^" << mu << " C = -[\\gamma^" << mu << "]^T ";
+    test(Cinv_gmu_C, mgmu_T);
+    std::cout << std::endl;
+  }
+  
+  //[C, \gamma^5] = 0
+  SpinMatrix Cg5 = C * testAlgebra[Gamma::Algebra::Gamma5];
+  SpinMatrix g5C = testAlgebra[Gamma::Algebra::Gamma5] * C;
+  std::cout << GridLogMessage << "C \\gamma^5 = \\gamma^5 C";
+  test(Cg5, g5C);
+  std::cout << std::endl;
+}
+
+
 int main(int argc, char *argv[])
 {
   Grid_init(&argc,&argv);
@@ -270,6 +323,13 @@ int main(int argc, char *argv[])
   {
     checkGammaL(i, sRNG);
   }
+
+  std::cout << GridLogMessage << "======== Charge conjugation matrix check" << std::endl;
+  checkChargeConjMatrix();
+  std::cout << GridLogMessage << std::endl;
+  
+
+
   
   Grid_finalize();
   
diff --git a/tests/core/Test_gparity.cc b/tests/core/Test_gparity.cc
index b2068901..5bf98ba6 100644
--- a/tests/core/Test_gparity.cc
+++ b/tests/core/Test_gparity.cc
@@ -55,13 +55,17 @@ static_assert(same_vComplex == 1, "Dirac Operators must have same underlying SIM
 int main (int argc, char ** argv)
 {
   int nu = 0;
-
+  int tbc_aprd = 0; //use antiperiodic BCs in the time direction?
+  
   Grid_init(&argc,&argv);
 
   for(int i=1;i<argc;i++){
     if(std::string(argv[i]) == "--Gparity-dir"){
       std::stringstream ss; ss << argv[i+1]; ss >> nu;
       std::cout << GridLogMessage << "Set Gparity direction to " << nu << std::endl;
+    }else if(std::string(argv[i]) == "--Tbc-APRD"){
+      tbc_aprd = 1;
+      std::cout << GridLogMessage << "Using antiperiodic BCs in the time direction" << std::endl;
     }
   }
 
@@ -155,13 +159,18 @@ int main (int argc, char ** argv)
 
   //Coordinate grid for reference
   LatticeInteger    xcoor_1f5(FGrid_1f);
-  LatticeCoordinate(xcoor_1f5,1+nu);
+  LatticeCoordinate(xcoor_1f5,1+nu); //note '1+nu'! This is because for 5D fields the s-direction is direction 0
   Replicate(src,src_1f);
   src_1f   = where( xcoor_1f5 >= Integer(L), 2.0*src_1f,src_1f );
 
   RealD mass=0.0;
   RealD M5=1.8;
-  StandardDiracOp Ddwf(Umu_1f,*FGrid_1f,*FrbGrid_1f,*UGrid_1f,*UrbGrid_1f,mass,M5 DOP_PARAMS);
+
+  //Standard Dirac op
+  AcceleratorVector<Complex,4> bc_std(Nd, 1.0);
+  if(tbc_aprd) bc_std[Nd-1] = -1.; //antiperiodic time BC
+  StandardDiracOp::ImplParams std_params(bc_std);
+  StandardDiracOp Ddwf(Umu_1f,*FGrid_1f,*FrbGrid_1f,*UGrid_1f,*UrbGrid_1f,mass,M5 DOP_PARAMS, std_params);
 
   StandardFermionField    src_o_1f(FrbGrid_1f);
   StandardFermionField result_o_1f(FrbGrid_1f);
@@ -172,9 +181,11 @@ int main (int argc, char ** argv)
   ConjugateGradient<StandardFermionField> CG(1.0e-8,10000);
   CG(HermOpEO,src_o_1f,result_o_1f);
   
-  //  const int nu = 3;
+  //Gparity Dirac op
   std::vector<int> twists(Nd,0);
   twists[nu] = 1;
+  if(tbc_aprd) twists[Nd-1] = 1;
+
   GparityDiracOp::ImplParams params;
   params.twists = twists;
   GparityDiracOp GPDdwf(Umu_2f,*FGrid_2f,*FrbGrid_2f,*UGrid_2f,*UrbGrid_2f,mass,M5 DOP_PARAMS,params);
@@ -271,8 +282,11 @@ int main (int argc, char ** argv)
   std::cout << "2f cb "<<result_o_2f.Checkerboard()<<std::endl;
   std::cout << "1f cb "<<result_o_1f.Checkerboard()<<std::endl;
 
-  std::cout << " result norms " <<norm2(result_o_2f)<<" " <<norm2(result_o_1f)<<std::endl;
+  //Compare norms
+  std::cout << " result norms 2f: " <<norm2(result_o_2f)<<" 1f: " <<norm2(result_o_1f)<<std::endl;
 
+
+  //Take the 2f solution and convert into the corresponding 1f solution (odd cb only)
   StandardFermionField    res0o  (FrbGrid_2f); 
   StandardFermionField    res1o  (FrbGrid_2f); 
   StandardFermionField    res0  (FGrid_2f); 
@@ -281,14 +295,15 @@ int main (int argc, char ** argv)
   res0=Zero();
   res1=Zero();
 
-  res0o = PeekIndex<0>(result_o_2f,0);
-  res1o = PeekIndex<0>(result_o_2f,1);
+  res0o = PeekIndex<0>(result_o_2f,0); //flavor 0, odd cb
+  res1o = PeekIndex<0>(result_o_2f,1); //flavor 1, odd cb
 
   std::cout << "res cb "<<res0o.Checkerboard()<<std::endl;
   std::cout << "res cb "<<res1o.Checkerboard()<<std::endl;
 
-  setCheckerboard(res0,res0o);
-  setCheckerboard(res1,res1o);
+  //poke odd onto non-cb field
+  setCheckerboard(res0,res0o); 
+  setCheckerboard(res1,res1o); 
 
   StandardFermionField replica (FGrid_1f);
   StandardFermionField replica0(FGrid_1f);
@@ -296,12 +311,13 @@ int main (int argc, char ** argv)
   Replicate(res0,replica0);
   Replicate(res1,replica1);
 
+  //2nd half of doubled lattice has f=1
   replica = where( xcoor_1f5 >= Integer(L), replica1,replica0 );
 
   replica0 = Zero();
   setCheckerboard(replica0,result_o_1f);
 
-  std::cout << "Norm2 solutions is " <<norm2(replica)<<" "<< norm2(replica0)<<std::endl;
+  std::cout << "Norm2 solutions 1f reconstructed from 2f: " <<norm2(replica)<<" Actual 1f: "<< norm2(replica0)<<std::endl;
 
   replica = replica - replica0;
   
diff --git a/tests/core/Test_gparity_flavour.cc b/tests/core/Test_gparity_flavour.cc
new file mode 100644
index 00000000..5b204e04
--- /dev/null
+++ b/tests/core/Test_gparity_flavour.cc
@@ -0,0 +1,177 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: ./tests/Test_gparity_flavour.cc
+
+Copyright (C) 2015-2017
+
+Author: Christopher Kelly <ckelly@bnl.gov>
+Author: Peter Boyle <paboyle@ph.ed.ac.uk>
+
+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>
+
+using namespace Grid;
+
+static constexpr double                      tolerance = 1.0e-6;
+static std::array<GparityFlavourMatrix, GparityFlavour::nSigma> testAlgebra;
+
+void print(const GparityFlavourMatrix &g)
+{
+  for(int i = 0; i < Ngp; i++)
+  {
+    std::cout << GridLogMessage << "(";
+    for(int j=0;j<Ngp;j++){
+      if ( abs( g(i,j)()() ) == 0 ) {
+        std::cout<< " 0";
+      } else if ( abs(g(i,j)()() - Complex(0,1)) == 0){
+        std::cout<< " i";
+      } else if ( abs(g(i,j)()() + Complex(0,1)) == 0){
+        std::cout<< "-i";
+      } else if ( abs(g(i,j)()() - Complex(1,0)) == 0){
+        std::cout<< " 1";
+      } else if ( abs(g(i,j)()() + Complex(1,0)) == 0){
+        std::cout<< "-1";
+      }
+      std::cout<<((j == Ngp-1) ? ")" : "," );
+    }
+    std::cout << std::endl;
+  }
+  std::cout << GridLogMessage << std::endl;
+}
+
+void createTestAlgebra(void)
+{
+  std::array<GparityFlavourMatrix, 3> testg;
+  const Complex             I(0., 1.), mI(0., -1.);
+
+  // 0 1
+  // 1 0
+  testg[0] = Zero();
+  testg[0](0, 1)()() = 1.;
+  testg[0](1, 0)()() = 1.;
+  std::cout << GridLogMessage << "test SigmaX= " << std::endl;
+  print(testg[0]);
+
+  // 0 -i
+  // i  0
+  testg[1] = Zero();
+  testg[1](0, 1)()() = mI;
+  testg[1](1, 0)()() = I;
+  std::cout << GridLogMessage << "test SigmaY= " << std::endl;
+  print(testg[1]);
+
+  // 1  0
+  // 0 -1
+  testg[2] = Zero();
+  testg[2](0, 0)()() = 1.0;
+  testg[2](1, 1)()() = -1.0;
+  std::cout << GridLogMessage << "test SigmaZ= " << std::endl;
+  print(testg[2]);
+
+  
+#define DEFINE_TEST_G(g, exp)\
+testAlgebra[GparityFlavour::Algebra::g]        = exp; \
+testAlgebra[GparityFlavour::Algebra::Minus##g] = -exp;
+  
+  DEFINE_TEST_G(SigmaX      , testg[0]);
+  DEFINE_TEST_G(SigmaY      , testg[1]);
+  DEFINE_TEST_G(SigmaZ      , testg[2]);
+  DEFINE_TEST_G(Identity    , 1.);
+
+  GparityFlavourMatrix pplus;
+  pplus = 1.0;
+  pplus = pplus + testg[1];
+  pplus = pplus * 0.5;
+
+  DEFINE_TEST_G(ProjPlus    , pplus);
+  
+  GparityFlavourMatrix pminus;
+  pminus = 1.0;
+  pminus = pminus - testg[1];
+  pminus = pminus * 0.5;
+
+  DEFINE_TEST_G(ProjMinus    , pminus);
+
+#undef DEFINE_TEST_G
+}
+
+template <typename Expr>
+void test(const Expr &a, const Expr &b)
+{
+  if (norm2(a - b) < tolerance)
+  {
+    std::cout << "[OK] ";
+  }
+  else
+  {
+    std::cout << "[fail]" << std::endl;
+    std::cout << GridLogError << "a= " << a << std::endl;
+    std::cout << GridLogError << "is different (tolerance= " << tolerance << ") from " << std::endl;
+    std::cout << GridLogError << "b= " << b << std::endl;
+    exit(EXIT_FAILURE);
+  }
+}
+
+void checkSigma(const GparityFlavour::Algebra a, GridSerialRNG &rng)
+{
+  GparityFlavourVector v;
+  GparityFlavourMatrix m, &testg = testAlgebra[a];
+  GparityFlavour      g(a);
+  
+  random(rng, v);
+  random(rng, m);
+  
+  std::cout << GridLogMessage << "Checking " << GparityFlavour::name[a] << ": ";
+  std::cout << "vecmul ";
+  test(g*v, testg*v);
+  std::cout << "matlmul ";
+  test(g*m, testg*m);
+  std::cout << "matrmul ";
+  test(m*g, m*testg);
+  std::cout << std::endl;
+}
+
+int main(int argc, char *argv[])
+{
+  Grid_init(&argc,&argv);
+  
+  Coordinate latt_size   = GridDefaultLatt();
+  Coordinate simd_layout = GridDefaultSimd(4,vComplex::Nsimd());
+  Coordinate mpi_layout  = GridDefaultMpi();
+  
+  GridCartesian Grid(latt_size,simd_layout,mpi_layout);
+  GridSerialRNG sRNG;
+  
+  sRNG.SeedFixedIntegers(std::vector<int>({45,12,81,9}));
+  
+  std::cout << GridLogMessage << "======== Test algebra" << std::endl;
+  createTestAlgebra();
+  std::cout << GridLogMessage << "======== Multiplication operators check" << std::endl;
+  for (int i = 0; i < GparityFlavour::nSigma; ++i)
+  {
+    checkSigma(i, sRNG);
+  }
+  std::cout << GridLogMessage << std::endl;
+  
+  Grid_finalize();
+  
+  return EXIT_SUCCESS;
+}
diff --git a/tests/forces/Test_dwf_gpforce.cc b/tests/forces/Test_dwf_gpforce.cc
index 1fa1c6e4..9db2c563 100644
--- a/tests/forces/Test_dwf_gpforce.cc
+++ b/tests/forces/Test_dwf_gpforce.cc
@@ -71,26 +71,14 @@ int main (int argc, char ** argv)
   ////////////////////////////////////
   RealD mass=0.2; //kills the diagonal term
   RealD M5=1.8;
-  //  const int nu = 3;
-  //  std::vector<int> twists(Nd,0); // twists[nu] = 1;
-  //  GparityDomainWallFermionR::ImplParams params;  params.twists = twists;
-  //  GparityDomainWallFermionR Ddwf(U,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,params);
 
-  //  DomainWallFermionR Dw     (U,     Grid,RBGrid,mass,M5);
-
-  const int nu = 3;
+  const int nu = 0; //gparity direction
   std::vector<int> twists(Nd,0);
   twists[nu] = 1;
+  twists[Nd-1] = 1; //antiperiodic in time
   GparityDomainWallFermionR::ImplParams params;
   params.twists = twists;
-
-  /*
-  params.boundary_phases[0] = 1.0;
-  params.boundary_phases[1] = 1.0;
-  params.boundary_phases[2] = 1.0;
-  params.boundary_phases[3] =- 1.0;
-  */
-  
+ 
   GparityDomainWallFermionR Dw(U,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,params);
 
   Dw.M   (phi,Mphi);
diff --git a/tests/forces/Test_gpdwf_force.cc b/tests/forces/Test_gpdwf_force.cc
index d6744080..af1ce82b 100644
--- a/tests/forces/Test_gpdwf_force.cc
+++ b/tests/forces/Test_gpdwf_force.cc
@@ -71,8 +71,10 @@ int main (int argc, char ** argv)
   RealD mass=0.01; 
   RealD M5=1.8; 
 
-  const int nu = 3;
-  std::vector<int> twists(Nd,0);  twists[nu] = 1;
+  const int nu = 1;
+  std::vector<int> twists(Nd,0);
+  twists[nu] = 1;
+  twists[3] = 1;
   GparityDomainWallFermionR::ImplParams params;  params.twists = twists;
   GparityDomainWallFermionR Ddwf(U,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,params);
   Ddwf.M   (phi,Mphi);
diff --git a/tests/forces/Test_gpwilson_force.cc b/tests/forces/Test_gpwilson_force.cc
index d731f27a..7ab2ddeb 100644
--- a/tests/forces/Test_gpwilson_force.cc
+++ b/tests/forces/Test_gpwilson_force.cc
@@ -64,8 +64,12 @@ int main (int argc, char ** argv)
   ////////////////////////////////////
   RealD mass=0.01; 
 
-  const int nu = 3;
-  std::vector<int> twists(Nd,0);  twists[nu] = 1;
+  const int nu = 1;
+  const int Lnu=latt_size[nu];
+
+  std::vector<int> twists(Nd,0);
+  twists[nu] = 1;
+  twists[3]=1;
   GparityWilsonFermionR::ImplParams params;  params.twists = twists;
   GparityWilsonFermionR Wil(U,*UGrid,*UrbGrid,mass,params);
   Wil.M   (phi,Mphi);
diff --git a/tests/lanczos/Test_dwf_lanczos.cc b/tests/lanczos/Test_dwf_lanczos.cc
index 00d29ec0..1fe29bb2 100644
--- a/tests/lanczos/Test_dwf_lanczos.cc
+++ b/tests/lanczos/Test_dwf_lanczos.cc
@@ -31,14 +31,38 @@ using namespace std;
 using namespace Grid;
  ;
 
-typedef typename GparityDomainWallFermionR::FermionField FermionField;
+template<typename Action>
+struct Setup{};
 
-RealD AllZero(RealD x){ return 0.;}
+template<>
+struct Setup<GparityMobiusFermionR>{
+  static GparityMobiusFermionR* getAction(LatticeGaugeField &Umu,
+					  GridCartesian* FGrid, GridRedBlackCartesian* FrbGrid, GridCartesian* UGrid, GridRedBlackCartesian* UrbGrid){
+    RealD mass=0.01;
+    RealD M5=1.8;
+    RealD mob_b=1.5;
+    GparityMobiusFermionD ::ImplParams params;
+    std::vector<int> twists({1,1,1,0});
+    params.twists = twists;
+    return new GparityMobiusFermionR(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,mob_b,mob_b-1.,params);
+  }
+};
 
-int main (int argc, char ** argv)
-{
-  Grid_init(&argc,&argv);
+template<>
+struct Setup<DomainWallFermionR>{
+  static DomainWallFermionR* getAction(LatticeGaugeField &Umu,
+					  GridCartesian* FGrid, GridRedBlackCartesian* FrbGrid, GridCartesian* UGrid, GridRedBlackCartesian* UrbGrid){
+    RealD mass=0.01;
+    RealD M5=1.8;
+    return new DomainWallFermionR(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5);
+  }
+};
 
+
+
+template<typename Action>
+void run(){
+  typedef typename Action::FermionField FermionField;
   const int Ls=8;
 
   GridCartesian         * UGrid   = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd,vComplex::Nsimd()),GridDefaultMpi());
@@ -56,24 +80,10 @@ int main (int argc, char ** argv)
   LatticeGaugeField Umu(UGrid); 
   SU<Nc>::HotConfiguration(RNG4, Umu);
 
-  std::vector<LatticeColourMatrix> U(4,UGrid);
-  for(int mu=0;mu<Nd;mu++){
-    U[mu] = PeekIndex<LorentzIndex>(Umu,mu);
-  }
-  
-  RealD mass=0.01;
-  RealD M5=1.8;
-  RealD mob_b=1.5;
-//  DomainWallFermionR Ddwf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5);
-  GparityMobiusFermionD ::ImplParams params;
-  std::vector<int> twists({1,1,1,0});
-  params.twists = twists;
-  GparityMobiusFermionR  Ddwf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,mob_b,mob_b-1.,params);
-
-//  MdagMLinearOperator<DomainWallFermionR,LatticeFermion> HermOp(Ddwf);
-//  SchurDiagTwoOperator<DomainWallFermionR,LatticeFermion> HermOp(Ddwf);
-  SchurDiagTwoOperator<GparityMobiusFermionR,FermionField> HermOp(Ddwf);
-//  SchurDiagMooeeOperator<DomainWallFermionR,LatticeFermion> HermOp(Ddwf);
+  Action *action = Setup<Action>::getAction(Umu,FGrid,FrbGrid,UGrid,UrbGrid);
+ 
+  //MdagMLinearOperator<Action,FermionField> HermOp(Ddwf);
+  SchurDiagTwoOperator<Action,FermionField> HermOp(*action);
 
   const int Nstop = 30;
   const int Nk = 40;
@@ -90,8 +100,7 @@ int main (int argc, char ** argv)
      PlainHermOp<FermionField> Op     (HermOp);
 
   ImplicitlyRestartedLanczos<FermionField> IRL(OpCheby,Op,Nstop,Nk,Nm,resid,MaxIt);
-
-  
+ 
   std::vector<RealD>          eval(Nm);
   FermionField    src(FrbGrid); 
   gaussian(RNG5rb,src);
@@ -103,6 +112,28 @@ int main (int argc, char ** argv)
   int Nconv;
   IRL.calc(eval,evec,src,Nconv);
 
+  delete action;
+}
+  
+int main (int argc, char ** argv)
+{
+  Grid_init(&argc,&argv);
 
+  std::string action = "GparityMobius";
+  for(int i=1;i<argc;i++){
+    if(std::string(argv[i]) == "-action"){
+      action = argv[i+1];
+    }
+  }
+
+  if(action == "GparityMobius"){
+    run<GparityMobiusFermionR>();
+  }else if(action == "DWF"){
+    run<DomainWallFermionR>();
+  }else{
+    std::cout << "Unknown action" << std::endl;
+    exit(1);
+  }
+  
   Grid_finalize();
 }
diff --git a/tests/solver/Test_dwf_multishift_mixedprec.cc b/tests/solver/Test_dwf_multishift_mixedprec.cc
new file mode 100644
index 00000000..bdede459
--- /dev/null
+++ b/tests/solver/Test_dwf_multishift_mixedprec.cc
@@ -0,0 +1,184 @@
+    /*************************************************************************************
+
+    Grid physics library, www.github.com/paboyle/Grid 
+
+    Source file: ./tests/Test_dwf_multishift_mixedprec.cc
+
+    Copyright (C) 2015
+
+Author: Christopher Kelly <ckelly@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>
+
+using namespace Grid;
+
+template<typename SpeciesD, typename SpeciesF, typename GaugeStatisticsType>
+void run_test(int argc, char ** argv, const typename SpeciesD::ImplParams &params){
+  const int Ls = 16;
+  GridCartesian* UGrid_d = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd, vComplexD::Nsimd()), GridDefaultMpi());
+  GridRedBlackCartesian* UrbGrid_d = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid_d);
+  GridCartesian* FGrid_d = SpaceTimeGrid::makeFiveDimGrid(Ls, UGrid_d);
+  GridRedBlackCartesian* FrbGrid_d = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls, UGrid_d);
+
+  GridCartesian* UGrid_f = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd, vComplexF::Nsimd()), GridDefaultMpi());
+  GridRedBlackCartesian* UrbGrid_f = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid_f);
+  GridCartesian* FGrid_f = SpaceTimeGrid::makeFiveDimGrid(Ls, UGrid_f);
+  GridRedBlackCartesian* FrbGrid_f = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls, UGrid_f);
+
+  typedef typename SpeciesD::FermionField FermionFieldD;
+  typedef typename SpeciesF::FermionField FermionFieldF;
+  
+  std::vector<int> seeds4({1, 2, 3, 4});
+  std::vector<int> seeds5({5, 6, 7, 8});
+  GridParallelRNG RNG5(FGrid_d);
+  RNG5.SeedFixedIntegers(seeds5);
+  GridParallelRNG RNG4(UGrid_d);
+  RNG4.SeedFixedIntegers(seeds4);
+
+  FermionFieldD src_d(FGrid_d);
+  random(RNG5, src_d);
+
+  LatticeGaugeFieldD Umu_d(UGrid_d);
+
+  //CPS-created G-parity ensembles have a factor of 2 error in the plaquette that causes the read to fail unless we workaround it
+  bool gparity_plaquette_fix = false;
+  for(int i=1;i<argc;i++){
+    if(std::string(argv[i]) == "--gparity_plaquette_fix"){
+      gparity_plaquette_fix=true;
+      break;
+    }
+  }
+
+  bool cfg_loaded=false;
+  for(int i=1;i<argc;i++){
+    if(std::string(argv[i]) == "--load_config"){
+      assert(i != argc-1);
+      std::string file = argv[i+1];
+      NerscIO io;
+      FieldMetaData metadata;
+
+      if(gparity_plaquette_fix) NerscIO::exitOnReadPlaquetteMismatch() = false;
+
+      io.readConfiguration<GaugeStatisticsType>(Umu_d, metadata, file);
+
+      if(gparity_plaquette_fix){
+	metadata.plaquette *= 2.; //correct header value
+
+	//Get the true plaquette
+	FieldMetaData tmp;
+	GaugeStatisticsType gs; gs(Umu_d, tmp);
+	
+	std::cout << "After correction: plaqs " << tmp.plaquette << " " << metadata.plaquette << std::endl;
+	assert(fabs(tmp.plaquette -metadata.plaquette ) < 1.0e-5 );
+      }
+
+      cfg_loaded=true;
+      break;
+    }
+  }
+
+  if(!cfg_loaded)
+    SU<Nc>::HotConfiguration(RNG4, Umu_d);
+
+  LatticeGaugeFieldF Umu_f(UGrid_f);
+  precisionChange(Umu_f, Umu_d);
+
+  std::cout << GridLogMessage << "Lattice dimensions: " << GridDefaultLatt() << "   Ls: " << Ls << std::endl;
+
+  RealD mass = 0.01;
+  RealD M5 = 1.8;
+  SpeciesD Ddwf_d(Umu_d, *FGrid_d, *FrbGrid_d, *UGrid_d, *UrbGrid_d, mass, M5, params);
+  SpeciesF Ddwf_f(Umu_f, *FGrid_f, *FrbGrid_f, *UGrid_f, *UrbGrid_f, mass, M5, params);
+
+  FermionFieldD src_o_d(FrbGrid_d);
+  pickCheckerboard(Odd, src_o_d, src_d);
+
+  SchurDiagMooeeOperator<SpeciesD, FermionFieldD> HermOpEO_d(Ddwf_d);
+  SchurDiagMooeeOperator<SpeciesF, FermionFieldF> HermOpEO_f(Ddwf_f);
+
+  AlgRemez remez(1e-4, 64, 50);
+  int order = 15;
+  remez.generateApprox(order, 1, 2); //sqrt
+
+  MultiShiftFunction shifts(remez, 1e-10, false);
+
+  int relup_freq = 50;
+  double t1=usecond();
+  ConjugateGradientMultiShiftMixedPrec<FermionFieldD,FermionFieldF> mcg(10000, shifts, FrbGrid_f, HermOpEO_f, relup_freq);
+
+  std::vector<FermionFieldD> results_o_d(order, FrbGrid_d);
+  mcg(HermOpEO_d, src_o_d, results_o_d);
+  double t2=usecond();
+
+  //Crosscheck double and mixed prec results
+  ConjugateGradientMultiShift<FermionFieldD> dmcg(10000, shifts);
+  std::vector<FermionFieldD> results_o_d_2(order, FrbGrid_d);
+  dmcg(HermOpEO_d, src_o_d, results_o_d_2);
+  double t3=usecond();
+
+  std::cout << GridLogMessage << "Comparison of mixed prec results to double prec results |mixed - double|^2 :" << std::endl;
+  FermionFieldD tmp(FrbGrid_d);
+  for(int i=0;i<order;i++){
+    RealD ndiff = axpy_norm(tmp, -1., results_o_d[i], results_o_d_2[i]);
+    std::cout << i << " " << ndiff << std::endl;
+  }
+
+  std::cout<<GridLogMessage << "Mixed precision algorithm: Total usec    =   "<< (t2-t1)<<std::endl;
+  std::cout<<GridLogMessage << "Double precision algorithm: Total usec    =   "<< (t3-t2)<<std::endl;
+}
+
+
+
+
+
+int main (int argc, char ** argv)
+{
+  Grid_init(&argc, &argv);
+
+  bool gparity = false;
+  int gpdir;
+
+  for(int i=1;i<argc;i++){
+    std::string arg(argv[i]);
+    if(arg == "--Gparity"){
+      assert(i!=argc-1);
+      gpdir = std::stoi(argv[i+1]);
+      assert(gpdir >= 0 && gpdir <= 2); //spatial!
+      gparity = true;
+    }
+  }
+  if(gparity){
+    std::cout << "Running test with G-parity BCs in " << gpdir << " direction" << std::endl;
+    GparityWilsonImplParams params;
+    params.twists[gpdir] = 1;
+    
+    std::vector<int> conj_dirs(Nd,0);
+    conj_dirs[gpdir] = 1;
+    ConjugateGimplD::setDirections(conj_dirs);
+
+    run_test<GparityDomainWallFermionD, GparityDomainWallFermionF, ConjugateGaugeStatistics>(argc,argv,params);
+  }else{
+    std::cout << "Running test with periodic BCs" << std::endl;
+    WilsonImplParams params;
+    run_test<DomainWallFermionD, DomainWallFermionF, PeriodicGaugeStatistics>(argc,argv,params);
+  }
+
+  Grid_finalize();
+}