Debugged HMC for Creutz relation

lib/qcd/action/ActionBase.h

@@ -150,11 +150,11 @@ struct ActionLevelHirep {
   // Loop on tuple for a callable function
   template <std::size_t I = 1, typename Callable, typename ...Args>
   inline typename std::enable_if<I == std::tuple_size<action_collection>::value, void>::type apply(
-      Callable, Repr& R,Args...) const {}
+      Callable, Repr& R,Args&...) const {}
 
   template <std::size_t I = 1, typename Callable, typename ...Args>
   inline typename std::enable_if<I < std::tuple_size<action_collection>::value, void>::type apply(
-      Callable fn, Repr& R, Args... arguments) const {
+      Callable fn, Repr& R, Args&... arguments) const {
     fn(std::get<I>(actions_hirep), std::get<I>(R.rep), arguments...);
     apply<I + 1>(fn, R, arguments...);
   }  

lib/qcd/action/fermion/WilsonKernels.h

@@ -54,7 +54,7 @@ class WilsonKernels : public FermionOperator<Impl>, public WilsonKernelsStatic {
 
  public:
   template <bool EnableBool = true>
-  typename std::enable_if<Impl::Dimension == 3 && EnableBool, void>::type
+  typename std::enable_if<Impl::Dimension == 3 && Nc == 3 &&EnableBool, void>::type
   DiracOptDhopSite(
       StencilImpl &st, LebesgueOrder &lo, DoubledGaugeField &U,
       std::vector<SiteHalfSpinor, alignedAllocator<SiteHalfSpinor> > &buf,
@@ -85,7 +85,7 @@ class WilsonKernels : public FermionOperator<Impl>, public WilsonKernelsStatic {
   }
 
   template <bool EnableBool = true>
-  typename std::enable_if<Impl::Dimension != 3 && EnableBool, void>::type
+    typename std::enable_if<(Impl::Dimension != 3 || (Impl::Dimension == 3 && Nc != 3)) && EnableBool, void>::type
   DiracOptDhopSite(
       StencilImpl &st, LebesgueOrder &lo, DoubledGaugeField &U,
       std::vector<SiteHalfSpinor, alignedAllocator<SiteHalfSpinor> > &buf,
@@ -102,7 +102,7 @@ class WilsonKernels : public FermionOperator<Impl>, public WilsonKernelsStatic {
   }
 
   template <bool EnableBool = true>
-  typename std::enable_if<Impl::Dimension == 3 && EnableBool, void>::type
+  typename std::enable_if<Impl::Dimension == 3 && Nc== 3 && EnableBool, void>::type
   DiracOptDhopSiteDag(
       StencilImpl &st, LebesgueOrder &lo, DoubledGaugeField &U,
       std::vector<SiteHalfSpinor, alignedAllocator<SiteHalfSpinor> > &buf,
@@ -126,7 +126,7 @@ class WilsonKernels : public FermionOperator<Impl>, public WilsonKernelsStatic {
   }
 
   template <bool EnableBool = true>
-  typename std::enable_if<Impl::Dimension != 3 && EnableBool, void>::type
+    typename std::enable_if<(Impl::Dimension != 3 || (Impl::Dimension == 3 && Nc != 3)) && EnableBool, void>::type
     DiracOptDhopSiteDag(
 			StencilImpl &st, LebesgueOrder &lo, DoubledGaugeField &U,
 			std::vector<SiteHalfSpinor, alignedAllocator<SiteHalfSpinor> > &buf,

lib/qcd/action/pseudofermion/TwoFlavour.h

@@ -133,6 +133,28 @@ class TwoFlavourPseudoFermionAction : public Action<typename Impl::GaugeField> {
     DerivativeSolver(MdagMOp, Phi, X);
     MdagMOp.Op(X, Y);
 
+    // Check hermiticity
+    /*
+    std::vector<int> seeds({1,2,3,4});
+    GridParallelRNG          RNG(U._grid);   RNG.SeedFixedIntegers(seeds);
+    FermionField RNGphi(FermOp.FermionGrid());
+    FermionField RNGchi(FermOp.FermionGrid());
+    FermionField Achi(FermOp.FermionGrid());
+    FermionField Aphi(FermOp.FermionGrid());
+
+    random(RNG, RNGphi);
+    random(RNG, RNGchi);
+    MdagMOp.Op(RNGchi, Achi);
+    MdagMOp.Op(RNGphi, Aphi);
+    ComplexD pAc = innerProduct(RNGphi, Achi);
+    ComplexD cAp = innerProduct(RNGchi, Aphi);
+    //these should be real
+    ComplexD pAp = innerProduct(RNGphi, Aphi);
+    ComplexD cAc = innerProduct(RNGchi, Achi);
+
+    std::cout<<GridLogMessage<< "pAc "<<pAc<<" cAp "<< cAp<< " diff "<<pAc-adj(cAp)<<std::endl;
+    std::cout << GridLogMessage << "pAp " << pAp << " cAc " << cAc << std::endl;
+*/
     // Our conventions really make this UdSdU; We do not differentiate wrt Udag
     // here.
     // So must take dSdU - adj(dSdU) and left multiply by mom to get dS/dt.

lib/qcd/hmc/HmcRunner.h

@@ -102,7 +102,8 @@ class NerscHmcRunnerTemplate {
     std::vector<int> ParSeed({6, 7, 8, 9, 10});
 
     // Create integrator, including the smearing policy
-    // 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
@@ -110,7 +111,9 @@ 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, SmearedConfiguration<Gimpl>, RepresentationsPolicy >
         IntegratorType;  // change here to change the algorithm
     IntegratorParameters MDpar(20, 1.0);
@@ -131,19 +134,19 @@ class NerscHmcRunnerTemplate {
       HMCpar.MetropolisTest = true;
       sRNG.SeedFixedIntegers(SerSeed);
       pRNG.SeedFixedIntegers(ParSeed);
-      SU3::HotConfiguration(pRNG, U);
+      SU<Nc>::HotConfiguration(pRNG, U);
     } else if (StartType == ColdStart) {
       // Cold start
       HMCpar.MetropolisTest = true;
       sRNG.SeedFixedIntegers(SerSeed);
       pRNG.SeedFixedIntegers(ParSeed);
-      SU3::ColdConfiguration(pRNG, U);
+      SU<Nc>::ColdConfiguration(pRNG, U);
     } else if (StartType == TepidStart) {
       // Tepid start
       HMCpar.MetropolisTest = true;
       sRNG.SeedFixedIntegers(SerSeed);
       pRNG.SeedFixedIntegers(ParSeed);
-      SU3::TepidConfiguration(pRNG, U);
+      SU<Nc>::TepidConfiguration(pRNG, U);
     } else if (StartType == CheckpointStart) {
       HMCpar.MetropolisTest = true;
       // CheckpointRestart

lib/qcd/hmc/integrators/Integrator.h

@@ -149,7 +149,7 @@ class Integrator {
     }
 
     // Force from the other representations
-    as[level].apply(update_P_hireps, Representations, Mom, U, ep);
+    //as[level].apply(update_P_hireps, Representations, Mom, U, ep);
   }
 
   void update_U(GaugeField& U, double ep) {
@@ -173,7 +173,7 @@ class Integrator {
     // Update the smeared fields, can be implemented as observer
     Smearer.set_GaugeField(U);
     // Update the higher representations fields
-    Representations.update(U);  // void functions if fundamental representation
+    //Representations.update(U);  // void functions if fundamental representation
   }
 
   virtual void step(GaugeField& U, int level, int first, int last) = 0;
@@ -203,6 +203,7 @@ class Integrator {
                     GridParallelRNG& pRNG) {
       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{};
 
@@ -216,7 +217,7 @@ class Integrator {
     // of the Metropolis
     Smearer.set_GaugeField(U);
     // Set the (eventual) representations gauge fields
-    // Representations.update(U);
+    //Representations.update(U);
 
     // The Smearer is attached to a pointer of the gauge field
     // automatically gets the correct field
@@ -230,7 +231,8 @@ class Integrator {
         as[level].actions.at(actionID)->refresh(Us, pRNG);
       }
 
-      as[level].apply(refresh_hireps, Representations, pRNG);
+      // Refresh the higher representation actions
+      //as[level].apply(refresh_hireps, Representations, pRNG);
     }
   }
 
@@ -240,12 +242,13 @@ class Integrator {
     template <class FieldType, class Repr>
     void operator()(std::vector<Action<FieldType>*> repr_set, Repr& Rep,
                     int level, RealD& H) {
-      RealD H_hirep = 0.0;
+      
       for (int a = 0; a < repr_set.size(); ++a) {
         RealD Hterm = repr_set.at(a)->S(Rep.U);
         std::cout << GridLogMessage << "S Level " << level << " term " << a
                   << " H Hirep = " << Hterm << std::endl;
         H += Hterm;
+
       }
     }
   } S_hireps{};
@@ -278,7 +281,7 @@ class Integrator {
                   << actionID << " H = " << Hterm << std::endl;
         H += Hterm;
       }
-      as[level].apply(S_hireps, Representations, level, H);
+      //as[level].apply(S_hireps, Representations, level, H);
     }
 
     return H;

lib/qcd/representations/adjoint.h

@@ -28,6 +28,7 @@ class AdjointRep {
 
   explicit AdjointRep(GridBase *grid) : U(grid) {}
   void update_representation(const LatticeGaugeField &Uin) {
+  	std::cout << GridLogDebug << "Updating adjoint representation\n" ;
     // Uin is in the fundamental representation
     // get the U in AdjointRep
     // (U_adj)_B = tr[e^a U e^b U^dag]
@@ -49,9 +50,18 @@ class AdjointRep {
       auto U_mu = peekLorentz(U, mu);
       for (int a = 0; a < Dimension; a++) {
         tmp = 2.0 * adj(Uin_mu) * ta[a] * Uin_mu;
-        for (int b = 0; b < (ncolour * ncolour - 1); b++)
+        for (int b = 0; b < Dimension; b++)
           pokeColour(U_mu, trace(tmp * ta[b]), a, b);
       }
+      // Check matrix U_mu, must be real orthogonal
+      //reality
+      LatticeMatrix Ucheck = U_mu - conjugate(U_mu);
+      std::cout << GridLogMessage << "Reality check: " << norm2(Ucheck) << std::endl;
+      LatticeMatrix uno(Uin._grid);
+      uno = 1.0;
+      Ucheck = U_mu * adj(U_mu) - uno;
+      std::cout << GridLogMessage << "orthogonality check: " << norm2(Ucheck) << std::endl;
+
       pokeLorentz(U, U_mu, mu);
     }
   }
@@ -59,6 +69,7 @@ class AdjointRep {
   LatticeGaugeField RtoFundamentalProject(const LatticeField &in,
                                           Real scale = 1.0) const {
     LatticeGaugeField out(in._grid);
+    out = zero;
 
     for (int mu = 0; mu < Nd; mu++) {
       LatticeColourMatrix out_mu(in._grid);  // fundamental representation
@@ -70,6 +81,8 @@ class AdjointRep {
       projectOnAlgebra(h, in_mu, scale);
       FundamentalLieAlgebraMatrix(h, out_mu, 1.0);  // apply scale only once
       pokeLorentz(out, out_mu, mu);
+	  // Returns traceless antihermitian matrix Nc * Nc.
+	  // Confirmed            
     }
     return out;
   }

lib/qcd/smearing/GaugeConfiguration.h

@@ -10,6 +10,29 @@ namespace Grid {
 
 namespace QCD {
 
+  //trivial class for no smearing
+  template< class Gimpl >
+class NoSmearing {
+public:
+  INHERIT_GIMPL_TYPES(Gimpl);
+
+  GaugeField*
+      ThinLinks;
+
+  NoSmearing(): ThinLinks(NULL) {}
+
+  void set_GaugeField(GaugeField& U) { ThinLinks = &U; }
+
+  void smeared_force(GaugeField& SigmaTilde) const {}
+
+  GaugeField& get_SmearedU() { return *ThinLinks; }
+
+  GaugeField& get_U(bool smeared = false) {
+    return *ThinLinks;
+  }
+
+};
+
 /*!
   @brief Smeared configuration container
 
@@ -201,6 +224,8 @@ class SmearedConfiguration {
   SmearedConfiguration()
       : smearingLevels(0), StoutSmearing(), SmearedSet(), ThinLinks(NULL) {}
 
+
+  
   // attach the smeared routines to the thin links U and fill the smeared set
   void set_GaugeField(GaugeField& U) { fill_smearedSet(U); }
 

lib/qcd/smearing/StoutSmearing.h

@@ -18,14 +18,12 @@ class Smear_Stout : public Smear<Gimpl> {
   INHERIT_GIMPL_TYPES(Gimpl)
 
   Smear_Stout(Smear<Gimpl>* base) : SmearBase(base) {
-    static_assert(Nc == 3,
+    assert(Nc == 3);//                  "Stout smearing currently implemented only for Nc==3");
-                  "Stout smearing currently implemented only for Nc==3");
   }
 
   /*! Default constructor */
   Smear_Stout(double rho = 1.0) : SmearBase(new Smear_APE<Gimpl>(rho)) {
-    static_assert(Nc == 3,
+    assert(Nc == 3);//                  "Stout smearing currently implemented only for Nc==3");
-                  "Stout smearing currently implemented only for Nc==3");
   }
 
   ~Smear_Stout() {}  // delete SmearBase...

lib/qcd/utils/SUn.h

@@ -48,7 +48,7 @@ class SU {
   using iSU2Matrix = iScalar<iScalar<iMatrix<vtype, 2> > >;
   template <typename vtype>
   using iSUnAlgebraVector =
-      iScalar<iScalar<iVector<vtype, (ncolour * ncolour - 1)> > >;
+      iScalar<iScalar<iVector<vtype, AdjointDimension> > >;
 
   //////////////////////////////////////////////////////////////////////////////////////////////////
   // Types can be accessed as SU<2>::Matrix , SU<2>::vSUnMatrix,
@@ -656,7 +656,7 @@ class SU {
     }
   }
 
-  static void FundamentalLieAlgebraMatrix(LatticeAlgebraVector &h,
+  static void FundamentalLieAlgebraMatrix(const LatticeAlgebraVector &h,
                                           LatticeMatrix &out,
                                           Real scale = 1.0) {
     conformable(h, out);
@@ -667,7 +667,7 @@ class SU {
     out = zero;
     for (int a = 0; a < AdjointDimension; a++) {
       generator(a, ta);
-      la = peekColour(h, a) * scale * ta;
+      la = peekColour(h, a) * timesI(ta) * scale;
       out += la;
     }
   }

lib/qcd/utils/SUnAdjoint.h

@@ -118,7 +118,7 @@ class SU_Adjoint : public SU<ncolour> {
 
     for (int a = 0; a < Dimension; a++) {
       generator(a, iTa);
-      auto tmp = real(trace(iTa * in)) * scale; 
+      auto tmp = - 2.0 * real(trace(iTa * in)) * scale; 
       pokeColour(h_out, tmp, a);
     }
   }
@@ -135,7 +135,7 @@ class SU_Adjoint : public SU<ncolour> {
     }
 
     for (int a = 0; a < Dimension; a++) {
-      auto tmp = real(trace(iTa[a] * in)) * scale; 
+      auto tmp = - 2.0*real(trace(iTa[a] * in)) * scale; 
       pokeColour(h_out, tmp, a);
     }
   }

lib/tensors/Tensor_exp.h

@@ -56,7 +56,7 @@ namespace Grid {
 	temp = unit + temp*arg;
       }
       
-      return ProjectOnGroup(temp);//maybe not strictly necessary
+      return temp;
       
     }
 

tests/Test_hmc_WilsonAdjointFermionGauge.cc

@@ -39,7 +39,8 @@ namespace Grid {
 namespace QCD {
 
 // Here change the allowed (higher) representations
-typedef Representations< FundamentalRepresentation, AdjointRepresentation > TheRepresentations;
+//typedef Representations< FundamentalRepresentation, AdjointRepresentation > TheRepresentations;
+typedef Representations< FundamentalRepresentation > TheRepresentations;
 
 
 class HmcRunner : public NerscHmcRunnerHirep< TheRepresentations > {
@@ -61,7 +62,7 @@ class HmcRunner : public NerscHmcRunnerHirep< TheRepresentations > {
 
     // temporarily need a gauge field
     LatticeGaugeField U(UGrid);
-    //AdjointRepresentation::LatticeField Ua(UGrid);
+    //AdjointRepresentation::LatticeField U(UGrid);
 
     // Gauge action
     WilsonGaugeActionR Waction(5.6);
@@ -69,7 +70,7 @@ class HmcRunner : public NerscHmcRunnerHirep< TheRepresentations > {
     Real mass = -0.77;
     FermionAction FermOp(U, *FGrid, *FrbGrid, mass);
 
-    ConjugateGradient<FermionField> CG(1.0e-6, 10000);
+    ConjugateGradient<FermionField> CG(1.0e-8, 10000);
 
     TwoFlavourPseudoFermionAction<ImplPolicy> Nf2(FermOp, CG, CG);