Refactor clover to align with other files and prepare for upcoming changes

Grid/qcd/action/fermion/WilsonCloverFermion.h

@@ -4,10 +4,11 @@
 
     Source file: ./lib/qcd/action/fermion/WilsonCloverFermion.h
 
-    Copyright (C) 2017
+    Copyright (C) 2017 - 2022
 
     Author: Guido Cossu <guido.cossu@ed.ac.uk>
     Author: David Preti <>
+    Author: Daniel Richtmann <daniel.richtmann@gmail.com>
 
     This program is free software; you can redistribute it and/or modify
     it under the terms of the GNU General Public License as published by
@@ -29,7 +30,8 @@
 
 #pragma once
 
-#include <Grid/Grid.h>
+#include <Grid/qcd/action/fermion/WilsonCloverTypes.h>
+#include <Grid/qcd/action/fermion/WilsonCloverHelpers.h>
 
 NAMESPACE_BEGIN(Grid);
 
@@ -50,18 +52,15 @@ NAMESPACE_BEGIN(Grid);
 //////////////////////////////////////////////////////////////////
 
 template <class Impl>
-class WilsonCloverFermion : public WilsonFermion<Impl>
+class WilsonCloverFermion : public WilsonFermion<Impl>,
+                            public WilsonCloverHelpers<Impl>
 {
 public:
-  // Types definitions
   INHERIT_IMPL_TYPES(Impl);
-  template <typename vtype>
+  INHERIT_CLOVER_TYPES(Impl);
-  using iImplClover = iScalar<iMatrix<iMatrix<vtype, Impl::Dimension>, Ns>>;
-  typedef iImplClover<Simd> SiteCloverType;
-  typedef Lattice<SiteCloverType> CloverFieldType;
 
-public:
   typedef WilsonFermion<Impl>       WilsonBase;
+  typedef WilsonCloverHelpers<Impl> Helpers;
 
   virtual int    ConstEE(void)     { return 0; };
   virtual void Instantiatable(void){};
@@ -72,42 +71,7 @@ public:
                       const RealD _csw_r = 0.0,
                       const RealD _csw_t = 0.0,
                       const WilsonAnisotropyCoefficients &clover_anisotropy = WilsonAnisotropyCoefficients(),
-                      const ImplParams &impl_p = ImplParams()) : WilsonFermion<Impl>(_Umu,
+                      const ImplParams &impl_p = ImplParams());
-                                                                                     Fgrid,
-                                                                                     Hgrid,
-                                                                                     _mass, impl_p, clover_anisotropy),
-                                                                 CloverTerm(&Fgrid),
-                                                                 CloverTermInv(&Fgrid),
-                                                                 CloverTermEven(&Hgrid),
-                                                                 CloverTermOdd(&Hgrid),
-                                                                 CloverTermInvEven(&Hgrid),
-                                                                 CloverTermInvOdd(&Hgrid),
-                                                                 CloverTermDagEven(&Hgrid),
-                                                                 CloverTermDagOdd(&Hgrid),
-                                                                 CloverTermInvDagEven(&Hgrid),
-                                                                 CloverTermInvDagOdd(&Hgrid)
-  {
-    assert(Nd == 4); // require 4 dimensions
-
-    if (clover_anisotropy.isAnisotropic)
-    {
-      csw_r = _csw_r * 0.5 / clover_anisotropy.xi_0;
-      diag_mass = _mass + 1.0 + (Nd - 1) * (clover_anisotropy.nu / clover_anisotropy.xi_0);
-    }
-    else
-    {
-      csw_r = _csw_r * 0.5;
-      diag_mass = 4.0 + _mass;
-    }
-    csw_t = _csw_t * 0.5;
-
-    if (csw_r == 0)
-      std::cout << GridLogWarning << "Initializing WilsonCloverFermion with csw_r = 0" << std::endl;
-    if (csw_t == 0)
-      std::cout << GridLogWarning << "Initializing WilsonCloverFermion with csw_t = 0" << std::endl;
-
-    ImportGauge(_Umu);
-  }
 
   virtual void M(const FermionField &in, FermionField &out);
   virtual void Mdag(const FermionField &in, FermionField &out);
@@ -124,250 +88,21 @@ public:
   void ImportGauge(const GaugeField &_Umu);
 
   // Derivative parts unpreconditioned pseudofermions
