Clover term force ok

lib/qcd/action/fermion/WilsonCloverFermion.cc

@@ -48,10 +48,7 @@ RealD WilsonCloverFermion<Impl>::M(const FermionField &in, FermionField &out)
   // Clover term
   Mooee(in, temp);
 
-  //hack
+  out= zero;
-  out = zero;
-
-
   out += temp;
   return norm2(out);
 }
@@ -68,9 +65,7 @@ RealD WilsonCloverFermion<Impl>::Mdag(const FermionField &in, FermionField &out)
   // Clover term
   MooeeDag(in, temp);
 
-   //hack
+  out=zero;
-  out = zero;
-
   out += temp;
   return norm2(out);
 }

lib/qcd/action/fermion/WilsonCloverFermion.h

@@ -122,6 +122,8 @@ public:
     // Clover term derivative
     ///////////////////////////////////////////////////////////
     Impl::outerProductImpl(Lambda, X, Y); 
+    //std::cout << "Lambda:" << Lambda << std::endl;
+
 
     Gamma::Algebra sigma[] = {
         Gamma::Algebra::SigmaXY,
@@ -153,16 +155,18 @@ public:
       for (int nu = 0; nu < 4; nu++)
       {
         if (mu == nu) continue;
-        PropagatorField Slambda = Gamma(sigma[count]) * Lambda;
+        PropagatorField Slambda = Gamma(sigma[count]) * Lambda; // sigma checked
-        Impl::TraceSpinImpl(lambda, Slambda); //traceSpin
+        Impl::TraceSpinImpl(lambda, Slambda); // traceSpin ok
-        force_mu += Cmunu(U, lambda, mu, nu);
+        force_mu -= Cmunu(U, lambda, mu, nu); // checked
         count++;
       }
 
       pokeLorentz(clover_force, U[mu] * force_mu, mu);
     }
-    clover_force *= csw / 8.;
+    clover_force *= csw;
     force += clover_force;
+
+
   }
 
   // Computing C_{\mu \nu}(x) as in Eq.(B.39) in Zbigniew Sroczynski's PhD thesis
