From 321f0f51b59109f9cb2b17d0e0f6a1883076be54 Mon Sep 17 00:00:00 2001
From: Peter Boyle <paboyle@ph.ed.ac.uk>
Date: Tue, 24 Nov 2020 21:46:10 -0500
Subject: [PATCH] Project to SU(N)

---
 Grid/qcd/action/gauge/GaugeImplTypes.h |  4 ++++
 Grid/qcd/action/scalar/ScalarImpl.h    |  8 ++++++++
 Grid/qcd/hmc/HMC.h                     |  2 +-
 Grid/qcd/hmc/integrators/Integrator.h  |  2 ++
 Grid/qcd/utils/SUn.h                   |  5 ++---
 tests/core/Test_reunitarise.cc         | 19 +++++++++++++------
 6 files changed, 30 insertions(+), 10 deletions(-)

diff --git a/Grid/qcd/action/gauge/GaugeImplTypes.h b/Grid/qcd/action/gauge/GaugeImplTypes.h
index 9b7d5a60..55a20eca 100644
--- a/Grid/qcd/action/gauge/GaugeImplTypes.h
+++ b/Grid/qcd/action/gauge/GaugeImplTypes.h
@@ -154,6 +154,10 @@ public:
     return Hsum.real();
   }
 
+  static inline void Project(Field &U) {
+    ProjectSUn(U);
+  }
+
   static inline void HotConfiguration(GridParallelRNG &pRNG, Field &U) {
     SU<Nc>::HotConfiguration(pRNG, U);
   }
diff --git a/Grid/qcd/action/scalar/ScalarImpl.h b/Grid/qcd/action/scalar/ScalarImpl.h
index 14675b11..403ea573 100644
--- a/Grid/qcd/action/scalar/ScalarImpl.h
+++ b/Grid/qcd/action/scalar/ScalarImpl.h
@@ -54,6 +54,10 @@ public:
   static inline void ColdConfiguration(GridParallelRNG &pRNG, Field &U) {
     U = 1.0;
   }
+
+  static inline void Project(Field &U) {
+    return;
+  }
     
   static void MomentumSpacePropagator(Field &out, RealD m)
   {
@@ -234,6 +238,10 @@ public:
 #endif //USE_FFT_ACCELERATION
   }
 
+  static inline void Project(Field &U) {
+    return;
+  }
+
   static inline void HotConfiguration(GridParallelRNG &pRNG, Field &U) {
     Group::GaussianFundamentalLieAlgebraMatrix(pRNG, U);
   }
diff --git a/Grid/qcd/hmc/HMC.h b/Grid/qcd/hmc/HMC.h
index 0f933204..f168b69a 100644
--- a/Grid/qcd/hmc/HMC.h
+++ b/Grid/qcd/hmc/HMC.h
@@ -95,7 +95,7 @@ private:
 
   typedef typename IntegratorType::Field Field;
   typedef std::vector< HmcObservable<Field> * > ObsListType;
-  
+
   //pass these from the resource manager
   GridSerialRNG &sRNG;   
   GridParallelRNG &pRNG; 
diff --git a/Grid/qcd/hmc/integrators/Integrator.h b/Grid/qcd/hmc/integrators/Integrator.h
index d5475704..70055754 100644
--- a/Grid/qcd/hmc/integrators/Integrator.h
+++ b/Grid/qcd/hmc/integrators/Integrator.h
@@ -313,6 +313,8 @@ public:
       std::cout << GridLogIntegrator << " times[" << level << "]= " << t_P[level] << " " << t_U << std::endl;
     }
 
+    FieldImplementation::Project(U);
+
     // and that we indeed got to the end of the trajectory
     assert(fabs(t_U - Params.trajL) < 1.0e-6);
 
diff --git a/Grid/qcd/utils/SUn.h b/Grid/qcd/utils/SUn.h
index 69ab4ebb..675493b3 100644
--- a/Grid/qcd/utils/SUn.h
+++ b/Grid/qcd/utils/SUn.h
@@ -820,7 +820,6 @@ LatticeComplexD Determinant(const Lattice<iScalar<iScalar<iMatrix<vComplexD, N>
       }}
     ComplexD det = EigenU.determinant();
     pokeLocalSite(det,ret_v,lcoor);
-    std::cout << " site " <<site<<" det " <<det <<std::endl;
   });
   return ret;
 }
@@ -830,8 +829,8 @@ static void ProjectSUn(Lattice<iScalar<iScalar<iMatrix<vComplexD, N> > > > &Umu)
   Umu      = ProjectOnGroup(Umu);
   auto det = Determinant(Umu);
 
-  det = pow(det,-1);
-  
+  det = conjugate(det);
+
   for(int i=0;i<N;i++){
     auto element = PeekIndex<ColourIndex>(Umu,N-1,i);
     element = element * det;
diff --git a/tests/core/Test_reunitarise.cc b/tests/core/Test_reunitarise.cc
index 3e78b961..9a6781f1 100644
--- a/tests/core/Test_reunitarise.cc
+++ b/tests/core/Test_reunitarise.cc
@@ -102,7 +102,8 @@ int main (int argc, char ** argv)
   LatticeComplexD detUU(grid);
 
   detU= Determinant(U) ;
-  std::cout << "Determinant before screw up " <<detU<<std::endl;
+  detU=detU-1.0;
+  std::cout << "Determinant before screw up " << norm2(detU)<<std::endl;
 
   std::cout << " Screwing up determinant " << std::endl;
 
@@ -116,18 +117,24 @@ int main (int argc, char ** argv)
   UU=U;
 
   detU= Determinant(U) ;
-  std::cout << "Determinant after screw up " <<detU<<std::endl;
-
+  detU=detU-1.0;
+  std::cout << "Determinant defect before projection " <<norm2(detU)<<std::endl;
+  tmp = U*adj(U) - ident;
+  std::cout << "Unitarity check before projection    " << norm2(tmp)<<std::endl; 
+  
   ProjectSU3(U);
   detU= Determinant(U) ;
-  std::cout << "Determinant ProjectSU3  " <<detU<<std::endl;
   detU= detU -1.0;
   std::cout << "Determinant ProjectSU3 defect " <<norm2(detU)<<std::endl;
+  tmp = U*adj(U) - ident;
+  std::cout << "Unitarity check after projection    " << norm2(tmp)<<std::endl; 
 
-  ProjectSUn<3>(UU);
+  ProjectSUn(UU);
   detUU= Determinant(UU);
-  std::cout << "Determinant ProjectSUn " <<detUU<<std::endl;
+  detUU= detUU -1.0;
   std::cout << "Determinant ProjectSUn defect " <<norm2(detUU)<<std::endl;
+  tmp = UU*adj(UU) - ident;
+  std::cout << "Unitarity check after projection    " << norm2(tmp)<<std::endl; 
   
   Grid_finalize();
 }