Adding force terms

lib/qcd/QCD.h

@@ -35,7 +35,11 @@ namespace Grid{
 
 namespace QCD {
 
-
+    static const int Xdir = 0;
+    static const int Ydir = 1;
+    static const int Zdir = 2;
+    static const int Tdir = 3;
+    
     static const int Xp = 0;
     static const int Yp = 1;
     static const int Zp = 2;

lib/qcd/action/fermion/WilsonCloverFermion.cc

@@ -76,6 +76,12 @@ namespace QCD {
   void WilsonCloverFermion<Impl>::ImportGauge(const GaugeField& _Umu) {
     this->ImportGauge(_Umu);
     // Compute the field strength terms
+    WilsonLoops<Impl>::FieldStrength(Bx, _Umu, Ydir, Zdir);
+    WilsonLoops<Impl>::FieldStrength(By, _Umu, Zdir, Xdir);
+    WilsonLoops<Impl>::FieldStrength(Bz, _Umu, Xdir, Ydir);
+    WilsonLoops<Impl>::FieldStrength(Ex, _Umu, Tdir, Xdir);
+    WilsonLoops<Impl>::FieldStrength(Ey, _Umu, Tdir, Ydir);
+    WilsonLoops<Impl>::FieldStrength(Ez, _Umu, Tdir, Zdir);
 
     // Invert the Moo, Mee terms (?)
   }
@@ -84,7 +90,7 @@ namespace QCD {
   template<class Impl>
   void WilsonCloverFermion<Impl>::Mooee(const FermionField &in, FermionField &out) {
     out.checkerboard = in.checkerboard;
-
+    assert(0); // to be completed
   }
 
   template<class Impl>
@@ -100,6 +106,27 @@ namespace QCD {
     assert(0); // not implemented yet
   }
 
+  // Derivative parts
+  template<class Impl>
+  void WilsonCloverFermion<Impl>::MDeriv(GaugeField&mat, const FermionField&U, const FermionField&V, int dag){
+    GaugeField tmp(mat._grid);
+    this->DhopDeriv(mat, U, V, dag);
+    MooDeriv(tmp, U, V, dag);
+    mat += tmp;
+  }
+
+ // Derivative parts
+  template<class Impl>
+  void WilsonCloverFermion<Impl>::MooDeriv(GaugeField&mat, const FermionField&U, const FermionField&V, int dag){
+    assert(0); // not implemented yet
+  }
+
+   // Derivative parts
+  template<class Impl>
+  void WilsonCloverFermion<Impl>::MeeDeriv(GaugeField&mat, const FermionField&U, const FermionField&V, int dag){
+    assert(0); // not implemented yet
+  }
+
   FermOpTemplateInstantiate(WilsonCloverFermion);
 
 }

lib/qcd/action/fermion/WilsonCloverFermion.h

@@ -61,7 +61,6 @@ public:
   {
     csw = _csw;
     assert(Nd == 4); // require 4 dimensions
-
   }
 
   virtual RealD M(const FermionField& in, FermionField& out);
@@ -72,6 +71,11 @@ public:
   virtual void MooeeInv(const FermionField &in, FermionField &out);
   virtual void MooeeInvDag(const FermionField &in, FermionField &out);
 
+  virtual void MDeriv(GaugeField&mat, const FermionField&U, const FermionField&V, int dag);
+  virtual void MooDeriv(GaugeField&mat, const FermionField&U, const FermionField&V, int dag);
+  virtual void MeeDeriv(GaugeField&mat, const FermionField&U, const FermionField&V, int dag);
+
+
   void ImportGauge(const GaugeField &_Umu);
 private:
   // here fixing the 4 dimensions, make it more general?

lib/qcd/utils/WilsonLoops.h

@@ -256,6 +256,56 @@ public:
     }
   }
 
+  ////////////////////////////////////////////////////////////////////////
+  // the sum over all staples on each site in direction mu,nu, lower part
+  ////////////////////////////////////////////////////////////////////////
+  static void StapleLower(GaugeMat &staple, const GaugeLorentz &Umu, int mu,
+                          int nu) {
+    if (nu != mu) {
+      GridBase *grid = Umu._grid;
+
+      std::vector<GaugeMat> U(Nd, grid);
+      for (int d = 0; d < Nd; d++) {
+        U[d] = PeekIndex<LorentzIndex>(Umu, d);// some redundant copies
+      }
+
+      // mu
+      // ^
+      // |__>  nu
+
+      //  __
+      // |
+      // |__
+      //
+      //
+      staple = Gimpl::ShiftStaple(
+          Gimpl::CovShiftBackward(U[nu], nu,
+                                  Gimpl::CovShiftBackward(U[mu], mu, U[nu])),
+          mu);
+    }
+  }
+
+  //////////////////////////////////////////////////////
+  //  Field Strength
+  //////////////////////////////////////////////////////
+  static void FieldStrength(GaugeMat &FS, const GaugeLorentz &Umu, int mu, int nu){
+      // Fmn +--<--+  Ut +--<--+
+      //     |     |     |     |
+      //  (x)+-->--+     +-->--+(x)
+      //     |     |     |     |
+      //     +--<--+     +--<--+
+
+      GaugeMat Vup(Umu._grid), Vdn(Umu._grid);
+      StapleUpper(Vup, Umu, mu, nu);// coalesce these two (up low)
+      StapleLower(Vdn, Umu, mu, nu);
+      GaugeMat v = adj(Vup) - adj(Vdn);
+      GaugeMat u = PeekIndex<LorentzIndex>(Umu, mu);  // some redundant copies
+      GaugeMat vu = v*u;
+      FS = 0.25*Ta(u*v - Cshift(vu, mu, +1));
+  }
+
+
+
   //////////////////////////////////////////////////////
   // Similar to above for rectangle is required
   //////////////////////////////////////////////////////