diff --git a/lib/Make.inc b/lib/Make.inc
index 5ab4d810..73d13eb5 100644
--- a/lib/Make.inc
+++ b/lib/Make.inc
@@ -1,4 +1,4 @@
 
-HFILES=./algorithms/approx/bigfloat.h ./algorithms/approx/bigfloat_double.h ./algorithms/approx/Chebyshev.h ./algorithms/approx/MultiShiftFunction.h ./algorithms/approx/Remez.h ./algorithms/approx/Zolotarev.h ./algorithms/CoarsenedMatrix.h ./algorithms/iterative/AdefGeneric.h ./algorithms/iterative/ConjugateGradient.h ./algorithms/iterative/ConjugateGradientMultiShift.h ./algorithms/iterative/ConjugateResidual.h ./algorithms/iterative/DenseMatrix.h ./algorithms/iterative/EigenSort.h ./algorithms/iterative/Francis.h ./algorithms/iterative/Householder.h ./algorithms/iterative/ImplicitlyRestartedLanczos.h ./algorithms/iterative/Matrix.h ./algorithms/iterative/MatrixUtils.h ./algorithms/iterative/NormalEquations.h ./algorithms/iterative/PrecConjugateResidual.h ./algorithms/iterative/PrecGeneralisedConjugateResidual.h ./algorithms/iterative/SchurRedBlack.h ./algorithms/LinearOperator.h ./algorithms/Preconditioner.h ./algorithms/SparseMatrix.h ./Algorithms.h ./AlignedAllocator.h ./cartesian/Cartesian_base.h ./cartesian/Cartesian_full.h ./cartesian/Cartesian_red_black.h ./Cartesian.h ./communicator/Communicator_base.h ./Communicator.h ./cshift/Cshift_common.h ./cshift/Cshift_mpi.h ./cshift/Cshift_none.h ./Cshift.h ./Grid.h ./Init.h ./lattice/Lattice_arith.h ./lattice/Lattice_base.h ./lattice/Lattice_comparison.h ./lattice/Lattice_comparison_utils.h ./lattice/Lattice_conformable.h ./lattice/Lattice_coordinate.h ./lattice/Lattice_ET.h ./lattice/Lattice_local.h ./lattice/Lattice_overload.h ./lattice/Lattice_peekpoke.h ./lattice/Lattice_reality.h ./lattice/Lattice_reduction.h ./lattice/Lattice_rng.h ./lattice/Lattice_trace.h ./lattice/Lattice_transfer.h ./lattice/Lattice_transpose.h ./lattice/Lattice_unary.h ./lattice/Lattice_where.h ./Lattice.h ./Log.h ./Old/Tensor_peek.h ./Old/Tensor_poke.h ./parallelIO/BinaryIO.h ./parallelIO/NerscIO.h ./PerfCount.h ./pugixml/pugixml.h ./qcd/action/ActionBase.h ./qcd/action/ActionParams.h ./qcd/action/Actions.h ./qcd/action/fermion/CayleyFermion5D.h ./qcd/action/fermion/ContinuedFractionFermion5D.h ./qcd/action/fermion/DomainWallFermion.h ./qcd/action/fermion/FermionOperator.h ./qcd/action/fermion/FermionOperatorImpl.h ./qcd/action/fermion/g5HermitianLinop.h ./qcd/action/fermion/MobiusFermion.h ./qcd/action/fermion/MobiusZolotarevFermion.h ./qcd/action/fermion/OverlapWilsonCayleyTanhFermion.h ./qcd/action/fermion/OverlapWilsonCayleyZolotarevFermion.h ./qcd/action/fermion/OverlapWilsonContfracTanhFermion.h ./qcd/action/fermion/OverlapWilsonContfracZolotarevFermion.h ./qcd/action/fermion/OverlapWilsonPartialFractionTanhFermion.h ./qcd/action/fermion/OverlapWilsonPartialFractionZolotarevFermion.h ./qcd/action/fermion/PartialFractionFermion5D.h ./qcd/action/fermion/ScaledShamirFermion.h ./qcd/action/fermion/ShamirZolotarevFermion.h ./qcd/action/fermion/WilsonCompressor.h ./qcd/action/fermion/WilsonFermion.h ./qcd/action/fermion/WilsonFermion5D.h ./qcd/action/fermion/WilsonKernels.h ./qcd/action/fermion/WilsonTMFermion.h ./qcd/action/gauge/PlaqPlusRectangleAction.h ./qcd/action/gauge/WilsonGaugeAction.h ./qcd/action/pseudofermion/EvenOddSchurDifferentiable.h ./qcd/action/pseudofermion/OneFlavourEvenOddRational.h ./qcd/action/pseudofermion/OneFlavourEvenOddRationalRatio.h ./qcd/action/pseudofermion/OneFlavourRational.h ./qcd/action/pseudofermion/OneFlavourRationalRatio.h ./qcd/action/pseudofermion/TwoFlavour.h ./qcd/action/pseudofermion/TwoFlavourEvenOdd.h ./qcd/action/pseudofermion/TwoFlavourEvenOddRatio.h ./qcd/action/pseudofermion/TwoFlavourRatio.h ./qcd/hmc/HMC.h ./qcd/hmc/HmcRunner.h ./qcd/hmc/integrators/Integrator.h ./qcd/hmc/integrators/Integrator_algorithm.h ./qcd/hmc/NerscCheckpointer.h ./qcd/QCD.h ./qcd/spin/Dirac.h ./qcd/spin/TwoSpinor.h ./qcd/utils/CovariantCshift.h ./qcd/utils/LinalgUtils.h ./qcd/utils/SpaceTimeGrid.h ./qcd/utils/SUn.h ./qcd/utils/WilsonLoops.h ./serialisation/BaseIO.h ./serialisation/BinaryIO.h ./serialisation/MacroMagic.h ./serialisation/Serialisation.h ./serialisation/TextIO.h ./serialisation/XmlIO.h ./simd/Avx512Asm.h ./simd/Grid_avx.h ./simd/Grid_avx512.h ./simd/Grid_empty.h ./simd/Grid_imci.h ./simd/Grid_neon.h ./simd/Grid_qpx.h ./simd/Grid_sse4.h ./simd/Grid_vector_types.h ./simd/Grid_vector_unops.h ./Simd.h ./stencil/Lebesgue.h ./Stencil.h ./tensors/Tensor_arith.h ./tensors/Tensor_arith_add.h ./tensors/Tensor_arith_mac.h ./tensors/Tensor_arith_mul.h ./tensors/Tensor_arith_scalar.h ./tensors/Tensor_arith_sub.h ./tensors/Tensor_class.h ./tensors/Tensor_determinant.h ./tensors/Tensor_exp.h ./tensors/Tensor_extract_merge.h ./tensors/Tensor_index.h ./tensors/Tensor_inner.h ./tensors/Tensor_logical.h ./tensors/Tensor_outer.h ./tensors/Tensor_reality.h ./tensors/Tensor_Ta.h ./tensors/Tensor_trace.h ./tensors/Tensor_traits.h ./tensors/Tensor_transpose.h ./tensors/Tensor_unary.h ./Tensors.h ./Threads.h ./Timer.h
+HFILES=./algorithms/approx/bigfloat.h ./algorithms/approx/bigfloat_double.h ./algorithms/approx/Chebyshev.h ./algorithms/approx/MultiShiftFunction.h ./algorithms/approx/Remez.h ./algorithms/approx/Zolotarev.h ./algorithms/CoarsenedMatrix.h ./algorithms/iterative/AdefGeneric.h ./algorithms/iterative/ConjugateGradient.h ./algorithms/iterative/ConjugateGradientMultiShift.h ./algorithms/iterative/ConjugateResidual.h ./algorithms/iterative/DenseMatrix.h ./algorithms/iterative/EigenSort.h ./algorithms/iterative/Francis.h ./algorithms/iterative/Householder.h ./algorithms/iterative/ImplicitlyRestartedLanczos.h ./algorithms/iterative/Matrix.h ./algorithms/iterative/MatrixUtils.h ./algorithms/iterative/NormalEquations.h ./algorithms/iterative/PrecConjugateResidual.h ./algorithms/iterative/PrecGeneralisedConjugateResidual.h ./algorithms/iterative/SchurRedBlack.h ./algorithms/LinearOperator.h ./algorithms/Preconditioner.h ./algorithms/SparseMatrix.h ./Algorithms.h ./AlignedAllocator.h ./cartesian/Cartesian_base.h ./cartesian/Cartesian_full.h ./cartesian/Cartesian_red_black.h ./Cartesian.h ./communicator/Communicator_base.h ./Communicator.h ./cshift/Cshift_common.h ./cshift/Cshift_mpi.h ./cshift/Cshift_none.h ./Cshift.h ./Grid.h ./Init.h ./lattice/Lattice_arith.h ./lattice/Lattice_base.h ./lattice/Lattice_comparison.h ./lattice/Lattice_comparison_utils.h ./lattice/Lattice_conformable.h ./lattice/Lattice_coordinate.h ./lattice/Lattice_ET.h ./lattice/Lattice_local.h ./lattice/Lattice_overload.h ./lattice/Lattice_peekpoke.h ./lattice/Lattice_reality.h ./lattice/Lattice_reduction.h ./lattice/Lattice_rng.h ./lattice/Lattice_trace.h ./lattice/Lattice_transfer.h ./lattice/Lattice_transpose.h ./lattice/Lattice_unary.h ./lattice/Lattice_where.h ./Lattice.h ./Log.h ./Old/Tensor_peek.h ./Old/Tensor_poke.h ./parallelIO/BinaryIO.h ./parallelIO/NerscIO.h ./PerfCount.h ./pugixml/pugixml.h ./qcd/action/ActionBase.h ./qcd/action/ActionParams.h ./qcd/action/Actions.h ./qcd/action/fermion/CayleyFermion5D.h ./qcd/action/fermion/ContinuedFractionFermion5D.h ./qcd/action/fermion/DomainWallFermion.h ./qcd/action/fermion/FermionOperator.h ./qcd/action/fermion/FermionOperatorImpl.h ./qcd/action/fermion/g5HermitianLinop.h ./qcd/action/fermion/MobiusFermion.h ./qcd/action/fermion/MobiusZolotarevFermion.h ./qcd/action/fermion/OverlapWilsonCayleyTanhFermion.h ./qcd/action/fermion/OverlapWilsonCayleyZolotarevFermion.h ./qcd/action/fermion/OverlapWilsonContfracTanhFermion.h ./qcd/action/fermion/OverlapWilsonContfracZolotarevFermion.h ./qcd/action/fermion/OverlapWilsonPartialFractionTanhFermion.h ./qcd/action/fermion/OverlapWilsonPartialFractionZolotarevFermion.h ./qcd/action/fermion/PartialFractionFermion5D.h ./qcd/action/fermion/ScaledShamirFermion.h ./qcd/action/fermion/ShamirZolotarevFermion.h ./qcd/action/fermion/WilsonCompressor.h ./qcd/action/fermion/WilsonFermion.h ./qcd/action/fermion/WilsonFermion5D.h ./qcd/action/fermion/WilsonKernels.h ./qcd/action/fermion/WilsonTMFermion.h ./qcd/action/gauge/GaugeImpl.h ./qcd/action/gauge/PlaqPlusRectangleAction.h ./qcd/action/gauge/WilsonGaugeAction.h ./qcd/action/pseudofermion/EvenOddSchurDifferentiable.h ./qcd/action/pseudofermion/OneFlavourEvenOddRational.h ./qcd/action/pseudofermion/OneFlavourEvenOddRationalRatio.h ./qcd/action/pseudofermion/OneFlavourRational.h ./qcd/action/pseudofermion/OneFlavourRationalRatio.h ./qcd/action/pseudofermion/TwoFlavour.h ./qcd/action/pseudofermion/TwoFlavourEvenOdd.h ./qcd/action/pseudofermion/TwoFlavourEvenOddRatio.h ./qcd/action/pseudofermion/TwoFlavourRatio.h ./qcd/hmc/HMC.h ./qcd/hmc/HmcRunner.h ./qcd/hmc/integrators/Integrator.h ./qcd/hmc/integrators/Integrator_algorithm.h ./qcd/hmc/NerscCheckpointer.h ./qcd/QCD.h ./qcd/spin/Dirac.h ./qcd/spin/TwoSpinor.h ./qcd/utils/CovariantCshift.h ./qcd/utils/LinalgUtils.h ./qcd/utils/SpaceTimeGrid.h ./qcd/utils/SUn.h ./qcd/utils/WilsonLoops.h ./serialisation/BaseIO.h ./serialisation/BinaryIO.h ./serialisation/MacroMagic.h ./serialisation/Serialisation.h ./serialisation/TextIO.h ./serialisation/XmlIO.h ./simd/Avx512Asm.h ./simd/Grid_avx.h ./simd/Grid_avx512.h ./simd/Grid_empty.h ./simd/Grid_imci.h ./simd/Grid_neon.h ./simd/Grid_qpx.h ./simd/Grid_sse4.h ./simd/Grid_vector_types.h ./simd/Grid_vector_unops.h ./Simd.h ./stencil/Lebesgue.h ./Stencil.h ./tensors/Tensor_arith.h ./tensors/Tensor_arith_add.h ./tensors/Tensor_arith_mac.h ./tensors/Tensor_arith_mul.h ./tensors/Tensor_arith_scalar.h ./tensors/Tensor_arith_sub.h ./tensors/Tensor_class.h ./tensors/Tensor_determinant.h ./tensors/Tensor_exp.h ./tensors/Tensor_extract_merge.h ./tensors/Tensor_index.h ./tensors/Tensor_inner.h ./tensors/Tensor_logical.h ./tensors/Tensor_outer.h ./tensors/Tensor_reality.h ./tensors/Tensor_Ta.h ./tensors/Tensor_trace.h ./tensors/Tensor_traits.h ./tensors/Tensor_transpose.h ./tensors/Tensor_unary.h ./Tensors.h ./Threads.h ./Timer.h
 
 CCFILES=./algorithms/approx/MultiShiftFunction.cc ./algorithms/approx/Remez.cc ./algorithms/approx/Zolotarev.cc ./Init.cc ./Log.cc ./PerfCount.cc ./pugixml/pugixml.cc ./qcd/action/fermion/CayleyFermion5D.cc ./qcd/action/fermion/ContinuedFractionFermion5D.cc ./qcd/action/fermion/PartialFractionFermion5D.cc ./qcd/action/fermion/WilsonFermion.cc ./qcd/action/fermion/WilsonFermion5D.cc ./qcd/action/fermion/WilsonKernels.cc ./qcd/action/fermion/WilsonKernelsAsm.cc ./qcd/action/fermion/WilsonKernelsHand.cc ./qcd/action/fermion/WilsonTMFermion.cc ./qcd/hmc/HMC.cc ./qcd/spin/Dirac.cc ./qcd/utils/SpaceTimeGrid.cc ./serialisation/BinaryIO.cc ./serialisation/TextIO.cc ./serialisation/XmlIO.cc ./stencil/Lebesgue.cc ./stencil/Stencil_common.cc