-  void MDeriv(GaugeField &force, const FermionField &X, const FermionField &Y, int dag)
+  void MDeriv(GaugeField &force, const FermionField &X, const FermionField &Y, int dag);
-  {
-    conformable(X.Grid(), Y.Grid());
-    conformable(X.Grid(), force.Grid());
-    GaugeLinkField force_mu(force.Grid()), lambda(force.Grid());
-    GaugeField clover_force(force.Grid());
-    PropagatorField Lambda(force.Grid());
 
-    // Guido: Here we are hitting some performance issues:
+public:
-    // need to extract the components of the DoubledGaugeField
-    // for each call
-    // Possible solution
-    // Create a vector object to store them? (cons: wasting space)
-    std::vector<GaugeLinkField> U(Nd, this->Umu.Grid());
-
-    Impl::extractLinkField(U, this->Umu);
-
-    force = Zero();
-    // Derivative of the Wilson hopping term
-    this->DhopDeriv(force, X, Y, dag);
-
-    ///////////////////////////////////////////////////////////
-    // Clover term derivative
-    ///////////////////////////////////////////////////////////
-    Impl::outerProductImpl(Lambda, X, Y);
-    //std::cout << "Lambda:" << Lambda << std::endl;
-
-    Gamma::Algebra sigma[] = {
-        Gamma::Algebra::SigmaXY,
-        Gamma::Algebra::SigmaXZ,
-        Gamma::Algebra::SigmaXT,
-        Gamma::Algebra::MinusSigmaXY,
-        Gamma::Algebra::SigmaYZ,
-        Gamma::Algebra::SigmaYT,
-        Gamma::Algebra::MinusSigmaXZ,
-        Gamma::Algebra::MinusSigmaYZ,
-        Gamma::Algebra::SigmaZT,
-        Gamma::Algebra::MinusSigmaXT,
-        Gamma::Algebra::MinusSigmaYT,
-        Gamma::Algebra::MinusSigmaZT};
-
-    /*
-      sigma_{\mu \nu}=
-      | 0         sigma[0]  sigma[1]  sigma[2] |
-      | sigma[3]    0       sigma[4]  sigma[5] |
-      | sigma[6]  sigma[7]     0      sigma[8] |
-      | sigma[9]  sigma[10] sigma[11]   0      |
-    */
-
-    int count = 0;
-    clover_force = Zero();
-    for (int mu = 0; mu < 4; mu++)
-    {
-      force_mu = Zero();
-      for (int nu = 0; nu < 4; nu++)
-      {
-        if (mu == nu)
-        continue;
-        
-        RealD factor;
-        if (nu == 4 || mu == 4)
-        {
-          factor = 2.0 * csw_t;
-        }
-        else
-        {
-          factor = 2.0 * csw_r;
-        }
-        PropagatorField Slambda = Gamma(sigma[count]) * Lambda; // sigma checked
-        Impl::TraceSpinImpl(lambda, Slambda);                   // traceSpin ok
-        force_mu -= factor*Cmunu(U, lambda, mu, nu);                   // checked
-        count++;
-      }
-
-      pokeLorentz(clover_force, U[mu] * force_mu, mu);
-    }
-    //clover_force *= csw;
-    force += clover_force;
-  }
-
-  // Computing C_{\mu \nu}(x) as in Eq.(B.39) in Zbigniew Sroczynski's PhD thesis
-  GaugeLinkField Cmunu(std::vector<GaugeLinkField> &U, GaugeLinkField &lambda, int mu, int nu)
-  {
-    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::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::ShiftStaple(adj(lambda), nu);
-    out += Impl::ShiftStaple(Impl::CovShiftForward(U[nu], nu, Impl::CovShiftBackward(U[mu], mu, Impl::CovShiftIdentityBackward(tmp, nu))), mu);
-
-    // C4+
-    out += Impl::ShiftStaple(Impl::CovShiftForward(U[nu], nu, Impl::CovShiftBackward(U[mu], mu, Impl::CovShiftIdentityBackward(U[nu], nu))), mu) * lambda;
-
-    // insertion in lower staple
-    // C1-
-    out -= Impl::ShiftStaple(lambda, mu) * Impl::ShiftStaple(Impl::CovShiftBackward(U[nu], nu, Impl::CovShiftBackward(U[mu], mu, U[nu])), mu);
-
-    // C2-
-    tmp = adj(lambda) * U[nu];
-    out -= Impl::ShiftStaple(Impl::CovShiftBackward(tmp, nu, Impl::CovShiftBackward(U[mu], mu, U[nu])), mu);
-
-    // C3-
-    tmp = lambda * U[nu];
-    out -= Impl::ShiftStaple(Impl::CovShiftBackward(U[nu], nu, Impl::CovShiftBackward(U[mu], mu, tmp)), mu);
-
-    // C4-
-    out -= Impl::ShiftStaple(Impl::CovShiftBackward(U[nu], nu, Impl::CovShiftBackward(U[mu], mu, U[nu])), mu) * lambda;
-
-    return out;
-  }
-
-protected:
   // here fixing the 4 dimensions, make it more general?
 
   RealD csw_r;                                               // Clover coefficient - spatial
   RealD csw_t;                                               // Clover coefficient - temporal
   RealD diag_mass;                                           // Mass term
