HMC for Adjoint fermions works

Accepts and reproduces known results

Check initial instability of inverters
when starting from hot configurations

lib/qcd/action/pseudofermion/TwoFlavour.h

@@ -133,29 +133,6 @@ class TwoFlavourPseudoFermionAction : public Action<typename Impl::GaugeField> {
     DerivativeSolver(MdagMOp, Phi, X); // X = (MdagM)^-1 phi    
     MdagMOp.Op(X, Y);                  // Y = M X = (Mdag)^-1 phi
 
-    // 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.HermOp(RNGchi, Achi);
-    MdagMOp.HermOp(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;
-    // These ones should be real
-    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/integrators/Integrator.h

@@ -120,7 +120,6 @@ 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: "
@@ -163,13 +162,14 @@ class Integrator {
               << " dt " << ep << " : t_U " << t_U << std::endl;
   }
   void update_U(GaugeField& Mom, GaugeField& U, double ep) {
-    // rewrite exponential to deal automatically  with the lorentz index?
+    // rewrite exponential to deal internally with the lorentz index?
     for (int mu = 0; mu < Nd; mu++) {
       auto Umu = PeekIndex<LorentzIndex>(U, mu);
       auto Pmu = PeekIndex<LorentzIndex>(Mom, mu);
       Umu = expMat(Pmu, ep, Params.Nexp) * Umu;
       PokeIndex<LorentzIndex>(U, ProjectOnGroup(Umu), mu);
     }
+    
     // Update the smeared fields, can be implemented as observer
     Smearer.set_GaugeField(U);
     // Update the higher representations fields

lib/qcd/representations/adjoint.h

@@ -57,7 +57,7 @@ class AdjointRep {
       for (int a = 0; a < Dimension; a++) {
         tmp = 2.0 * adj(Uin_mu) * ta[a] * Uin_mu;
         for (int b = 0; b < Dimension; b++)
-          pokeColour(U_mu, trace(tmp * ta[b]), b, a);
+          pokeColour(U_mu, trace(tmp * ta[b]), a, b);
       }
       pokeLorentz(U, U_mu, mu);
       // Check matrix U_mu, must be real orthogonal
@@ -71,7 +71,6 @@ class AdjointRep {
       std::cout << GridLogMessage << "orthogonality check: " << norm2(Ucheck) <<
       std::endl;
       */
-
     }
   }
 
@@ -87,8 +86,8 @@ class AdjointRep {
       out_mu = zero;
 
       typename SU<ncolour>::LatticeAlgebraVector h(in._grid);
-      projectOnAlgebra(h, in_mu, scale);
+      projectOnAlgebra(h, in_mu, double(Nc) * 2.0);  // factor C(r)/C(fund)
-      FundamentalLieAlgebraMatrix(h, out_mu, 1.0);  // apply scale only once
+      FundamentalLieAlgebraMatrix(h, out_mu);   // apply scale only once
       pokeLorentz(out, out_mu, mu);
       // Returns traceless antihermitian matrix Nc * Nc.
       // Confirmed

lib/qcd/utils/SUn.h

@@ -160,8 +160,9 @@ class SU {
     else
       generatorSigmaX(su2Index, ta);
   }
+  
   template <class cplx>
-  static void generatorSigmaX(int su2Index, iSUnMatrix<cplx> &ta) {
+  static void generatorSigmaY(int su2Index, iSUnMatrix<cplx> &ta) {
     ta = zero;
     int i1, i2;
     su2SubGroupIndex(i1, i2, su2Index);
@@ -170,15 +171,16 @@ class SU {
     ta = ta * 0.5;
   }
   template <class cplx>
-  static void generatorSigmaY(int su2Index, iSUnMatrix<cplx> &ta) {
+  static void generatorSigmaX(int su2Index, iSUnMatrix<cplx> &ta) {
     ta = zero;
     cplx i(0.0, 1.0);
     int i1, i2;
     su2SubGroupIndex(i1, i2, su2Index);
-    ta()()(i1, i2) = -i;
+    ta()()(i1, i2) = i;
-    ta()()(i2, i1) = i;
+    ta()()(i2, i1) = -i;
     ta = ta * 0.5;
   }
+
   template <class cplx>
   static void generatorDiagonal(int diagIndex, iSUnMatrix<cplx> &ta) {
     // diag ({1, 1, ..., 1}(k-times), -k, 0, 0, ...)
@@ -651,9 +653,10 @@ class SU {
     for (int a = 0; a < AdjointDimension; a++) {
       gaussian(pRNG, ca);
       generator(a, ta);
-      la = toComplex(ca) * ci * ta;
+      la = toComplex(ca) * ta;
       out += la;
     }
+    out *= ci;
   }
 
   static void FundamentalLieAlgebraMatrix(const LatticeAlgebraVector &h,
@@ -684,7 +687,6 @@ class SU {
       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;
   }
 
 

lib/qcd/utils/SUnAdjoint.h

@@ -75,7 +75,8 @@ class SU_Adjoint : public SU<ncolour> {
         typename SU<ncolour>::template iSUnMatrix<cplx> tmp1 =
             2.0 * tmp * ta[b];  // 2.0 from the normalization
         Complex iTr = TensorRemove(timesI(trace(tmp1)));
-        iAdjTa()()(b, a) = iTr;
+        //iAdjTa()()(b, a) = iTr;
+        iAdjTa()()(a, b) = iTr;
       }
     }
   }
@@ -121,9 +122,10 @@ class SU_Adjoint : public SU<ncolour> {
     out = zero;
     for (int a = 0; a < Dimension; a++) {
       generator(a, iTa);
-      la = peekColour(h, a) * iTa * scale;
+      la = peekColour(h, a) * iTa;
       out += la;
     }
+    out *= scale;
   }
 
   // Projects the algebra components a lattice matrix (of dimension ncol*ncol -1 )
@@ -131,16 +133,16 @@ class SU_Adjoint : public SU<ncolour> {
     conformable(h_out, in);
     h_out = zero;
     AMatrix iTa;
+    Real coefficient = - 1.0/(ncolour) * scale;// 1/Nc for the normalization of the trace in the adj rep
 
     for (int a = 0; a < Dimension; a++) {
       generator(a, iTa);
-      auto tmp = - 1.0/(ncolour) * (trace(iTa * in)) * scale;// 1/Nc for the normalization of the trace in the adj rep
+      auto tmp = real(trace(iTa * in)) * coefficient;
       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. 
+  // a projector that keeps the generators stored to avoid the overhead of recomputing them 
   static void projector(typename SU<ncolour>::LatticeAlgebraVector &h_out, const LatticeAdjMatrix &in, Real scale = 1.0) {
     conformable(h_out, in);
     static std::vector<AMatrix> iTa(Dimension);  // to store the generators
@@ -151,8 +153,11 @@ class SU_Adjoint : public SU<ncolour> {
         for (int a = 0; a < Dimension; a++) generator(a, iTa[a]);
     }
 
+    Real coefficient = -1.0 / (ncolour)*scale;  // 1/Nc for the normalization of
+                                                // the trace in the adj rep
+
     for (int a = 0; a < Dimension; a++) {
-      auto tmp = - 1.0/(ncolour) * (trace(iTa[a] * in)) * scale; 
+      auto tmp = real(trace(iTa[a] * in)) * coefficient; 
       pokeColour(h_out, tmp, a);
     }
   }

tests/Test_hmc_WilsonAdjointFermionGauge.cc

@@ -65,14 +65,15 @@ class HmcRunner : public NerscHmcRunnerHirep< TheRepresentations > {
     AdjointRepresentation::LatticeField U(UGrid);
 
     // Gauge action
-    WilsonGaugeActionR Waction(2.0);
+    WilsonGaugeActionR Waction(2.25);
 
-    Real mass = -1.0;
+    Real mass = -0.95;
     FermionAction FermOp(U, *FGrid, *FrbGrid, mass);
 
     ConjugateGradient<FermionField> CG(1.0e-8, 10000);
     ConjugateResidual<FermionField> CR(1.0e-8, 10000);
 
+    // Pass two solvers: one for the force computation and one for the action
     TwoFlavourPseudoFermionAction<ImplPolicy> Nf2(FermOp, CG, CG);
 
     // Set smearing (true/false), default: false

tests/Test_lie_generators.cc

@@ -121,6 +121,7 @@ int main(int argc, char** argv) {
   // Test group structure
   // (U_f * V_f)_r = U_r * V_r
   LatticeGaugeField UV(grid);
+  UV = zero;
   for (int mu = 0; mu < Nd; mu++) {
     SU<Nc>::LatticeMatrix Umu = peekLorentz(U,mu);
     SU<Nc>::LatticeMatrix Vmu = peekLorentz(V,mu);
@@ -138,6 +139,7 @@ int main(int argc, char** argv) {
   typename AdjointRep<Nc>::LatticeField Vr = AdjRep.U;  // V_r
 
   typename AdjointRep<Nc>::LatticeField UrVr(grid);
+  UrVr = zero;
   for (int mu = 0; mu < Nd; mu++) {
     typename AdjointRep<Nc>::LatticeMatrix Urmu = peekLorentz(Ur,mu);
     typename AdjointRep<Nc>::LatticeMatrix Vrmu = peekLorentz(Vr,mu);
@@ -145,7 +147,7 @@ int main(int argc, char** argv) {
   }
 
   typename AdjointRep<Nc>::LatticeField Diff_check = UVr - UrVr;
-  std::cout << GridLogMessage << "Group structure check difference : " << norm2(Diff_check) << std::endl;
+  std::cout << GridLogMessage << "Group structure SU("<<Nc<<") check difference : " << norm2(Diff_check) << std::endl;
 
   // Check correspondence of algebra and group transformations
   // Create a random vector

tests/Test_wilson_force.cc

@@ -58,8 +58,8 @@ int main (int argc, char ** argv)
 
   LatticeGaugeField U(&Grid);
 
-  SU2::HotConfiguration(pRNG,U);
+  //SU2::HotConfiguration(pRNG,U);
-  //  SU3::ColdConfiguration(pRNG,U);
+  SU3::ColdConfiguration(pRNG,U);
   
   ////////////////////////////////////
   // Unmodified matrix element
@@ -95,7 +95,7 @@ int main (int argc, char ** argv)
 
   for(int mu=0;mu<Nd;mu++){
 
-    SU2::GaussianFundamentalLieAlgebraMatrix(pRNG, mommu); // Traceless antihermitian momentum; gaussian in lie alg
+    SU3::GaussianFundamentalLieAlgebraMatrix(pRNG, mommu); // Traceless antihermitian momentum; gaussian in lie alg
 
     Hmom -= real(sum(trace(mommu*mommu)));
 
@@ -142,12 +142,14 @@ int main (int argc, char ** argv)
   LatticeComplex dS(&Grid); dS = zero;
   LatticeComplex dSmom(&Grid); dSmom = zero;
   LatticeComplex dSmom2(&Grid); dSmom2 = zero;
+  
   for(int mu=0;mu<Nd;mu++){
     mommu   = PeekIndex<LorentzIndex>(UdSdU,mu);
     mommu=Ta(mommu)*2.0;
     PokeIndex<LorentzIndex>(UdSdU,mommu,mu);
   }
   
+
   for(int mu=0;mu<Nd;mu++){
     mommu   = PeekIndex<LorentzIndex>(mom,mu);
     std::cout << GridLogMessage<< " Mommu  " << norm2(mommu)<<std::endl;