diff --git a/lib/parallelIO/NerscIO.h b/lib/parallelIO/NerscIO.h
index 4af7cad7..a45cc387 100644
--- a/lib/parallelIO/NerscIO.h
+++ b/lib/parallelIO/NerscIO.h
@@ -51,7 +51,6 @@ class NerscField {
     std::string floating_point;
 };
 
-
 //////////////////////////////////////////////////////////////////////
 // Bit and Physical Checksumming and QA of data
 //////////////////////////////////////////////////////////////////////
@@ -69,8 +68,9 @@ inline void NerscGrid(GridBase *grid,NerscField &header)
 template<class GaugeField>
 inline void NerscStatistics(GaugeField & data,NerscField &header)
 {
-  header.link_trace=Grid::QCD::WilsonLoops<GaugeField>::linkTrace(data);
-  header.plaquette =Grid::QCD::WilsonLoops<GaugeField>::avgPlaquette(data);
+  // How to convert data precision etc...
+  header.link_trace=Grid::QCD::WilsonLoops<PeriodicGimplR>::linkTrace(data);
+  header.plaquette =Grid::QCD::WilsonLoops<PeriodicGimplR>::avgPlaquette(data);
 }
 
 inline void NerscMachineCharacteristics(NerscField &header)
diff --git a/lib/qcd/QCD.h b/lib/qcd/QCD.h
index 60a00243..d4185ae6 100644
--- a/lib/qcd/QCD.h
+++ b/lib/qcd/QCD.h
@@ -425,7 +425,6 @@ namespace QCD {
 #include <qcd/spin/TwoSpinor.h>
 #include <qcd/utils/LinalgUtils.h>
 #include <qcd/utils/CovariantCshift.h>
-#include <qcd/utils/WilsonLoops.h>
 #include <qcd/utils/SUn.h>
 #include <qcd/action/Actions.h>
 #include <qcd/hmc/integrators/Integrator.h>
diff --git a/lib/qcd/action/Actions.h b/lib/qcd/action/Actions.h
index ab5f25a5..a1fd9e1c 100644
--- a/lib/qcd/action/Actions.h
+++ b/lib/qcd/action/Actions.h
@@ -1,11 +1,6 @@
 #ifndef GRID_QCD_ACTIONS_H
 #define GRID_QCD_ACTIONS_H
 
-
-// Some reorganisation likely required as both Chroma and IroIro
-// are separating the concept of the operator from that of action.
-//
-// The FermAction contains methods to create 
 // * Linear operators             (Hermitian and non-hermitian)  .. my LinearOperator
 // * System solvers               (Hermitian and non-hermitian)  .. my OperatorFunction
 // * MultiShift System solvers    (Hermitian and non-hermitian)  .. my OperatorFunction
@@ -19,6 +14,9 @@
 ////////////////////////////////////////////
 // Utility functions
 ////////////////////////////////////////////
+#include <qcd/action/gauge/GaugeImpl.h>
+#include <qcd/utils/WilsonLoops.h>
+
 #include <qcd/action/fermion/WilsonCompressor.h>     //used by all wilson type fermions
 #include <qcd/action/fermion/FermionOperatorImpl.h>
 #include <qcd/action/fermion/FermionOperator.h>
@@ -29,20 +27,37 @@
 ////////////////////////////////////////////
 #include <qcd/action/gauge/WilsonGaugeAction.h>
 #include <qcd/action/gauge/PlaqPlusRectangleAction.h>
+
 namespace Grid {
 namespace QCD {
-typedef WilsonGaugeAction<LatticeGaugeField>     WilsonGaugeActionR;
-typedef WilsonGaugeAction<LatticeGaugeFieldF>    WilsonGaugeActionF;
-typedef WilsonGaugeAction<LatticeGaugeFieldD>    WilsonGaugeActionD;
-typedef PlaqPlusRectangleAction<LatticeGaugeField>     PlaqPlusRectangleActionR;
-typedef PlaqPlusRectangleAction<LatticeGaugeFieldF>    PlaqPlusRectangleActionF;
-typedef PlaqPlusRectangleAction<LatticeGaugeFieldD>    PlaqPlusRectangleActionD;
-typedef IwasakiGaugeAction<LatticeGaugeField>     IwasakiGaugeActionR;
-typedef IwasakiGaugeAction<LatticeGaugeFieldF>    IwasakiGaugeActionF;
-typedef IwasakiGaugeAction<LatticeGaugeFieldD>    IwasakiGaugeActionD;
-typedef SymanzikGaugeAction<LatticeGaugeField>     SymanzikGaugeActionR;
-typedef SymanzikGaugeAction<LatticeGaugeFieldF>    SymanzikGaugeActionF;
-typedef SymanzikGaugeAction<LatticeGaugeFieldD>    SymanzikGaugeActionD;
+
+typedef WilsonGaugeAction<PeriodicGimplR>          WilsonGaugeActionR;
+typedef WilsonGaugeAction<PeriodicGimplF>          WilsonGaugeActionF;
+typedef WilsonGaugeAction<PeriodicGimplD>          WilsonGaugeActionD;
+typedef PlaqPlusRectangleAction<PeriodicGimplR>    PlaqPlusRectangleActionR;
+typedef PlaqPlusRectangleAction<PeriodicGimplF>    PlaqPlusRectangleActionF;
+typedef PlaqPlusRectangleAction<PeriodicGimplD>    PlaqPlusRectangleActionD;
+typedef IwasakiGaugeAction<PeriodicGimplR>         IwasakiGaugeActionR;
+typedef IwasakiGaugeAction<PeriodicGimplF>         IwasakiGaugeActionF;
+typedef IwasakiGaugeAction<PeriodicGimplD>         IwasakiGaugeActionD;
+typedef SymanzikGaugeAction<PeriodicGimplR>        SymanzikGaugeActionR;
+typedef SymanzikGaugeAction<PeriodicGimplF>        SymanzikGaugeActionF;
+typedef SymanzikGaugeAction<PeriodicGimplD>        SymanzikGaugeActionD;
+
+
+typedef WilsonGaugeAction<ConjugateGimplR>          ConjugateWilsonGaugeActionR;
+typedef WilsonGaugeAction<ConjugateGimplF>          ConjugateWilsonGaugeActionF;
+typedef WilsonGaugeAction<ConjugateGimplD>          ConjugateWilsonGaugeActionD;
+typedef PlaqPlusRectangleAction<ConjugateGimplR>    ConjugatePlaqPlusRectangleActionR;
+typedef PlaqPlusRectangleAction<ConjugateGimplF>    ConjugatePlaqPlusRectangleActionF;
+typedef PlaqPlusRectangleAction<ConjugateGimplD>    ConjugatePlaqPlusRectangleActionD;
+typedef IwasakiGaugeAction<ConjugateGimplR>         ConjugateIwasakiGaugeActionR;
+typedef IwasakiGaugeAction<ConjugateGimplF>         ConjugateIwasakiGaugeActionF;
+typedef IwasakiGaugeAction<ConjugateGimplD>         ConjugateIwasakiGaugeActionD;
+typedef SymanzikGaugeAction<ConjugateGimplR>        ConjugateSymanzikGaugeActionR;
+typedef SymanzikGaugeAction<ConjugateGimplF>        ConjugateSymanzikGaugeActionF;
+typedef SymanzikGaugeAction<ConjugateGimplD>        ConjugateSymanzikGaugeActionD;
+
 }}
 
 ////////////////////////////////////////////////////////////////////////////////////////////////////
@@ -177,11 +192,6 @@ typedef DomainWallFermion<GparityWilsonImplD> GparityDomainWallFermionD;
 #include <qcd/action/pseudofermion/TwoFlavourEvenOdd.h>
 #include <qcd/action/pseudofermion/TwoFlavourEvenOddRatio.h>
 
-//IroIro inserted general "Nf" param; could also be done,
-//but not clear why unless into large Nf BSM studies
-//Even there, don't want the explicit (2) on power denominator
-//if even number of flavours, so further generalised interface
-//would be required but easy.
 #include <qcd/action/pseudofermion/OneFlavourRational.h>
 #include <qcd/action/pseudofermion/OneFlavourRationalRatio.h>
 #include <qcd/action/pseudofermion/OneFlavourEvenOddRational.h>
diff --git a/lib/qcd/action/fermion/FermionOperatorImpl.h b/lib/qcd/action/fermion/FermionOperatorImpl.h
index 21187870..9f18580e 100644
--- a/lib/qcd/action/fermion/FermionOperatorImpl.h
+++ b/lib/qcd/action/fermion/FermionOperatorImpl.h
@@ -58,38 +58,6 @@ namespace Grid {
     //  }
     //////////////////////////////////////////////
 
-    
-    ////////////////////////////////////////////////////////////////////////
-    // Implementation dependent gauge types
-    ////////////////////////////////////////////////////////////////////////
-
-#define INHERIT_IMPL_TYPES(Base) \
-    INHERIT_GIMPL_TYPES(Base)\
-    INHERIT_FIMPL_TYPES(Base)
-
-#define INHERIT_GIMPL_TYPES(GImpl) \
-    typedef typename GImpl::Simd                           Simd;\
-    typedef typename GImpl::GaugeLinkField       GaugeLinkField;\
-    typedef typename GImpl::GaugeField               GaugeField;	
-    
-    // Composition with smeared link, bc's etc.. probably need multiple inheritance
-    // Variable precision "S" and variable Nc
-    template<class S,int Nrepresentation=Nc>
-    class ImplGauge { 
-    public:
-    
-      typedef S Simd;
-    
-      template<typename vtype> using iImplGaugeLink          = iScalar<iScalar<iMatrix<vtype, Nrepresentation> > >;
-      template<typename vtype> using iImplGaugeField         = iVector<iScalar<iMatrix<vtype, Nrepresentation> >, Nd  >;
-    
-      typedef iImplGaugeLink    <Simd>           SiteGaugeLink;
-      typedef iImplGaugeField   <Simd>           SiteGaugeField;
-    
-      typedef Lattice<SiteGaugeLink>                GaugeLinkField; // bit ugly naming; polarised gauge field, lorentz... all ugly
-      typedef Lattice<SiteGaugeField>                   GaugeField;
-
-    };
 
     ////////////////////////////////////////////////////////////////////////
     // Implementation dependent fermion types
@@ -104,14 +72,18 @@ namespace Grid {
     typedef typename Impl::StencilImpl              StencilImpl;	\
     typedef typename Impl::ImplParams ImplParams;
 
+#define INHERIT_IMPL_TYPES(Base) \
+    INHERIT_GIMPL_TYPES(Base)\
+    INHERIT_FIMPL_TYPES(Base)
+
     ///////
     // Single flavour four spinors with colour index
     ///////
     template<class S,int Nrepresentation=Nc>
-    class WilsonImpl :  public ImplGauge<S,Nrepresentation> { 
+    class WilsonImpl :  public PeriodicGaugeImpl<S,Nrepresentation> { 
     public:
 
-      typedef ImplGauge<S,Nrepresentation> Gimpl;
+      typedef PeriodicGaugeImpl<S,Nrepresentation> Gimpl;
 
       INHERIT_GIMPL_TYPES(Gimpl);
 
@@ -180,10 +152,10 @@ PARALLEL_FOR_LOOP
     ////////////////////////////////////////////////////////////////////////////////////////
 
     template<class S,int Nrepresentation>
-    class GparityWilsonImpl : public ImplGauge<S,Nrepresentation> { 
+    class GparityWilsonImpl : public ConjugateGaugeImpl<S,Nrepresentation> { 
     public:
 
-      typedef ImplGauge<S,Nrepresentation> Gimpl;
+      typedef ConjugateGaugeImpl<S,Nrepresentation> Gimpl;
 
       INHERIT_GIMPL_TYPES(Gimpl);
 
diff --git a/lib/qcd/action/gauge/GaugeImpl.h b/lib/qcd/action/gauge/GaugeImpl.h
new file mode 100644
index 00000000..8f9a55dc
--- /dev/null
+++ b/lib/qcd/action/gauge/GaugeImpl.h
@@ -0,0 +1,150 @@
+#ifndef  GRID_QCD_GAUGE_IMPL_H
+#define  GRID_QCD_GAUGE_IMPL_H
+
+namespace Grid {
+
+  namespace QCD {
+
+    
+    ////////////////////////////////////////////////////////////////////////
+    // Implementation dependent gauge types
+    ////////////////////////////////////////////////////////////////////////
+
+template<class Gimpl> class WilsonLoops;
+
+#define INHERIT_GIMPL_TYPES(GImpl) \
+    typedef typename GImpl::Simd                           Simd;\
+    typedef typename GImpl::GaugeLinkField       GaugeLinkField;\
+    typedef typename GImpl::GaugeField               GaugeField;	
+
+    
+    // Composition with smeared link, bc's etc.. probably need multiple inheritance
+    // Variable precision "S" and variable Nc
+    template<class S,int Nrepresentation=Nc>
+    class PeriodicGaugeImpl { 
+    public:
+    
+      typedef S Simd;
+    
+      template<typename vtype> using iImplGaugeLink          = iScalar<iScalar<iMatrix<vtype, Nrepresentation> > >;
+      template<typename vtype> using iImplGaugeField         = iVector<iScalar<iMatrix<vtype, Nrepresentation> >, Nd  >;
+    
+      typedef iImplGaugeLink    <Simd>           SiteGaugeLink;
+      typedef iImplGaugeField   <Simd>           SiteGaugeField;
+    
+      typedef Lattice<SiteGaugeLink>                GaugeLinkField; // bit ugly naming; polarised gauge field, lorentz... all ugly
+      typedef Lattice<SiteGaugeField>                   GaugeField;
+     
+    ////////////////////////////////////////////////////////////////////////////////////////////////////////////
+    // Support needed for the assembly of loops including all boundary condition effects such as conjugate bcs
+    ////////////////////////////////////////////////////////////////////////////////////////////////////////////
+    
+      template<class covariant>  static inline
+      Lattice<covariant> CovShiftForward (const GaugeLinkField &Link, int mu, const Lattice<covariant> &field) {
+	return PeriodicBC::CovShiftForward(Link,mu,field);
+      }
+
+      template<class covariant> static inline
+      Lattice<covariant> CovShiftBackward(const GaugeLinkField &Link, int mu,const Lattice<covariant> &field) {
+	return PeriodicBC::CovShiftBackward(Link,mu,field);
+      }
+      static inline
+      GaugeLinkField CovShiftIdentityBackward(const GaugeLinkField &Link, int mu) {
+	return Cshift(adj(Link),mu,-1);
+      }
+      static inline
+      GaugeLinkField CovShiftIdentityForward(const GaugeLinkField &Link, int mu) {
+	return Link;
+      }
+      static inline
+      GaugeLinkField ShiftStaple(const GaugeLinkField &Link, int mu) {
+	return Cshift(Link,mu,1);
+      }
+
+      static inline bool isPeriodicGaugeField(void) {
+	return true;
+      }
+
+
+    };
+
+    
+    // Composition with smeared link, bc's etc.. probably need multiple inheritance
+    // Variable precision "S" and variable Nc
+    template<class S,int Nrepresentation=Nc>
+    class ConjugateGaugeImpl { 
+    public:
+    
+      typedef S Simd;
+    
+      template<typename vtype> using iImplGaugeLink          = iScalar<iScalar<iMatrix<vtype, Nrepresentation> > >;
+      template<typename vtype> using iImplGaugeField         = iVector<iScalar<iMatrix<vtype, Nrepresentation> >, Nd  >;
+    
+      typedef iImplGaugeLink    <Simd>           SiteGaugeLink;
+      typedef iImplGaugeField   <Simd>           SiteGaugeField;
+    
+      typedef Lattice<SiteGaugeLink>                GaugeLinkField; // bit ugly naming; polarised gauge field, lorentz... all ugly
+      typedef Lattice<SiteGaugeField>                   GaugeField;
+
+    ////////////////////////////////////////////////////////////////////////////////////////////////////////////
+    // Support needed for the assembly of loops including all boundary condition effects such as Gparity.
+    ////////////////////////////////////////////////////////////////////////////////////////////////////////////
+    template<class covariant>  static
+    Lattice<covariant> CovShiftForward (const GaugeLinkField &Link, int mu, const Lattice<covariant> &field) {
+      return GparityBC::CovShiftForward(Link,mu,field);
+    }
+
+    template<class covariant> static
+    Lattice<covariant> CovShiftBackward(const GaugeLinkField &Link, int mu,const Lattice<covariant> &field) {
+      return GparityBC::CovShiftBackward(Link,mu,field);
+    }
+
+    static inline
+    GaugeLinkField CovShiftIdentityBackward(const GaugeLinkField &Link, int mu) {
+      GridBase *grid = Link._grid;
+      int Lmu = grid->GlobalDimensions()[mu]-1;
+      
+      Lattice<iScalar<vInteger> > coor(grid);    LatticeCoordinate(coor,mu);
+
+      GaugeLinkField tmp (grid);
+      tmp=adj(Link);
+      tmp = where(coor==Lmu,conjugate(tmp),tmp);
+      return Cshift(tmp,mu,-1);// moves towards positive mu
+    }
+    static inline
+    GaugeLinkField CovShiftIdentityForward(const GaugeLinkField &Link, int mu) {
+      return Link;
+    }
+
+    static inline
+    GaugeLinkField ShiftStaple(const GaugeLinkField &Link, int mu) {
+      GridBase *grid = Link._grid;
+      int Lmu = grid->GlobalDimensions()[mu]-1;
+      
+      Lattice<iScalar<vInteger> > coor(grid);    LatticeCoordinate(coor,mu);
+
+      GaugeLinkField tmp (grid);
+      tmp=Cshift(Link,mu,1);
+      tmp=where(coor==Lmu,conjugate(tmp),tmp);
+      return tmp;
+    }
+
+    static inline bool isPeriodicGaugeField(void) {
+      return false;
+    }
+    
+    };
+
+
+    typedef PeriodicGaugeImpl<vComplex ,Nc> PeriodicGimplR; // Real.. whichever prec
+    typedef PeriodicGaugeImpl<vComplexF,Nc> PeriodicGimplF; // Float
+    typedef PeriodicGaugeImpl<vComplexD,Nc> PeriodicGimplD; // Double
+
+    typedef ConjugateGaugeImpl<vComplex ,Nc> ConjugateGimplR; // Real.. whichever prec
+    typedef ConjugateGaugeImpl<vComplexF,Nc> ConjugateGimplF; // Float
+    typedef ConjugateGaugeImpl<vComplexD,Nc> ConjugateGimplD; // Double
+
+  }
+}
+
+#endif
diff --git a/lib/qcd/action/gauge/PlaqPlusRectangleAction.h b/lib/qcd/action/gauge/PlaqPlusRectangleAction.h
index 0f39b8f3..4b19cafd 100644
--- a/lib/qcd/action/gauge/PlaqPlusRectangleAction.h
+++ b/lib/qcd/action/gauge/PlaqPlusRectangleAction.h
@@ -7,11 +7,11 @@ namespace Grid{
     ////////////////////////////////////////////////////////////////////////
     // PlaqPlusRectangleActoin
     ////////////////////////////////////////////////////////////////////////
-    template<class GaugeField>
-    class PlaqPlusRectangleAction : public Action<GaugeField> {
+    template<class Gimpl>
+    class PlaqPlusRectangleAction : public Action<typename Gimpl::GaugeField> {
     public:
 
-      typedef LorentzScalar<GaugeField> GaugeLinkField;
+      INHERIT_GIMPL_TYPES(Gimpl);
 
     private:
       RealD c_plaq;
@@ -25,8 +25,8 @@ namespace Grid{
       virtual RealD S(const GaugeField &U) {
 	RealD vol = U._grid->gSites();
 
-	RealD plaq = WilsonLoops<GaugeField>::avgPlaquette(U);
-	RealD rect = WilsonLoops<GaugeField>::avgRectangle(U);
+	RealD plaq = WilsonLoops<Gimpl>::avgPlaquette(U);
+	RealD rect = WilsonLoops<Gimpl>::avgRectangle(U);
 
 	RealD action=c_plaq*(1.0 -plaq)*(Nd*(Nd-1.0))*vol*0.5
 	            +c_rect*(1.0 -rect)*(Nd*(Nd-1.0))*vol;
@@ -46,7 +46,7 @@ namespace Grid{
 
 	for(int mu=0;mu<Nd;mu++){
 	  U[mu] = PeekIndex<LorentzIndex>(Umu,mu);
-	  WilsonLoops<GaugeField>::RectStapleDouble(U2[mu],U[mu],mu);
+	  WilsonLoops<Gimpl>::RectStapleDouble(U2[mu],U[mu],mu);
 	}
 
 	GaugeLinkField dSdU_mu(grid);
@@ -56,13 +56,11 @@ namespace Grid{
 
 	  // Staple in direction mu
 
-	  WilsonLoops<GaugeField>::Staple(staple,Umu,mu);
+	  WilsonLoops<Gimpl>::Staple(staple,Umu,mu);
 
 	  dSdU_mu = Ta(U[mu]*staple)*factor_p;
 
-	  //	  WilsonLoops<GaugeField>::RectStaple(staple,Umu,mu);
-
-	  WilsonLoops<GaugeField>::RectStapleOptimised(staple,U2,U,mu);
+	  WilsonLoops<Gimpl>::RectStaple(Umu,staple,U2,U,mu);
 
 	  dSdU_mu = dSdU_mu + Ta(U[mu]*staple)*factor_r;
 	  
@@ -78,31 +76,35 @@ namespace Grid{
     // RBC c1 parameterisation is not really RBC but don't have good
     // reference and we are happy to change name if prior use of this plaq coeff
     // parameterisation is made known to us. 
-    template<class GaugeField>
-    class RBCGaugeAction : public PlaqPlusRectangleAction<GaugeField> {
+    template<class Gimpl>
+    class RBCGaugeAction : public PlaqPlusRectangleAction<Gimpl> {
     public:
-      RBCGaugeAction(RealD beta,RealD c1) : PlaqPlusRectangleAction<GaugeField>(beta*(1.0-8.0*c1), beta*c1) {
+      INHERIT_GIMPL_TYPES(Gimpl);
+      RBCGaugeAction(RealD beta,RealD c1) : PlaqPlusRectangleAction<Gimpl>(beta*(1.0-8.0*c1), beta*c1) {
       };
     };
 
-    template<class GaugeField>
-    class IwasakiGaugeAction : public RBCGaugeAction<GaugeField> {
+    template<class Gimpl>
+    class IwasakiGaugeAction : public RBCGaugeAction<Gimpl> {
     public:
-      IwasakiGaugeAction(RealD beta) : RBCGaugeAction<GaugeField>(beta,-0.331) {
+      INHERIT_GIMPL_TYPES(Gimpl);
+      IwasakiGaugeAction(RealD beta) : RBCGaugeAction<Gimpl>(beta,-0.331) {
       };
     };
 
-    template<class GaugeField>
-    class SymanzikGaugeAction : public RBCGaugeAction<GaugeField> {
+    template<class Gimpl>
+    class SymanzikGaugeAction : public RBCGaugeAction<Gimpl> {
     public:
-      SymanzikGaugeAction(RealD beta) : RBCGaugeAction<GaugeField>(beta,-1.0/12.0) {
+      INHERIT_GIMPL_TYPES(Gimpl);
+      SymanzikGaugeAction(RealD beta) : RBCGaugeAction<Gimpl>(beta,-1.0/12.0) {
       };
     };
 
-    template<class GaugeField>
-    class DBW2GaugeAction : public RBCGaugeAction<GaugeField> {
+    template<class Gimpl>
+    class DBW2GaugeAction : public RBCGaugeAction<Gimpl> {
     public:
-      DBW2GaugeAction(RealD beta) : RBCGaugeAction<GaugeField>(beta,-1.4067) {
+      INHERIT_GIMPL_TYPES(Gimpl);
+      DBW2GaugeAction(RealD beta) : RBCGaugeAction<Gimpl>(beta,-1.4067) {
       };
     };
 
diff --git a/lib/qcd/action/gauge/WilsonGaugeAction.h b/lib/qcd/action/gauge/WilsonGaugeAction.h
index 2f8f355e..0f15245e 100644
--- a/lib/qcd/action/gauge/WilsonGaugeAction.h
+++ b/lib/qcd/action/gauge/WilsonGaugeAction.h
@@ -7,11 +7,13 @@ namespace Grid{
     ////////////////////////////////////////////////////////////////////////
     // Wilson Gauge Action .. should I template the Nc etc..
     ////////////////////////////////////////////////////////////////////////
-    template<class GaugeField>
-    class WilsonGaugeAction : public Action<GaugeField> {
+    template<class Gimpl>
+    class WilsonGaugeAction : public Action<typename Gimpl::GaugeField> {
     public:
 
-      typedef LorentzScalar<GaugeField> GaugeLinkField;
+      INHERIT_GIMPL_TYPES(Gimpl);
+
+      //      typedef LorentzScalar<GaugeField> GaugeLinkField;
 
     private:
       RealD beta;
@@ -21,7 +23,7 @@ namespace Grid{
       virtual void refresh(const GaugeField &U, GridParallelRNG& pRNG) {}; // noop as no pseudoferms
       
       virtual RealD S(const GaugeField &U) {
-	RealD plaq = WilsonLoops<GaugeField>::avgPlaquette(U);
+	RealD plaq = WilsonLoops<Gimpl>::avgPlaquette(U);
 	RealD vol = U._grid->gSites();
 	RealD action=beta*(1.0 -plaq)*(Nd*(Nd-1.0))*vol*0.5;
 	return action;
@@ -41,7 +43,7 @@ namespace Grid{
 	  Umu = PeekIndex<LorentzIndex>(U,mu);
 
 	  // Staple in direction mu
-	  WilsonLoops<GaugeField>::Staple(dSdU_mu,U,mu);
+	  WilsonLoops<Gimpl>::Staple(dSdU_mu,U,mu);
 	  dSdU_mu = Ta(Umu*dSdU_mu)*factor;
 	  PokeIndex<LorentzIndex>(dSdU, dSdU_mu, mu);
 	}
diff --git a/lib/qcd/hmc/HMC.h b/lib/qcd/hmc/HMC.h
index f8699892..ebf6fdd5 100644
--- a/lib/qcd/hmc/HMC.h
+++ b/lib/qcd/hmc/HMC.h
@@ -39,11 +39,12 @@ namespace Grid{
       virtual void TrajectoryComplete (int traj, GaugeField &U, GridSerialRNG &sRNG, GridParallelRNG & pRNG )=0;
     };
 
-    template<class GaugeField> 
-    class PlaquetteLogger : public HmcObservable<GaugeField> {
+    template<class Gimpl> 
+    class PlaquetteLogger : public HmcObservable<typename Gimpl::GaugeField> {
     private:
       std::string Stem;
     public:
+      INHERIT_GIMPL_TYPES(Gimpl);
       PlaquetteLogger(std::string cf) {
         Stem  = cf;
       };
@@ -52,11 +53,16 @@ namespace Grid{
       {
 	  std::string file;   { std::ostringstream os; os << Stem     <<"."<< traj; file = os.str(); }
 	  std::ofstream of(file);
-	  RealD plaq = WilsonLoops<GaugeField>::avgPlaquette(U);
-	  RealD rect = WilsonLoops<GaugeField>::avgRectangle(U);
-	  of << plaq << " " << rect << std::endl;
-	  std::cout<< GridLogMessage<< "Plaquette for trajectory "<< traj << " is " << plaq <<std::endl;
-	  std::cout<< GridLogMessage<< "Rectangle for trajectory "<< traj << " is " << rect <<std::endl;
+
+	  RealD peri_plaq = WilsonLoops<PeriodicGimplR>::avgPlaquette(U);
+	  RealD peri_rect = WilsonLoops<PeriodicGimplR>::avgRectangle(U);
+
+	  RealD impl_plaq = WilsonLoops<Gimpl>::avgPlaquette(U);
+	  RealD impl_rect = WilsonLoops<Gimpl>::avgRectangle(U);
+
+	  of << traj<<" "<< impl_plaq << " " << impl_rect << "  "<< peri_plaq<<" "<<peri_rect<<std::endl;
+	  std::cout<< GridLogMessage<< "traj"<<" "<< "plaq " << " " << " rect  " << "  "<< "peri_plaq" <<" "<<"peri_rect"<<std::endl;
+	  std::cout<< GridLogMessage<< traj<<" "<< impl_plaq << " " << impl_rect << "  "<< peri_plaq<<" "<<peri_rect<<std::endl;
       }
     };
 
diff --git a/lib/qcd/hmc/HmcRunner.h b/lib/qcd/hmc/HmcRunner.h
index 4f3be698..184ce52b 100644
--- a/lib/qcd/hmc/HmcRunner.h
+++ b/lib/qcd/hmc/HmcRunner.h
@@ -5,12 +5,15 @@ namespace Grid{
   namespace QCD{
 
 
-class NerscHmcRunner {
+template<class Gimpl>
+class NerscHmcRunnerTemplate {
 public:
 
+  INHERIT_GIMPL_TYPES(Gimpl);
+
   enum StartType_t { ColdStart, HotStart, TepidStart, CheckpointStart };
 
-  ActionSet<LatticeGaugeField> TheAction;
+  ActionSet<GaugeField> TheAction;
 
   GridCartesian         * UGrid   ;
   GridCartesian         * FGrid   ;
@@ -52,13 +55,13 @@ public:
     }
 
     // Create integrator
-    typedef MinimumNorm2<LatticeGaugeField>  IntegratorType;// change here to change the algorithm
+    typedef MinimumNorm2<GaugeField>  IntegratorType;// change here to change the algorithm
     IntegratorParameters MDpar(20);
     IntegratorType MDynamics(UGrid,MDpar, TheAction);
 
     // Checkpoint strategy
-    NerscHmcCheckpointer<LatticeGaugeField> Checkpoint(std::string("ckpoint_lat"),std::string("ckpoint_rng"),1);
-    PlaquetteLogger<LatticeGaugeField> PlaqLog(std::string("plaq"));
+    NerscHmcCheckpointer<Gimpl> Checkpoint(std::string("ckpoint_lat"),std::string("ckpoint_rng"),1);
+    PlaquetteLogger<Gimpl>      PlaqLog(std::string("plaq"));
 
     HMCparameters HMCpar;
     HMCpar.StartTrajectory = StartTraj;
@@ -99,7 +102,7 @@ public:
       Checkpoint.CheckpointRestore(StartTraj, U, sRNG, pRNG);
     }
 
-    HybridMonteCarlo<LatticeGaugeField,IntegratorType>  HMC(HMCpar, MDynamics,sRNG,pRNG,U);
+    HybridMonteCarlo<GaugeField,IntegratorType>  HMC(HMCpar, MDynamics,sRNG,pRNG,U);
     HMC.AddObservable(&Checkpoint);
     HMC.AddObservable(&PlaqLog);
     
@@ -109,5 +112,18 @@ public:
   }
   
 };
+
+ typedef NerscHmcRunnerTemplate<PeriodicGimplR> NerscHmcRunner;
+ typedef NerscHmcRunnerTemplate<PeriodicGimplF> NerscHmcRunnerF;
+ typedef NerscHmcRunnerTemplate<PeriodicGimplD> NerscHmcRunnerD;
+
+ typedef NerscHmcRunnerTemplate<PeriodicGimplR> PeriodicNerscHmcRunner;
+ typedef NerscHmcRunnerTemplate<PeriodicGimplF> PeriodicNerscHmcRunnerF;
+ typedef NerscHmcRunnerTemplate<PeriodicGimplD> PeriodicNerscHmcRunnerD;
+
+ typedef NerscHmcRunnerTemplate<ConjugateGimplR> ConjugateNerscHmcRunner;
+ typedef NerscHmcRunnerTemplate<ConjugateGimplF> ConjugateNerscHmcRunnerF;
+ typedef NerscHmcRunnerTemplate<ConjugateGimplD> ConjugateNerscHmcRunnerD;
+
 }}
 #endif
diff --git a/lib/qcd/hmc/NerscCheckpointer.h b/lib/qcd/hmc/NerscCheckpointer.h
index 2309eb0e..1f49ec94 100644
--- a/lib/qcd/hmc/NerscCheckpointer.h
+++ b/lib/qcd/hmc/NerscCheckpointer.h
@@ -10,13 +10,15 @@ namespace Grid{
   namespace QCD{
     
 
-    template<class GaugeField> 
-    class NerscHmcCheckpointer : public HmcObservable<GaugeField> {
+    template<class Gimpl> 
+    class NerscHmcCheckpointer : public HmcObservable<typename Gimpl::GaugeField> {
     private:
       std::string configStem;
       std::string rngStem;
       int SaveInterval;
     public:
+      INHERIT_GIMPL_TYPES(Gimpl);
+
       NerscHmcCheckpointer(std::string cf, std::string rn,int savemodulo) {
         configStem  = cf;
         rngStem     = rn;
diff --git a/lib/qcd/utils/CovariantCshift.h b/lib/qcd/utils/CovariantCshift.h
index 8bc0986f..aa3c3789 100644
--- a/lib/qcd/utils/CovariantCshift.h
+++ b/lib/qcd/utils/CovariantCshift.h
@@ -1,25 +1,102 @@
 #ifndef QCD_UTILS_COVARIANT_CSHIFT_H
 #define QCD_UTILS_COVARIANT_CSHIFT_H
+
 namespace Grid {
 namespace QCD {
 ////////////////////////////////////////////////////////////////////////
 // Low performance implementation of CovariantCshift API
 ////////////////////////////////////////////////////////////////////////
 // Make these members of an Impl class for BC's.
-template<class covariant,class gauge> Lattice<covariant> CovShiftForward(const Lattice<gauge> &Link, 
+
+namespace PeriodicBC { 
+
+  template<class covariant,class gauge> Lattice<covariant> CovShiftForward(const Lattice<gauge> &Link, 
 									    int mu,
 									    const Lattice<covariant> &field)
-{
-  return Link*Cshift(field,mu,1);// moves towards negative mu
+  {
+    return Link*Cshift(field,mu,1);// moves towards negative mu
+  }
+  template<class covariant,class gauge> Lattice<covariant> CovShiftBackward(const Lattice<gauge> &Link, 
+									    int mu,
+									    const Lattice<covariant> &field)
+  {
+    Lattice<covariant> tmp(field._grid);
+    tmp = adj(Link)*field;
+    return Cshift(tmp,mu,-1);// moves towards positive mu
+  }
 }
-template<class covariant,class gauge> Lattice<covariant> CovShiftBackward(const Lattice<gauge> &Link, 
-									     int mu,
-									     const Lattice<covariant> &field)
-{
-  Lattice<covariant> tmp(field._grid);
-  tmp = adj(Link)*field;
-  return Cshift(tmp,mu,-1);// moves towards positive mu
+
+
+namespace GparityBC { 
+
+  // Must give right answers across boundary
+  //     <----
+  //       --    
+  //      |  |              
+  // xxxxxxxxxxxxxxxxxxxx
+  //      |  | 
+  //  
+  //     Stap= Cshift(GImpl::CovShiftForward(U[nu],nu,
+  //	    	      GImpl::CovShiftForward(U[nu],nu,
+  //                  GImpl::CovShiftBackward(U[mu],mu,
+  //                  GImpl::CovShiftBackward(U[nu],nu,
+  //		      GImpl::CovShiftIdentityBackward(U[nu],nu,-1))))) , mu, 1);
+  //  
+  //                      U  U^* U^* U^T U^adj =  U  (U U U^dag U^T )^*
+  //                                           =  U  (U U U^dag)^* ( U^T )^*
+  //
+  // So covariant shift rule: conjugate inward shifted plane when crossing boundary applies.
+  //
+  // This conjugate should be applied to BOTH the link and the covariant field on backward shift
+  // boundary wrap.
+  // 
+  //      |  |              
+  // xxxxxxxxxxxxxxxxx         
+  //      |  | <---- this link is conjugated, and the path leading into it. Segment crossing in and out is double conjugated.
+  //       -- 
+  //    ------->
+  template<class covariant,class gauge> Lattice<covariant> CovShiftForward(const Lattice<gauge> &Link, 
+									    int mu,
+									    const Lattice<covariant> &field)
+  {
+    GridBase * grid = Link._grid;
+
+    int Lmu = grid->GlobalDimensions()[mu]-1;
+
+    conformable(field,Link);
+
+    Lattice<iScalar<vInteger> > coor(grid);    LatticeCoordinate(coor,mu);
+
+    Lattice<covariant> field_bc = Cshift(field,mu,1);// moves towards negative mu;
+
+    field_bc = where(coor==Lmu,conjugate(field_bc),field_bc);
+    //    std::cout<<"Gparity::CovCshiftForward mu="<<mu<<std::endl;
+    return Link*field_bc;
+  }
+
+  template<class covariant,class gauge> Lattice<covariant> CovShiftBackward(const Lattice<gauge> &Link, 
+									    int mu,
+									    const Lattice<covariant> &field)
+  {
+    GridBase * grid = field._grid;
+
+    int Lmu = grid->GlobalDimensions()[mu]-1;
+
+    conformable(field,Link);
+
+    Lattice<iScalar<vInteger> > coor(grid);    LatticeCoordinate(coor,mu);
+
+    Lattice<covariant> tmp(grid);
+
+    tmp = adj(Link)*field;
+    tmp = where(coor==Lmu,conjugate(tmp),tmp);
+    //    std::cout<<"Gparity::CovCshiftBackward mu="<<mu<<std::endl;
+    return Cshift(tmp,mu,-1);// moves towards positive mu
+  }
+
+
 }
 
+
 }}
 #endif
diff --git a/lib/qcd/utils/WilsonLoops.h b/lib/qcd/utils/WilsonLoops.h
index 36fbe51c..c3277196 100644
--- a/lib/qcd/utils/WilsonLoops.h
+++ b/lib/qcd/utils/WilsonLoops.h
@@ -4,18 +4,27 @@ namespace Grid {
 namespace QCD {
 
 // Common wilson loop observables
-template<class GaugeLorentz>
-class WilsonLoops {
+template<class Gimpl>
+class WilsonLoops : public Gimpl {
 public:
+  
+  INHERIT_GIMPL_TYPES(Gimpl);
 
-  typedef LorentzScalar<GaugeLorentz> GaugeMat;
+  typedef typename Gimpl::GaugeLinkField GaugeMat;
+  typedef typename Gimpl::GaugeField     GaugeLorentz;
 
   //////////////////////////////////////////////////
   // directed plaquette oriented in mu,nu plane
   //////////////////////////////////////////////////
   static void dirPlaquette(GaugeMat &plaq,const std::vector<GaugeMat> &U, const int mu, const int nu)
   {
-    plaq=CovShiftForward(U[mu],mu,U[nu])*adj(CovShiftForward(U[nu],nu,U[mu]));
+    // Annoyingly, must use either scope resolution to find dependent base class, 
+    // or this-> ; there is no "this" in a static method. This forces explicit Gimpl scope
+    // resolution throughout the usage in this file, and rather defeats the purpose of deriving
+    // from Gimpl.
+    plaq= Gimpl::CovShiftBackward(U[mu],mu,
+	  Gimpl::CovShiftBackward(U[nu],nu,
+          Gimpl::CovShiftForward (U[mu],mu,U[nu])));
   }
   //////////////////////////////////////////////////
   // trace of directed plaquette oriented in mu,nu plane
@@ -98,10 +107,10 @@ public:
     for(int d=0;d<Nd;d++){
       U[d] = PeekIndex<LorentzIndex>(Umu,d);
     }
-
     staple = zero;
     GaugeMat tmp(grid);
 
+    
     for(int nu=0;nu<Nd;nu++){
 
       if(nu != mu) {
@@ -110,35 +119,24 @@ public:
       // ^
       // |__>  nu
 
-      //    __                                         __
-      //      |         |                          
-      //    __|  =    __|                         *
+      //    __ 
+      //      |
+      //    __|
       //
 
-	//	staple   += CovShiftForward(U[nu],nu,U[mu])*Cshift(adj(U[nu]),mu,+1);
-	staple   += 
-	  Cshift(CovShiftForward (U[nu],nu, 
-		 CovShiftBackward(U[mu],mu,Cshift(adj(U[nu]),nu,-1))),mu,1);
+	staple+=Gimpl::ShiftStaple(
+	        Gimpl::CovShiftForward (U[nu],nu, 
+		Gimpl::CovShiftBackward(U[mu],mu,
+		Gimpl::CovShiftIdentityBackward(U[nu],nu))),mu);
 
-
-	// Unu(x) Umu(x+nu) Unu^dag(x+mu) ; left mult by Umu^dag(x) to close ring
-
-      //
-      //  __        __                         
-      // |         |                          
-      // |__     = |                              *       __
+      //  __ 
+      // |   
+      // |__ 
       //
       //
-	staple   += 
-          Cshift(
-	  CovShiftBackward(U[nu],nu,		  		  
-	  CovShiftBackward(U[mu],mu,U[nu])),mu,1);
-
-      //	tmp    = CovShiftForward (U[mu],mu,U[nu]);
-      //	staple+= CovShiftBackward(U[nu],nu,tmp);
-
-	// Unu^dag(x-nu) Umu(x-nu) Unu(x+mu-nu) ; left mult by Umu^dag(x) to close ring.
-      
+	staple+=Gimpl::ShiftStaple(  
+                Gimpl::CovShiftBackward(U[nu],nu,		  		  
+		Gimpl::CovShiftBackward(U[mu],mu,U[nu])),mu);
       }
     }
   }
@@ -148,11 +146,11 @@ public:
   //////////////////////////////////////////////////////
   static void dirRectangle(GaugeMat &rect,const std::vector<GaugeMat> &U, const int mu, const int nu)
   {
-    rect = CovShiftForward(U[mu],mu,CovShiftForward(U[mu],mu,U[nu]))* // ->->|
-	adj(CovShiftForward(U[nu],nu,CovShiftForward(U[mu],mu,U[mu]))) ;
+    rect =  Gimpl::CovShiftForward(U[mu],mu,Gimpl::CovShiftForward(U[mu],mu,U[nu]))* // ->->|
+	adj(Gimpl::CovShiftForward(U[nu],nu,Gimpl::CovShiftForward(U[mu],mu,U[mu]))) ;
     rect = rect + 
-      CovShiftForward(U[mu],mu,CovShiftForward(U[nu],nu,U[nu]))* // ->||
-      adj(CovShiftForward(U[nu],nu,CovShiftForward(U[nu],nu,U[mu]))) ;
+          Gimpl::CovShiftForward(U[mu],mu,Gimpl::CovShiftForward(U[nu],nu,U[nu]))* // ->||
+      adj(Gimpl::CovShiftForward(U[nu],nu,Gimpl::CovShiftForward(U[nu],nu,U[mu]))) ;
   }
   static void traceDirRectangle(LatticeComplex &rect, const std::vector<GaugeMat> &U, const int mu, const int nu)
   {
@@ -209,6 +207,12 @@ public:
   static void RectStapleDouble(GaugeMat &U2,const GaugeMat & U,int mu){
     U2 = U * Cshift(U,mu,1);
   }
+
+  ////////////////////////////////////////////////////////////////////////////
+  // Hop by two optimisation strategy does not work nicely with Gparity. (could do,
+  // but need to track two deep where cross boundary and apply a conjugation).
+  // Must differentiate this in Gimpl, and use Gimpl::isPeriodicGaugeField to do so .
+  ////////////////////////////////////////////////////////////////////////////
   static void RectStapleOptimised(GaugeMat &Stap,std::vector<GaugeMat> &U2,std::vector<GaugeMat> &U,int mu){
 
     Stap = zero;
@@ -227,14 +231,14 @@ public:
 	tmp = adj(U2[mu])*tmp;
 	tmp = Cshift(tmp,mu,-2);
 
-	Staple2x1 = CovShiftForward (U[nu],nu,tmp);
+	Staple2x1 = Gimpl::CovShiftForward (U[nu],nu,tmp);
 
 
 	// Down staple
 	//             |___ ___|
 	//
 	tmp = adj(U2[mu])*U[nu];
-	Staple2x1+= CovShiftBackward(U[nu],nu,Cshift(tmp,mu,-2));
+	Staple2x1+= Gimpl::CovShiftBackward(U[nu],nu,Cshift(tmp,mu,-2));
 
 
 	//              ___ ___
@@ -242,7 +246,7 @@ public:
 	//             |___ ___|
 	//
 
-	Stap+= Cshift(CovShiftForward (U[mu],mu,Staple2x1),mu,1);
+	Stap+= Cshift(Gimpl::CovShiftForward (U[mu],mu,Staple2x1),mu,1);
 
 	//              ___ ___
 	//             |___    |
@@ -259,7 +263,7 @@ public:
 	//      |  | 
 	
 	tmp = Cshift(adj(U2[nu]),nu,-2);
-	tmp = CovShiftBackward(U[mu],mu,tmp);
+	tmp = Gimpl::CovShiftBackward(U[mu],mu,tmp);
 	tmp = U2[nu]*Cshift(tmp,nu,2);
 	Stap+= Cshift(tmp, mu, 1);
 
@@ -268,7 +272,7 @@ public:
 	//      |  | 
 	//       -- 
 	
-	tmp = CovShiftBackward(U[mu],mu,U2[nu]);
+	tmp = Gimpl::CovShiftBackward(U[mu],mu,U2[nu]);
 	tmp = adj(U2[nu])*tmp;
 	tmp = Cshift(tmp,nu,-2);
 	Stap+=Cshift(tmp, mu, 1);
@@ -276,23 +280,20 @@ public:
 
 
   }
-  static void RectStaple(GaugeMat &Stap,const GaugeLorentz &Umu,int mu){
 
-#if 1
+  static void RectStaple(GaugeMat &Stap,const GaugeLorentz & Umu,int mu)
+  {
     RectStapleUnoptimised(Stap,Umu,mu);
-#else
-    GridBase *grid = Umu._grid;
-
-    std::vector<GaugeMat> U (Nd,grid);
-    std::vector<GaugeMat> U2(Nd,grid);
-
-    for(int mu=0;mu<Nd;mu++){
-      U[mu] = PeekIndex<LorentzIndex>(Umu,mu);
-      RectStapleDouble(U2[mu],U[mu],mu);
+  }
+  static void RectStaple(const GaugeLorentz & Umu,GaugeMat &Stap,
+			 std::vector<GaugeMat> &U2,
+			 std::vector<GaugeMat> &U, int mu)
+  {
+    if ( Gimpl::isPeriodicGaugeField() ){ 
+      RectStapleOptimised(Stap,U2,U,mu);
+    } else {
+      RectStapleUnoptimised(Stap,Umu,mu);
     }
-    
-    RectStapleOptimised(Stap,U2,U,mu);
-#endif
   }
 
   static void RectStapleUnoptimised(GaugeMat &Stap,const GaugeLorentz &Umu,int mu){
@@ -310,46 +311,51 @@ public:
     //           __ ___ 
     //          |    __ |
     //
-    Stap+= Cshift(CovShiftForward (U[mu],mu,
-		  CovShiftForward (U[nu],nu,
-		  CovShiftBackward(U[mu],mu,
-                  CovShiftBackward(U[mu],mu,
-	          Cshift(adj(U[nu]),nu,-1))))) , mu, 1);
+    Stap+= Gimpl::ShiftStaple(
+		  Gimpl::CovShiftForward (U[mu],mu,
+		  Gimpl::CovShiftForward (U[nu],nu,
+		  Gimpl::CovShiftBackward(U[mu],mu,
+                  Gimpl::CovShiftBackward(U[mu],mu,
+		  Gimpl::CovShiftIdentityBackward(U[nu],nu))))) , mu);
 
     //              __ 
     //          |__ __ |
 
-    Stap+= Cshift(CovShiftForward (U[mu],mu,
-		  CovShiftBackward(U[nu],nu,
-		  CovShiftBackward(U[mu],mu,
-                  CovShiftBackward(U[mu],mu, U[nu])))) , mu, 1);
+    Stap+= Gimpl::ShiftStaple(
+                  Gimpl::CovShiftForward (U[mu],mu,
+		  Gimpl::CovShiftBackward(U[nu],nu,
+		  Gimpl::CovShiftBackward(U[mu],mu,
+                  Gimpl::CovShiftBackward(U[mu],mu, U[nu])))) , mu);
 
     //           __ 
     //          |__ __ |
 
-    Stap+= Cshift(CovShiftBackward(U[nu],nu,
-		  CovShiftBackward(U[mu],mu,
-		  CovShiftBackward(U[mu],mu,
-		  CovShiftForward(U[nu],nu,U[mu])))) , mu, 1);
+    Stap+= Gimpl::ShiftStaple(
+		  Gimpl::CovShiftBackward(U[nu],nu,
+		  Gimpl::CovShiftBackward(U[mu],mu,
+		  Gimpl::CovShiftBackward(U[mu],mu,
+		  Gimpl::CovShiftForward(U[nu],nu,U[mu])))) , mu);
 
     //           __ ___ 
     //          |__    |
 
-     Stap+= Cshift(CovShiftForward (U[nu],nu,
-	           CovShiftBackward(U[mu],mu,
-                   CovShiftBackward(U[mu],mu,
-                   CovShiftBackward(U[nu],nu,U[mu])))) , mu, 1);
+    Stap+= Gimpl::ShiftStaple(
+		   Gimpl::CovShiftForward (U[nu],nu,
+	           Gimpl::CovShiftBackward(U[mu],mu,
+                   Gimpl::CovShiftBackward(U[mu],mu,
+                   Gimpl::CovShiftBackward(U[nu],nu,U[mu])))) , mu);
 
      //       --    
      //      |  |              
      //          
      //      |  | 
-
-     Stap+= Cshift(CovShiftForward(U[nu],nu,
-		   CovShiftForward(U[nu],nu,
-                   CovShiftBackward(U[mu],mu,
-                   CovShiftBackward(U[nu],nu,
-		   Cshift(adj(U[nu]),nu,-1))))) , mu, 1);
+     
+    Stap+= Gimpl::ShiftStaple(
+		   Gimpl::CovShiftForward(U[nu],nu,
+		   Gimpl::CovShiftForward(U[nu],nu,
+                   Gimpl::CovShiftBackward(U[mu],mu,
+                   Gimpl::CovShiftBackward(U[nu],nu,
+		   Gimpl::CovShiftIdentityBackward(U[nu],nu))))) , mu);
 
 
      //      |  |              
@@ -357,10 +363,11 @@ public:
      //      |  | 
      //       -- 
      
-     Stap+= Cshift(CovShiftBackward(U[nu],nu,
-		   CovShiftBackward(U[nu],nu,
-                   CovShiftBackward(U[mu],mu,
-                   CovShiftForward (U[nu],nu,U[nu])))) , mu, 1);
+    Stap+= Gimpl::ShiftStaple(
+		   Gimpl::CovShiftBackward(U[nu],nu,
+		   Gimpl::CovShiftBackward(U[nu],nu,
+                   Gimpl::CovShiftBackward(U[mu],mu,
+                   Gimpl::CovShiftForward (U[nu],nu,U[nu])))) , mu);
     }}
   }
 
@@ -368,10 +375,10 @@ public:
 };
 
 
- typedef WilsonLoops<LatticeGaugeField> ColourWilsonLoops;
- typedef WilsonLoops<LatticeGaugeField> U1WilsonLoops;
- typedef WilsonLoops<LatticeGaugeField> SU2WilsonLoops;
- typedef WilsonLoops<LatticeGaugeField> SU3WilsonLoops;
+ typedef WilsonLoops<PeriodicGimplR> ColourWilsonLoops;
+ typedef WilsonLoops<PeriodicGimplR> U1WilsonLoops;
+ typedef WilsonLoops<PeriodicGimplR> SU2WilsonLoops;
+ typedef WilsonLoops<PeriodicGimplR> SU3WilsonLoops;
 
 }}
 
diff --git a/tests/Make.inc b/tests/Make.inc
index f69b9e5b..f7c83671 100644
--- a/tests/Make.inc
+++ b/tests/Make.inc
@@ -1,5 +1,5 @@
 
-bin_PROGRAMS = Test_GaugeAction Test_RectPlaq Test_cayley_cg Test_cayley_coarsen_support Test_cayley_even_odd Test_cayley_ldop_cr Test_cf_coarsen_support Test_cf_cr_unprec Test_cheby Test_contfrac_cg Test_contfrac_even_odd Test_contfrac_force Test_cshift Test_cshift_red_black Test_dwf_cg_prec Test_dwf_cg_schur Test_dwf_cg_unprec Test_dwf_cr_unprec Test_dwf_even_odd Test_dwf_force Test_dwf_fpgcr Test_dwf_gpforce Test_dwf_hdcr Test_dwf_lanczos Test_gamma Test_gparity Test_gpdwf_force Test_gpwilson_even_odd Test_hmc_EODWFRatio Test_hmc_EODWFRatio_Gparity Test_hmc_EOWilsonFermionGauge Test_hmc_EOWilsonRatio Test_hmc_IwasakiGauge Test_hmc_RectGauge Test_hmc_WilsonFermionGauge Test_hmc_WilsonGauge Test_hmc_WilsonRatio Test_lie_generators Test_main Test_multishift_sqrt Test_nersc_io Test_partfrac_force Test_quenched_update Test_rect_force Test_remez Test_rhmc_EOWilson1p1 Test_rhmc_EOWilsonRatio Test_rhmc_Wilson1p1 Test_rhmc_WilsonRatio Test_rng Test_rng_fixed Test_serialisation Test_simd Test_stencil Test_synthetic_lanczos Test_wilson_cg_prec Test_wilson_cg_schur Test_wilson_cg_unprec Test_wilson_cr_unprec Test_wilson_even_odd Test_wilson_force Test_wilson_force_phiMdagMphi Test_wilson_force_phiMphi Test_wilson_tm_even_odd 
+bin_PROGRAMS = Test_GaugeAction Test_RectPlaq Test_cayley_cg Test_cayley_coarsen_support Test_cayley_even_odd Test_cayley_ldop_cr Test_cf_coarsen_support Test_cf_cr_unprec Test_cheby Test_contfrac_cg Test_contfrac_even_odd Test_contfrac_force Test_cshift Test_cshift_red_black Test_dwf_cg_prec Test_dwf_cg_schur Test_dwf_cg_unprec Test_dwf_cr_unprec Test_dwf_even_odd Test_dwf_force Test_dwf_fpgcr Test_dwf_gpforce Test_dwf_hdcr Test_dwf_lanczos Test_gamma Test_gp_rect_force Test_gparity Test_gpdwf_force Test_gpwilson_even_odd Test_hmc_EODWFRatio Test_hmc_EODWFRatio_Gparity Test_hmc_EOWilsonFermionGauge Test_hmc_EOWilsonRatio Test_hmc_GparityIwasakiGauge Test_hmc_GparityWilsonGauge Test_hmc_IwasakiGauge Test_hmc_RectGauge Test_hmc_WilsonFermionGauge Test_hmc_WilsonGauge Test_hmc_WilsonRatio Test_lie_generators Test_main Test_multishift_sqrt Test_nersc_io Test_partfrac_force Test_quenched_update Test_rect_force Test_remez Test_rhmc_EOWilson1p1 Test_rhmc_EOWilsonRatio Test_rhmc_Wilson1p1 Test_rhmc_WilsonRatio Test_rng Test_rng_fixed Test_serialisation Test_simd Test_stencil Test_synthetic_lanczos Test_wilson_cg_prec Test_wilson_cg_schur Test_wilson_cg_unprec Test_wilson_cr_unprec Test_wilson_even_odd Test_wilson_force Test_wilson_force_phiMdagMphi Test_wilson_force_phiMphi Test_wilson_tm_even_odd 
 
 
 Test_GaugeAction_SOURCES=Test_GaugeAction.cc
@@ -102,6 +102,10 @@ Test_gamma_SOURCES=Test_gamma.cc
 Test_gamma_LDADD=-lGrid
 
 
+Test_gp_rect_force_SOURCES=Test_gp_rect_force.cc
+Test_gp_rect_force_LDADD=-lGrid
+
+
 Test_gparity_SOURCES=Test_gparity.cc
 Test_gparity_LDADD=-lGrid
 
@@ -130,6 +134,14 @@ Test_hmc_EOWilsonRatio_SOURCES=Test_hmc_EOWilsonRatio.cc
 Test_hmc_EOWilsonRatio_LDADD=-lGrid
 
 
+Test_hmc_GparityIwasakiGauge_SOURCES=Test_hmc_GparityIwasakiGauge.cc
+Test_hmc_GparityIwasakiGauge_LDADD=-lGrid
+
+
+Test_hmc_GparityWilsonGauge_SOURCES=Test_hmc_GparityWilsonGauge.cc
+Test_hmc_GparityWilsonGauge_LDADD=-lGrid
+
+
 Test_hmc_IwasakiGauge_SOURCES=Test_hmc_IwasakiGauge.cc
 Test_hmc_IwasakiGauge_LDADD=-lGrid
 
diff --git a/tests/Test_GaugeAction.cc b/tests/Test_GaugeAction.cc
index cc593c71..150af081 100644
--- a/tests/Test_GaugeAction.cc
+++ b/tests/Test_GaugeAction.cc
@@ -87,7 +87,7 @@ int main (int argc, char ** argv)
   Plaq = zero;
   for(int mu=1;mu<Nd;mu++){
     for(int nu=0;nu<mu;nu++){
-      Plaq = Plaq + trace(CovShiftForward(U[mu],mu,U[nu])*adj(CovShiftForward(U[nu],nu,U[mu])));
+      Plaq = Plaq + trace(PeriodicBC::CovShiftForward(U[mu],mu,U[nu])*adj(PeriodicBC::CovShiftForward(U[nu],nu,U[mu])));
     }
   }
   
diff --git a/tests/Test_RectPlaq.cc b/tests/Test_RectPlaq.cc
index 57c85211..cf22ba21 100644
--- a/tests/Test_RectPlaq.cc
+++ b/tests/Test_RectPlaq.cc
@@ -34,11 +34,11 @@ void RectPlaq(const std::vector<LatticeColourMatrix> &U, LatticeComplex &RectPla
   for(int mu=1;mu<Nd;mu++){
     for(int nu=0;nu<mu;nu++){
       RectPlaqValue = RectPlaqValue + trace(
-             CovShiftForward(U[mu],mu,CovShiftForward(U[mu],mu,U[nu]))* // ->->|
-         adj(CovShiftForward(U[nu],nu,CovShiftForward(U[mu],mu,U[mu]))) );
+             PeriodicBC::CovShiftForward(U[mu],mu,PeriodicBC::CovShiftForward(U[mu],mu,U[nu]))* // ->->|
+         adj(PeriodicBC::CovShiftForward(U[nu],nu,PeriodicBC::CovShiftForward(U[mu],mu,U[mu]))) );
       RectPlaqValue = RectPlaqValue + trace(
-             CovShiftForward(U[mu],mu,CovShiftForward(U[nu],nu,U[nu]))* // ->||
-         adj(CovShiftForward(U[nu],nu,CovShiftForward(U[nu],nu,U[mu]))) );
+             PeriodicBC::CovShiftForward(U[mu],mu,PeriodicBC::CovShiftForward(U[nu],nu,U[nu]))* // ->||
+         adj(PeriodicBC::CovShiftForward(U[nu],nu,PeriodicBC::CovShiftForward(U[nu],nu,U[mu]))) );
     }
   }
 }
@@ -87,7 +87,7 @@ int main (int argc, char ** argv)
   Plaq = zero;
   for(int mu=1;mu<Nd;mu++){
     for(int nu=0;nu<mu;nu++){
-      Plaq = Plaq + trace(CovShiftForward(U[mu],mu,U[nu])*adj(CovShiftForward(U[nu],nu,U[mu])));
+      Plaq = Plaq + trace(PeriodicBC::CovShiftForward(U[mu],mu,U[nu])*adj(PeriodicBC::CovShiftForward(U[nu],nu,U[mu])));
     }
   }
 
@@ -134,7 +134,7 @@ int main (int argc, char ** argv)
 /*     
            (x)  ---> --->  : U(x,mu)*U(x+mu, mu)
 */
-           left_2=  CovShiftForward(U[mu],mu,U[mu]);
+           left_2=  PeriodicBC::CovShiftForward(U[mu],mu,U[mu]);
 /*
   upper_l =   <---- <---
                          ^
@@ -144,7 +144,7 @@ int main (int argc, char ** argv)
 */
 
           tmp=Cshift(U[nu],mu,2);
-          upper_l=  CovShiftForward(tmp,nu,adj(left_2)); // i.e. upper_l
+          upper_l=  PeriodicBC::CovShiftForward(tmp,nu,adj(left_2)); // i.e. upper_l
 /*
  upper_staple=               <---- <--- ^
                              |          |
@@ -170,7 +170,7 @@ int main (int argc, char ** argv)
                             |          |
                             <-----<--- V 
 */          
-          tmp=CovShiftBackward(U[mu],nu,down_staple);
+          tmp=PeriodicBC::CovShiftBackward(U[mu],nu,down_staple);
           ds_U+=Cshift(tmp,mu,-1); 
 /*
      ds_U+=                 <----<---- ^ 
@@ -193,7 +193,7 @@ int main (int argc, char ** argv)
                      | 
                    (x)
 */    
-     LatticeColourMatrix up2= CovShiftForward(U[nu],nu,U[nu]);
+     LatticeColourMatrix up2= PeriodicBC::CovShiftForward(U[nu],nu,U[nu]);
 /*
                   <----^
                        | 
@@ -202,7 +202,7 @@ int main (int argc, char ** argv)
                  (x)
 */
     // Unu(x+mu)Unu(x+mu+nu) UmuDag(x+nu+nu) lives at X
-    upper_l= CovShiftForward(Cshift(up2,mu,1),nu,Cshift(adj(U[mu]),nu,1));
+    upper_l= PeriodicBC::CovShiftForward(Cshift(up2,mu,1),nu,Cshift(adj(U[mu]),nu,1));
 /*
                      |<----^
    upper_staple =    V     | 
@@ -218,7 +218,7 @@ int main (int argc, char ** argv)
    downer_l=           |  
                (x)<----V                 
 */    
-    upper_l= adj(CovShiftForward(U[mu],mu,up2)); //downer_l
+    upper_l= adj(PeriodicBC::CovShiftForward(U[mu],mu,up2)); //downer_l
 /*
                      ^     |
    down_staple  =    |     V 
diff --git a/tests/Test_gp_rect_force.cc b/tests/Test_gp_rect_force.cc
new file mode 100644
index 00000000..ade8a894
--- /dev/null
+++ b/tests/Test_gp_rect_force.cc
@@ -0,0 +1,98 @@
+#include <Grid.h>
+
+using namespace std;
+using namespace Grid;
+using namespace Grid::QCD;
+
+#define parallel_for PARALLEL_FOR_LOOP for
+
+int main (int argc, char ** argv)
+{
+  Grid_init(&argc,&argv);
+
+  std::vector<int> latt_size   = GridDefaultLatt();
+  std::vector<int> simd_layout = GridDefaultSimd(Nd,vComplex::Nsimd());
+  std::vector<int> mpi_layout  = GridDefaultMpi();
+
+  GridCartesian               Grid(latt_size,simd_layout,mpi_layout);
+  GridRedBlackCartesian     RBGrid(latt_size,simd_layout,mpi_layout);
+
+  int threads = GridThread::GetThreads();
+  std::cout<<GridLogMessage << "Grid is setup to use "<<threads<<" threads"<<std::endl;
+
+  std::vector<int> seeds({1,2,3,4});
+
+  GridParallelRNG          pRNG(&Grid);
+  pRNG.SeedRandomDevice();
+
+  LatticeGaugeField U(&Grid);
+
+  SU3::HotConfiguration(pRNG,U);
+  
+  double beta = 1.0;
+  double c1   = 0.331;
+
+  //GparityPlaqPlusRectangleActionR Action(beta,c1);
+  ConjugateWilsonGaugeActionR Action(beta);
+  //WilsonGaugeActionR Action(beta);
+
+  ComplexD S    = Action.S(U);
+
+  // get the deriv of phidag MdagM phi with respect to "U"
+  LatticeGaugeField UdSdU(&Grid);
+
+  Action.deriv(U,UdSdU);
+
+  ////////////////////////////////////
+  // Modify the gauge field a little 
+  ////////////////////////////////////
+  RealD dt = 0.0001;
+
+  LatticeColourMatrix mommu(&Grid); 
+  LatticeColourMatrix forcemu(&Grid); 
+  LatticeGaugeField mom(&Grid); 
+  LatticeGaugeField Uprime(&Grid); 
+
+  for(int mu=0;mu<Nd;mu++){
+
+    SU3::GaussianLieAlgebraMatrix(pRNG, mommu); // Traceless antihermitian momentum; gaussian in lie alg
+
+    PokeIndex<LorentzIndex>(mom,mommu,mu);
+
+    // fourth order exponential approx
+    parallel_for(auto i=mom.begin();i<mom.end();i++){ // exp(pmu dt) * Umu
+      Uprime[i](mu) = U[i](mu) + mom[i](mu)*U[i](mu)*dt ;
+    }
+  }
+
+  ComplexD Sprime    = Action.S(Uprime);
+
+  //////////////////////////////////////////////
+  // Use derivative to estimate dS
+  //////////////////////////////////////////////
+
+  LatticeComplex dS(&Grid); dS = zero;
+
+  for(int mu=0;mu<Nd;mu++){
+
+    auto UdSdUmu = PeekIndex<LorentzIndex>(UdSdU,mu);
+         mommu   = PeekIndex<LorentzIndex>(mom,mu);
+
+    // Update gauge action density
+    // U = exp(p dt) U
+    // dU/dt = p U
+    // so dSdt = trace( dUdt dSdU) = trace( p UdSdUmu ) 
+
+    dS = dS - trace(mommu*UdSdUmu)*dt*2.0;
+
+  }
+  Complex dSpred    = sum(dS);
+
+  std::cout << GridLogMessage << " S      "<<S<<std::endl;
+  std::cout << GridLogMessage << " Sprime "<<Sprime<<std::endl;
+  std::cout << GridLogMessage << "dS      "<<Sprime-S<<std::endl;
+  std::cout << GridLogMessage << "pred dS "<< dSpred <<std::endl;
+
+  std::cout<< GridLogMessage << "Done" <<std::endl;
+  Grid_finalize();
+}
diff --git a/tests/Test_hmc_EODWFRatio_Gparity.cc b/tests/Test_hmc_EODWFRatio_Gparity.cc
index 0ea2ca3f..521c496c 100644
--- a/tests/Test_hmc_EODWFRatio_Gparity.cc
+++ b/tests/Test_hmc_EODWFRatio_Gparity.cc
@@ -8,7 +8,7 @@ namespace Grid {
   namespace QCD { 
 
 
-class HmcRunner : public NerscHmcRunner {
+class HmcRunner : public ConjugateNerscHmcRunner {
 public:
 
   void BuildTheAction (int argc, char **argv)
@@ -31,7 +31,7 @@ public:
     LatticeGaugeField  U(UGrid);
 
     // Gauge action
-    WilsonGaugeActionR Waction(5.6);
+    ConjugateWilsonGaugeActionR Waction(5.6);
 
     Real mass=0.04;
     Real pv  =1.0;
diff --git a/tests/Test_hmc_GparityIwasakiGauge.cc b/tests/Test_hmc_GparityIwasakiGauge.cc
new file mode 100644
index 00000000..ffe572c4
--- /dev/null
+++ b/tests/Test_hmc_GparityIwasakiGauge.cc
@@ -0,0 +1,53 @@
+#include "Grid.h"
+
+using namespace std;
+using namespace Grid;
+using namespace Grid::QCD;
+
+namespace Grid { 
+  namespace QCD { 
+
+
+class HmcRunner : public ConjugateNerscHmcRunner {
+public:
+
+  void BuildTheAction (int argc, char **argv)
+
+  {
+    UGrid   = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd,vComplex::Nsimd()),GridDefaultMpi());
+    UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid);
+  
+    FGrid   = UGrid;
+    FrbGrid = UrbGrid;
+
+    // temporarily need a gauge field
+    LatticeGaugeField  U(UGrid);
+
+    // Gauge action
+    ConjugateIwasakiGaugeActionR Gaction(2.6);
+
+    //Collect actions
+    ActionLevel<LatticeGaugeField> Level1(1);
+    Level1.push_back(&Gaction);
+    TheAction.push_back(Level1);
+
+    Run(argc,argv);
+  };
+
+};
+
+}}
+
+int main (int argc, char ** argv)
+{
+  Grid_init(&argc,&argv);
+
+  int threads = GridThread::GetThreads();
+  std::cout<<GridLogMessage << "Grid is setup to use "<<threads<<" threads"<<std::endl;
+
+  HmcRunner TheHMC;
+  
+  TheHMC.BuildTheAction(argc,argv);
+
+}
+
diff --git a/tests/Test_hmc_GparityWilsonGauge.cc b/tests/Test_hmc_GparityWilsonGauge.cc
new file mode 100644
index 00000000..9dead670
--- /dev/null
+++ b/tests/Test_hmc_GparityWilsonGauge.cc
@@ -0,0 +1,53 @@
+#include "Grid.h"
+
+using namespace std;
+using namespace Grid;
+using namespace Grid::QCD;
+
+namespace Grid { 
+  namespace QCD { 
+
+
+class HmcRunner : public ConjugateNerscHmcRunner {
+public:
+
+  void BuildTheAction (int argc, char **argv)
+
+  {
+    UGrid   = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd,vComplex::Nsimd()),GridDefaultMpi());
+    UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid);
+  
+    FGrid   = UGrid;
+    FrbGrid = UrbGrid;
+
+    // temporarily need a gauge field
+    LatticeGaugeField  U(UGrid);
+
+    // Gauge action
+    ConjugateWilsonGaugeActionR Waction(5.6);
+
+    //Collect actions
+    ActionLevel<LatticeGaugeField> Level1(1);
+    Level1.push_back(&Waction);
+    TheAction.push_back(Level1);
+
+    Run(argc,argv);
+  };
+
+};
+
+}}
+
+int main (int argc, char ** argv)
+{
+  Grid_init(&argc,&argv);
+
+  int threads = GridThread::GetThreads();
+  std::cout<<GridLogMessage << "Grid is setup to use "<<threads<<" threads"<<std::endl;
+
+  HmcRunner TheHMC;
+  
+  TheHMC.BuildTheAction(argc,argv);
+
+}
+