-  CloverFieldType CloverTerm, CloverTermInv;                 // Clover term
+  CloverField CloverTerm, CloverTermInv;                     // Clover term
-  CloverFieldType CloverTermEven, CloverTermOdd;             // Clover term EO
+  CloverField CloverTermEven, CloverTermOdd;                 // Clover term EO
-  CloverFieldType CloverTermInvEven, CloverTermInvOdd;       // Clover term Inv EO
+  CloverField CloverTermInvEven, CloverTermInvOdd;           // Clover term Inv EO
-  CloverFieldType CloverTermDagEven, CloverTermDagOdd;       // Clover term Dag EO
+  CloverField CloverTermDagEven, CloverTermDagOdd;           // Clover term Dag EO
-  CloverFieldType CloverTermInvDagEven, CloverTermInvDagOdd; // Clover term Inv Dag EO
+  CloverField CloverTermInvDagEven, CloverTermInvDagOdd;     // Clover term Inv Dag EO
-
- public:
-  // eventually these can be compressed into 6x6 blocks instead of the 12x12
-  // using the DeGrand-Rossi basis for the gamma matrices
-  CloverFieldType fillCloverYZ(const GaugeLinkField &F)
-  {
-    CloverFieldType T(F.Grid());
-    T = Zero();
-    autoView(T_v,T,AcceleratorWrite);
-    autoView(F_v,F,AcceleratorRead);
-    accelerator_for(i, CloverTerm.Grid()->oSites(),1,
-    {
-      T_v[i]()(0, 1) = timesMinusI(F_v[i]()());
-      T_v[i]()(1, 0) = timesMinusI(F_v[i]()());
-      T_v[i]()(2, 3) = timesMinusI(F_v[i]()());
-      T_v[i]()(3, 2) = timesMinusI(F_v[i]()());
-    });
-
-    return T;
-  }
-
-  CloverFieldType fillCloverXZ(const GaugeLinkField &F)
-  {
-    CloverFieldType T(F.Grid());
-    T = Zero();
-    
-    autoView(T_v, T,AcceleratorWrite);
-    autoView(F_v, F,AcceleratorRead);
-    accelerator_for(i, CloverTerm.Grid()->oSites(),1,
-    {
-      T_v[i]()(0, 1) = -F_v[i]()();
-      T_v[i]()(1, 0) = F_v[i]()();
-      T_v[i]()(2, 3) = -F_v[i]()();
-      T_v[i]()(3, 2) = F_v[i]()();
-    });
-
-    return T;
-  }
-
-  CloverFieldType fillCloverXY(const GaugeLinkField &F)
-  {
-    CloverFieldType T(F.Grid());
-    T = Zero();
-
-    autoView(T_v,T,AcceleratorWrite);
-    autoView(F_v,F,AcceleratorRead);
-    accelerator_for(i, CloverTerm.Grid()->oSites(),1,
-    {
-      T_v[i]()(0, 0) = timesMinusI(F_v[i]()());
-      T_v[i]()(1, 1) = timesI(F_v[i]()());
-      T_v[i]()(2, 2) = timesMinusI(F_v[i]()());
-      T_v[i]()(3, 3) = timesI(F_v[i]()());
-    });
-
-    return T;
-  }
-
-  CloverFieldType fillCloverXT(const GaugeLinkField &F)
-  {
-    CloverFieldType T(F.Grid());
-    T = Zero();
-
-    autoView( T_v , T, AcceleratorWrite);
-    autoView( F_v , F, AcceleratorRead);
-    accelerator_for(i, CloverTerm.Grid()->oSites(),1,
-    {
-      T_v[i]()(0, 1) = timesI(F_v[i]()());
-      T_v[i]()(1, 0) = timesI(F_v[i]()());
-      T_v[i]()(2, 3) = timesMinusI(F_v[i]()());
-      T_v[i]()(3, 2) = timesMinusI(F_v[i]()());
-    });
-
-    return T;
-  }
-
-  CloverFieldType fillCloverYT(const GaugeLinkField &F)
-  {
-    CloverFieldType T(F.Grid());
-    T = Zero();
-    
-    autoView( T_v ,T,AcceleratorWrite);
-    autoView( F_v ,F,AcceleratorRead);
-    accelerator_for(i, CloverTerm.Grid()->oSites(),1,
-    {
-      T_v[i]()(0, 1) = -(F_v[i]()());
-      T_v[i]()(1, 0) = (F_v[i]()());
-      T_v[i]()(2, 3) = (F_v[i]()());
-      T_v[i]()(3, 2) = -(F_v[i]()());
-    });
-
-    return T;
-  }
-
-  CloverFieldType fillCloverZT(const GaugeLinkField &F)
-  {
-    CloverFieldType T(F.Grid());
-
-    T = Zero();
-
-    autoView( T_v , T,AcceleratorWrite);
-    autoView( F_v , F,AcceleratorRead);
-    accelerator_for(i, CloverTerm.Grid()->oSites(),1,
-    {
-      T_v[i]()(0, 0) = timesI(F_v[i]()());
-      T_v[i]()(1, 1) = timesMinusI(F_v[i]()());
-      T_v[i]()(2, 2) = timesMinusI(F_v[i]()());
-      T_v[i]()(3, 3) = timesI(F_v[i]()());
-    });
-
-    return T;
-  }
 };
