Added unit tests on the representation transformations

Status: Passing all tests

lib/qcd/QCD.h

@@ -45,7 +45,7 @@ namespace QCD {
     static const int Zm = 6;
     static const int Tm = 7;
 
-    static const int Nc=3;
+    static const int Nc=2;
     static const int Ns=4;
     static const int Nd=4;
     static const int Nhs=2; // half spinor

lib/qcd/action/pseudofermion/TwoFlavour.h

@@ -130,8 +130,8 @@ class TwoFlavourPseudoFermionAction : public Action<typename Impl::GaugeField> {
     MdagMLinearOperator<FermionOperator<Impl>, FermionField> MdagMOp(FermOp);
 
     X = zero;
-    DerivativeSolver(MdagMOp, Phi, X);
-    MdagMOp.Op(X, Y);
+    DerivativeSolver(MdagMOp, Phi, X); // X = (MdagM)^-1 phi
+    MdagMOp.Op(X, Y);                  // Y = M X = (Mdag)^-1 phi
 
     // Check hermiticity
     

lib/qcd/hmc/HmcRunner.h

@@ -103,7 +103,7 @@ class NerscHmcRunnerTemplate {
 
     // Create integrator, including the smearing policy
     // Smearing policy, only defined for Nc=3
-    
+    /*
     std::cout << GridLogDebug << " Creating the Stout class\n";
     double rho = 0.1;  // smearing parameter, now hardcoded
     int Nsmear = 1;    // number of smearing levels
@@ -111,7 +111,7 @@ class NerscHmcRunnerTemplate {
     std::cout << GridLogDebug << " Creating the SmearedConfiguration class\n";
     //SmearedConfiguration<Gimpl> SmearingPolicy(UGrid, Nsmear, Stout);
     std::cout << GridLogDebug << " done\n";
-    
+    */
     //////////////
     NoSmearing<Gimpl> SmearingPolicy;
     typedef MinimumNorm2<GaugeField, NoSmearing<Gimpl>, RepresentationsPolicy >

lib/qcd/hmc/integrators/Integrator.h

@@ -120,6 +120,7 @@ class Integrator {
         FieldType forceR(U._grid);
         // Implement smearing only for the fundamental representation now
         repr_set.at(a)->deriv(Rep.U, forceR);
+        forceR -= adj(forceR);
         GF force =
             Rep.RtoFundamentalProject(forceR);  // Ta for the fundamental rep
         std::cout << GridLogIntegrator << "Hirep Force average: "
@@ -200,10 +201,12 @@ class Integrator {
     template <class FieldType, class Repr>
     void operator()(std::vector<Action<FieldType>*> repr_set, Repr& Rep,
                     GridParallelRNG& pRNG) {
-      for (int a = 0; a < repr_set.size(); ++a)
+      for (int a = 0; a < repr_set.size(); ++a){
         repr_set.at(a)->refresh(Rep.U, pRNG);
+      
       std::cout << GridLogDebug << "Hirep refreshing pseudofermions" << std::endl;
     }
+    }
   } refresh_hireps{};
 
   // Initialization of momenta and actions

lib/qcd/utils/SUn.h

@@ -672,6 +672,20 @@ class SU {
     }
   }
 
+  // Projects the algebra components a lattice matrix (of dimension ncol*ncol -1 )
+  // inverse operation: FundamentalLieAlgebraMatrix
+  static void projectOnAlgebra(LatticeAlgebraVector &h_out, const LatticeMatrix &in, Real scale = 1.0) {
+    conformable(h_out, in);
+    h_out = zero;
+    Matrix Ta;
+
+    for (int a = 0; a < AdjointDimension; a++) {
+      generator(a, Ta);
+      auto tmp = - 2.0 * (trace(timesI(Ta) * in)) * scale;// 2.0 for the normalization of the trace in the fundamental rep
+      pokeColour(h_out, tmp, a);
+    }
+    std::cout << "h_out " << h_out << std::endl;
+  }
 
 
   template <typename GaugeField>

lib/qcd/utils/SUnAdjoint.h

@@ -110,6 +110,22 @@ class SU_Adjoint : public SU<ncolour> {
     std::cout << GridLogMessage << std::endl;
   }
 
+  static void AdjointLieAlgebraMatrix(
+      const typename SU<ncolour>::LatticeAlgebraVector &h,
+      LatticeAdjMatrix &out, Real scale = 1.0) {
+    conformable(h, out);
+    GridBase *grid = out._grid;
+    LatticeAdjMatrix la(grid);
+    AMatrix iTa;
+
+    out = zero;
+    for (int a = 0; a < Dimension; a++) {
+      generator(a, iTa);
+      la = peekColour(h, a) * iTa * scale;
+      out += la;
+    }
+  }
+
   // Projects the algebra components a lattice matrix (of dimension ncol*ncol -1 )
   static void projectOnAlgebra(typename SU<ncolour>::LatticeAlgebraVector &h_out, const LatticeAdjMatrix &in, Real scale = 1.0) {
     conformable(h_out, in);
@@ -118,9 +134,10 @@ class SU_Adjoint : public SU<ncolour> {
 
     for (int a = 0; a < Dimension; a++) {
       generator(a, iTa);
-      auto tmp = - 2.0 * real(trace(iTa * in)) * scale; 
+      auto tmp = - 1.0/(ncolour) * (trace(iTa * in)) * scale;// 1/Nc for the normalization of the trace in the adj rep
       pokeColour(h_out, tmp, a);
     }
+    //std::cout << "h_out " << h_out << std::endl;
   }
 
   // a projector that keeps the generators stored to avoid the overhead of recomputing. 
@@ -130,12 +147,12 @@ class SU_Adjoint : public SU<ncolour> {
     h_out = zero;
     static bool precalculated = false; 
     if (!precalculated){
-    	precalculated = true;
+      precalculated = true;
         for (int a = 0; a < Dimension; a++) generator(a, iTa[a]);
     }
 
     for (int a = 0; a < Dimension; a++) {
-      auto tmp = - 2.0*real(trace(iTa[a] * in)) * scale; 
+      auto tmp = - 1.0/(ncolour) * (trace(iTa[a] * in)) * scale; 
       pokeColour(h_out, tmp, a);
     }
   }

tests/Test_hmc_WilsonAdjointFermionGauge.cc

@@ -65,12 +65,13 @@ class HmcRunner : public NerscHmcRunnerHirep< TheRepresentations > {
     AdjointRepresentation::LatticeField U(UGrid);
 
     // Gauge action
-    WilsonGaugeActionR Waction(5.6);
+    WilsonGaugeActionR Waction(2.0);
 
-    Real mass = -0.77;
+    Real mass = -1.0;
     FermionAction FermOp(U, *FGrid, *FrbGrid, mass);
 
     ConjugateGradient<FermionField> CG(1.0e-8, 10000);
+    ConjugateResidual<FermionField> CR(1.0e-8, 10000);
 
     TwoFlavourPseudoFermionAction<ImplPolicy> Nf2(FermOp, CG, CG);
 

tests/Test_lie_generators.cc

@@ -109,10 +109,104 @@ int main(int argc, char** argv) {
 
 
   // Testing HMC representation classes
-  AdjointRep<3> AdjRep(grid);
+  AdjointRep<Nc> AdjRep(grid);
 
   // AdjointRepresentation has the predefined number of colours Nc
   Representations<FundamentalRepresentation, AdjointRepresentation> RepresentationTypes(grid);  
 
+  LatticeGaugeField U(grid), V(grid);
+  SU<Nc>::HotConfiguration<LatticeGaugeField>(gridRNG, U);
+  SU<Nc>::HotConfiguration<LatticeGaugeField>(gridRNG, V);
+
+  // Test group structure
+  // (U_f * V_f)_r = U_r * V_r
+  LatticeGaugeField UV(grid);
+  for (int mu = 0; mu < Nd; mu++) {
+    SU<Nc>::LatticeMatrix Umu = peekLorentz(U,mu);
+    SU<Nc>::LatticeMatrix Vmu = peekLorentz(V,mu);
+    pokeLorentz(UV,Umu*Vmu, mu);
+  }
+
+  AdjRep.update_representation(UV);
+  typename AdjointRep<Nc>::LatticeField UVr = AdjRep.U;  // (U_f * V_f)_r
+
+
+  AdjRep.update_representation(U);
+  typename AdjointRep<Nc>::LatticeField Ur = AdjRep.U;  // U_r
+
+  AdjRep.update_representation(V);
+  typename AdjointRep<Nc>::LatticeField Vr = AdjRep.U;  // V_r
+
+  typename AdjointRep<Nc>::LatticeField UrVr(grid);
+  for (int mu = 0; mu < Nd; mu++) {
+    typename AdjointRep<Nc>::LatticeMatrix Urmu = peekLorentz(Ur,mu);
+    typename AdjointRep<Nc>::LatticeMatrix Vrmu = peekLorentz(Vr,mu);
+    pokeLorentz(UrVr,Urmu*Vrmu, mu);
+  }
+
+  typename AdjointRep<Nc>::LatticeField Diff_check = UVr - UrVr;
+  std::cout << GridLogMessage << "Group structure check difference : " << norm2(Diff_check) << std::endl;
+
+  //  Check correspondence of algebra and group transformations
+  // Create a random vector
+  SU<Nc>::LatticeAlgebraVector h_adj(grid);
+  typename AdjointRep<Nc>::LatticeMatrix Ar(grid);
+  random(gridRNG,h_adj);
+  h_adj = real(h_adj);
+  SU_Adjoint<Nc>::AdjointLieAlgebraMatrix(h_adj,Ar);
+
+  // Re-extract h_adj
+  SU<Nc>::LatticeAlgebraVector h_adj2(grid);
+  SU_Adjoint<Nc>::projectOnAlgebra(h_adj2, Ar);
+  SU<Nc>::LatticeAlgebraVector h_diff = h_adj - h_adj2;
+  std::cout << GridLogMessage << "Projections structure check vector difference : " << norm2(h_diff) << std::endl;
+
+  // Exponentiate
+  typename AdjointRep<Nc>::LatticeMatrix Uadj(grid);
+  Uadj  = expMat(Ar, 1.0, 16);
+
+  typename AdjointRep<Nc>::LatticeMatrix uno(grid);
+  uno = 1.0;
+  // Check matrix Uadj, must be real orthogonal
+  typename AdjointRep<Nc>::LatticeMatrix Ucheck = Uadj - conjugate(Uadj);
+  std::cout << GridLogMessage << "Reality check: " << norm2(Ucheck)
+            << std::endl;
+
+  Ucheck = Uadj * adj(Uadj) - uno;
+  std::cout << GridLogMessage << "orthogonality check 1: " << norm2(Ucheck)
+            << std::endl;
+  Ucheck = adj(Uadj) * Uadj - uno;
+  std::cout << GridLogMessage << "orthogonality check 2: " << norm2(Ucheck)
+            << std::endl;
+      
+
+  // Construct the fundamental matrix in the group
+  SU<Nc>::LatticeMatrix Af(grid);
+  SU<Nc>::FundamentalLieAlgebraMatrix(h_adj,Af);
+  SU<Nc>::LatticeMatrix Ufund(grid);
+  Ufund  = expMat(Af, 1.0, 16);
+  // Check unitarity
+  SU<Nc>::LatticeMatrix uno_f(grid);
+  uno_f = 1.0;
+  SU<Nc>::LatticeMatrix UnitCheck(grid);
+  UnitCheck = Ufund * adj(Ufund) - uno_f;
+  std::cout << GridLogMessage << "unitarity check 1: " << norm2(UnitCheck)
+            << std::endl;
+  UnitCheck = adj(Ufund) * Ufund - uno_f;
+  std::cout << GridLogMessage << "unitarity check 2: " << norm2(UnitCheck)
+            << std::endl;
+
+
+  // Tranform to the adjoint representation
+  U = zero; // fill this with only one direction
+  pokeLorentz(U,Ufund,0); // the representation transf acts on full gauge fields
+
+  AdjRep.update_representation(U);
+  Ur = AdjRep.U;  // U_r  
+  typename AdjointRep<Nc>::LatticeMatrix Ur0 = peekLorentz(Ur,0); // this should be the same as Uadj
+
+  typename AdjointRep<Nc>::LatticeMatrix Diff_check_mat = Ur0 - Uadj;
+  std::cout << GridLogMessage << "Projections structure check group difference : " << norm2(Diff_check_mat) << std::endl;
+
   Grid_finalize();
 }

tests/Test_wilson_force.cc

@@ -58,7 +58,7 @@ int main (int argc, char ** argv)
 
   LatticeGaugeField U(&Grid);
 
-  SU3::HotConfiguration(pRNG,U);
+  SU2::HotConfiguration(pRNG,U);
   //  SU3::ColdConfiguration(pRNG,U);
   
   ////////////////////////////////////
@@ -95,7 +95,7 @@ int main (int argc, char ** argv)
 
   for(int mu=0;mu<Nd;mu++){
 
-    SU3::GaussianFundamentalLieAlgebraMatrix(pRNG, mommu); // Traceless antihermitian momentum; gaussian in lie alg
+    SU2::GaussianFundamentalLieAlgebraMatrix(pRNG, mommu); // Traceless antihermitian momentum; gaussian in lie alg
 
     Hmom -= real(sum(trace(mommu*mommu)));