Added calculation of 5Li topological charge

WilsonFlow code now calls topological charge calculation with correct gauge implementation rather than assuming periodic
Added version of WilsonFlow::flowMeasureEnergyDensityPlaquette that outputs the smeared gauge field at the end

Grid/qcd/smearing/WilsonFlow.h

@@ -82,6 +82,9 @@ public:
   RealD energyDensityPlaquette(const GaugeField& U) const;
 
   //Evolve the gauge field by Nstep steps of epsilon and return the energy density computed every interval steps
+  //The smeared field is output as V
+  std::vector<RealD> flowMeasureEnergyDensityPlaquette(GaugeField &V, const GaugeField& U) const;
+
   std::vector<RealD> flowMeasureEnergyDensityPlaquette(const GaugeField& U) const;
 };
 
@@ -162,9 +165,9 @@ RealD WilsonFlow<Gimpl>::energyDensityPlaquette(const GaugeField& U) const {
 }
 
 template <class Gimpl>
-std::vector<RealD> WilsonFlow<Gimpl>::flowMeasureEnergyDensityPlaquette(const GaugeField& U) const{
+std::vector<RealD> WilsonFlow<Gimpl>::flowMeasureEnergyDensityPlaquette(GaugeField &V, const GaugeField& U) const{
   std::vector<RealD> out;
-  GaugeField V(U);
+  V = U;
   for (unsigned int step = 0; step < Nstep; step++) { //bn  tau = epsilon*(step+1) so tau after performing step=0 is epsilon
     std::cout << GridLogMessage << "[WilsonFlow] Evolving step " << step << std::endl;
     evolve_step(V);
@@ -175,7 +178,14 @@ std::vector<RealD> WilsonFlow<Gimpl>::flowMeasureEnergyDensityPlaquette(const Ga
   }
   return out;
 }
-  
+
+template <class Gimpl>
+std::vector<RealD> WilsonFlow<Gimpl>::flowMeasureEnergyDensityPlaquette(const GaugeField& U) const{
+  GaugeField V(U);
+  return flowMeasureEnergyDensityPlaquette(V,U);
+}
+
+
 //#define WF_TIMING 
 template <class Gimpl>
 void WilsonFlow<Gimpl>::smear(GaugeField& out, const GaugeField& in) const {
@@ -194,7 +204,7 @@ void WilsonFlow<Gimpl>::smear(GaugeField& out, const GaugeField& in) const {
     if( step % measure_interval == 0){
       std::cout << GridLogMessage << "[WilsonFlow] Top. charge           : "
 		<< step << "  " 
-		<< WilsonLoops<PeriodicGimplR>::TopologicalCharge(out) << std::endl;
+		<< WilsonLoops<Gimpl>::TopologicalCharge(out) << std::endl;
     }
   }
 }
@@ -214,7 +224,7 @@ void WilsonFlow<Gimpl>::smear_adaptive(GaugeField& out, const GaugeField& in, Re
     if( step % measure_interval == 0){
       std::cout << GridLogMessage << "[WilsonFlow] Top. charge           : "
 		<< step << "  " 
-		<< WilsonLoops<PeriodicGimplR>::TopologicalCharge(out) << std::endl;
+		<< WilsonLoops<Gimpl>::TopologicalCharge(out) << std::endl;
     }
   } while (taus < maxTau);
 

Grid/qcd/utils/WilsonLoops.h

@@ -390,6 +390,171 @@ public:
   }
 
 
+  //Clover-leaf Wilson loop combination for arbitrary mu-extent M and nu extent N,  mu >= nu
+  //cf  https://arxiv.org/pdf/hep-lat/9701012.pdf Eq 7  for 1x2 Wilson loop    
+  //Clockwise ordering
+  static void CloverleafMxN(GaugeMat &FS, const GaugeMat &Umu, const GaugeMat &Unu, int mu, int nu, int M, int N){  
+#define Fmu(A) Gimpl::CovShiftForward(Umu, mu, A)
+#define Bmu(A) Gimpl::CovShiftBackward(Umu, mu, A)
+#define Fnu(A) Gimpl::CovShiftForward(Unu, nu, A)
+#define Bnu(A) Gimpl::CovShiftBackward(Unu, nu, A)
+#define FmuI Gimpl::CovShiftIdentityForward(Umu, mu)
+#define BmuI Gimpl::CovShiftIdentityBackward(Umu, mu)
+#define FnuI Gimpl::CovShiftIdentityForward(Unu, nu)
+#define BnuI Gimpl::CovShiftIdentityBackward(Unu, nu)
+
+    //Upper right loop
+    GaugeMat tmp = BmuI;
+    for(int i=1;i<M;i++)
+      tmp = Bmu(tmp);
+    for(int j=0;j<N;j++)
+      tmp = Bnu(tmp);
+    for(int i=0;i<M;i++)
+      tmp = Fmu(tmp);
+    for(int j=0;j<N;j++)
+      tmp = Fnu(tmp);
+      
+    FS = tmp;
+
+    //Upper left loop
+    tmp = BnuI;
+    for(int j=1;j<N;j++)
+      tmp = Bnu(tmp);
+    for(int i=0;i<M;i++)
+      tmp = Fmu(tmp);
+    for(int j=0;j<N;j++)
+      tmp = Fnu(tmp);
+    for(int i=0;i<M;i++)
+      tmp = Bmu(tmp);
+      
+    FS = FS + tmp;
+
+    //Lower right loop
+    tmp = FnuI;
+    for(int j=1;j<N;j++)
+      tmp = Fnu(tmp);
+    for(int i=0;i<M;i++)
+      tmp = Bmu(tmp);
+    for(int j=0;j<N;j++)
+      tmp = Bnu(tmp);
+    for(int i=0;i<M;i++)
+      tmp = Fmu(tmp);
+      
+    FS = FS + tmp;
+
+    //Lower left loop
+    tmp = FmuI;
+    for(int i=1;i<M;i++)
+      tmp = Fmu(tmp);
+    for(int j=0;j<N;j++)
+      tmp = Fnu(tmp);
+    for(int i=0;i<M;i++)
+      tmp = Bmu(tmp);
+    for(int j=0;j<N;j++)
+      tmp = Bnu(tmp);
+
+    FS = FS + tmp;
+
+#undef Fmu
+#undef Bmu
+#undef Fnu
+#undef Bnu
+#undef FmuI
+#undef BmuI
+#undef FnuI
+#undef BnuI
+  }
+
+  //Field strength from MxN Wilson loop
+  //Note F_numu = - F_munu
+  static void FieldStrengthMxN(GaugeMat &FS, const GaugeLorentz &U, int mu, int nu, int M, int N){  
+    GaugeMat Umu = PeekIndex<LorentzIndex>(U, mu);
+    GaugeMat Unu = PeekIndex<LorentzIndex>(U, nu);
+    if(M == N){
+      GaugeMat F(Umu.Grid());
+      CloverleafMxN(F, Umu, Unu, mu, nu, M, N);
+      FS = 0.125 * ( F - adj(F) );
+    }else{
+      //Average over both orientations
+      GaugeMat horizontal(Umu.Grid()), vertical(Umu.Grid());
+      CloverleafMxN(horizontal, Umu, Unu, mu, nu, M, N);
+      CloverleafMxN(vertical, Umu, Unu, mu, nu, N, M);
+      FS = 0.0625 * ( horizontal - adj(horizontal) + vertical - adj(vertical) );
+    }
+  }
+
+  //Topological charge contribution from MxN Wilson loops
+  //cf  https://arxiv.org/pdf/hep-lat/9701012.pdf  Eq 6
+  static Real TopologicalChargeMxN(const GaugeLorentz &U, int M, int N){
+    assert(Nd == 4);
+    std::vector<std::vector<GaugeMat*> > F(Nd,std::vector<GaugeMat*>(Nd,nullptr));
+    //Note F_numu = - F_munu
+    //hence we only need to loop over mu,nu,rho,sigma that aren't related by permuting mu,nu  or rho,sigma
+    //Use nu > mu
+    for(int mu=0;mu<Nd-1;mu++){
+      for(int nu=mu+1; nu<Nd; nu++){
+	F[mu][nu] = new GaugeMat(U.Grid());
+	FieldStrengthMxN(*F[mu][nu], U, mu, nu, M, N);
+      }
+    }
+    Real coeff = -1./(32 * M_PI*M_PI * M*M * N*N); //overall sign to match CPS and Grid conventions, possibly related to time direction = 3 vs 0
+
+    static const int combs[3][4] = { {0,1,2,3}, {0,2,1,3}, {0,3,1,2} };
+    static const int signs[3] = { 1, -1, 1 }; //epsilon_{mu nu rho sigma}
+
+    ComplexField fsum(U.Grid());
+    fsum = Zero();
+    for(int c=0;c<3;c++){
+      int mu = combs[c][0], nu = combs[c][1], rho = combs[c][2], sigma = combs[c][3];
+      int eps = signs[c];
+      fsum = fsum + (8. * coeff * eps) * trace( (*F[mu][nu]) * (*F[rho][sigma]) ); 
+    }
+
+    for(int mu=0;mu<Nd-1;mu++)
+      for(int nu=mu+1; nu<Nd; nu++)
+	delete F[mu][nu];
+
+    auto Tq = sum(fsum);
+    return TensorRemove(Tq).real();
+  }
+
+  //Generate the contributions to the 5Li topological charge from Wilson loops of the following sizes
+  //Use coefficients from hep-lat/9701012
+  //1x1 : c1=(19.-55.*c5)/9.
+  //2x2 : c2=(1-64.*c5)/9.
+  //1x2 : c3=(-64.+640.*c5)/45.
+  //1x3 : c4=1./5.-2.*c5
+  //3x3 : c5=1./20.
+  //Output array contains the loops in the above order
+  static std::vector<Real> TopologicalCharge5LiContributions(const GaugeLorentz &U){
+    static const int exts[5][2] = { {1,1}, {2,2}, {1,2}, {1,3}, {3,3} };       
+    std::vector<Real> out(5);
+    std::cout << GridLogMessage << "Computing topological charge" << std::endl;
+    for(int i=0;i<5;i++){	
+      out[i] = TopologicalChargeMxN(U,exts[i][0],exts[i][1]);
+      std::cout << GridLogMessage << exts[i][0] << "x" << exts[i][1] << " Wilson loop contribution " << out[i] << std::endl;
+    }
+    return out;
+  }
+
+  //Compute the 5Li topological charge
+  static Real TopologicalCharge5Li(const GaugeLorentz &U){
+    std::vector<Real> loops = TopologicalCharge5LiContributions(U);
+
+    double c5=1./20.;
+    double c4=1./5.-2.*c5;
+    double c3=(-64.+640.*c5)/45.;
+    double c2=(1-64.*c5)/9.;
+    double c1=(19.-55.*c5)/9.;
+
+    double Q = c1*loops[0] + c2*loops[1] + c3*loops[2] + c4*loops[3] + c5*loops[4];
+
+    std::cout << GridLogMessage << "5Li Topological charge: " << Q << std::endl;
+    return Q;
+  }
+
+
+
   //////////////////////////////////////////////////////
   // Similar to above for rectangle is required
   //////////////////////////////////////////////////////