+
 NAMESPACE_END(Grid);
 
 

Grid/qcd/action/fermion/WilsonCloverHelpers.h

@@ -0,0 +1,191 @@
+/*************************************************************************************
+
+    Grid physics library, www.github.com/paboyle/Grid
+
+    Source file: ./lib/qcd/action/fermion/WilsonCloverHelpers.h
+
+    Copyright (C) 2021 - 2022
+
+    Author: Daniel Richtmann <daniel.richtmann@gmail.com>
+
+    This program is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation; either version 2 of the License, or
+    (at your option) any later version.
+
+    This program is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License along
+    with this program; if not, write to the Free Software Foundation, Inc.,
+    51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+
+    See the full license in the file "LICENSE" in the top level distribution directory
+*************************************************************************************/
+/*  END LEGAL */
+
+#pragma once
+
+// Helper routines that implement common clover functionality
+
+NAMESPACE_BEGIN(Grid);
+
+template<class Impl> class WilsonCloverHelpers {
+public:
+
+  INHERIT_IMPL_TYPES(Impl);
+  INHERIT_CLOVER_TYPES(Impl);
+
+  // Computing C_{\mu \nu}(x) as in Eq.(B.39) in Zbigniew Sroczynski's PhD thesis
+  static GaugeLinkField Cmunu(std::vector<GaugeLinkField> &U, GaugeLinkField &lambda, int mu, int nu)
+  {
+    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::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::ShiftStaple(adj(lambda), nu);
+    out += Impl::ShiftStaple(Impl::CovShiftForward(U[nu], nu, Impl::CovShiftBackward(U[mu], mu, Impl::CovShiftIdentityBackward(tmp, nu))), mu);
+
+    // C4+
+    out += Impl::ShiftStaple(Impl::CovShiftForward(U[nu], nu, Impl::CovShiftBackward(U[mu], mu, Impl::CovShiftIdentityBackward(U[nu], nu))), mu) * lambda;
+
+    // insertion in lower staple
+    // C1-
+    out -= Impl::ShiftStaple(lambda, mu) * Impl::ShiftStaple(Impl::CovShiftBackward(U[nu], nu, Impl::CovShiftBackward(U[mu], mu, U[nu])), mu);
+
+    // C2-
+    tmp = adj(lambda) * U[nu];
+    out -= Impl::ShiftStaple(Impl::CovShiftBackward(tmp, nu, Impl::CovShiftBackward(U[mu], mu, U[nu])), mu);
+
+    // C3-
+    tmp = lambda * U[nu];
+    out -= Impl::ShiftStaple(Impl::CovShiftBackward(U[nu], nu, Impl::CovShiftBackward(U[mu], mu, tmp)), mu);
+
+    // C4-
+    out -= Impl::ShiftStaple(Impl::CovShiftBackward(U[nu], nu, Impl::CovShiftBackward(U[mu], mu, U[nu])), mu) * lambda;
+
+    return out;
+  }
+
+  static CloverField fillCloverYZ(const GaugeLinkField &F)
+  {
+    CloverField T(F.Grid());
+    T = Zero();
+    autoView(T_v,T,AcceleratorWrite);
+    autoView(F_v,F,AcceleratorRead);
+    accelerator_for(i, T.Grid()->oSites(),1,
+    {
+      T_v[i]()(0, 1) = timesMinusI(F_v[i]()());
+      T_v[i]()(1, 0) = timesMinusI(F_v[i]()());
+      T_v[i]()(2, 3) = timesMinusI(F_v[i]()());
+      T_v[i]()(3, 2) = timesMinusI(F_v[i]()());
+    });
+
+    return T;
+  }
+
+  static CloverField fillCloverXZ(const GaugeLinkField &F)
+  {
+    CloverField T(F.Grid());
+    T = Zero();
+    
+    autoView(T_v, T,AcceleratorWrite);
+    autoView(F_v, F,AcceleratorRead);
+    accelerator_for(i, T.Grid()->oSites(),1,
+    {
+      T_v[i]()(0, 1) = -F_v[i]()();
+      T_v[i]()(1, 0) = F_v[i]()();
+      T_v[i]()(2, 3) = -F_v[i]()();
+      T_v[i]()(3, 2) = F_v[i]()();
+    });
+
+    return T;
+  }
+
+  static CloverField fillCloverXY(const GaugeLinkField &F)
+  {
+    CloverField T(F.Grid());
+    T = Zero();
+
+    autoView(T_v,T,AcceleratorWrite);
+    autoView(F_v,F,AcceleratorRead);
+    accelerator_for(i, T.Grid()->oSites(),1,
+    {
+      T_v[i]()(0, 0) = timesMinusI(F_v[i]()());
+      T_v[i]()(1, 1) = timesI(F_v[i]()());
+      T_v[i]()(2, 2) = timesMinusI(F_v[i]()());
+      T_v[i]()(3, 3) = timesI(F_v[i]()());
+    });
+
+    return T;
+  }
+
+  static CloverField fillCloverXT(const GaugeLinkField &F)
+  {
+    CloverField T(F.Grid());
+    T = Zero();
+
+    autoView( T_v , T, AcceleratorWrite);
+    autoView( F_v , F, AcceleratorRead);
+    accelerator_for(i, T.Grid()->oSites(),1,
+    {
+      T_v[i]()(0, 1) = timesI(F_v[i]()());
+      T_v[i]()(1, 0) = timesI(F_v[i]()());
+      T_v[i]()(2, 3) = timesMinusI(F_v[i]()());
+      T_v[i]()(3, 2) = timesMinusI(F_v[i]()());
+    });
+
+    return T;
+  }
+
+  static CloverField fillCloverYT(const GaugeLinkField &F)
+  {
+    CloverField T(F.Grid());
+    T = Zero();
+    
+    autoView( T_v ,T,AcceleratorWrite);
+    autoView( F_v ,F,AcceleratorRead);
+    accelerator_for(i, T.Grid()->oSites(),1,
+    {
+      T_v[i]()(0, 1) = -(F_v[i]()());
+      T_v[i]()(1, 0) = (F_v[i]()());
+      T_v[i]()(2, 3) = (F_v[i]()());
+      T_v[i]()(3, 2) = -(F_v[i]()());
+    });
+
+    return T;
+  }
+
+  static CloverField fillCloverZT(const GaugeLinkField &F)
+  {
+    CloverField T(F.Grid());
+
+    T = Zero();
+
+    autoView( T_v , T,AcceleratorWrite);
+    autoView( F_v , F,AcceleratorRead);
+    accelerator_for(i, T.Grid()->oSites(),1,
+    {
+      T_v[i]()(0, 0) = timesI(F_v[i]()());
+      T_v[i]()(1, 1) = timesMinusI(F_v[i]()());
+      T_v[i]()(2, 2) = timesMinusI(F_v[i]()());
+      T_v[i]()(3, 3) = timesI(F_v[i]()());
+    });
+
+    return T;
+  }
+};
+
+NAMESPACE_END(Grid);

