Full clover force correct

lib/qcd/action/fermion/WilsonCloverFermion.cc

@@ -43,12 +43,11 @@ RealD WilsonCloverFermion<Impl>::M(const FermionField &in, FermionField &out)
 
   // Wilson term
   out.checkerboard = in.checkerboard;
-  //this->Dhop(in, out, DaggerNo);
+  this->Dhop(in, out, DaggerNo);
 
   // Clover term
   Mooee(in, temp);
 
-  out= zero;
   out += temp;
   return norm2(out);
 }
@@ -60,12 +59,11 @@ RealD WilsonCloverFermion<Impl>::Mdag(const FermionField &in, FermionField &out)
 
   // Wilson term
   out.checkerboard = in.checkerboard;
-  //this->Dhop(in, out, DaggerYes);
+  this->Dhop(in, out, DaggerYes);
 
   // Clover term
   MooeeDag(in, temp);
 
-  out=zero;
   out += temp;
   return norm2(out);
 }

lib/qcd/action/fermion/WilsonCloverFermion.h

@@ -116,7 +116,7 @@ public:
     
     force = zero;
     // Derivative of the Wilson hopping term
-    //this->DhopDeriv(force, X, Y, dag);
+    this->DhopDeriv(force, X, Y, dag);
 
     ///////////////////////////////////////////////////////////
     // Clover term derivative

tests/forces/Test_wilsonclover_force.cc

@@ -48,13 +48,12 @@ int main(int argc, char **argv)
   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);
@@ -64,50 +63,14 @@ int main(int argc, char **argv)
 
   LatticeGaugeField U(&Grid);
 
-/*
-  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);
+  //SU3::ColdConfiguration(pRNG, U);// Clover term zero
-
 
   ////////////////////////////////////
   // Unmodified matrix element
   ////////////////////////////////////
-  RealD mass = -4.0; //kills the diagonal term
+  RealD mass = 0.1;
   Real csw = 1.0;
   WilsonCloverFermionR Dw(U, Grid, RBGrid, mass, csw);
   Dw.ImportGauge(U);
@@ -125,103 +88,42 @@ int main(int argc, char **argv)
   UdSdU += tmp;
   /////////////////////////////////////////////
 
-  // Take the traceless antihermitian component
-  //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
   ////////////////////////////////////
-  RealD dt = 0.0001;
+  RealD dt = 0.00005;
   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++)
-  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++)
-  parallel_for(int ss=0;ss<mom._grid->oSites();ss++){
-    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[ss]._internal[mu] = ProjectOnGroup(Exponentiate(mom[ss]._internal[mu], dt, 12) * U[ss]._internal[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);
 
   //////////////////////////////////////////////
   // Use derivative to estimate dS
   //////////////////////////////////////////////
 
-  ///////////////////////////////////////////////////////
-  std::cout << GridLogMessage << "Antihermiticity tests - 1 " << std::endl;
-  for (int mu = 0; mu < Nd; mu++)
-  {
-    mommu = PeekIndex<LorentzIndex>(mom, mu);
-    std::cout << GridLogMessage << " Mommu  " << norm2(mommu) << std::endl;
-    mommu = mommu + adj(mommu);
-    std::cout << GridLogMessage << " Test: Mommu + Mommudag " << norm2(mommu) << std::endl;
-    mommu = PeekIndex<LorentzIndex>(UdSdU, mu);
-    std::cout << GridLogMessage << " dsdumu  " << norm2(mommu) << std::endl;
-    mommu = mommu + adj(mommu);
-    std::cout << GridLogMessage << " Test: dsdumu + dag  " << norm2(mommu) << std::endl;
-    std::cout << "" << std::endl;
-  }
-  ////////////////////////////////////////////////////////
-
   LatticeComplex dS(&Grid);
   dS = zero;
   LatticeComplex dSmom(&Grid);
@@ -229,7 +131,6 @@ int main(int argc, char **argv)
   LatticeComplex dSmom2(&Grid);
   dSmom2 = zero;
 
-
   for (int mu = 0; mu < Nd; mu++)
   {
     mommu = PeekIndex<LorentzIndex>(UdSdU, mu); // P_mu =
@@ -237,7 +138,7 @@ int main(int argc, char **argv)
     PokeIndex<LorentzIndex>(UdSdU, mommu, mu);  // UdSdU_mu = Mom
   }
 
-  std::cout << GridLogMessage << "Antihermiticity tests - 2 " << std::endl;
+  std::cout << GridLogMessage << "Antihermiticity tests" << std::endl;
   for (int mu = 0; mu < Nd; mu++)
   {
     mommu = PeekIndex<LorentzIndex>(mom, mu);
@@ -279,8 +180,8 @@ int main(int argc, char **argv)
 
   std::cout << GridLogMessage << " S      " << S << std::endl;
   std::cout << GridLogMessage << " Sprime " << Sprime << std::endl;
-  std::cout << GridLogMessage << "dS      " << Sprime - S << std::endl;
+  std::cout << GridLogMessage << "dS (S' - S)          :" << Sprime - S << std::endl;
-  std::cout << GridLogMessage << "predict dS    " << dSpred << std::endl;
+  std::cout << GridLogMessage << "predict dS (force)   :" << dSpred << std::endl;
   std::cout << GridLogMessage << "dSm " << dSm << std::endl;
   std::cout << GridLogMessage << "dSm2" << dSm2 << std::endl;