Debugging the Clover term

.gitignore

@@ -123,3 +123,5 @@ make-bin-BUCK.sh
 lib/qcd/spin/gamma-gen/*.h
 lib/qcd/spin/gamma-gen/*.cc
 
+.vscode/settings.json
+settings.json

.vscode/settings.json

@@ -36,6 +36,15 @@
         "tuple": "cpp",
         "type_traits": "cpp",
         "typeinfo": "cpp",
-        "utility": "cpp"
+        "utility": "cpp",
+        "iostream": "cpp",
+        "strstream": "cpp",
+        "complex": "cpp",
+        "fstream": "cpp",
+        "iomanip": "cpp",
+        "istream": "cpp",
+        "ostream": "cpp",
+        "sstream": "cpp",
+        "streambuf": "cpp"
     }
 }

lib/cartesian/Cartesian_base.h

@@ -69,6 +69,8 @@ public:
     std::vector<int> _lstart;     // local start of array in gcoors _processor_coor[d]*_ldimensions[d]
     std::vector<int> _lend  ;     // local end of array in gcoors   _processor_coor[d]*_ldimensions[d]+_ldimensions_[d]-1
 
+    bool _isCheckerBoarded; 
+
 public:
 
     ////////////////////////////////////////////////////////////////

lib/cartesian/Cartesian_full.h

@@ -69,6 +69,7 @@ public:
         ///////////////////////
         // Grid information
         ///////////////////////
+        _isCheckerBoarded = false;
         _ndimension = dimensions.size();
             
         _fdimensions.resize(_ndimension);

lib/cartesian/Cartesian_red_black.h

@@ -139,6 +139,7 @@ public:
     ///////////////////////
     // Grid information
     ///////////////////////
+      _isCheckerBoarded = true;
       _checker_dim = checker_dim;
       assert(checker_dim_mask[checker_dim]==1);
       _ndimension = dimensions.size();

lib/qcd/action/fermion/WilsonCloverFermion.cc

@@ -35,27 +35,6 @@ namespace Grid
 namespace QCD
 {
 
-//WilsonLoop::CloverPlaquette
-/////////////////////////////////////////////////////
-//// Clover plaquette combination in mu,nu plane with Double Stored U
-////////////////////////////////////////////////////
-//static void CloverPlaquette(GaugeMat &Q, const std::vector<GaugeMat> &U,
-//                             const int mu, const int nu){
-//   Q = zero;
-//   Q += Gimpl::CovShiftBackward(
-//    U[mu], mu, Gimpl::CovShiftBackward(
-//             U[nu], nu, Gimpl::CovShiftForward(U[mu], mu, U[nu] )));
-//   Q += Gimpl::CovShiftForward(
-//    U[mu], mu, Gimpl::CovShiftForward(
-//         U[nu], nu, Gimpl::CovShiftBackward(U[mu], mu, U[nu+Nd] )));
-//   Q += Gimpl::CovShiftBackward(
-//    U[nu], nu, Gimpl::CovShiftForward(
-//         U[mu], mu, Gimpl::CovShiftForward(U[nu], nu, U[mu+Nd] )));
-//   Q += Gimpl::CovShiftForward(
-//    U[mu], mu, Gimpl::CovShiftBackward(
-//             U[nu], nu, Gimpl::CovShiftBackward(U[mu], mu, U[nu] )));
-// }
-
 // *NOT* EO
 template <class Impl>
 RealD WilsonCloverFermion<Impl>::M(const FermionField &in, FermionField &out)
@@ -89,7 +68,7 @@ RealD WilsonCloverFermion<Impl>::Mdag(const FermionField &in, FermionField &out)
 template <class Impl>
 void WilsonCloverFermion<Impl>::ImportGauge(const GaugeField &_Umu)
 {
-  this->ImportGauge(_Umu);
+  WilsonFermion<Impl>::ImportGauge(_Umu);
   GridBase *grid = _Umu._grid;
   typename Impl::GaugeLinkField Bx(grid), By(grid), Bz(grid), Ex(grid), Ey(grid), Ez(grid);
 
@@ -102,12 +81,12 @@ void WilsonCloverFermion<Impl>::ImportGauge(const GaugeField &_Umu)
   WilsonLoops<Impl>::FieldStrength(Ez, _Umu, Tdir, Zdir);
 
   // Compute the Clover Operator acting on Colour and Spin
-  CloverTerm  = fillClover(Bx) * (Gamma(Gamma::Algebra::SigmaYZ));
+  CloverTerm  = fillCloverYZ(Bx);
-  CloverTerm += fillClover(By) * (Gamma(Gamma::Algebra::MinusSigmaXZ));
+  CloverTerm += fillCloverXZ(By);
-  CloverTerm += fillClover(Bz) * (Gamma(Gamma::Algebra::SigmaXY));
+  CloverTerm += fillCloverXY(Bz);
-  CloverTerm += fillClover(Ex) * (Gamma(Gamma::Algebra::MinusSigmaXT));
+  CloverTerm += fillCloverXT(Ex);
-  CloverTerm += fillClover(Ey) * (Gamma(Gamma::Algebra::MinusSigmaYT));
+  CloverTerm += fillCloverYT(Ey);
-  CloverTerm += fillClover(Ez) * (Gamma(Gamma::Algebra::MinusSigmaZT));
+  CloverTerm += fillCloverZT(Ez) ;
   CloverTerm *= csw;
 
   int lvol = _Umu._grid->lSites();
@@ -130,8 +109,11 @@ void WilsonCloverFermion<Impl>::ImportGauge(const GaugeField &_Umu)
         for (int a = 0; a < DimRep; a++)
           for (int b = 0; b < DimRep; b++)
             EigenCloverOp(a + j * DimRep, b + k * DimRep) = Qx()(j, k)(a, b);
+    //std::cout << EigenCloverOp << std::endl;
+
 
     EigenInvCloverOp = EigenCloverOp.inverse();
+    //std::cout << EigenInvCloverOp << std::endl;
     for (int j = 0; j < Ns; j++)
       for (int k = 0; k < Ns; k++)
         for (int a = 0; a < DimRep; a++)
@@ -139,17 +121,21 @@ void WilsonCloverFermion<Impl>::ImportGauge(const GaugeField &_Umu)
             Qxinv()(j, k)(a, b) = EigenInvCloverOp(a + j * DimRep, b + k * DimRep);
 
     pokeLocalSite(Qxinv, CloverTermInv, lcoor);
+  }
+
   // Separate the even and odd parts.
   pickCheckerboard(Even, CloverTermEven, CloverTerm);
   pickCheckerboard( Odd,  CloverTermOdd, CloverTerm);
+
   pickCheckerboard(Even, CloverTermInvEven, CloverTermInv);
   pickCheckerboard( Odd,  CloverTermInvOdd, CloverTermInv);
-  }
+
 }
 
 template <class Impl>
 void WilsonCloverFermion<Impl>::Mooee(const FermionField &in, FermionField &out)
 {
+  conformable(in,out);
   this->MooeeInternal(in, out, DaggerNo, InverseNo);
 }
 
@@ -176,15 +162,27 @@ void WilsonCloverFermion<Impl>::MooeeInternal(const FermionField &in, FermionFie
 {
   out.checkerboard = in.checkerboard;
   CloverFieldType *Clover;
-  if (in.checkerboard == Odd){
+  assert(in.checkerboard == Odd || in.checkerboard == Even);
+
+  if (in._grid->_isCheckerBoarded)
+  {
+    if (in.checkerboard == Odd)
+    {
       std::cout << "Calling clover term Odd" << std::endl;
       Clover = (inv) ? &CloverTermInvOdd : &CloverTermOdd;
     }
-  if (in.checkerboard == Even){
+    else
+    {
       std::cout << "Calling clover term Even" << std::endl;
       Clover = (inv) ? &CloverTermInvEven : &CloverTermEven;
     }
+  }
+  else
+  {
+    Clover = (inv) ? &CloverTermInv : &CloverTerm;
+  }
 
+  std::cout << GridLogMessage << "*Clover.checkerboard "  << (*Clover).checkerboard << std::endl;
   if (dag){ out = adj(*Clover) * in;} else { out = *Clover * in;}
 } // MooeeInternal
 

lib/qcd/action/fermion/WilsonCloverFermion.h

@@ -94,17 +94,102 @@ private:
   // eventually these two can be compressed into 6x6 blocks instead of the 12x12
   // using the DeGrand-Rossi basis for the gamma matrices
 
-  CloverFieldType fillClover(const GaugeLinkField& F){
+  CloverFieldType fillCloverYZ(const GaugeLinkField &F)
+  {
     CloverFieldType T(F._grid);
+    T = zero;
     PARALLEL_FOR_LOOP
-    for (int i = 0; i < CloverTerm._grid->oSites(); i++){
+    for (int i = 0; i < CloverTerm._grid->oSites(); i++)
-      for (int s1 = 0; s1 < Nc; s1++)
+    {
-      for (int s2 = 0; s2 < Nc; s2++)
+      T._odata[i]()(0, 1) = timesMinusI(F._odata[i]()());
-      T._odata[i]()(s1,s2) = F._odata[i]()();
+      T._odata[i]()(1, 0) = timesMinusI(F._odata[i]()());
+      T._odata[i]()(2, 3) = timesMinusI(F._odata[i]()());
+      T._odata[i]()(3, 2) = timesMinusI(F._odata[i]()());
     }
+  
   return T;
+}
+
+  CloverFieldType fillCloverXZ(const GaugeLinkField &F)
+  {
+    CloverFieldType T(F._grid);
+    T = zero;
+    PARALLEL_FOR_LOOP
+    for (int i = 0; i < CloverTerm._grid->oSites(); i++)
+    {
+      T._odata[i]()(0, 1) = -F._odata[i]()();
+      T._odata[i]()(1, 0) = F._odata[i]()();
+      T._odata[i]()(2, 3) = -F._odata[i]()();
+      T._odata[i]()(3, 2) = F._odata[i]()();
     }
   
+  return T;
+}
+
+  CloverFieldType fillCloverXY(const GaugeLinkField &F)
+  {
+    CloverFieldType T(F._grid);
+    T = zero;
+    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]()());
+      T._odata[i]()(3, 3) = timesI(F._odata[i]()());
+    }
+  
+  return T;
+}
+
+  CloverFieldType fillCloverXT(const GaugeLinkField &F)
+  {
+    CloverFieldType T(F._grid);
+    T = zero;
+    PARALLEL_FOR_LOOP
+    for (int i = 0; i < CloverTerm._grid->oSites(); i++)
+    {
+      T._odata[i]()(0, 1) = timesMinusI(F._odata[i]()());
+      T._odata[i]()(1, 0) = timesMinusI(F._odata[i]()());
+      T._odata[i]()(2, 3) = timesI(F._odata[i]()());
+      T._odata[i]()(3, 2) = timesI(F._odata[i]()());
+    }
+  
+  return T;
+}
+
+  CloverFieldType fillCloverYT(const GaugeLinkField &F)
+  {
+    CloverFieldType T(F._grid);
+    T = zero;
+    PARALLEL_FOR_LOOP
+    for (int i = 0; i < CloverTerm._grid->oSites(); i++)
+    {
+      T._odata[i]()(0, 1) = (F._odata[i]()());
+      T._odata[i]()(1, 0) = -(F._odata[i]()());
+      T._odata[i]()(2, 3) = -(F._odata[i]()());
+      T._odata[i]()(3, 2) = (F._odata[i]()());
+    }
+  
+  return T;
+}
+
+  CloverFieldType fillCloverZT(const GaugeLinkField &F)
+  {
+    CloverFieldType T(F._grid);
+    T = zero;
+    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) = timesI(F._odata[i]()());
+      T._odata[i]()(3, 3) = timesMinusI(F._odata[i]()());
+    }
+  
+  return T;
+}
+
 };
 }
 }

tests/core/Test_wilson_clover.cc

@@ -91,7 +91,7 @@ int main (int argc, char ** argv)
   }
 
   WilsonCloverFermionR Dwc(Umu,Grid,RBGrid,mass,csw,params);
-  
+  Dwc.ImportGauge(Umu); 
 
   std::cout<<GridLogMessage<<"=========================================================="<<std::endl;
   std::cout<<GridLogMessage<<"= Testing Dhop against cshift implementation         "<<std::endl;
@@ -191,7 +191,6 @@ int main (int argc, char ** argv)
   Dwc.MooeeInv(src_e,phi_e);
 
   Dwc.Mooee(chi_o,src_o);
-exit(1);
   Dwc.MooeeInv(src_o,phi_o);
 
   setCheckerboard(phi,phi_e);