Grid/qcd/action/fermion/WilsonCloverTypes.h

@@ -0,0 +1,49 @@
+/*************************************************************************************
+
+    Grid physics library, www.github.com/paboyle/Grid
+
+    Source file: ./lib/qcd/action/fermion/WilsonCloverTypes.h
+
+    Copyright (C) 2021 - 2022
+
+    Author: Daniel Richtmann <daniel.richtmann@gmail.com>
+
+    This program is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation; either version 2 of the License, or
+    (at your option) any later version.
+
+    This program is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License along
+    with this program; if not, write to the Free Software Foundation, Inc.,
+    51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+
+    See the full license in the file "LICENSE" in the top level distribution directory
+*************************************************************************************/
+/*  END LEGAL */
+
+#pragma once
+
+NAMESPACE_BEGIN(Grid);
+
+template<class Impl>
+class WilsonCloverTypes {
+public:
+  INHERIT_IMPL_TYPES(Impl);
+
+  template <typename vtype> using iImplClover = iScalar<iMatrix<iMatrix<vtype, Impl::Dimension>, Ns>>;
+
+  typedef iImplClover<Simd> SiteClover;
+
+  typedef Lattice<SiteClover> CloverField;
+};
+
+#define INHERIT_CLOVER_TYPES(Impl)                                 \
+  typedef typename WilsonCloverTypes<Impl>::SiteClover SiteClover; \
+  typedef typename WilsonCloverTypes<Impl>::CloverField CloverField;
+
+NAMESPACE_END(Grid);