@@ -170,20 +174,19 @@ public:
   {
     conformable(lambda._grid, U[0]._grid);
     GaugeLinkField out(lambda._grid), tmp(lambda._grid);
-
     // insertion in upper staple
     // please check redundancy of shift operations
-
+    
     // C1+
     tmp = lambda * U[nu];
     out = Impl::ShiftStaple(Impl::CovShiftForward(tmp, nu, Impl::CovShiftBackward(U[mu], mu, Impl::CovShiftIdentityBackward(U[nu], nu))), mu);
-
+    
     // C2+
-    tmp = U[mu] * Impl::CovShiftIdentityForward(adj(lambda), mu);
+    tmp = U[mu] * Impl::ShiftStaple(adj(lambda), mu);
     out += Impl::ShiftStaple(Impl::CovShiftForward(U[nu], nu, Impl::CovShiftBackward(tmp, mu, Impl::CovShiftIdentityBackward(U[nu], nu))), mu);
-
+    
     // C3+
-    tmp = U[nu] * Impl::CovShiftIdentityForward(adj(lambda), nu);
+    tmp = U[nu] * Impl::ShiftStaple(adj(lambda), nu);
     out += Impl::ShiftStaple(Impl::CovShiftForward(U[nu], nu, Impl::CovShiftBackward(U[mu], mu, Impl::CovShiftIdentityBackward(tmp, nu))), mu);
 
     // C4+
@@ -259,6 +262,7 @@ private:
     PARALLEL_FOR_LOOP
     for (int i = 0; i < CloverTerm._grid->oSites(); i++)
     {
+
       T._odata[i]()(0, 0) = timesMinusI(F._odata[i]()());
       T._odata[i]()(1, 1) = timesI(F._odata[i]()());
       T._odata[i]()(2, 2) = timesMinusI(F._odata[i]()());

lib/qcd/utils/WilsonLoops.h

@@ -327,7 +327,7 @@ static void StapleMult(GaugeMat &staple, const GaugeLorentz &Umu, int mu) {
   static void FieldStrength(GaugeMat &FS, const GaugeLorentz &Umu, int mu, int nu){
       // Fmn +--<--+  Ut +--<--+
       //     |     |     |     |
-      //  (x)+-->--+     +-->--+(x)
+      //  (x)+-->--+     +-->--+(x)  - h.c.
       //     |     |     |     |
       //     +--<--+     +--<--+
 

tests/forces/Test_wilsonclover_force.cc

@@ -1,6 +1,6 @@
 /*************************************************************************************
 
-    Grid physics library, www.github.com/paboyle/Grid 
+    Grid physics library, www.github.com/paboyle/Grid
 
     Source file: ./tests/Test_wilson_force.cc
 
@@ -45,14 +45,17 @@ int main(int argc, char **argv)
   int threads = GridThread::GetThreads();
   std::cout << GridLogMessage << "Grid is setup to use " << threads << " threads" << std::endl;
 
-  std::vector<int> seeds({1, 2, 3, 4});
+  std::vector<int> seeds({1, 2, 30, 50});
 
   GridParallelRNG pRNG(&Grid);
+  
   std::vector<int> vrand(4);
   std::srand(std::time(0));
   std::generate(vrand.begin(), vrand.end(), std::rand);
   std::cout << GridLogMessage << vrand << std::endl;
   pRNG.SeedFixedIntegers(vrand);
+  
+  //pRNG.SeedFixedIntegers(seeds);
 
   LatticeFermion phi(&Grid);
   gaussian(pRNG, phi);
@@ -61,16 +64,53 @@ int main(int argc, char **argv)
 
   LatticeGaugeField U(&Grid);
 
-  SU3::HotConfiguration(pRNG,U);
+/*
+  std::vector<int> x(4); // 4d fermions
+  std::vector<int> gd = Grid.GlobalDimensions();
+  Grid::QCD::SpinColourVector F;
+  Grid::Complex c;
+
+  phi = zero;
+  for (x[0] = 0; x[0] < 1; x[0]++)
+  {
+    for (x[1] = 0; x[1] < 1; x[1]++)
+    {
+      for (x[2] = 0; x[2] < 1; x[2]++)
+      {
+        for (x[3] = 0; x[3] < 1; x[3]++)
+        {
+          for (int sp = 0; sp < 4; sp++)
+          {
+            for (int j = 0; j < 3; j++) // colours
+            {
+              F()(sp)(j) = Grid::Complex(0.0,0.0);
+              if (((sp == 0) && (j==0)))
+              {
+                c = Grid::Complex(1.0, 0.0);
+                F()(sp)(j) = c;
+              }
+            }
+          }
+          Grid::pokeSite(F, phi, x);
+
+        }
+      }
+    }
+  }
+*/
+
+  std::vector<int> site = {0, 0, 0, 0};
+  SU3::HotConfiguration(pRNG, U);
   //SU3::ColdConfiguration(pRNG, U);
 
+
   ////////////////////////////////////
   // Unmodified matrix element
   ////////////////////////////////////
   RealD mass = -4.0; //kills the diagonal term
   Real csw = 1.0;
   WilsonCloverFermionR Dw(U, Grid, RBGrid, mass, csw);
-
+  Dw.ImportGauge(U);
   Dw.M(phi, Mphi);
   ComplexD S = innerProduct(Mphi, Mphi); // Action : pdag MdagM p
 
@@ -78,11 +118,23 @@ int main(int argc, char **argv)
   LatticeGaugeField UdSdU(&Grid);
   LatticeGaugeField tmp(&Grid);
 
-  Dw.MDeriv(tmp, Mphi, phi, DaggerNo);  UdSdU  = tmp;
+  ////////////////////////////////////////////
-  Dw.MDeriv(tmp, phi, Mphi, DaggerYes); UdSdU += tmp;
+  Dw.MDeriv(tmp, Mphi, phi, DaggerNo);
+  UdSdU = tmp;
+  Dw.MDeriv(tmp, phi, Mphi, DaggerYes);
+  UdSdU += tmp;
+  /////////////////////////////////////////////
+
   // Take the traceless antihermitian component
-  UdSdU = Ta(UdSdU);
+  //UdSdU = Ta(UdSdU);
-  
+
+  //std::cout << UdSdU << std::endl;
+  //SU3::LatticeAlgebraVector hforce(&Grid);
+  LatticeColourMatrix mommu(&Grid);
+  //mommu = PeekIndex<LorentzIndex>(UdSdU, 0);
+  //SU3::projectOnAlgebra(hforce, mommu);
+  //std::cout << hforce << std::endl;
+
   ////////////////////////////////////
   // Modify the gauge field a little
   ////////////////////////////////////
@@ -90,28 +142,63 @@ int main(int argc, char **argv)
   RealD Hmom = 0.0;
   RealD Hmomprime = 0.0;
   RealD Hmompp = 0.0;
-  LatticeColourMatrix mommu(&Grid);
   LatticeColourMatrix forcemu(&Grid);
   LatticeGaugeField mom(&Grid);
   LatticeGaugeField Uprime(&Grid);
 
+  
   for (int mu = 0; mu < Nd; mu++) {
     // Traceless antihermitian momentum; gaussian in lie alg
     SU3::GaussianFundamentalLieAlgebraMatrix(pRNG, mommu);
     Hmom -= real(sum(trace(mommu * mommu)));
     PokeIndex<LorentzIndex>(mom, mommu, mu);
   }
+  /*
+  SU3::AlgebraVector h;
+  SU3::LatticeAlgebraVector hl(&Grid);
+  h()()(0) = 1.0;
+  hl = zero;
+  pokeSite(h, hl, site);
+  SU3::FundamentalLieAlgebraMatrix(hl, mommu);
+  mom = zero;
+  PokeIndex<LorentzIndex>(mom, mommu, 0);
+  Hmom -= real(sum(trace(mommu * mommu)));
+  */
 
+  /*
   parallel_for(int ss=0;ss<mom._grid->oSites();ss++){
-    for (int mu = 0; mu < Nd; mu++) 
+    for (int mu = 0; mu < Nd; mu++)
       Uprime[ss]._internal[mu] = ProjectOnGroup(Exponentiate(mom[ss]._internal[mu], dt, 12) * U[ss]._internal[mu]);
-  } 
+  }
+*/
+
+  for (int mu = 0; mu < Nd; mu++)
+  {
+    parallel_for(auto i = mom.begin(); i < mom.end(); i++)
+    {
+      Uprime[i](mu) = U[i](mu);
+      Uprime[i](mu) += mom[i](mu) * U[i](mu) * dt;
+      Uprime[i](mu) += mom[i](mu) * mom[i](mu) * U[i](mu) * (dt * dt / 2.0);
+      Uprime[i](mu) += mom[i](mu) * mom[i](mu) * mom[i](mu) * U[i](mu) * (dt * dt * dt / 6.0);
+      Uprime[i](mu) += mom[i](mu) * mom[i](mu) * mom[i](mu) * mom[i](mu) * U[i](mu) * (dt * dt * dt * dt / 24.0);
+      Uprime[i](mu) += mom[i](mu) * mom[i](mu) * mom[i](mu) * mom[i](mu) * mom[i](mu) * U[i](mu) * (dt * dt * dt * dt * dt / 120.0);
+      Uprime[i](mu) += mom[i](mu) * mom[i](mu) * mom[i](mu) * mom[i](mu) * mom[i](mu) * mom[i](mu) * U[i](mu) * (dt * dt * dt * dt * dt * dt / 720.0);
+    }
+  }
 
   std::cout << GridLogMessage << "Initial mom hamiltonian is " << Hmom << std::endl;
-  
+
   // New action
+  LatticeGaugeField diff(&Grid);
+  diff = Uprime - U;
+  //std::cout << "Diff:" << diff << std::endl;
   Dw.ImportGauge(Uprime);
   Dw.M(phi, MphiPrime);
+  LatticeFermion DiffFermion(&Grid);
+  DiffFermion = MphiPrime - Mphi;
+  //std::cout << "DiffFermion:" << DiffFermion << std::endl;
+  //std::cout << "Mphi:" << Mphi << std::endl;
+  //std::cout << "MphiPrime:" << MphiPrime << std::endl;
 
   ComplexD Sprime = innerProduct(MphiPrime, MphiPrime);
 
@@ -143,16 +230,14 @@ int main(int argc, char **argv)
   dSmom2 = zero;
 
 
-  // need for this??? 
-  // ultimately it is just a 2.0 factor in UdSdU
   for (int mu = 0; mu < Nd; mu++)
   {
-    mommu = PeekIndex<LorentzIndex>(UdSdU, mu); // P_mu = 
+    mommu = PeekIndex<LorentzIndex>(UdSdU, mu); // P_mu =
     mommu = Ta(mommu) * 2.0;                    // Mom = (P_mu - P_mu^dag) - trace(P_mu - P_mu^dag)
     PokeIndex<LorentzIndex>(UdSdU, mommu, mu);  // UdSdU_mu = Mom
   }
 
-  std::cout << GridLogMessage<< "Antihermiticity tests - 2 " << std::endl;
+  std::cout << GridLogMessage << "Antihermiticity tests - 2 " << std::endl;
   for (int mu = 0; mu < Nd; mu++)
   {
     mommu = PeekIndex<LorentzIndex>(mom, mu);
@@ -167,7 +252,6 @@ int main(int argc, char **argv)
   }
   /////////////////////////////////////////////////////
 
-
   for (int mu = 0; mu < Nd; mu++)
   {
     forcemu = PeekIndex<LorentzIndex>(UdSdU, mu);