Merge branch 'feature/conjugate-bc-dirs' into develop

Grid/parallelIO/IldgIO.h

@@ -123,7 +123,7 @@ assert(GRID_FIELD_NORM_CALC(FieldNormMetaData_, n2ck) < 1.0e-5);
  ////////////////////////////////////////////////////////////
  // Helper to fill out metadata
  ////////////////////////////////////////////////////////////
- template<class vobj> void ScidacMetaData(Lattice<vobj> & field,
+template<class vobj> void ScidacMetaData(Lattice<vobj> & field,
 					  FieldMetaData &header,
 					  scidacRecord & _scidacRecord,
 					  scidacFile   & _scidacFile) 
@@ -619,12 +619,12 @@ class IldgWriter : public ScidacWriter {
   // Don't require scidac records EXCEPT checksum
   // Use Grid MetaData object if present.
   ////////////////////////////////////////////////////////////////
-  template <class vsimd>
+  template <class stats = PeriodicGaugeStatistics>
-  void writeConfiguration(Lattice<iLorentzColourMatrix<vsimd> > &Umu,int sequence,std::string LFN,std::string description) 
+  void writeConfiguration(Lattice<vLorentzColourMatrixD > &Umu,int sequence,std::string LFN,std::string description) 
   {
     GridBase * grid = Umu.Grid();
-    typedef Lattice<iLorentzColourMatrix<vsimd> > GaugeField;
+    typedef Lattice<vLorentzColourMatrixD> GaugeField;
-    typedef iLorentzColourMatrix<vsimd> vobj;
+    typedef vLorentzColourMatrixD vobj;
     typedef typename vobj::scalar_object sobj;
 
     ////////////////////////////////////////
@@ -636,6 +636,9 @@ class IldgWriter : public ScidacWriter {
 
     ScidacMetaData(Umu,header,_scidacRecord,_scidacFile);
 
+    stats Stats;
+    Stats(Umu,header);
+    
     std::string format = header.floating_point;
     header.ensemble_id    = description;
     header.ensemble_label = description;
@@ -705,10 +708,10 @@ class IldgReader : public GridLimeReader {
   // Else use ILDG MetaData object if present.
   // Else use SciDAC MetaData object if present.
   ////////////////////////////////////////////////////////////////
-  template <class vsimd>
+  template <class stats = PeriodicGaugeStatistics>
-  void readConfiguration(Lattice<iLorentzColourMatrix<vsimd> > &Umu, FieldMetaData &FieldMetaData_) {
+  void readConfiguration(Lattice<vLorentzColourMatrixD> &Umu, FieldMetaData &FieldMetaData_) {
 
-    typedef Lattice<iLorentzColourMatrix<vsimd> > GaugeField;
+    typedef Lattice<vLorentzColourMatrixD > GaugeField;
     typedef typename GaugeField::vector_object  vobj;
     typedef typename vobj::scalar_object sobj;
 
@@ -921,7 +924,8 @@ class IldgReader : public GridLimeReader {
 
     if ( found_FieldMetaData || found_usqcdInfo ) {
       FieldMetaData checker;
-      GaugeStatistics(Umu,checker);
+      stats Stats;
+      Stats(Umu,checker);
       assert(fabs(checker.plaquette  - FieldMetaData_.plaquette )<1.0e-5);
       assert(fabs(checker.link_trace - FieldMetaData_.link_trace)<1.0e-5);
       std::cout << GridLogMessage<<"Plaquette and link trace match " << std::endl;

Grid/parallelIO/MetaData.h

@@ -176,29 +176,18 @@ template<class vobj> inline void PrepareMetaData(Lattice<vobj> & field, FieldMet
   GridMetaData(grid,header); 
   MachineCharacteristics(header);
 }
-inline void GaugeStatistics(Lattice<vLorentzColourMatrixF> & data,FieldMetaData &header)
+template<class Impl>
+class GaugeStatistics
 {
-  // How to convert data precision etc...
+public:
-  header.link_trace=WilsonLoops<PeriodicGimplF>::linkTrace(data);
+  void operator()(Lattice<vLorentzColourMatrixD> & data,FieldMetaData &header)
-  header.plaquette =WilsonLoops<PeriodicGimplF>::avgPlaquette(data);
+  {
-}
+    header.link_trace=WilsonLoops<Impl>::linkTrace(data);
-inline void GaugeStatistics(Lattice<vLorentzColourMatrixD> & data,FieldMetaData &header)
+    header.plaquette =WilsonLoops<Impl>::avgPlaquette(data);
-{
+  }
-  // How to convert data precision etc...
+};
-  header.link_trace=WilsonLoops<PeriodicGimplD>::linkTrace(data);
+typedef GaugeStatistics<PeriodicGimplD> PeriodicGaugeStatistics;
-  header.plaquette =WilsonLoops<PeriodicGimplD>::avgPlaquette(data);
+typedef GaugeStatistics<ConjugateGimplD> ConjugateGaugeStatistics;
-}
-template<> inline void PrepareMetaData<vLorentzColourMatrixF>(Lattice<vLorentzColourMatrixF> & field, FieldMetaData &header)
-{
-   
-  GridBase *grid = field.Grid();
-  std::string format = getFormatString<vLorentzColourMatrixF>();
-  header.floating_point = format;
-  header.checksum = 0x0; // Nersc checksum unused in ILDG, Scidac
-  GridMetaData(grid,header); 
-  GaugeStatistics(field,header);
-  MachineCharacteristics(header);
-}
 template<> inline void PrepareMetaData<vLorentzColourMatrixD>(Lattice<vLorentzColourMatrixD> & field, FieldMetaData &header)
 {
   GridBase *grid = field.Grid();
@@ -206,7 +195,6 @@ template<> inline void PrepareMetaData<vLorentzColourMatrixD>(Lattice<vLorentzCo
   header.floating_point = format;
   header.checksum = 0x0; // Nersc checksum unused in ILDG, Scidac
   GridMetaData(grid,header); 
-  GaugeStatistics(field,header);
   MachineCharacteristics(header);
 }
 

Grid/parallelIO/NerscIO.h

@@ -40,6 +40,8 @@ using namespace Grid;
 class NerscIO : public BinaryIO { 
 public:
 
+  typedef Lattice<vLorentzColourMatrixD> GaugeField;
+
   static inline void truncate(std::string file){
     std::ofstream fout(file,std::ios::out);
   }
@@ -129,12 +131,12 @@ public:
   // Now the meat: the object readers
   /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
 
-  template<class vsimd>
+  template<class GaugeStats=PeriodicGaugeStatistics>
-  static inline void readConfiguration(Lattice<iLorentzColourMatrix<vsimd> > &Umu,
+  static inline void readConfiguration(GaugeField &Umu,
 				       FieldMetaData& header,
-				       std::string file)
+				       std::string file,
+				       GaugeStats GaugeStatisticsCalculator=GaugeStats())
   {
-    typedef Lattice<iLorentzColourMatrix<vsimd> > GaugeField;
 
     GridBase *grid = Umu.Grid();
     uint64_t offset = readHeader(file,Umu.Grid(),header);
@@ -153,23 +155,23 @@ public:
     // munger is a function of <floating point, Real, data_type>
     if ( header.data_type == std::string("4D_SU3_GAUGE") ) {
       if ( ieee32 || ieee32big ) {
-	BinaryIO::readLatticeObject<iLorentzColourMatrix<vsimd>, LorentzColour2x3F> 
+	BinaryIO::readLatticeObject<vLorentzColourMatrixD, LorentzColour2x3F> 
 	  (Umu,file,Gauge3x2munger<LorentzColour2x3F,LorentzColourMatrix>(), offset,format,
 	   nersc_csum,scidac_csuma,scidac_csumb);
       }
       if ( ieee64 || ieee64big ) {
-	BinaryIO::readLatticeObject<iLorentzColourMatrix<vsimd>, LorentzColour2x3D> 
+	BinaryIO::readLatticeObject<vLorentzColourMatrixD, LorentzColour2x3D> 
 	  (Umu,file,Gauge3x2munger<LorentzColour2x3D,LorentzColourMatrix>(),offset,format,
 	   nersc_csum,scidac_csuma,scidac_csumb);
       }
     } else if ( header.data_type == std::string("4D_SU3_GAUGE_3x3") ) {
       if ( ieee32 || ieee32big ) {
-	BinaryIO::readLatticeObject<iLorentzColourMatrix<vsimd>,LorentzColourMatrixF>
+	BinaryIO::readLatticeObject<vLorentzColourMatrixD,LorentzColourMatrixF>
 	  (Umu,file,GaugeSimpleMunger<LorentzColourMatrixF,LorentzColourMatrix>(),offset,format,
 	   nersc_csum,scidac_csuma,scidac_csumb);
       }
       if ( ieee64 || ieee64big ) {
-	BinaryIO::readLatticeObject<iLorentzColourMatrix<vsimd>,LorentzColourMatrixD>
+	BinaryIO::readLatticeObject<vLorentzColourMatrixD,LorentzColourMatrixD>
 	  (Umu,file,GaugeSimpleMunger<LorentzColourMatrixD,LorentzColourMatrix>(),offset,format,
 	   nersc_csum,scidac_csuma,scidac_csumb);
       }
@@ -177,7 +179,7 @@ public:
       assert(0);
     }
 
-    GaugeStatistics(Umu,clone);
+    GaugeStats Stats; Stats(Umu,clone);
 
     std::cout<<GridLogMessage <<"NERSC Configuration "<<file<<" checksum "<<std::hex<<nersc_csum<< std::dec
 	     <<" header   "<<std::hex<<header.checksum<<std::dec <<std::endl;
@@ -203,15 +205,13 @@ public:
     std::cout<<GridLogMessage <<"NERSC Configuration "<<file<< " and plaquette, link trace, and checksum agree"<<std::endl;
   }
 
-  template<class vsimd>
+  template<class GaugeStats=PeriodicGaugeStatistics>
-  static inline void writeConfiguration(Lattice<iLorentzColourMatrix<vsimd> > &Umu,
+  static inline void writeConfiguration(Lattice<vLorentzColourMatrixD > &Umu,
 					std::string file, 
 					int two_row,
 					int bits32)
   {
-    typedef Lattice<iLorentzColourMatrix<vsimd> > GaugeField;
+    typedef vLorentzColourMatrixD vobj;
-
-    typedef iLorentzColourMatrix<vsimd> vobj;
     typedef typename vobj::scalar_object sobj;
 
     FieldMetaData header;
@@ -229,7 +229,7 @@ public:
 
     GridMetaData(grid,header);
     assert(header.nd==4);
-    GaugeStatistics(Umu,header);
+    GaugeStats Stats; Stats(Umu,header);
     MachineCharacteristics(header);
 
 	uint64_t offset;
@@ -238,19 +238,19 @@ public:
     header.floating_point = std::string("IEEE64BIG");
     header.data_type      = std::string("4D_SU3_GAUGE_3x3");
     GaugeSimpleUnmunger<fobj3D,sobj> munge;
-	if ( grid->IsBoss() ) { 
+    if ( grid->IsBoss() ) { 
-	  truncate(file);
+      truncate(file);
-    offset = writeHeader(header,file);
+      offset = writeHeader(header,file);
-	}
+    }
-	grid->Broadcast(0,(void *)&offset,sizeof(offset));
+    grid->Broadcast(0,(void *)&offset,sizeof(offset));
 
     uint32_t nersc_csum,scidac_csuma,scidac_csumb;
     BinaryIO::writeLatticeObject<vobj,fobj3D>(Umu,file,munge,offset,header.floating_point,
 					      nersc_csum,scidac_csuma,scidac_csumb);
     header.checksum = nersc_csum;
-	if ( grid->IsBoss() ) { 
+    if ( grid->IsBoss() ) { 
-    writeHeader(header,file);
+      writeHeader(header,file);
-	}
+    }
 
     std::cout<<GridLogMessage <<"Written NERSC Configuration on "<< file << " checksum "
 	     <<std::hex<<header.checksum

Grid/parallelIO/OpenQcdIO.h

@@ -154,7 +154,7 @@ public:
     grid->Barrier(); timer.Stop();
     std::cout << Grid::GridLogMessage << "OpenQcdIO::readConfiguration: redistribute overhead " << timer.Elapsed() << std::endl;
 
-    GaugeStatistics(Umu, clone);
+    PeriodicGaugeStatistics Stats; Stats(Umu, clone);
 
     RealD plaq_diff = fabs(clone.plaquette - header.plaquette);
 

Grid/parallelIO/OpenQcdIOChromaReference.h

@@ -208,7 +208,7 @@ public:
 
     FieldMetaData clone(header);
 
-    GaugeStatistics(Umu, clone);
+    PeriodicGaugeStatistics Stats; Stats(Umu, clone);
 
     RealD plaq_diff = fabs(clone.plaquette - header.plaquette);
 

Grid/qcd/action/gauge/Gauge.cc

@@ -0,0 +1,38 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid
+
+Source file: ./lib/qcd/action/gauge/Gauge.cc
+
+Copyright (C) 2020
+
+Author: Peter Boyle <paboyle@ph.ed.ac.uk>
+Author: Peter Boyle <peterboyle@Peters-MacBook-Pro-2.local>
+Author: paboyle <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/qcd/action/fermion/FermionCore.h>
+
+NAMESPACE_BEGIN(Grid);
+
+std::vector<int> ConjugateGaugeImplBase::_conjDirs;
+
+NAMESPACE_END(Grid);
+

Grid/qcd/action/gauge/GaugeImplementations.h

@@ -59,14 +59,14 @@ public:
   }
   static inline GaugeLinkField
   CovShiftIdentityBackward(const GaugeLinkField &Link, int mu) {
-    return Cshift(adj(Link), mu, -1);
+    return PeriodicBC::CovShiftIdentityBackward(Link, mu);
   }
   static inline GaugeLinkField
   CovShiftIdentityForward(const GaugeLinkField &Link, int mu) {
-    return Link;
+    return PeriodicBC::CovShiftIdentityForward(Link,mu);
   }
   static inline GaugeLinkField ShiftStaple(const GaugeLinkField &Link, int mu) {
-    return Cshift(Link, mu, 1);
+    return PeriodicBC::ShiftStaple(Link,mu);
   }
 
   static inline bool isPeriodicGaugeField(void) { return true; }
@@ -74,7 +74,13 @@ public:
 
 // Composition with smeared link, bc's etc.. probably need multiple inheritance
 // Variable precision "S" and variable Nc
-template <class GimplTypes> class ConjugateGaugeImpl : public GimplTypes {
+class ConjugateGaugeImplBase {
+protected:
+  static std::vector<int> _conjDirs;
+};
+
+  template <class GimplTypes> class ConjugateGaugeImpl : public GimplTypes, ConjugateGaugeImplBase {
+private:
 public:
   INHERIT_GIMPL_TYPES(GimplTypes);
 
@@ -84,47 +90,56 @@ public:
   ////////////////////////////////////////////////////////////////////////////////////////////////////////////
   template <class covariant>
   static Lattice<covariant> CovShiftForward(const GaugeLinkField &Link, int mu,
-                                            const Lattice<covariant> &field) {
+                                            const Lattice<covariant> &field)
-    return ConjugateBC::CovShiftForward(Link, mu, field);
+  {
+    assert(_conjDirs.size() == Nd);
+    if(_conjDirs[mu]) 
+      return ConjugateBC::CovShiftForward(Link, mu, field);
+    else
+      return PeriodicBC::CovShiftForward(Link, mu, field);
   }
 
   template <class covariant>
   static Lattice<covariant> CovShiftBackward(const GaugeLinkField &Link, int mu,
-                                             const Lattice<covariant> &field) {
+                                             const Lattice<covariant> &field)
-    return ConjugateBC::CovShiftBackward(Link, mu, field);
+  {
+    assert(_conjDirs.size() == Nd);
+    if(_conjDirs[mu]) 
+      return ConjugateBC::CovShiftBackward(Link, mu, field);
+    else 
+      return PeriodicBC::CovShiftBackward(Link, mu, field);
   }
 
   static inline GaugeLinkField
-  CovShiftIdentityBackward(const GaugeLinkField &Link, int mu) {
+  CovShiftIdentityBackward(const GaugeLinkField &Link, int mu)
-    GridBase *grid = Link.Grid();
+  {
-    int Lmu = grid->GlobalDimensions()[mu] - 1;
+    assert(_conjDirs.size() == Nd);
-
+    if(_conjDirs[mu]) 
-    Lattice<iScalar<vInteger>> coor(grid);
+      return ConjugateBC::CovShiftIdentityBackward(Link, mu);
-    LatticeCoordinate(coor, mu);
+    else 
-
+      return PeriodicBC::CovShiftIdentityBackward(Link, 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) {
+  CovShiftIdentityForward(const GaugeLinkField &Link, int mu)
-    return Link;
+  {
+    assert(_conjDirs.size() == Nd);
+    if(_conjDirs[mu]) 
+      return ConjugateBC::CovShiftIdentityForward(Link,mu);
+    else
+      return PeriodicBC::CovShiftIdentityForward(Link,mu);
   }
 
-  static inline GaugeLinkField ShiftStaple(const GaugeLinkField &Link, int mu) {
+  static inline GaugeLinkField ShiftStaple(const GaugeLinkField &Link, int mu)
-    GridBase *grid = Link.Grid();
+  {
-    int Lmu = grid->GlobalDimensions()[mu] - 1;
+    assert(_conjDirs.size() == Nd);
-
+    if(_conjDirs[mu]) 
-    Lattice<iScalar<vInteger>> coor(grid);
+      return ConjugateBC::ShiftStaple(Link,mu);
-    LatticeCoordinate(coor, mu);
+    else     
-
+      return PeriodicBC::ShiftStaple(Link,mu);
-    GaugeLinkField tmp(grid);
-    tmp = Cshift(Link, mu, 1);
-    tmp = where(coor == Lmu, conjugate(tmp), tmp);
-    return tmp;
   }
 
+  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; }
 };
 

Grid/qcd/hmc/checkpointers/BaseCheckpointer.h

@@ -74,7 +74,7 @@ public:
       conf_file = os.str();
     }
   } 
-
+  virtual ~BaseHmcCheckpointer(){};
   void check_filename(const std::string &filename){
     std::ifstream f(filename.c_str());
     if(!f.good()){
@@ -82,7 +82,6 @@ public:
       abort();
     };
   }
-
   virtual void initialize(const CheckpointerParameters &Params) = 0;
 
   virtual void CheckpointRestore(int traj, typename Impl::Field &U,

Grid/qcd/hmc/checkpointers/ILDGCheckpointer.h

@@ -45,6 +45,7 @@ private:
 
 public:
   INHERIT_GIMPL_TYPES(Implementation);
+  typedef GaugeStatistics<Implementation> GaugeStats;
 
   ILDGHmcCheckpointer(const CheckpointerParameters &Params_) { initialize(Params_); }
 
@@ -78,7 +79,7 @@ public:
       BinaryIO::writeRNG(sRNG, pRNG, rng, 0,nersc_csum,scidac_csuma,scidac_csumb);
       IldgWriter _IldgWriter(grid->IsBoss());
       _IldgWriter.open(config);
-      _IldgWriter.writeConfiguration(U, traj, config, config);
+      _IldgWriter.writeConfiguration<GaugeStats>(U, traj, config, config);
       _IldgWriter.close();
 
       std::cout << GridLogMessage << "Written ILDG Configuration on " << config
@@ -105,7 +106,7 @@ public:
     FieldMetaData header;
     IldgReader _IldgReader;
     _IldgReader.open(config);
-    _IldgReader.readConfiguration(U,header);  // format from the header
+    _IldgReader.readConfiguration<GaugeStats>(U,header);  // format from the header
     _IldgReader.close();
 
     std::cout << GridLogMessage << "Read ILDG Configuration from " << config

Grid/qcd/hmc/checkpointers/NerscCheckpointer.h

@@ -43,7 +43,8 @@ private:
 
 public:
   INHERIT_GIMPL_TYPES(Gimpl);  // only for gauge configurations
-
+  typedef GaugeStatistics<Gimpl> GaugeStats;
+  
   NerscHmcCheckpointer(const CheckpointerParameters &Params_) { initialize(Params_); }
 
   void initialize(const CheckpointerParameters &Params_) {
@@ -60,7 +61,7 @@ public:
       int precision32 = 1;
       int tworow = 0;
       NerscIO::writeRNGState(sRNG, pRNG, rng);
-      NerscIO::writeConfiguration(U, config, tworow, precision32);
+      NerscIO::writeConfiguration<GaugeStats>(U, config, tworow, precision32);
     }
   };
 
@@ -74,7 +75,7 @@ public:
 
     FieldMetaData header;
     NerscIO::readRNGState(sRNG, pRNG, header, rng);
-    NerscIO::readConfiguration(U, header, config);
+    NerscIO::readConfiguration<GaugeStats>(U, header, config);
   };
 };
 

Grid/qcd/modules/Modules.h

@@ -99,7 +99,7 @@ public:
   virtual Prod* getPtr() = 0;
 
   // add a getReference? 
-  
+  virtual ~HMCModuleBase(){};
   virtual void print_parameters(){};  // default to nothing
 };
 

Grid/qcd/utils/CovariantCshift.h

@@ -53,6 +53,24 @@ namespace PeriodicBC {
     return Cshift(tmp,mu,-1);// moves towards positive mu
   }
 
+  template<class gauge> Lattice<gauge>
+  CovShiftIdentityBackward(const Lattice<gauge> &Link, int mu) 
+  {
+    return Cshift(adj(Link), mu, -1);
+  }
+
+  template<class gauge> Lattice<gauge>
+  CovShiftIdentityForward(const Lattice<gauge> &Link, int mu)
+  {
+    return Link;
+  }
+
+  template<class gauge> Lattice<gauge>
+  ShiftStaple(const Lattice<gauge> &Link, int mu)
+  {
+    return Cshift(Link, mu, 1);
+  }
+  
   template<class gauge,class Expr,typename std::enable_if<is_lattice_expr<Expr>::value,void>::type * = nullptr>
     auto  CovShiftForward(const Lattice<gauge> &Link, 
 			  int mu,
@@ -70,6 +88,7 @@ namespace PeriodicBC {
     return CovShiftBackward(Link,mu,arg);
   }
 
+
 }
 
 
@@ -139,6 +158,38 @@ namespace ConjugateBC {
     //    std::cout<<"Gparity::CovCshiftBackward mu="<<mu<<std::endl;
     return Cshift(tmp,mu,-1);// moves towards positive mu
   }
+  template<class gauge> Lattice<gauge>
+  CovShiftIdentityBackward(const Lattice<gauge> &Link, int mu) {
+    GridBase *grid = Link.Grid();
+    int Lmu = grid->GlobalDimensions()[mu] - 1;
+
+    Lattice<iScalar<vInteger>> coor(grid);
+    LatticeCoordinate(coor, mu);
+
+    Lattice<gauge> tmp(grid);
+    tmp = adj(Link);
+    tmp = where(coor == Lmu, conjugate(tmp), tmp);
+    return Cshift(tmp, mu, -1); // moves towards positive mu
+  }
+  template<class gauge> Lattice<gauge>
+  CovShiftIdentityForward(const Lattice<gauge> &Link, int mu) {
+    return Link;
+  }
+
+  template<class gauge> Lattice<gauge>
+  ShiftStaple(const Lattice<gauge> &Link, int mu)
+  {
+    GridBase *grid = Link.Grid();
+    int Lmu = grid->GlobalDimensions()[mu] - 1;
+
+    Lattice<iScalar<vInteger>> coor(grid);
+    LatticeCoordinate(coor, mu);
+
+    Lattice<gauge> tmp(grid);
+    tmp = Cshift(Link, mu, 1);
+    tmp = where(coor == Lmu, conjugate(tmp), tmp);
+    return tmp;
+  }
 
   template<class gauge,class Expr,typename std::enable_if<is_lattice_expr<Expr>::value,void>::type * = nullptr>
     auto  CovShiftForward(const Lattice<gauge> &Link, 

Grid/tensors/Tensor_Ta.h

@@ -117,7 +117,19 @@ accelerator_inline iMatrix<vtype,N> ProjectOnGroup(const iMatrix<vtype,N> &arg)
 	ret._internal[b][c] -= pr * ret._internal[c1][c];
       }
     }
-	  
+  }
+
+  // Normalise last row
+  {
+    int c1 = N-1;
+    zeroit(inner);	
+    for(int c2=0;c2<N;c2++)
+      inner += innerProduct(ret._internal[c1][c2],ret._internal[c1][c2]);
+
+    nrm = sqrt(inner);
+    nrm = 1.0/nrm;
+    for(int c2=0;c2<N;c2++)
+      ret._internal[c1][c2]*= nrm;
   }
   // assuming the determinant is ok
   return ret;

tests/core/Test_reunitarise.cc

@@ -103,7 +103,7 @@ int main (int argc, char ** argv)
 
   detU= Determinant(U) ;
   detU=detU-1.0;
-  std::cout << "Determinant before screw up " << norm2(detU)<<std::endl;
+  std::cout << "Determinant defect before screw up " << norm2(detU)<<std::endl;
 
   std::cout << " Screwing up determinant " << std::endl;
 
@@ -113,7 +113,8 @@ int main (int argc, char ** argv)
     auto element = PeekIndex<ColourIndex>(U,Nc-1,i);
     element = element * phase;
     PokeIndex<ColourIndex>(U,element,Nc-1,i);
-  }  
+  }
+  U=U*0.1;
   UU=U;
 
   detU= Determinant(U) ;

tests/hmc/Test_hmc_EODWFRatio_Gparity.cc

@@ -81,6 +81,10 @@ int main(int argc, char **argv) {
   // that have a complex construction
   // standard
   RealD beta = 5.6 ;
+  const int nu = 3;
+  std::vector<int> twists(Nd,0);
+  twists[nu] = 1;
+  ConjugateGimplD::setDirections(twists);
   ConjugateWilsonGaugeActionR Waction(beta);
 
   const int Ls = 8;
@@ -93,9 +97,6 @@ int main(int argc, char **argv) {
   // temporarily need a gauge field
   LatticeGaugeField U(GridPtr);
 
-  const int nu = 3;
-  std::vector<int> twists(Nd,0);
-  twists[nu] = 1;
   FermionAction::ImplParams params;
   params.twists = twists;
   Real mass=0.04;

tests/hmc/Test_hmc_GparityIwasakiGauge.cc

@@ -79,6 +79,10 @@ int main(int argc, char **argv) {
   // that have a complex construction
   // standard
   RealD beta = 2.6 ;
+  const int nu = 3;
+  std::vector<int> twists(Nd,0);
+  twists[nu] = 1;
+  ConjugateGimplD::setDirections(twists);
   ConjugateIwasakiGaugeActionR Waction(beta);
 
  

tests/hmc/Test_hmc_GparityWilsonGauge.cc

@@ -80,6 +80,9 @@ int main(int argc, char **argv) {
   // that have a complex construction
   // standard
   RealD beta = 5.6 ;
+  std::vector<int> twists(Nd,0);
+  twists[3] = 1;
+  ConjugateGimplD::setDirections(twists);
   ConjugateWilsonGaugeActionR Waction(beta);