Grid/qcd/action/fermion/implementation/WilsonCloverFermionImplementation.h

@@ -2,12 +2,13 @@
 
     Grid physics library, www.github.com/paboyle/Grid
 
-    Source file: ./lib/qcd/action/fermion/WilsonCloverFermion.cc
+    Source file: ./lib/qcd/action/fermion/WilsonCloverFermionImplementation.h
 
-    Copyright (C) 2017
+    Copyright (C) 2017 - 2022
 
     Author: paboyle <paboyle@ph.ed.ac.uk>
     Author: Guido Cossu <guido.cossu@ed.ac.uk>
+    Author: Daniel Richtmann <daniel.richtmann@gmail.com>
 
     This program is free software; you can redistribute it and/or modify
     it under the terms of the GNU General Public License as published by
@@ -33,6 +34,45 @@
 
 NAMESPACE_BEGIN(Grid);
 
+template<class Impl>
+WilsonCloverFermion<Impl>::WilsonCloverFermion(GaugeField&                         _Umu,
+                                               GridCartesian&                      Fgrid,
+                                               GridRedBlackCartesian&              Hgrid,
+                                               const RealD                         _mass,
+                                               const RealD                         _csw_r,
+                                               const RealD                         _csw_t,
+                                               const WilsonAnisotropyCoefficients& clover_anisotropy,
+                                               const ImplParams&                   impl_p)
+  : WilsonFermion<Impl>(_Umu, Fgrid, Hgrid, _mass, impl_p, clover_anisotropy)
+  , CloverTerm(&Fgrid)
+  , CloverTermInv(&Fgrid)
+  , CloverTermEven(&Hgrid)
+  , CloverTermOdd(&Hgrid)
+  , CloverTermInvEven(&Hgrid)
+  , CloverTermInvOdd(&Hgrid)
+  , CloverTermDagEven(&Hgrid)
+  , CloverTermDagOdd(&Hgrid)
+  , CloverTermInvDagEven(&Hgrid)
+  , CloverTermInvDagOdd(&Hgrid) {
+  assert(Nd == 4); // require 4 dimensions
+
+  if(clover_anisotropy.isAnisotropic) {
+    csw_r     = _csw_r * 0.5 / clover_anisotropy.xi_0;
+    diag_mass = _mass + 1.0 + (Nd - 1) * (clover_anisotropy.nu / clover_anisotropy.xi_0);
+  } else {
+    csw_r     = _csw_r * 0.5;
+    diag_mass = 4.0 + _mass;
+  }
+  csw_t = _csw_t * 0.5;
+
+  if(csw_r == 0)
+    std::cout << GridLogWarning << "Initializing WilsonCloverFermion with csw_r = 0" << std::endl;
+  if(csw_t == 0)
+    std::cout << GridLogWarning << "Initializing WilsonCloverFermion with csw_t = 0" << std::endl;
+
+  ImportGauge(_Umu);
+}
+
 // *NOT* EO
 template <class Impl>
 void WilsonCloverFermion<Impl>::M(const FermionField &in, FermionField &out)
@@ -67,10 +107,13 @@ void WilsonCloverFermion<Impl>::Mdag(const FermionField &in, FermionField &out)
 template <class Impl>
 void WilsonCloverFermion<Impl>::ImportGauge(const GaugeField &_Umu)
 {
+  double t0 = usecond();
   WilsonFermion<Impl>::ImportGauge(_Umu);
+  double t1 = usecond();
   GridBase *grid = _Umu.Grid();
   typename Impl::GaugeLinkField Bx(grid), By(grid), Bz(grid), Ex(grid), Ey(grid), Ez(grid);
 
+  double t2 = usecond();
   // Compute the field strength terms mu>nu
   WilsonLoops<Impl>::FieldStrength(Bx, _Umu, Zdir, Ydir);
   WilsonLoops<Impl>::FieldStrength(By, _Umu, Zdir, Xdir);
@@ -79,19 +122,22 @@ void WilsonCloverFermion<Impl>::ImportGauge(const GaugeField &_Umu)
   WilsonLoops<Impl>::FieldStrength(Ey, _Umu, Tdir, Ydir);
   WilsonLoops<Impl>::FieldStrength(Ez, _Umu, Tdir, Zdir);
 
+  double t3 = usecond();
   // Compute the Clover Operator acting on Colour and Spin
   // multiply here by the clover coefficients for the anisotropy
-  CloverTerm  = fillCloverYZ(Bx) * csw_r;
+  CloverTerm  = Helpers::fillCloverYZ(Bx) * csw_r;
-  CloverTerm += fillCloverXZ(By) * csw_r;
+  CloverTerm += Helpers::fillCloverXZ(By) * csw_r;
-  CloverTerm += fillCloverXY(Bz) * csw_r;
+  CloverTerm += Helpers::fillCloverXY(Bz) * csw_r;
-  CloverTerm += fillCloverXT(Ex) * csw_t;
+  CloverTerm += Helpers::fillCloverXT(Ex) * csw_t;
-  CloverTerm += fillCloverYT(Ey) * csw_t;
+  CloverTerm += Helpers::fillCloverYT(Ey) * csw_t;
-  CloverTerm += fillCloverZT(Ez) * csw_t;
+  CloverTerm += Helpers::fillCloverZT(Ez) * csw_t;
   CloverTerm += diag_mass;
 
+  double t4 = usecond();
   int lvol = _Umu.Grid()->lSites();
   int DimRep = Impl::Dimension;
 
+  double t5 = usecond();
   {
     autoView(CTv,CloverTerm,CpuRead);
     autoView(CTIv,CloverTermInv,CpuWrite);
@@ -100,7 +146,7 @@ void WilsonCloverFermion<Impl>::ImportGauge(const GaugeField &_Umu)
       grid->LocalIndexToLocalCoor(site, lcoor);
       Eigen::MatrixXcd EigenCloverOp = Eigen::MatrixXcd::Zero(Ns * DimRep, Ns * DimRep);
       Eigen::MatrixXcd EigenInvCloverOp = Eigen::MatrixXcd::Zero(Ns * DimRep, Ns * DimRep);
-      typename SiteCloverType::scalar_object Qx = Zero(), Qxinv = Zero();
+      typename SiteClover::scalar_object Qx = Zero(), Qxinv = Zero();
       peekLocalSite(Qx, CTv, lcoor);
       //if (csw!=0){
       for (int j = 0; j < Ns; j++)
@@ -125,6 +171,7 @@ void WilsonCloverFermion<Impl>::ImportGauge(const GaugeField &_Umu)
     });
   }
 
+  double t6 = usecond();
   // Separate the even and odd parts
   pickCheckerboard(Even, CloverTermEven, CloverTerm);
   pickCheckerboard(Odd, CloverTermOdd, CloverTerm);
@@ -137,6 +184,20 @@ void WilsonCloverFermion<Impl>::ImportGauge(const GaugeField &_Umu)
 
   pickCheckerboard(Even, CloverTermInvDagEven, adj(CloverTermInv));
   pickCheckerboard(Odd, CloverTermInvDagOdd, adj(CloverTermInv));
+  double t7 = usecond();
+
+#if 0
+  std::cout << GridLogMessage << "WilsonCloverFermion::ImportGauge timings:"
+            << " WilsonFermion::Importgauge = " << (t1 - t0) / 1e6
+            << ", allocations = "               << (t2 - t1) / 1e6
+            << ", field strength = "            << (t3 - t2) / 1e6
+            << ", fill clover = "               << (t4 - t3) / 1e6
+            << ", misc = "                      << (t5 - t4) / 1e6
+            << ", inversions = "                << (t6 - t5) / 1e6
+            << ", pick cbs = "                  << (t7 - t6) / 1e6
+            << ", total = "                     << (t7 - t0) / 1e6
+            << std::endl;
+#endif
 }
 
 template <class Impl>
@@ -167,7 +228,7 @@ template <class Impl>
 void WilsonCloverFermion<Impl>::MooeeInternal(const FermionField &in, FermionField &out, int dag, int inv)
 {
   out.Checkerboard() = in.Checkerboard();
-  CloverFieldType *Clover;
+  CloverField *Clover;
   assert(in.Checkerboard() == Odd || in.Checkerboard() == Even);
 
   if (dag)
@@ -214,9 +275,89 @@ void WilsonCloverFermion<Impl>::MooeeInternal(const FermionField &in, FermionFie
       out = *Clover * in;
     }
   }
-
 } // MooeeInternal
 
+// Derivative parts unpreconditioned pseudofermions
+template <class Impl>
+void WilsonCloverFermion<Impl>::MDeriv(GaugeField &force, const FermionField &X, const FermionField &Y, int dag)
+{
+  conformable(X.Grid(), Y.Grid());
+  conformable(X.Grid(), force.Grid());
+  GaugeLinkField force_mu(force.Grid()), lambda(force.Grid());
+  GaugeField clover_force(force.Grid());
+  PropagatorField Lambda(force.Grid());
+
+  // Guido: Here we are hitting some performance issues:
+  // need to extract the components of the DoubledGaugeField
+  // for each call
+  // Possible solution
+  // Create a vector object to store them? (cons: wasting space)
+  std::vector<GaugeLinkField> U(Nd, this->Umu.Grid());
+
+  Impl::extractLinkField(U, this->Umu);
+
+  force = Zero();
+  // Derivative of the Wilson hopping term
+  this->DhopDeriv(force, X, Y, dag);
+
+  ///////////////////////////////////////////////////////////
+  // Clover term derivative
+  ///////////////////////////////////////////////////////////
+  Impl::outerProductImpl(Lambda, X, Y);
+  //std::cout << "Lambda:" << Lambda << std::endl;
+
+  Gamma::Algebra sigma[] = {
+      Gamma::Algebra::SigmaXY,
+      Gamma::Algebra::SigmaXZ,
+      Gamma::Algebra::SigmaXT,
+      Gamma::Algebra::MinusSigmaXY,
+      Gamma::Algebra::SigmaYZ,
+      Gamma::Algebra::SigmaYT,
+      Gamma::Algebra::MinusSigmaXZ,
+      Gamma::Algebra::MinusSigmaYZ,
+      Gamma::Algebra::SigmaZT,
+      Gamma::Algebra::MinusSigmaXT,
+      Gamma::Algebra::MinusSigmaYT,
+      Gamma::Algebra::MinusSigmaZT};
+
+  /*
+    sigma_{\mu \nu}=
+    | 0         sigma[0]  sigma[1]  sigma[2] |
+    | sigma[3]    0       sigma[4]  sigma[5] |
+    | sigma[6]  sigma[7]     0      sigma[8] |
+    | sigma[9]  sigma[10] sigma[11]   0      |
+  */
+
+  int count = 0;
+  clover_force = Zero();
+  for (int mu = 0; mu < 4; mu++)
+  {
+    force_mu = Zero();
+    for (int nu = 0; nu < 4; nu++)
+    {
+      if (mu == nu)
+      continue;
+
+      RealD factor;
+      if (nu == 4 || mu == 4)
+      {
+        factor = 2.0 * csw_t;
+      }
+      else
+      {
+        factor = 2.0 * csw_r;
+      }
+      PropagatorField Slambda = Gamma(sigma[count]) * Lambda; // sigma checked
+      Impl::TraceSpinImpl(lambda, Slambda);                   // traceSpin ok
+      force_mu -= factor*Helpers::Cmunu(U, lambda, mu, nu);                   // checked
+      count++;
+    }
+
+    pokeLorentz(clover_force, U[mu] * force_mu, mu);
+  }
+  //clover_force *= csw;
+  force += clover_force;
+}
 
 // Derivative parts
 template <class Impl>