Merge branch 'develop' into feature/hadrons

# Conflicts:
#	extras/Hadrons/Modules.hpp
#	extras/Hadrons/modules.inc

.gitignore

@@ -123,3 +123,9 @@ make-bin-BUCK.sh
 #####################
 lib/qcd/spin/gamma-gen/*.h
 lib/qcd/spin/gamma-gen/*.cc
+
+# vs code editor files #
+########################
+.vscode/
+.vscode/settings.json
+settings.json

.travis.yml

@@ -44,3 +44,4 @@ script:
     - make -j4
     - ./benchmarks/Benchmark_dwf --threads 1 --debug-signals
     - make check
+

extras/Hadrons/Modules.hpp

@@ -47,6 +47,7 @@ See the full license in the file "LICENSE" in the top level distribution directo
 #include <Grid/Hadrons/Modules/MGauge/Unit.hpp>
 #include <Grid/Hadrons/Modules/MGauge/Random.hpp>
 #include <Grid/Hadrons/Modules/MGauge/StochEm.hpp>
+#include <Grid/Hadrons/Modules/MGauge/FundtoHirep.hpp>
 #include <Grid/Hadrons/Modules/MUtilities/TestSeqGamma.hpp>
 #include <Grid/Hadrons/Modules/MUtilities/TestSeqConserved.hpp>
 #include <Grid/Hadrons/Modules/MLoop/NoiseLoop.hpp>
@@ -55,6 +56,7 @@ See the full license in the file "LICENSE" in the top level distribution directo
 #include <Grid/Hadrons/Modules/MScalar/ChargedProp.hpp>
 #include <Grid/Hadrons/Modules/MAction/DWF.hpp>
 #include <Grid/Hadrons/Modules/MAction/Wilson.hpp>
+#include <Grid/Hadrons/Modules/MAction/WilsonClover.hpp>
 #include <Grid/Hadrons/Modules/MScalarSUN/Div.hpp>
 #include <Grid/Hadrons/Modules/MScalarSUN/TrMag.hpp>
 #include <Grid/Hadrons/Modules/MScalarSUN/TwoPoint.hpp>

extras/Hadrons/Modules/MAction/WilsonClover.hpp

@@ -0,0 +1,153 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: extras/Hadrons/Modules/MAction/Wilson.hpp
+
+Copyright (C) 2015
+Copyright (C) 2016
+
+Author: Antonin Portelli <antonin.portelli@me.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 */
+
+#ifndef Hadrons_MAction_WilsonClover_hpp_
+#define Hadrons_MAction_WilsonClover_hpp_
+
+#include <Grid/Hadrons/Global.hpp>
+#include <Grid/Hadrons/Module.hpp>
+#include <Grid/Hadrons/ModuleFactory.hpp>
+
+BEGIN_HADRONS_NAMESPACE
+
+/******************************************************************************
+ *                            TWilson quark action                            *
+ ******************************************************************************/
+BEGIN_MODULE_NAMESPACE(MAction)
+
+class WilsonCloverPar: Serializable
+{
+public:
+    GRID_SERIALIZABLE_CLASS_MEMBERS(WilsonCloverPar,
+                                    std::string, gauge,
+                                    double     , mass,
+				                    double     , csw_r,
+				                    double     , csw_t,
+				                    WilsonAnisotropyCoefficients ,clover_anisotropy,
+                                    std::string, boundary
+				    );
+};
+
+template <typename FImpl>
+class TWilsonClover: public Module<WilsonCloverPar>
+{
+public:
+    FGS_TYPE_ALIASES(FImpl,);
+public:
+    // constructor
+    TWilsonClover(const std::string name);
+    // destructor
+    virtual ~TWilsonClover(void) = default;
+    // dependencies/products
+    virtual std::vector<std::string> getInput(void);
+    virtual std::vector<std::string> getOutput(void);
+    // setup
+    virtual void setup(void);
+    // execution
+    virtual void execute(void);
+};
+
+MODULE_REGISTER_NS(WilsonClover, TWilsonClover<FIMPL>, MAction);
+
+/******************************************************************************
+ *                     TWilsonClover template implementation                        *
+ ******************************************************************************/
+// constructor /////////////////////////////////////////////////////////////////
+template <typename FImpl>
+TWilsonClover<FImpl>::TWilsonClover(const std::string name)
+: Module<WilsonCloverPar>(name)
+{}
+
+// dependencies/products ///////////////////////////////////////////////////////
+template <typename FImpl>
+std::vector<std::string> TWilsonClover<FImpl>::getInput(void)
+{
+    std::vector<std::string> in = {par().gauge};
+
+    return in;
+}
+
+template <typename FImpl>
+std::vector<std::string> TWilsonClover<FImpl>::getOutput(void)
+{
+    std::vector<std::string> out = {getName()};
+
+    return out;
+}
+
+// setup ///////////////////////////////////////////////////////////////////////
+template <typename FImpl>
+void TWilsonClover<FImpl>::setup(void)
+{
+    //unsigned int size;
+
+    // size = 2*env().template lattice4dSize<typename FImpl::DoubledGaugeField>();
+    // env().registerObject(getName(), size);
+
+
+    LOG(Message) << "Setting up TWilsonClover fermion matrix with m= " << par().mass
+                 << " using gauge field '" << par().gauge << "'" << std::endl;
+    LOG(Message) << "Fermion boundary conditions: " << par().boundary 
+                 << std::endl;
+    LOG(Message) << "Clover term csw_r: " << par().csw_r
+                 << " csw_t: " << par().csw_t
+                 << std::endl;
+    auto &U      = envGet(LatticeGaugeField, par().gauge);
+    auto &grid   = *env().getGrid();
+    auto &gridRb = *env().getRbGrid();
+    std::vector<Complex> boundary = strToVec<Complex>(par().boundary);
+    typename WilsonCloverFermion<FImpl>::ImplParams implParams(boundary);
+    envCreateDerived(FMat, WilsonCloverFermion<FImpl>, getName(), 1, U, grid, gridRb, par().mass,
+						  par().csw_r,
+						  par().csw_t,
+					      par().clover_anisotropy,
+						  implParams); 
+
+
+    //FMat *fMatPt = new WilsonCloverFermion<FImpl>(U, grid, gridRb, par().mass,
+	//					  par().csw_r,
+	//					  par().csw_t,
+	//				      par().clover_anisotropy,
+	//					  implParams);
+    //env().setObject(getName(), fMatPt);
+
+}
+
+// execution ///////////////////////////////////////////////////////////////////
+template <typename FImpl>
+void TWilsonClover<FImpl>::execute()
+{
+
+}
+
+END_MODULE_NAMESPACE
+
+END_HADRONS_NAMESPACE
+
+#endif // Hadrons_WilsonClover_hpp_

extras/Hadrons/Modules/MFermion/GaugeProp.hpp

@@ -94,7 +94,6 @@ private:
 };
 
 MODULE_REGISTER_NS(GaugeProp, TGaugeProp<FIMPL>, MFermion);
-
 /******************************************************************************
  *                      TGaugeProp implementation                             *
  ******************************************************************************/
@@ -154,7 +153,7 @@ void TGaugeProp<FImpl>::execute(void)
     LOG(Message) << "Inverting using solver '" << par().solver
                  << "' on source '" << par().source << "'" << std::endl;
     for (unsigned int s = 0; s < Ns; ++s)
-    for (unsigned int c = 0; c < Nc; ++c)
+      for (unsigned int c = 0; c < FImpl::Dimension; ++c)
     {
         LOG(Message) << "Inversion for spin= " << s << ", color= " << c
                      << std::endl;
@@ -163,11 +162,11 @@ void TGaugeProp<FImpl>::execute(void)
         {
             if (Ls_ == 1)
             {
-                PropToFerm(source, fullSrc, s, c);
+               PropToFerm<FImpl>(source, fullSrc, s, c);
             }
             else
             {
-                PropToFerm(tmp, fullSrc, s, c);
+                PropToFerm<FImpl>(tmp, fullSrc, s, c);
                 make_5D(tmp, source, Ls_);
             }
         }
@@ -180,18 +179,18 @@ void TGaugeProp<FImpl>::execute(void)
             }
             else
             {
-                PropToFerm(source, fullSrc, s, c);
+                PropToFerm<FImpl>(source, fullSrc, s, c);
             }
         }
         sol = zero;
         solver(sol, source);
-        FermToProp(prop, sol, s, c);
+        FermToProp<FImpl>(prop, sol, s, c);
         // create 4D propagators from 5D one if necessary
         if (Ls_ > 1)
         {
             PropagatorField &p4d = envGet(PropagatorField, getName());
             make_4D(sol, tmp, Ls_);
-            FermToProp(p4d, tmp, s, c);
+            FermToProp<FImpl>(p4d, tmp, s, c);
         }
     }
 }

extras/Hadrons/Modules/MGauge/FundtoHirep.cc

@@ -0,0 +1,75 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: extras/Hadrons/Modules/MGauge/FundtoHirep.cc
+
+Copyright (C) 2015
+Copyright (C) 2016
+
+
+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 */
+
+#include <Grid/Hadrons/Modules/MGauge/FundtoHirep.hpp>
+
+using namespace Grid;
+using namespace Hadrons;
+using namespace MGauge;
+
+// constructor /////////////////////////////////////////////////////////////////
+template <class Rep>
+TFundtoHirep<Rep>::TFundtoHirep(const std::string name)
+: Module<FundtoHirepPar>(name)
+{}
+
+// dependencies/products ///////////////////////////////////////////////////////
+template <class Rep>
+std::vector<std::string> TFundtoHirep<Rep>::getInput(void)
+{
+    std::vector<std::string> in;
+    return in;
+}
+
+template <class Rep>
+std::vector<std::string> TFundtoHirep<Rep>::getOutput(void)
+{
+    std::vector<std::string> out = {getName()};
+    return out;
+}
+
+// setup ///////////////////////////////////////////////////////////////////////
+template <typename Rep>
+void TFundtoHirep<Rep>::setup(void)
+{
+    env().template registerLattice<typename Rep::LatticeField>(getName());
+}
+
+// execution ///////////////////////////////////////////////////////////////////
+template <class Rep>
+void TFundtoHirep<Rep>::execute(void)
+{
+    auto &U      = *env().template getObject<LatticeGaugeField>(par().gaugeconf);
+    LOG(Message) << "Transforming Representation" << std::endl;
+
+    Rep TargetRepresentation(U._grid);
+    TargetRepresentation.update_representation(U);
+
+   typename Rep::LatticeField &URep = *env().template createLattice<typename Rep::LatticeField>(getName());
+    URep = TargetRepresentation.U;
+}

extras/Hadrons/Modules/MGauge/FundtoHirep.hpp

@@ -0,0 +1,77 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: extras/Hadrons/Modules/MGauge/FundtoHirep.hpp
+
+Copyright (C) 2015
+Copyright (C) 2016
+
+Author: David Preti <david.preti@to.infn.it>
+	Guido Cossu <guido.cossu@ed.ac.uk>
+
+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 */
+
+#ifndef Hadrons_MGauge_FundtoHirep_hpp_
+#define Hadrons_MGauge_FundtoHirep_hpp_
+
+#include <Grid/Hadrons/Global.hpp>
+#include <Grid/Hadrons/Module.hpp>
+#include <Grid/Hadrons/ModuleFactory.hpp>
+
+BEGIN_HADRONS_NAMESPACE
+
+/******************************************************************************
+ *                         Load a NERSC configuration                         *
+ ******************************************************************************/
+BEGIN_MODULE_NAMESPACE(MGauge)
+
+class FundtoHirepPar: Serializable
+{
+public:
+    GRID_SERIALIZABLE_CLASS_MEMBERS(FundtoHirepPar,
+                                    std::string, gaugeconf);
+};
+
+template <class Rep>
+class TFundtoHirep: public Module<FundtoHirepPar>
+{
+public:
+    // constructor
+    TFundtoHirep(const std::string name);
+    // destructor
+    virtual ~TFundtoHirep(void) = default;
+    // dependency relation
+    virtual std::vector<std::string> getInput(void);
+    virtual std::vector<std::string> getOutput(void);
+    // setup
+    void setup(void);
+    // execution
+    void execute(void);
+};
+
+//MODULE_REGISTER_NS(FundtoAdjoint,   TFundtoHirep<AdjointRepresentation>, MGauge);
+//MODULE_REGISTER_NS(FundtoTwoIndexSym, TFundtoHirep<TwoIndexSymmetricRepresentation>, MGauge);
+//MODULE_REGISTER_NS(FundtoTwoIndexAsym, TFundtoHirep<TwoIndexAntiSymmetricRepresentation>, MGauge);
+
+END_MODULE_NAMESPACE
+
+END_HADRONS_NAMESPACE
+
+#endif // Hadrons_MGauge_FundtoHirep_hpp_

extras/Hadrons/modules.inc

@@ -1,10 +1,13 @@
 modules_cc =\
+  Modules/MScalar/ChargedProp.cc \
+  Modules/MScalar/FreeProp.cc \
   Modules/MContraction/WeakHamiltonianEye.cc \
   Modules/MContraction/WeakNeutral4ptDisc.cc \
   Modules/MContraction/WeakHamiltonianNonEye.cc \
   Modules/MGauge/Unit.cc \
   Modules/MGauge/StochEm.cc \
   Modules/MGauge/Random.cc \
+  Modules/MGauge/FundtoHirep.cc \
   Modules/MScalar/FreeProp.cc \
   Modules/MScalar/ChargedProp.cc \
   Modules/MIO/LoadNersc.cc
@@ -31,6 +34,7 @@ modules_hpp =\
   Modules/MGauge/Unit.hpp \
   Modules/MGauge/Random.hpp \
   Modules/MGauge/StochEm.hpp \
+  Modules/MGauge/FundtoHirep.hpp \
   Modules/MUtilities/TestSeqGamma.hpp \
   Modules/MUtilities/TestSeqConserved.hpp \
   Modules/MLoop/NoiseLoop.hpp \
@@ -39,6 +43,7 @@ modules_hpp =\
   Modules/MScalar/ChargedProp.hpp \
   Modules/MAction/DWF.hpp \
   Modules/MAction/Wilson.hpp \
+  Modules/MAction/WilsonClover.hpp \
   Modules/MScalarSUN/Div.hpp \
   Modules/MScalarSUN/TrMag.hpp \
   Modules/MScalarSUN/TwoPoint.hpp \

lib/algorithms/LinearOperator.h

@@ -184,10 +184,12 @@ namespace Grid {
       virtual  RealD MpcDag   (const Field &in, Field &out) =0;
       virtual void MpcDagMpc(const Field &in, Field &out,RealD &ni,RealD &no) {
       Field tmp(in._grid);
+      tmp.checkerboard = in.checkerboard;
 	ni=Mpc(in,tmp);
 	no=MpcDag(tmp,out);
       }
       virtual void HermOpAndNorm(const Field &in, Field &out,RealD &n1,RealD &n2){
+      out.checkerboard = in.checkerboard;
 	MpcDagMpc(in,out,n1,n2);
       }
       virtual void HermOp(const Field &in, Field &out){
@@ -216,12 +218,14 @@ namespace Grid {
       SchurDiagMooeeOperator (Matrix &Mat): _Mat(Mat){};
       virtual  RealD Mpc      (const Field &in, Field &out) {
       Field tmp(in._grid);
+      tmp.checkerboard = !in.checkerboard;
 	//std::cout <<"grid pointers: in._grid="<< in._grid << " out._grid=" << out._grid << "  _Mat.Grid=" << _Mat.Grid() << " _Mat.RedBlackGrid=" << _Mat.RedBlackGrid() << std::endl;
 
 	_Mat.Meooe(in,tmp);
 	_Mat.MooeeInv(tmp,out);
 	_Mat.Meooe(out,tmp);
 
+      //std::cout << "cb in " << in.checkerboard << "  cb out " << out.checkerboard << std::endl;
 	_Mat.Mooee(in,out);
 	return axpy_norm(out,-1.0,tmp,out);
       }

lib/cartesian/Cartesian_base.h

@@ -79,6 +79,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

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

lib/cartesian/Cartesian_red_black.h

@@ -171,9 +171,8 @@ public:
               const std::vector<int> &checker_dim_mask,
               int checker_dim)
     {
-      ///////////////////////
+
-      // Grid information
+      _isCheckerBoarded = true;
-      ///////////////////////
       _checker_dim = checker_dim;
       assert(checker_dim_mask[checker_dim] == 1);
       _ndimension = dimensions.size();

lib/qcd/QCD.h

@@ -39,6 +39,7 @@ namespace QCD {
     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;
@@ -420,15 +421,16 @@ namespace QCD {
     //////////////////////////////////////////////
     // Fermion <-> propagator assignements
     //////////////////////////////////////////////
-    template <class Prop, class Ferm>
+    //template <class Prop, class Ferm>
-    void FermToProp(Prop &p, const Ferm &f, const int s, const int c)
+    template <class Fimpl>
+      void FermToProp(typename Fimpl::PropagatorField &p, const typename Fimpl::FermionField &f, const int s, const int c)
     {
       for(int j = 0; j < Ns; ++j)
         {
             auto pjs = peekSpin(p, j, s);
             auto fj  = peekSpin(f, j);
             
-            for(int i = 0; i < Nc; ++i)
+            for(int i = 0; i < Fimpl::Dimension; ++i)
             {
                 pokeColour(pjs, peekColour(fj, i), i, c);
             }
@@ -436,15 +438,16 @@ namespace QCD {
         }
     }
     
-    template <class Prop, class Ferm>
+    //template <class Prop, class Ferm>
-    void PropToFerm(Ferm &f, const Prop &p, const int s, const int c)
+    template <class Fimpl>
+      void PropToFerm(typename Fimpl::FermionField &f, const typename Fimpl::PropagatorField &p, const int s, const int c)
     {
         for(int j = 0; j < Ns; ++j)
         {
             auto pjs = peekSpin(p, j, s);
             auto fj  = peekSpin(f, j);
             
-            for(int i = 0; i < Nc; ++i)
+            for(int i = 0; i < Fimpl::Dimension; ++i)
             {
                 pokeColour(fj, peekColour(pjs, i, c), i);
             }
@@ -503,38 +506,6 @@ namespace QCD {
 }   //namespace QCD
 } // Grid
 
-/*
-<<<<<<< HEAD
-#include <Grid/qcd/utils/SpaceTimeGrid.h>
-#include <Grid/qcd/spin/Dirac.h>
-#include <Grid/qcd/spin/TwoSpinor.h>
-#include <Grid/qcd/utils/LinalgUtils.h>
-#include <Grid/qcd/utils/CovariantCshift.h>
-
-// Include representations  
-#include <Grid/qcd/utils/SUn.h>
-#include <Grid/qcd/utils/SUnAdjoint.h>
-#include <Grid/qcd/utils/SUnTwoIndex.h>
-#include <Grid/qcd/representations/hmc_types.h>
-
-// Scalar field
-#include <Grid/qcd/utils/ScalarObjs.h>
-
-#include <Grid/qcd/action/Actions.h>
-
-#include <Grid/qcd/smearing/Smearing.h>
-
-#include <Grid/qcd/hmc/integrators/Integrator.h>
-#include <Grid/qcd/hmc/integrators/Integrator_algorithm.h>
-#include <Grid/qcd/observables/hmc_observable.h>
-#include <Grid/qcd/hmc/HMC.h>
-
-
-//#include <Grid/qcd/modules/mods.h>
-=======
-
->>>>>>> develop
-*/
 
 
 #endif

lib/qcd/action/fermion/Fermion.h

@@ -51,10 +51,12 @@ Author: Peter Boyle <pabobyle@ph.ed.ac.uk>
 
 #include <Grid/qcd/action/fermion/WilsonFermion.h>       // 4d wilson like
 #include <Grid/qcd/action/fermion/WilsonTMFermion.h>     // 4d wilson like
+#include <Grid/qcd/action/fermion/WilsonCloverFermion.h> // 4d wilson clover fermions
 #include <Grid/qcd/action/fermion/WilsonFermion5D.h>     // 5d base used by all 5d overlap types
-//#include <Grid/qcd/action/fermion/CloverFermion.h>
+
 #include <Grid/qcd/action/fermion/ImprovedStaggeredFermion.h>
 #include <Grid/qcd/action/fermion/ImprovedStaggeredFermion5D.h>
+
 #include <Grid/qcd/action/fermion/CayleyFermion5D.h>     // Cayley types
 #include <Grid/qcd/action/fermion/DomainWallFermion.h>
 #include <Grid/qcd/action/fermion/DomainWallEOFAFermion.h>
@@ -104,10 +106,33 @@ typedef WilsonFermion<WilsonTwoIndexSymmetricImplR> WilsonTwoIndexSymmetricFermi
 typedef WilsonFermion<WilsonTwoIndexSymmetricImplF> WilsonTwoIndexSymmetricFermionF;
 typedef WilsonFermion<WilsonTwoIndexSymmetricImplD> WilsonTwoIndexSymmetricFermionD;
 
+typedef WilsonFermion<WilsonTwoIndexAntiSymmetricImplR> WilsonTwoIndexAntiSymmetricFermionR;
+typedef WilsonFermion<WilsonTwoIndexAntiSymmetricImplF> WilsonTwoIndexAntiSymmetricFermionF;
+typedef WilsonFermion<WilsonTwoIndexAntiSymmetricImplD> WilsonTwoIndexAntiSymmetricFermionD;
+
+// Twisted mass fermion
 typedef WilsonTMFermion<WilsonImplR> WilsonTMFermionR;
 typedef WilsonTMFermion<WilsonImplF> WilsonTMFermionF;
 typedef WilsonTMFermion<WilsonImplD> WilsonTMFermionD;
 
+// Clover fermions
+typedef WilsonCloverFermion<WilsonImplR> WilsonCloverFermionR;
+typedef WilsonCloverFermion<WilsonImplF> WilsonCloverFermionF;
+typedef WilsonCloverFermion<WilsonImplD> WilsonCloverFermionD;
+
+typedef WilsonCloverFermion<WilsonAdjImplR> WilsonCloverAdjFermionR;
+typedef WilsonCloverFermion<WilsonAdjImplF> WilsonCloverAdjFermionF;
+typedef WilsonCloverFermion<WilsonAdjImplD> WilsonCloverAdjFermionD;
+
+typedef WilsonCloverFermion<WilsonTwoIndexSymmetricImplR> WilsonCloverTwoIndexSymmetricFermionR;
+typedef WilsonCloverFermion<WilsonTwoIndexSymmetricImplF> WilsonCloverTwoIndexSymmetricFermionF;
+typedef WilsonCloverFermion<WilsonTwoIndexSymmetricImplD> WilsonCloverTwoIndexSymmetricFermionD;
+
+typedef WilsonCloverFermion<WilsonTwoIndexAntiSymmetricImplR> WilsonCloverTwoIndexAntiSymmetricFermionR;
+typedef WilsonCloverFermion<WilsonTwoIndexAntiSymmetricImplF> WilsonCloverTwoIndexAntiSymmetricFermionF;
+typedef WilsonCloverFermion<WilsonTwoIndexAntiSymmetricImplD> WilsonCloverTwoIndexAntiSymmetricFermionD;
+
+// Domain Wall fermions
 typedef DomainWallFermion<WilsonImplR> DomainWallFermionR;
 typedef DomainWallFermion<WilsonImplF> DomainWallFermionF;
 typedef DomainWallFermion<WilsonImplD> DomainWallFermionD;

lib/qcd/action/fermion/FermionCore.h

@@ -70,7 +70,9 @@ Author: Peter Boyle <pabobyle@ph.ed.ac.uk>
 
 #define TwoIndexFermOpTemplateInstantiate(A) \
   template class A<WilsonTwoIndexSymmetricImplF>; \
-  template class A<WilsonTwoIndexSymmetricImplD>; 
+  template class A<WilsonTwoIndexSymmetricImplD>; \
+  template class A<WilsonTwoIndexAntiSymmetricImplF>; \
+  template class A<WilsonTwoIndexAntiSymmetricImplD>;
 
 #define FermOp5dVecTemplateInstantiate(A) \
   template class A<DomainWallVec5dImplF>;	\

lib/qcd/action/fermion/FermionOperatorImpl.h

@@ -164,6 +164,7 @@ namespace QCD {
     public:
 
     static const int Dimension = Representation::Dimension;
+    static const bool isFundamental = Representation::isFundamental;
     static const bool LsVectorised=false;
     static const int Nhcs = Options::Nhcs;
 
@@ -263,6 +264,20 @@ namespace QCD {
       PokeIndex<LorentzIndex>(mat,link,mu);
     }  
     
+    inline void outerProductImpl(PropagatorField &mat, const FermionField &B, const FermionField &A){
+      mat = outerProduct(B,A); 
+    }  
+
+    inline void TraceSpinImpl(GaugeLinkField &mat, PropagatorField&P) {
+      mat = TraceIndex<SpinIndex>(P); 
+    }
+      
+    inline void extractLinkField(std::vector<GaugeLinkField> &mat, DoubledGaugeField &Uds){
+      for (int mu = 0; mu < Nd; mu++)
+      mat[mu] = PeekIndex<LorentzIndex>(Uds, mu);
+    }
+
+
     inline void InsertForce5D(GaugeField &mat, FermionField &Btilde, FermionField &Atilde,int mu){
       
       int Ls=Btilde._grid->_fdimensions[0];
@@ -284,27 +299,28 @@ namespace QCD {
   ////////////////////////////////////////////////////////////////////////////////////
   // Single flavour four spinors with colour index, 5d redblack
   ////////////////////////////////////////////////////////////////////////////////////
-template<class S,int Nrepresentation=Nc, class Options=CoeffReal>
+template<class S,class Representation = FundamentalRepresentation, class Options=CoeffReal>
-class DomainWallVec5dImpl :  public PeriodicGaugeImpl< GaugeImplTypes< S,Nrepresentation> > { 
+class DomainWallVec5dImpl :  public PeriodicGaugeImpl< GaugeImplTypes< S,Representation::Dimension> > { 
   public:
 
-  typedef PeriodicGaugeImpl<GaugeImplTypes<S, Nrepresentation> > Gimpl;
+  typedef PeriodicGaugeImpl<GaugeImplTypes<S, Representation::Dimension> > Gimpl;
   INHERIT_GIMPL_TYPES(Gimpl);
 
-  static const int Dimension = Nrepresentation;
+  static const int Dimension = Representation::Dimension;
+  static const bool isFundamental = Representation::isFundamental;
   static const bool LsVectorised=true;
   static const int Nhcs = Options::Nhcs;
       
   typedef typename Options::_Coeff_t Coeff_t;      
   typedef typename Options::template PrecisionMapper<Simd>::LowerPrecVector SimdL;
   
-  template <typename vtype> using iImplSpinor            = iScalar<iVector<iVector<vtype, Nrepresentation>, Ns> >;
+  template <typename vtype> using iImplSpinor            = iScalar<iVector<iVector<vtype, Dimension>, Ns> >;
-  template <typename vtype> using iImplPropagator        = iScalar<iMatrix<iMatrix<vtype, Nrepresentation>, Ns> >;
+  template <typename vtype> using iImplPropagator        = iScalar<iMatrix<iMatrix<vtype, Dimension>, Ns> >;
-  template <typename vtype> using iImplHalfSpinor        = iScalar<iVector<iVector<vtype, Nrepresentation>, Nhs> >;
+  template <typename vtype> using iImplHalfSpinor        = iScalar<iVector<iVector<vtype, Dimension>, Nhs> >;
-  template <typename vtype> using iImplHalfCommSpinor    = iScalar<iVector<iVector<vtype, Nrepresentation>, Nhcs> >;
+  template <typename vtype> using iImplHalfCommSpinor    = iScalar<iVector<iVector<vtype, Dimension>, Nhcs> >;
-  template <typename vtype> using iImplDoubledGaugeField = iVector<iScalar<iMatrix<vtype, Nrepresentation> >, Nds>;
+  template <typename vtype> using iImplDoubledGaugeField = iVector<iScalar<iMatrix<vtype, Dimension> >, Nds>;
-  template <typename vtype> using iImplGaugeField        = iVector<iScalar<iMatrix<vtype, Nrepresentation> >, Nd>;
+  template <typename vtype> using iImplGaugeField        = iVector<iScalar<iMatrix<vtype, Dimension> >, Nd>;
-  template <typename vtype> using iImplGaugeLink         = iScalar<iScalar<iMatrix<vtype, Nrepresentation> > >;
+  template <typename vtype> using iImplGaugeLink         = iScalar<iScalar<iMatrix<vtype, Dimension> > >;
   
   typedef iImplSpinor<Simd>            SiteSpinor;
   typedef iImplPropagator<Simd>        SitePropagator;
@@ -340,8 +356,8 @@ class DomainWallVec5dImpl :  public PeriodicGaugeImpl< GaugeImplTypes< S,Nrepres
                        const SiteHalfSpinor &chi, int mu, StencilEntry *SE,
                        StencilImpl &St) {
     SiteGaugeLink UU;
-    for (int i = 0; i < Nrepresentation; i++) {
+    for (int i = 0; i < Dimension; i++) {
-      for (int j = 0; j < Nrepresentation; j++) {
+      for (int j = 0; j < Dimension; j++) {
         vsplat(UU()()(i, j), U(mu)()(i, j));
       }
     }
@@ -353,8 +369,8 @@ class DomainWallVec5dImpl :  public PeriodicGaugeImpl< GaugeImplTypes< S,Nrepres
                            const SitePropagator &chi,
                            int mu) {
     SiteGaugeLink UU;
-    for (int i = 0; i < Nrepresentation; i++) {
+    for (int i = 0; i < Dimension; i++) {
-      for (int j = 0; j < Nrepresentation; j++) {
+      for (int j = 0; j < Dimension; j++) {
         vsplat(UU()()(i, j), U(mu)()(i, j));
       }
     }
@@ -393,6 +409,19 @@ class DomainWallVec5dImpl :  public PeriodicGaugeImpl< GaugeImplTypes< S,Nrepres
     assert(0);
   }
 
+  inline void outerProductImpl(PropagatorField &mat, const FermionField &Btilde, const FermionField &A){
+    assert(0);
+  } 
+
+  inline void TraceSpinImpl(GaugeLinkField &mat, PropagatorField&P) {
+    assert(0);
+  }
+
+  inline void extractLinkField(std::vector<GaugeLinkField> &mat, DoubledGaugeField &Uds){
+    assert(0);
+  }
+
+
   inline void InsertForce5D(GaugeField &mat, FermionField &Btilde, FermionField &Atilde, int mu) {
 
     assert(0);
@@ -445,25 +474,26 @@ class DomainWallVec5dImpl :  public PeriodicGaugeImpl< GaugeImplTypes< S,Nrepres
     ////////////////////////////////////////////////////////////////////////////////////////
     // Flavour doubled spinors; is Gparity the only? what about C*?
     ////////////////////////////////////////////////////////////////////////////////////////
-template <class S, int Nrepresentation, class Options=CoeffReal>
+template <class S, class Representation = FundamentalRepresentation, class Options=CoeffReal>
-class GparityWilsonImpl : public ConjugateGaugeImpl<GaugeImplTypes<S, Nrepresentation> > {
+class GparityWilsonImpl : public ConjugateGaugeImpl<GaugeImplTypes<S, Representation::Dimension> > {
  public:
 
- static const int Dimension = Nrepresentation;
+ static const int Dimension = Representation::Dimension;
+ static const bool isFundamental = Representation::isFundamental;
  static const int Nhcs = Options::Nhcs;
  static const bool LsVectorised=false;
 
- typedef ConjugateGaugeImpl< GaugeImplTypes<S,Nrepresentation> > Gimpl;
+ typedef ConjugateGaugeImpl< GaugeImplTypes<S,Dimension> > Gimpl;
  INHERIT_GIMPL_TYPES(Gimpl);
 
  typedef typename Options::_Coeff_t Coeff_t;
  typedef typename Options::template PrecisionMapper<Simd>::LowerPrecVector SimdL;
       
- template <typename vtype> using iImplSpinor            = iVector<iVector<iVector<vtype, Nrepresentation>, Ns>,   Ngp>;
+ template <typename vtype> using iImplSpinor            = iVector<iVector<iVector<vtype, Dimension>, Ns>,   Ngp>;
- template <typename vtype> using iImplPropagator        = iVector<iMatrix<iMatrix<vtype, Nrepresentation>, Ns>,   Ngp>;
+ template <typename vtype> using iImplPropagator        = iVector<iMatrix<iMatrix<vtype, Dimension>, Ns>,   Ngp>;
- template <typename vtype> using iImplHalfSpinor        = iVector<iVector<iVector<vtype, Nrepresentation>, Nhs>,  Ngp>;
+ template <typename vtype> using iImplHalfSpinor        = iVector<iVector<iVector<vtype, Dimension>, Nhs>,  Ngp>;
- template <typename vtype> using iImplHalfCommSpinor    = iVector<iVector<iVector<vtype, Nrepresentation>, Nhcs>, Ngp>;
+ template <typename vtype> using iImplHalfCommSpinor    = iVector<iVector<iVector<vtype, Dimension>, Nhcs>, Ngp>;
- template <typename vtype> using iImplDoubledGaugeField = iVector<iVector<iScalar<iMatrix<vtype, Nrepresentation> >, Nds>, Ngp>;
+ template <typename vtype> using iImplDoubledGaugeField = iVector<iVector<iScalar<iMatrix<vtype, Dimension> >, Nds>, Ngp>;
 
  typedef iImplSpinor<Simd>            SiteSpinor;
  typedef iImplPropagator<Simd>        SitePropagator;
@@ -636,6 +666,25 @@ class GparityWilsonImpl : public ConjugateGaugeImpl<GaugeImplTypes<S, Nrepresent
    return;
  }
       
+ inline void outerProductImpl(PropagatorField &mat, const FermionField &Btilde, const FermionField &A){
+   //mat = outerProduct(Btilde, A);
+   assert(0);
+  }
+
+  inline void TraceSpinImpl(GaugeLinkField &mat, PropagatorField&P) {
+    assert(0);
+    /*
+    auto tmp = TraceIndex<SpinIndex>(P);
+    parallel_for(auto ss = tmp.begin(); ss < tmp.end(); ss++) {
+      mat[ss]() = tmp[ss](0, 0) + conjugate(tmp[ss](1, 1));
+    }
+    */
+  }
+
+  inline void extractLinkField(std::vector<GaugeLinkField> &mat, DoubledGaugeField &Uds){
+    assert(0);
+  }
+  
  inline void InsertForce5D(GaugeField &mat, FermionField &Btilde, FermionField &Atilde, int mu) {
 
    int Ls = Btilde._grid->_fdimensions[0];
@@ -665,6 +714,7 @@ class StaggeredImpl : public PeriodicGaugeImpl<GaugeImplTypes<S, Representation:
 
     typedef RealD  _Coeff_t ;
     static const int Dimension = Representation::Dimension;
+    static const bool isFundamental = Representation::isFundamental;
     static const bool LsVectorised=false;
     typedef PeriodicGaugeImpl<GaugeImplTypes<S, Dimension > > Gimpl;
       
@@ -793,6 +843,7 @@ class StaggeredImpl : public PeriodicGaugeImpl<GaugeImplTypes<S, Representation:
     public:
 
     static const int Dimension = Representation::Dimension;
+    static const bool isFundamental = Representation::isFundamental;
     static const bool LsVectorised=true;
     typedef RealD   Coeff_t ;
     typedef PeriodicGaugeImpl<GaugeImplTypes<S, Dimension > > Gimpl;
@@ -983,29 +1034,33 @@ typedef WilsonImpl<vComplex,  TwoIndexSymmetricRepresentation, CoeffReal > Wilso
 typedef WilsonImpl<vComplexF, TwoIndexSymmetricRepresentation, CoeffReal > WilsonTwoIndexSymmetricImplF;  // Float
 typedef WilsonImpl<vComplexD, TwoIndexSymmetricRepresentation, CoeffReal > WilsonTwoIndexSymmetricImplD;  // Double
  
-typedef DomainWallVec5dImpl<vComplex ,Nc, CoeffReal> DomainWallVec5dImplR; // Real.. whichever prec
+typedef WilsonImpl<vComplex,  TwoIndexAntiSymmetricRepresentation, CoeffReal > WilsonTwoIndexAntiSymmetricImplR;   // Real.. whichever prec
-typedef DomainWallVec5dImpl<vComplexF,Nc, CoeffReal> DomainWallVec5dImplF; // Float
+typedef WilsonImpl<vComplexF, TwoIndexAntiSymmetricRepresentation, CoeffReal > WilsonTwoIndexAntiSymmetricImplF;  // Float
-typedef DomainWallVec5dImpl<vComplexD,Nc, CoeffReal> DomainWallVec5dImplD; // Double
+typedef WilsonImpl<vComplexD, TwoIndexAntiSymmetricRepresentation, CoeffReal > WilsonTwoIndexAntiSymmetricImplD;  // Double
 
-typedef DomainWallVec5dImpl<vComplex ,Nc, CoeffRealHalfComms> DomainWallVec5dImplRL; // Real.. whichever prec
+typedef DomainWallVec5dImpl<vComplex ,FundamentalRepresentation, CoeffReal> DomainWallVec5dImplR; // Real.. whichever prec
-typedef DomainWallVec5dImpl<vComplexF,Nc, CoeffRealHalfComms> DomainWallVec5dImplFH; // Float
+typedef DomainWallVec5dImpl<vComplexF,FundamentalRepresentation, CoeffReal> DomainWallVec5dImplF; // Float
-typedef DomainWallVec5dImpl<vComplexD,Nc, CoeffRealHalfComms> DomainWallVec5dImplDF; // Double
+typedef DomainWallVec5dImpl<vComplexD,FundamentalRepresentation, CoeffReal> DomainWallVec5dImplD; // Double
  
-typedef DomainWallVec5dImpl<vComplex ,Nc,CoeffComplex> ZDomainWallVec5dImplR; // Real.. whichever prec
+typedef DomainWallVec5dImpl<vComplex ,FundamentalRepresentation, CoeffRealHalfComms> DomainWallVec5dImplRL; // Real.. whichever prec
-typedef DomainWallVec5dImpl<vComplexF,Nc,CoeffComplex> ZDomainWallVec5dImplF; // Float
+typedef DomainWallVec5dImpl<vComplexF,FundamentalRepresentation, CoeffRealHalfComms> DomainWallVec5dImplFH; // Float
-typedef DomainWallVec5dImpl<vComplexD,Nc,CoeffComplex> ZDomainWallVec5dImplD; // Double
+typedef DomainWallVec5dImpl<vComplexD,FundamentalRepresentation, CoeffRealHalfComms> DomainWallVec5dImplDF; // Double
  
-typedef DomainWallVec5dImpl<vComplex ,Nc,CoeffComplexHalfComms> ZDomainWallVec5dImplRL; // Real.. whichever prec
+typedef DomainWallVec5dImpl<vComplex ,FundamentalRepresentation,CoeffComplex> ZDomainWallVec5dImplR; // Real.. whichever prec
-typedef DomainWallVec5dImpl<vComplexF,Nc,CoeffComplexHalfComms> ZDomainWallVec5dImplFH; // Float
+typedef DomainWallVec5dImpl<vComplexF,FundamentalRepresentation,CoeffComplex> ZDomainWallVec5dImplF; // Float
-typedef DomainWallVec5dImpl<vComplexD,Nc,CoeffComplexHalfComms> ZDomainWallVec5dImplDF; // Double
+typedef DomainWallVec5dImpl<vComplexD,FundamentalRepresentation,CoeffComplex> ZDomainWallVec5dImplD; // Double
  
-typedef GparityWilsonImpl<vComplex , Nc,CoeffReal> GparityWilsonImplR;  // Real.. whichever prec
+typedef DomainWallVec5dImpl<vComplex ,FundamentalRepresentation,CoeffComplexHalfComms> ZDomainWallVec5dImplRL; // Real.. whichever prec
-typedef GparityWilsonImpl<vComplexF, Nc,CoeffReal> GparityWilsonImplF;  // Float
+typedef DomainWallVec5dImpl<vComplexF,FundamentalRepresentation,CoeffComplexHalfComms> ZDomainWallVec5dImplFH; // Float
-typedef GparityWilsonImpl<vComplexD, Nc,CoeffReal> GparityWilsonImplD;  // Double
+typedef DomainWallVec5dImpl<vComplexD,FundamentalRepresentation,CoeffComplexHalfComms> ZDomainWallVec5dImplDF; // Double
  
-typedef GparityWilsonImpl<vComplex , Nc,CoeffRealHalfComms> GparityWilsonImplRL;  // Real.. whichever prec
+typedef GparityWilsonImpl<vComplex , FundamentalRepresentation,CoeffReal> GparityWilsonImplR;  // Real.. whichever prec
-typedef GparityWilsonImpl<vComplexF, Nc,CoeffRealHalfComms> GparityWilsonImplFH;  // Float
+typedef GparityWilsonImpl<vComplexF, FundamentalRepresentation,CoeffReal> GparityWilsonImplF;  // Float
-typedef GparityWilsonImpl<vComplexD, Nc,CoeffRealHalfComms> GparityWilsonImplDF;  // Double
+typedef GparityWilsonImpl<vComplexD, FundamentalRepresentation,CoeffReal> GparityWilsonImplD;  // Double
+ 
+typedef GparityWilsonImpl<vComplex , FundamentalRepresentation,CoeffRealHalfComms> GparityWilsonImplRL;  // Real.. whichever prec
+typedef GparityWilsonImpl<vComplexF, FundamentalRepresentation,CoeffRealHalfComms> GparityWilsonImplFH;  // Float
+typedef GparityWilsonImpl<vComplexD, FundamentalRepresentation,CoeffRealHalfComms> GparityWilsonImplDF;  // Double
 
 typedef StaggeredImpl<vComplex,  FundamentalRepresentation > StaggeredImplR;   // Real.. whichever prec
 typedef StaggeredImpl<vComplexF, FundamentalRepresentation > StaggeredImplF;  // Float

lib/qcd/action/fermion/WilsonCloverFermion.cc

@@ -0,0 +1,243 @@
+/*************************************************************************************
+
+    Grid physics library, www.github.com/paboyle/Grid
+
+    Source file: ./lib/qcd/action/fermion/WilsonCloverFermion.cc
+
+    Copyright (C) 2017
+
+    Author: paboyle <paboyle@ph.ed.ac.uk>
+    Author: Guido Cossu <guido.cossu@ed.ac.uk>
+
+    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 */
+#include <Grid/Grid.h>
+#include <Grid/Eigen/Dense>
+#include <Grid/qcd/spin/Dirac.h>
+
+namespace Grid
+{
+namespace QCD
+{
+
+// *NOT* EO
+template <class Impl>
+RealD WilsonCloverFermion<Impl>::M(const FermionField &in, FermionField &out)
+{
+  FermionField temp(out._grid);
+
+  // Wilson term
+  out.checkerboard = in.checkerboard;
+  this->Dhop(in, out, DaggerNo);
+
+  // Clover term
+  Mooee(in, temp);
+
+  out += temp;
+  return norm2(out);
+}
+
+template <class Impl>
+RealD WilsonCloverFermion<Impl>::Mdag(const FermionField &in, FermionField &out)
+{
+  FermionField temp(out._grid);
+
+  // Wilson term
+  out.checkerboard = in.checkerboard;
+  this->Dhop(in, out, DaggerYes);
+
+  // Clover term
+  MooeeDag(in, temp);
+
+  out += temp;
+  return norm2(out);
+}
+
+template <class Impl>
+void WilsonCloverFermion<Impl>::ImportGauge(const GaugeField &_Umu)
+{
+  WilsonFermion<Impl>::ImportGauge(_Umu);
+  GridBase *grid = _Umu._grid;
+  typename Impl::GaugeLinkField Bx(grid), By(grid), Bz(grid), Ex(grid), Ey(grid), Ez(grid);
+
+  // Compute the field strength terms mu>nu
+  WilsonLoops<Impl>::FieldStrength(Bx, _Umu, Zdir, Ydir);
+  WilsonLoops<Impl>::FieldStrength(By, _Umu, Zdir, Xdir);
+  WilsonLoops<Impl>::FieldStrength(Bz, _Umu, Ydir, Xdir);
+  WilsonLoops<Impl>::FieldStrength(Ex, _Umu, Tdir, Xdir);
+  WilsonLoops<Impl>::FieldStrength(Ey, _Umu, Tdir, Ydir);
+  WilsonLoops<Impl>::FieldStrength(Ez, _Umu, Tdir, Zdir);
+
+  // Compute the Clover Operator acting on Colour and Spin
+  // multiply here by the clover coefficients for the anisotropy
+  CloverTerm  = fillCloverYZ(Bx) * csw_r;
+  CloverTerm += fillCloverXZ(By) * csw_r;
+  CloverTerm += fillCloverXY(Bz) * csw_r;
+  CloverTerm += fillCloverXT(Ex) * csw_t;
+  CloverTerm += fillCloverYT(Ey) * csw_t;
+  CloverTerm += fillCloverZT(Ez) * csw_t;
+  CloverTerm += diag_mass;
+
+  int lvol = _Umu._grid->lSites();
+  int DimRep = Impl::Dimension;
+
+  Eigen::MatrixXcd EigenCloverOp = Eigen::MatrixXcd::Zero(Ns * DimRep, Ns * DimRep);
+  Eigen::MatrixXcd EigenInvCloverOp = Eigen::MatrixXcd::Zero(Ns * DimRep, Ns * DimRep);
+
+  std::vector<int> lcoor;
+  typename SiteCloverType::scalar_object Qx = zero, Qxinv = zero;
+
+  for (int site = 0; site < lvol; site++)
+  {
+    grid->LocalIndexToLocalCoor(site, lcoor);
+    EigenCloverOp = Eigen::MatrixXcd::Zero(Ns * DimRep, Ns * DimRep);
+    peekLocalSite(Qx, CloverTerm, lcoor);
+    Qxinv = zero;
+    //if (csw!=0){
+    for (int j = 0; j < Ns; j++)
+      for (int k = 0; k < Ns; k++)
+        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);
+    //   if (site==0) std::cout << "site =" << site << "\n" << 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++)
+          for (int b = 0; b < DimRep; b++)
+            Qxinv()(j, k)(a, b) = EigenInvCloverOp(a + j * DimRep, b + k * DimRep);
+    //    if (site==0) std::cout << "site =" << site << "\n" << EigenInvCloverOp << std::endl;
+    //  }
+    pokeLocalSite(Qxinv, CloverTermInv, lcoor);
+  }
+
+  // Separate the even and odd parts
+  pickCheckerboard(Even, CloverTermEven, CloverTerm);
+  pickCheckerboard(Odd, CloverTermOdd, CloverTerm);
+
+  pickCheckerboard(Even, CloverTermDagEven, adj(CloverTerm));
+  pickCheckerboard(Odd, CloverTermDagOdd, adj(CloverTerm));
+
+  pickCheckerboard(Even, CloverTermInvEven, CloverTermInv);
+  pickCheckerboard(Odd, CloverTermInvOdd, CloverTermInv);
+
+  pickCheckerboard(Even, CloverTermInvDagEven, adj(CloverTermInv));
+  pickCheckerboard(Odd, CloverTermInvDagOdd, adj(CloverTermInv));
+}
+
+template <class Impl>
+void WilsonCloverFermion<Impl>::Mooee(const FermionField &in, FermionField &out)
+{
+  this->MooeeInternal(in, out, DaggerNo, InverseNo);
+}
+
+template <class Impl>
+void WilsonCloverFermion<Impl>::MooeeDag(const FermionField &in, FermionField &out)
+{
+  this->MooeeInternal(in, out, DaggerYes, InverseNo);
+}
+
+template <class Impl>
+void WilsonCloverFermion<Impl>::MooeeInv(const FermionField &in, FermionField &out)
+{
+  this->MooeeInternal(in, out, DaggerNo, InverseYes);
+}
+
+template <class Impl>
+void WilsonCloverFermion<Impl>::MooeeInvDag(const FermionField &in, FermionField &out)
+{
+  this->MooeeInternal(in, out, DaggerYes, InverseYes);
+}
+
+template <class Impl>
+void WilsonCloverFermion<Impl>::MooeeInternal(const FermionField &in, FermionField &out, int dag, int inv)
+{
+  out.checkerboard = in.checkerboard;
+  CloverFieldType *Clover;
+  assert(in.checkerboard == Odd || in.checkerboard == Even);
+
+  if (dag)
+  {
+    if (in._grid->_isCheckerBoarded)
+    {
+      if (in.checkerboard == Odd)
+      {
+        Clover = (inv) ? &CloverTermInvDagOdd : &CloverTermDagOdd;
+      }
+      else
+      {
+        Clover = (inv) ? &CloverTermInvDagEven : &CloverTermDagEven;
+      }
+      out = *Clover * in;
+    }
+    else
+    {
+      Clover = (inv) ? &CloverTermInv : &CloverTerm;
+      out = adj(*Clover) * in;
+    }
+  }
+  else
+  {
+    if (in._grid->_isCheckerBoarded)
+    {
+
+      if (in.checkerboard == Odd)
+      {
+        //  std::cout << "Calling clover term Odd" << std::endl;
+        Clover = (inv) ? &CloverTermInvOdd : &CloverTermOdd;
+      }
+      else
+      {
+        //  std::cout << "Calling clover term Even" << std::endl;
+        Clover = (inv) ? &CloverTermInvEven : &CloverTermEven;
+      }
+      out = *Clover * in;
+      //  std::cout << GridLogMessage << "*Clover.checkerboard "  << (*Clover).checkerboard << std::endl;
+    }
+    else
+    {
+      Clover = (inv) ? &CloverTermInv : &CloverTerm;
+      out = *Clover * in;
+    }
+  }
+
+} // MooeeInternal
+
+
+// Derivative parts
+template <class Impl>
+void WilsonCloverFermion<Impl>::MooDeriv(GaugeField &mat, const FermionField &X, const FermionField &Y, int dag)
+{
+  assert(0);
+}
+
+// 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);
+AdjointFermOpTemplateInstantiate(WilsonCloverFermion);
+TwoIndexFermOpTemplateInstantiate(WilsonCloverFermion);
+//GparityFermOpTemplateInstantiate(WilsonCloverFermion);
+}
+}

lib/qcd/action/fermion/WilsonCloverFermion.h

@@ -0,0 +1,366 @@
+/*************************************************************************************
+
+    Grid physics library, www.github.com/paboyle/Grid
+
+    Source file: ./lib/qcd/action/fermion/WilsonCloverFermion.h
+
+    Copyright (C) 2017
+
+    Author: Guido Cossu <guido.cossu@ed.ac.uk>
+    Author: David Preti <>
+
+    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 */
+
+#ifndef GRID_QCD_WILSON_CLOVER_FERMION_H
+#define GRID_QCD_WILSON_CLOVER_FERMION_H
+
+#include <Grid/Grid.h>
+
+namespace Grid
+{
+namespace QCD
+{
+
+///////////////////////////////////////////////////////////////////
+// Wilson Clover
+//
+// Operator ( with anisotropy coefficients):
+//
+// Q =   1 + (Nd-1)/xi_0 + m
+//     + W_t + (nu/xi_0) * W_s
+//     - 1/2*[ csw_t * sum_s (sigma_ts F_ts) + (csw_s/xi_0) * sum_ss (sigma_ss F_ss)  ]
+//
+// s spatial, t temporal directions.
+// where W_t and W_s are the temporal and spatial components of the
+// Wilson Dirac operator
+//
+// csw_r = csw_t to recover the isotropic version
+//////////////////////////////////////////////////////////////////
+
+template <class Impl>
+class WilsonCloverFermion : public WilsonFermion<Impl>
+{
+public:
+  // Types definitions
+  INHERIT_IMPL_TYPES(Impl);
+  template <typename vtype>
+  using iImplClover = iScalar<iMatrix<iMatrix<vtype, Impl::Dimension>, Ns>>;
+  typedef iImplClover<Simd> SiteCloverType;
+  typedef Lattice<SiteCloverType> CloverFieldType;
+
+public:
+  typedef WilsonFermion<Impl> WilsonBase;
+
+  virtual void Instantiatable(void){};
+  // Constructors
+  WilsonCloverFermion(GaugeField &_Umu, GridCartesian &Fgrid,
+                      GridRedBlackCartesian &Hgrid,
+                      const RealD _mass,
+                      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,
+                                                                                     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 RealD M(const FermionField &in, FermionField &out);
+  virtual RealD Mdag(const FermionField &in, FermionField &out);
+
+  virtual void Mooee(const FermionField &in, FermionField &out);
+  virtual void MooeeDag(const FermionField &in, FermionField &out);
+  virtual void MooeeInv(const FermionField &in, FermionField &out);
+  virtual void MooeeInvDag(const FermionField &in, FermionField &out);
+  virtual void MooeeInternal(const FermionField &in, FermionField &out, int dag, int inv);
+
+  //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);
+
+  // Derivative parts unpreconditioned pseudofermions
+  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:
+    // 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;
+  }
+
+private:
+  // 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
+  CloverFieldType CloverTermEven, CloverTermOdd;             // Clover term EO
+  CloverFieldType CloverTermInvEven, CloverTermInvOdd;       // Clover term Inv EO
+  CloverFieldType CloverTermDagEven, CloverTermDagOdd;       // Clover term Dag EO
+  CloverFieldType CloverTermInvDagEven, CloverTermInvDagOdd; // Clover term Inv Dag EO
+
+  // 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;
+    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) = 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) = timesI(F._odata[i]()());
+      T._odata[i]()(1, 0) = timesI(F._odata[i]()());
+      T._odata[i]()(2, 3) = timesMinusI(F._odata[i]()());
+      T._odata[i]()(3, 2) = timesMinusI(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) = timesI(F._odata[i]()());
+      T._odata[i]()(1, 1) = timesMinusI(F._odata[i]()());
+      T._odata[i]()(2, 2) = timesMinusI(F._odata[i]()());
+      T._odata[i]()(3, 3) = timesI(F._odata[i]()());
+    }
+
+    return T;
+  }
+};
+}
+}
+
+#endif // GRID_QCD_WILSON_CLOVER_FERMION_H

lib/qcd/action/fermion/WilsonFermion.cc

@@ -47,7 +47,8 @@ int WilsonFermionStatic::HandOptDslash;
 template <class Impl>
 WilsonFermion<Impl>::WilsonFermion(GaugeField &_Umu, GridCartesian &Fgrid,
                                    GridRedBlackCartesian &Hgrid, RealD _mass,
-                                   const ImplParams &p)
+                                   const ImplParams &p,
+                                   const WilsonAnisotropyCoefficients &anis)
     : Kernels(p),
       _grid(&Fgrid),
       _cbgrid(&Hgrid),
@@ -60,16 +61,41 @@ WilsonFermion<Impl>::WilsonFermion(GaugeField &_Umu, GridCartesian &Fgrid,
       Umu(&Fgrid),
       UmuEven(&Hgrid),
       UmuOdd(&Hgrid),
-      _tmp(&Hgrid)
+      _tmp(&Hgrid),
+      anisotropyCoeff(anis)
 {
   // Allocate the required comms buffer
   ImportGauge(_Umu);
+  if  (anisotropyCoeff.isAnisotropic){
+    diag_mass = mass + 1.0 + (Nd-1)*(anisotropyCoeff.nu / anisotropyCoeff.xi_0);
+  } else {
+    diag_mass = 4.0 + mass;
+  }
+
+
 }
 
 template <class Impl>
 void WilsonFermion<Impl>::ImportGauge(const GaugeField &_Umu) {
   GaugeField HUmu(_Umu._grid);
+
+  //Here multiply the anisotropy coefficients
+  if (anisotropyCoeff.isAnisotropic)
+  {
+
+    for (int mu = 0; mu < Nd; mu++)
+    {
+      GaugeLinkField U_dir = (-0.5)*PeekIndex<LorentzIndex>(_Umu, mu);
+      if (mu != anisotropyCoeff.t_direction)
+        U_dir *= (anisotropyCoeff.nu / anisotropyCoeff.xi_0);
+
+      PokeIndex<LorentzIndex>(HUmu, U_dir, mu);
+    }
+  }
+  else
+  {
     HUmu = _Umu * (-0.5);
+  }
   Impl::DoubleStore(GaugeGrid(), Umu, HUmu);
   pickCheckerboard(Even, UmuEven, Umu);
   pickCheckerboard(Odd, UmuOdd, Umu);
@@ -83,14 +109,14 @@ template <class Impl>
 RealD WilsonFermion<Impl>::M(const FermionField &in, FermionField &out) {
   out.checkerboard = in.checkerboard;
   Dhop(in, out, DaggerNo);
-  return axpy_norm(out, 4 + mass, in, out);
+  return axpy_norm(out, diag_mass, in, out);
 }
 
 template <class Impl>
 RealD WilsonFermion<Impl>::Mdag(const FermionField &in, FermionField &out) {
   out.checkerboard = in.checkerboard;
   Dhop(in, out, DaggerYes);
-  return axpy_norm(out, 4 + mass, in, out);
+  return axpy_norm(out, diag_mass, in, out);
 }
 
 template <class Impl>
@@ -114,7 +140,7 @@ void WilsonFermion<Impl>::MeooeDag(const FermionField &in, FermionField &out) {
 template <class Impl>
 void WilsonFermion<Impl>::Mooee(const FermionField &in, FermionField &out) {
   out.checkerboard = in.checkerboard;
-  typename FermionField::scalar_type scal(4.0 + mass);
+  typename FermionField::scalar_type scal(diag_mass);
   out = scal * in;
 }
 
@@ -127,7 +153,7 @@ void WilsonFermion<Impl>::MooeeDag(const FermionField &in, FermionField &out) {
 template<class Impl>
 void WilsonFermion<Impl>::MooeeInv(const FermionField &in, FermionField &out) {
   out.checkerboard = in.checkerboard;
-  out = (1.0/(4.0+mass))*in;
+  out = (1.0/(diag_mass))*in;
 }
   
 template<class Impl>
@@ -204,7 +230,7 @@ void WilsonFermion<Impl>::DerivInternal(StencilImpl &st, DoubledGaugeField &U,
 
   FermionField Btilde(B._grid);
   FermionField Atilde(B._grid);
-  Atilde = A;
+  Atilde = A;//redundant
 
   st.HaloExchange(B, compressor);
 
@@ -393,7 +419,7 @@ void WilsonFermion<Impl>::SeqConservedCurrent(PropagatorField &q_in,
     conformable(_grid, q_in._grid);
     conformable(_grid, q_out._grid);
     Lattice<iSinglet<Simd>> ph(_grid), coor(_grid);
-    Complex i(0.0,1.0);
+    ComplexD i(0.0,1.0);
     PropagatorField tmpFwd(_grid), tmpBwd(_grid), tmp(_grid);
     unsigned int tshift = (mu == Tp) ? 1 : 0;
     unsigned int LLt    = GridDefaultLatt()[Tp];
@@ -405,7 +431,7 @@ void WilsonFermion<Impl>::SeqConservedCurrent(PropagatorField &q_in,
         LatticeCoordinate(coor, mu);
         ph = ph + mom[mu]*coor*((1./(_grid->_fdimensions[mu])));
     }
-    ph = exp((Real)(2*M_PI)*i*ph);
+    ph = exp((RealD)(2*M_PI)*i*ph);
 
     q_out = zero;
     LatticeInteger coords(_grid);

lib/qcd/action/fermion/WilsonFermion.h

@@ -44,6 +44,21 @@ class WilsonFermionStatic {
   static const int npoint = 8;
 };
 
+ struct WilsonAnisotropyCoefficients: Serializable
+ {
+  GRID_SERIALIZABLE_CLASS_MEMBERS(WilsonAnisotropyCoefficients,
+  bool, isAnisotropic,
+  int, t_direction,
+  double, xi_0,
+  double, nu);
+
+  WilsonAnisotropyCoefficients():
+    isAnisotropic(false), 
+    t_direction(Nd-1), 
+    xi_0(1.0), 
+    nu(1.0){}
+};
+
 template <class Impl>
 class WilsonFermion : public WilsonKernels<Impl>, public WilsonFermionStatic {
  public:
@@ -65,8 +80,8 @@ class WilsonFermion : public WilsonKernels<Impl>, public WilsonFermionStatic {
   // override multiply; cut number routines if pass dagger argument
   // and also make interface more uniformly consistent
   //////////////////////////////////////////////////////////////////
-  RealD M(const FermionField &in, FermionField &out);
+  virtual RealD M(const FermionField &in, FermionField &out);
-  RealD Mdag(const FermionField &in, FermionField &out);
+  virtual RealD Mdag(const FermionField &in, FermionField &out);
 
   /////////////////////////////////////////////////////////
   // half checkerboard operations
@@ -118,7 +133,8 @@ class WilsonFermion : public WilsonKernels<Impl>, public WilsonFermionStatic {
   // Constructor
   WilsonFermion(GaugeField &_Umu, GridCartesian &Fgrid,
                 GridRedBlackCartesian &Hgrid, RealD _mass, 
-                const ImplParams &p = ImplParams());
+                const ImplParams &p = ImplParams(), 
+                const WilsonAnisotropyCoefficients &anis = WilsonAnisotropyCoefficients() );
 
   // DoubleStore impl dependent
   void ImportGauge(const GaugeField &_Umu);
@@ -130,6 +146,7 @@ class WilsonFermion : public WilsonKernels<Impl>, public WilsonFermionStatic {
   //    protected:
  public:
   RealD mass;
+  RealD diag_mass;
 
   GridBase *_grid;
   GridBase *_cbgrid;
@@ -147,6 +164,8 @@ class WilsonFermion : public WilsonKernels<Impl>, public WilsonFermionStatic {
   LebesgueOrder Lebesgue;
   LebesgueOrder LebesgueEvenOdd;
 
+  WilsonAnisotropyCoefficients anisotropyCoeff;
+  
   ///////////////////////////////////////////////////////////////
   // Conserved current utilities
   ///////////////////////////////////////////////////////////////

lib/qcd/action/fermion/WilsonFermion5D.cc

@@ -793,7 +793,7 @@ void WilsonFermion5D<Impl>::SeqConservedCurrent(PropagatorField &q_in,
     Lattice<iSinglet<Simd>> ph(FermionGrid()), coor(FermionGrid());
     PropagatorField tmpFwd(FermionGrid()), tmpBwd(FermionGrid()),
                     tmp(FermionGrid());
-    Complex i(0.0, 1.0);
+    ComplexD i(0.0, 1.0);
     unsigned int tshift = (mu == Tp) ? 1 : 0;
     unsigned int LLs = q_in._grid->_rdimensions[0];
     unsigned int LLt    = GridDefaultLatt()[Tp];
@@ -806,7 +806,7 @@ void WilsonFermion5D<Impl>::SeqConservedCurrent(PropagatorField &q_in,
         LatticeCoordinate(coor, nu + 1);
         ph = ph + mom[nu]*coor*((1./(_FourDimGrid->_fdimensions[nu])));
     }
-    ph = exp((Real)(2*M_PI)*i*ph);
+    ph = exp((RealD)(2*M_PI)*i*ph);
 
     q_out = zero;
     LatticeInteger coords(_FourDimGrid);

lib/qcd/action/fermion/WilsonKernels.h

@@ -55,7 +55,7 @@ template<class Impl> class WilsonKernels : public FermionOperator<Impl> , public
 public:
    
   template <bool EnableBool = true>
-  typename std::enable_if<Impl::Dimension == 3 && Nc == 3 &&EnableBool, void>::type
+  typename std::enable_if<Impl::isFundamental==true && Nc == 3 &&EnableBool, void>::type
   DhopSite(StencilImpl &st, LebesgueOrder &lo, DoubledGaugeField &U, SiteHalfSpinor * buf,
 		   int sF, int sU, int Ls, int Ns, const FermionField &in, FermionField &out,int interior=1,int exterior=1) 
   {
@@ -99,7 +99,7 @@ public:
   }
      
   template <bool EnableBool = true>
-  typename std::enable_if<(Impl::Dimension != 3 || (Impl::Dimension == 3 && Nc != 3)) && EnableBool, void>::type
+  typename std::enable_if<(Impl::isFundamental==false || (Impl::isFundamental==true && Nc != 3)) && EnableBool, void>::type
   DhopSite(StencilImpl &st, LebesgueOrder &lo, DoubledGaugeField &U, SiteHalfSpinor * buf,
 	   int sF, int sU, int Ls, int Ns, const FermionField &in, FermionField &out,int interior=1,int exterior=1 ) {
     // no kernel choice  
@@ -116,7 +116,7 @@ public:
   }
      
   template <bool EnableBool = true>
-  typename std::enable_if<Impl::Dimension == 3 && Nc == 3 && EnableBool,void>::type
+  typename std::enable_if<Impl::isFundamental==true && Nc == 3 && EnableBool,void>::type
   DhopSiteDag(StencilImpl &st, LebesgueOrder &lo, DoubledGaugeField &U, SiteHalfSpinor * buf,
 	      int sF, int sU, int Ls, int Ns, const FermionField &in, FermionField &out,int interior=1,int exterior=1) 
 {
@@ -161,7 +161,7 @@ public:
   }
 
   template <bool EnableBool = true>
-  typename std::enable_if<(Impl::Dimension != 3 || (Impl::Dimension == 3 && Nc != 3)) && EnableBool,void>::type
+  typename std::enable_if<(Impl::isFundamental==false || (Impl::isFundamental==true && Nc != 3)) && EnableBool,void>::type
   DhopSiteDag(StencilImpl &st, LebesgueOrder &lo, DoubledGaugeField &U,SiteHalfSpinor * buf,
 		      int sF, int sU, int Ls, int Ns, const FermionField &in, FermionField &out,int interior=1,int exterior=1) {
 

lib/qcd/action/fermion/WilsonKernelsHand.cc

@@ -946,5 +946,6 @@ INSTANTIATE_THEM(DomainWallVec5dImplFH);
 INSTANTIATE_THEM(DomainWallVec5dImplDF);
 INSTANTIATE_THEM(ZDomainWallVec5dImplFH);
 INSTANTIATE_THEM(ZDomainWallVec5dImplDF);
-
+INSTANTIATE_THEM(WilsonTwoIndexAntiSymmetricImplF);
+INSTANTIATE_THEM(WilsonTwoIndexAntiSymmetricImplD);
 }}

lib/qcd/modules/ObservableModules.h

@@ -92,6 +92,19 @@ class PlaquetteMod: public ObservableModule<PlaquetteLogger<Impl>, NoParameters>
   PlaquetteMod(): ObsBase(NoParameters()){}
 };
 
+template < class Impl >
+class PolyakovMod: public ObservableModule<PolyakovLogger<Impl>, NoParameters>{
+  typedef ObservableModule<PolyakovLogger<Impl>, NoParameters> ObsBase;
+  using ObsBase::ObsBase; // for constructors
+
+  // acquire resource
+  virtual void initialize(){
+    this->ObservablePtr.reset(new PolyakovLogger<Impl>());
+  }
+  public:
+  PolyakovMod(): ObsBase(NoParameters()){}
+};
+
 
 template < class Impl >
 class TopologicalChargeMod: public ObservableModule<TopologicalCharge<Impl>, TopologyObsParameters>{

lib/qcd/observables/hmc_observable.h

@@ -45,5 +45,7 @@ class HmcObservable {
 
 #include "plaquette.h"
 #include "topological_charge.h"
+#include "polyakov_loop.h"
+
 
 #endif  //  HMC_OBSERVABLE_H

lib/qcd/observables/polyakov_loop.h

@@ -0,0 +1,68 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid
+
+Source file: ./lib/qcd/modules/polyakov_line.h
+
+Copyright (C) 2017
+
+Author: David Preti <david.preti@csic.es>
+
+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 */
+
+#ifndef HMC_POLYAKOV_H
+#define HMC_POLYAKOV_H
+
+namespace Grid {
+namespace QCD {
+
+// this is only defined for a gauge theory
+template <class Impl>
+class PolyakovLogger : public HmcObservable<typename Impl::Field> {
+ public:
+  // here forces the Impl to be of gauge fields
+  // if not the compiler will complain
+  INHERIT_GIMPL_TYPES(Impl);
+
+  // necessary for HmcObservable compatibility
+  typedef typename Impl::Field Field;
+
+  void TrajectoryComplete(int traj,
+                          Field &U,
+                          GridSerialRNG &sRNG,
+                          GridParallelRNG &pRNG) {
+
+    ComplexD polyakov = WilsonLoops<Impl>::avgPolyakovLoop(U);
+
+    int def_prec = std::cout.precision();
+
+    std::cout << GridLogMessage
+        << std::setprecision(std::numeric_limits<Real>::digits10 + 1)
+        << "Polyakov Loop: [ " << traj << " ] "<< polyakov << std::endl;
+
+    std::cout.precision(def_prec);
+
+  }
+};
+
+}  // namespace QCD
+}  // namespace Grid
+
+#endif  // HMC_POLYAKOV_H

lib/qcd/representations/adjoint.h

@@ -23,6 +23,7 @@ class AdjointRep {
   typedef typename SU_Adjoint<ncolour>::LatticeAdjMatrix LatticeMatrix;
   typedef typename SU_Adjoint<ncolour>::LatticeAdjField LatticeField;
   static const int Dimension = ncolour * ncolour - 1;
+  static const bool isFundamental = false;
 
   LatticeField U;
 

lib/qcd/representations/fundamental.h

@@ -19,6 +19,7 @@ template <int ncolour>
 class FundamentalRep {
  public:
   static const int Dimension = ncolour;
+  static const bool isFundamental = true;
 
   // typdef to be used by the Representations class in HMC to get the
   // types for the higher representation fields

lib/qcd/representations/two_index.h

@@ -29,6 +29,7 @@ class TwoIndexRep {
   typedef typename SU_TwoIndex<ncolour, S>::LatticeTwoIndexMatrix LatticeMatrix;
   typedef typename SU_TwoIndex<ncolour, S>::LatticeTwoIndexField LatticeField;
   static const int Dimension = ncolour * (ncolour + S) / 2;
+  static const bool isFundamental = false;
 
   LatticeField U;
 

lib/qcd/utils/WilsonLoops.h

@@ -123,6 +123,28 @@ public:
     return sumplaq / vol / faces / Nc; // Nd , Nc dependent... FIXME
   }
 
+
+  //////////////////////////////////////////////////
+  // average over all x,y,z the temporal loop
+  //////////////////////////////////////////////////
+  static ComplexD avgPolyakovLoop(const GaugeField &Umu) {  //assume Nd=4
+    GaugeMat Ut(Umu._grid), P(Umu._grid);
+    ComplexD out;
+    int T = Umu._grid->GlobalDimensions()[3];
+    int X = Umu._grid->GlobalDimensions()[0];
+    int Y = Umu._grid->GlobalDimensions()[1];
+    int Z = Umu._grid->GlobalDimensions()[2];
+
+    Ut = peekLorentz(Umu,3); //Select temporal direction
+    P = Ut;
+    for (int t=1;t<T;t++){ 
+      P = Gimpl::CovShiftForward(Ut,3,P);
+    }
+   RealD norm = 1.0/(Nc*X*Y*Z*T);
+   out = sum(trace(P))*norm;
+   return out;   
+}
+
   //////////////////////////////////////////////////
   // average over traced single links
   //////////////////////////////////////////////////
@@ -291,9 +313,9 @@ static void StapleMult(GaugeMat &staple, const GaugeLorentz &Umu, int mu) {
     }
   }
 
-  //////////////////////////////////////////////////
+  ////////////////////////////////////////////////////////////////////////
   // 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) {
@@ -315,7 +337,9 @@ static void StapleMult(GaugeMat &staple, const GaugeLorentz &Umu, int mu) {
       //
       staple = Gimpl::ShiftStaple(
           Gimpl::CovShiftBackward(U[nu], nu,
-                                  Gimpl::CovShiftBackward(U[mu], mu, U[nu])), mu);
+                                  Gimpl::CovShiftBackward(U[mu], mu, U[nu])),
+          mu);
+
     }
   }
 
@@ -325,7 +349,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.
       //     |     |     |     |
       //     +--<--+     +--<--+
 
@@ -335,7 +359,9 @@ static void StapleMult(GaugeMat &staple, const GaugeLorentz &Umu, int mu) {
       GaugeMat v = Vup - Vdn;
       GaugeMat u = PeekIndex<LorentzIndex>(Umu, mu);  // some redundant copies
       GaugeMat vu = v*u;
-      FS = 0.25*Ta(u*v + Cshift(vu, mu, -1));
+      //FS = 0.25*Ta(u*v + Cshift(vu, mu, -1));
+      FS = (u*v + Cshift(vu, mu, -1));
+      FS = 0.125*(FS - adj(FS));
   }
 
   static Real TopologicalCharge(GaugeLorentz &U){
@@ -360,6 +386,7 @@ static void StapleMult(GaugeMat &staple, const GaugeLorentz &Umu, int mu) {
     return TensorRemove(Tq).real();
   }
 
+
   //////////////////////////////////////////////////////
   // Similar to above for rectangle is required
   //////////////////////////////////////////////////////

tests/core/Test_wilson_clover.cc

@@ -0,0 +1,357 @@
+/*************************************************************************************
+
+    Grid physics library, www.github.com/paboyle/Grid 
+
+    Source file: ./benchmarks/Benchmark_wilson.cc
+
+    Copyright (C) 2015
+
+    Author: Guido Cossu <guido.cossu@ed.ac.uk>
+
+    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 */
+#include <Grid/Grid.h>
+
+using namespace std;
+using namespace Grid;
+using namespace Grid::QCD;
+
+int main(int argc, char **argv)
+{
+  Grid_init(&argc, &argv);
+
+  std::vector<int> latt_size = GridDefaultLatt();
+  std::vector<int> simd_layout = GridDefaultSimd(Nd, vComplex::Nsimd());
+  std::vector<int> mpi_layout = GridDefaultMpi();
+  GridCartesian Grid(latt_size, simd_layout, mpi_layout);
+  GridRedBlackCartesian RBGrid(&Grid);
+
+  int threads = GridThread::GetThreads();
+  std::cout << GridLogMessage << "Grid is setup to use " << threads << " threads" << std::endl;
+  std::cout << GridLogMessage << "Grid floating point word size is REALF" << sizeof(RealF) << std::endl;
+  std::cout << GridLogMessage << "Grid floating point word size is REALD" << sizeof(RealD) << std::endl;
+  std::cout << GridLogMessage << "Grid floating point word size is REAL" << sizeof(Real) << std::endl;
+
+  std::vector<int> seeds({1, 2, 3, 4});
+  GridParallelRNG pRNG(&Grid);
+  pRNG.SeedFixedIntegers(seeds);
+  //  pRNG.SeedFixedIntegers(std::vector<int>({45,12,81,9});
+
+  typedef typename WilsonCloverFermionR::FermionField FermionField;
+  typename WilsonCloverFermionR::ImplParams params;
+  WilsonAnisotropyCoefficients anis;
+
+  FermionField src(&Grid);
+  random(pRNG, src);
+  FermionField result(&Grid);
+  result = zero;
+  FermionField result2(&Grid);
+  result2 = zero;
+  FermionField ref(&Grid);
+  ref = zero;
+  FermionField tmp(&Grid);
+  tmp = zero;
+  FermionField err(&Grid);
+  err = zero;
+  FermionField err2(&Grid);
+  err2 = zero;
+  FermionField phi(&Grid);
+  random(pRNG, phi);
+  FermionField chi(&Grid);
+  random(pRNG, chi);
+  LatticeGaugeField Umu(&Grid);
+  SU3::HotConfiguration(pRNG, Umu);
+  std::vector<LatticeColourMatrix> U(4, &Grid);
+
+  double volume = 1;
+  for (int mu = 0; mu < Nd; mu++)
+  {
+    volume = volume * latt_size[mu];
+  }
+
+  RealD mass = 0.1;
+  RealD csw_r = 1.0;
+  RealD csw_t = 1.0;
+
+  WilsonCloverFermionR Dwc(Umu, Grid, RBGrid, mass, csw_r, csw_t, anis, params);
+  //Dwc.ImportGauge(Umu); // not necessary, included in the constructor
+
+  std::cout << GridLogMessage << "==========================================================" << std::endl;
+  std::cout << GridLogMessage << "= Testing that Deo + Doe = Dunprec " << std::endl;
+  std::cout << GridLogMessage << "==========================================================" << std::endl;
+
+  FermionField src_e(&RBGrid);
+  FermionField src_o(&RBGrid);
+  FermionField r_e(&RBGrid);
+  FermionField r_o(&RBGrid);
+  FermionField r_eo(&Grid);
+  pickCheckerboard(Even, src_e, src);
+  pickCheckerboard(Odd, src_o, src);
+
+  Dwc.Meooe(src_e, r_o);
+  std::cout << GridLogMessage << "Applied Meo" << std::endl;
+  Dwc.Meooe(src_o, r_e);
+  std::cout << GridLogMessage << "Applied Moe" << std::endl;
+  Dwc.Dhop(src, ref, DaggerNo);
+
+  setCheckerboard(r_eo, r_o);
+  setCheckerboard(r_eo, r_e);
+
+  err = ref - r_eo;
+  std::cout << GridLogMessage << "EO norm diff   " << norm2(err) << " " << norm2(ref) << " " << norm2(r_eo) << std::endl;
+
+  std::cout << GridLogMessage << "==============================================================" << std::endl;
+  std::cout << GridLogMessage << "= Test Ddagger is the dagger of D by requiring                " << std::endl;
+  std::cout << GridLogMessage << "=  < phi | Deo | chi > * = < chi | Deo^dag| phi>  " << std::endl;
+  std::cout << GridLogMessage << "==============================================================" << std::endl;
+
+  FermionField chi_e(&RBGrid);
+  FermionField chi_o(&RBGrid);
+
+  FermionField dchi_e(&RBGrid);
+  FermionField dchi_o(&RBGrid);
+
+  FermionField phi_e(&RBGrid);
+  FermionField phi_o(&RBGrid);
+
+  FermionField dphi_e(&RBGrid);
+  FermionField dphi_o(&RBGrid);
+
+  pickCheckerboard(Even, chi_e, chi);
+  pickCheckerboard(Odd, chi_o, chi);
+  pickCheckerboard(Even, phi_e, phi);
+  pickCheckerboard(Odd, phi_o, phi);
+
+  Dwc.Meooe(chi_e, dchi_o);
+  Dwc.Meooe(chi_o, dchi_e);
+  Dwc.MeooeDag(phi_e, dphi_o);
+  Dwc.MeooeDag(phi_o, dphi_e);
+
+  ComplexD pDce = innerProduct(phi_e, dchi_e);
+  ComplexD pDco = innerProduct(phi_o, dchi_o);
+  ComplexD cDpe = innerProduct(chi_e, dphi_e);
+  ComplexD cDpo = innerProduct(chi_o, dphi_o);
+
+  std::cout << GridLogMessage << "e " << pDce << " " << cDpe << std::endl;
+  std::cout << GridLogMessage << "o " << pDco << " " << cDpo << std::endl;
+
+  std::cout << GridLogMessage << "pDce - conj(cDpo) " << pDce - conj(cDpo) << std::endl;
+  std::cout << GridLogMessage << "pDco - conj(cDpe) " << pDco - conj(cDpe) << std::endl;
+
+  std::cout << GridLogMessage << "==============================================================" << std::endl;
+  std::cout << GridLogMessage << "= Test MeeInv Mee = 1   (if csw!=0)                           " << std::endl;
+  std::cout << GridLogMessage << "==============================================================" << std::endl;
+
+  pickCheckerboard(Even, chi_e, chi);
+  pickCheckerboard(Odd, chi_o, chi);
+
+  Dwc.Mooee(chi_e, src_e);
+  Dwc.MooeeInv(src_e, phi_e);
+
+  Dwc.Mooee(chi_o, src_o);
+  Dwc.MooeeInv(src_o, phi_o);
+
+  setCheckerboard(phi, phi_e);
+  setCheckerboard(phi, phi_o);
+
+  err = phi - chi;
+  std::cout << GridLogMessage << "norm diff   " << norm2(err) << std::endl;
+
+  std::cout << GridLogMessage << "==============================================================" << std::endl;
+  std::cout << GridLogMessage << "= Test MeeDag MeeInvDag = 1    (if csw!=0)                    " << std::endl;
+  std::cout << GridLogMessage << "==============================================================" << std::endl;
+
+  pickCheckerboard(Even, chi_e, chi);
+  pickCheckerboard(Odd, chi_o, chi);
+
+  Dwc.MooeeDag(chi_e, src_e);
+  Dwc.MooeeInvDag(src_e, phi_e);
+
+  Dwc.MooeeDag(chi_o, src_o);
+  Dwc.MooeeInvDag(src_o, phi_o);
+
+  setCheckerboard(phi, phi_e);
+  setCheckerboard(phi, phi_o);
+
+  err = phi - chi;
+  std::cout << GridLogMessage << "norm diff   " << norm2(err) << std::endl;
+
+  std::cout << GridLogMessage << "==============================================================" << std::endl;
+  std::cout << GridLogMessage << "= Test MeeInv MeeDag = 1      (if csw!=0)                     " << std::endl;
+  std::cout << GridLogMessage << "==============================================================" << std::endl;
+
+  pickCheckerboard(Even, chi_e, chi);
+  pickCheckerboard(Odd, chi_o, chi);
+
+  Dwc.MooeeDag(chi_e, src_e);
+  Dwc.MooeeInv(src_e, phi_e);
+
+  Dwc.MooeeDag(chi_o, src_o);
+  Dwc.MooeeInv(src_o, phi_o);
+
+  setCheckerboard(phi, phi_e);
+  setCheckerboard(phi, phi_o);
+
+  err = phi - chi;
+  std::cout << GridLogMessage << "norm diff   " << norm2(err) << std::endl;
+
+  std::cout << GridLogMessage << "================================================================" << std::endl;
+  std::cout << GridLogMessage << "= Testing gauge covariance Clover term with EO preconditioning  " << std::endl;
+  std::cout << GridLogMessage << "================================================================" << std::endl;
+
+  chi = zero;
+  phi = zero;
+  tmp = zero;
+  pickCheckerboard(Even, chi_e, chi);
+  pickCheckerboard(Odd, chi_o, chi);
+  pickCheckerboard(Even, phi_e, phi);
+  pickCheckerboard(Odd, phi_o, phi);
+
+  Dwc.Mooee(src_e, chi_e);
+  Dwc.Mooee(src_o, chi_o);
+  setCheckerboard(chi, chi_e);
+  setCheckerboard(chi, chi_o);
+  setCheckerboard(src, src_e);
+  setCheckerboard(src, src_o);
+
+  ////////////////////// Gauge Transformation
+  std::vector<int> seeds2({5, 6, 7, 8});
+  GridParallelRNG pRNG2(&Grid);
+  pRNG2.SeedFixedIntegers(seeds2);
+  LatticeColourMatrix Omega(&Grid);
+  LatticeColourMatrix ShiftedOmega(&Grid);
+  LatticeGaugeField U_prime(&Grid);
+  U_prime = zero;
+  LatticeColourMatrix U_prime_mu(&Grid);
+  U_prime_mu = zero;
+  SU<Nc>::LieRandomize(pRNG2, Omega, 1.0);
+  for (int mu = 0; mu < Nd; mu++)
+  {
+    U[mu] = peekLorentz(Umu, mu);
+    ShiftedOmega = Cshift(Omega, mu, 1);
+    U_prime_mu = Omega * U[mu] * adj(ShiftedOmega);
+    pokeLorentz(U_prime, U_prime_mu, mu);
+  }
+  /////////////////
+
+  WilsonCloverFermionR Dwc_prime(U_prime, Grid, RBGrid, mass, csw_r, csw_t, anis, params);
+  Dwc_prime.ImportGauge(U_prime);
+
+  tmp = Omega * src;
+  pickCheckerboard(Even, src_e, tmp);
+  pickCheckerboard(Odd, src_o, tmp);
+
+  Dwc_prime.Mooee(src_e, phi_e);
+  Dwc_prime.Mooee(src_o, phi_o);
+
+  setCheckerboard(phi, phi_e);
+  setCheckerboard(phi, phi_o);
+
+  err = chi - adj(Omega) * phi;
+  std::cout << GridLogMessage << "norm diff   " << norm2(err) << std::endl;
+
+  std::cout << GridLogMessage << "=================================================================" << std::endl;
+  std::cout << GridLogMessage << "= Testing gauge covariance Clover term w/o EO preconditioning  " << std::endl;
+  std::cout << GridLogMessage << "================================================================" << std::endl;
+
+  chi = zero;
+  phi = zero;
+
+  WilsonFermionR Dw(Umu, Grid, RBGrid, mass, params);
+  Dw.ImportGauge(Umu);
+
+  Dw.M(src, result);
+  Dwc.M(src, chi);
+
+  Dwc_prime.M(Omega * src, phi);
+
+  WilsonFermionR Dw_prime(U_prime, Grid, RBGrid, mass, params);
+  Dw_prime.ImportGauge(U_prime);
+  Dw_prime.M(Omega * src, result2);
+
+  err = chi - adj(Omega) * phi;
+  err2 = result - adj(Omega) * result2;
+  std::cout << GridLogMessage << "norm diff Wilson   " << norm2(err) << std::endl;
+  std::cout << GridLogMessage << "norm diff WilsonClover  " << norm2(err2) << std::endl;
+
+  std::cout << GridLogMessage << "==========================================================" << std::endl;
+  std::cout << GridLogMessage << "= Testing Mooee(csw=0) Clover to reproduce Mooee Wilson   " << std::endl;
+  std::cout << GridLogMessage << "==========================================================" << std::endl;
+
+  chi = zero;
+  phi = zero;
+  err = zero;
+  WilsonCloverFermionR Dwc_csw0(Umu, Grid, RBGrid, mass, 0.0, 0.0, anis, params); //  <-- Notice: csw=0
+  Dwc_csw0.ImportGauge(Umu);
+
+  pickCheckerboard(Even, phi_e, phi);
+  pickCheckerboard(Odd, phi_o, phi);
+  pickCheckerboard(Even, chi_e, chi);
+  pickCheckerboard(Odd, chi_o, chi);
+
+  Dw.Mooee(src_e, chi_e);
+  Dw.Mooee(src_o, chi_o);
+  Dwc_csw0.Mooee(src_e, phi_e);
+  Dwc_csw0.Mooee(src_o, phi_o);
+
+  setCheckerboard(chi, chi_e);
+  setCheckerboard(chi, chi_o);
+  setCheckerboard(phi, phi_e);
+  setCheckerboard(phi, phi_o);
+  setCheckerboard(src, src_e);
+  setCheckerboard(src, src_o);
+
+  err = chi - phi;
+  std::cout << GridLogMessage << "norm diff  " << norm2(err) << std::endl;
+
+  std::cout << GridLogMessage << "==========================================================" << std::endl;
+  std::cout << GridLogMessage << "= Testing EO operator is equal to the unprec              " << std::endl;
+  std::cout << GridLogMessage << "==========================================================" << std::endl;
+
+  chi = zero;
+  phi = zero;
+  err = zero;
+
+  pickCheckerboard(Even, phi_e, phi);
+  pickCheckerboard(Odd, phi_o, phi);
+  pickCheckerboard(Even, chi_e, chi);
+  pickCheckerboard(Odd, chi_o, chi);
+
+  // M phi = (Mooee src_e + Meooe src_o , Meooe src_e + Mooee src_o)
+
+  Dwc.M(src, ref); // Reference result from the unpreconditioned operator
+
+  // EO matrix
+  Dwc.Mooee(src_e, chi_e); 
+  Dwc.Mooee(src_o, chi_o);
+  Dwc.Meooe(src_o, phi_e);
+  Dwc.Meooe(src_e, phi_o);
+
+  phi_o += chi_o;
+  phi_e += chi_e;
+
+  setCheckerboard(phi, phi_e);
+  setCheckerboard(phi, phi_o);
+
+  err = ref - phi;
+  std::cout << GridLogMessage << "ref (unpreconditioned operator) diff  :" << norm2(ref) << std::endl;
+  std::cout << GridLogMessage << "phi (EO decomposition)          diff  :" << norm2(phi) << std::endl;
+  std::cout << GridLogMessage << "norm diff                             :" << norm2(err) << std::endl;
+
+  Grid_finalize();
+}

tests/forces/Test_wilson_force.cc

@@ -50,7 +50,12 @@ int main (int argc, char ** argv)
   std::vector<int> seeds({1,2,3,4});
 
   GridParallelRNG          pRNG(&Grid);
-  pRNG.SeedFixedIntegers(std::vector<int>({45,12,81,9}));
+  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(std::vector<int>({45,12,81,9}));
 
   LatticeFermion phi        (&Grid); gaussian(pRNG,phi);
   LatticeFermion Mphi       (&Grid); 

tests/forces/Test_wilsonclover_force.cc

@@ -0,0 +1,194 @@
+/*************************************************************************************
+
+    Grid physics library, www.github.com/paboyle/Grid
+
+    Source file: ./tests/Test_wilson_force.cc
+
+    Copyright (C) 2015
+
+Author: Peter Boyle <paboyle@ph.ed.ac.uk>
+
+    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 */
+#include <Grid/Grid.h>
+
+using namespace std;
+using namespace Grid;
+using namespace Grid::QCD;
+
+int main(int argc, char **argv)
+{
+  Grid_init(&argc, &argv);
+
+  std::vector<int> latt_size = GridDefaultLatt();
+  std::vector<int> simd_layout = GridDefaultSimd(Nd, vComplex::Nsimd());
+  std::vector<int> mpi_layout = GridDefaultMpi();
+
+  GridCartesian Grid(latt_size, simd_layout, mpi_layout);
+  GridRedBlackCartesian RBGrid(&Grid);
+
+  int threads = GridThread::GetThreads();
+  std::cout << GridLogMessage << "Grid is setup to use " << threads << " threads" << std::endl;
+
+  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);
+  LatticeFermion Mphi(&Grid);
+  LatticeFermion MphiPrime(&Grid);
+
+  LatticeGaugeField U(&Grid);
+
+  std::vector<int> site = {0, 0, 0, 0};
+  SU3::HotConfiguration(pRNG, U);
+  //SU3::ColdConfiguration(pRNG, U);// Clover term zero
+
+  ////////////////////////////////////
+  // Unmodified matrix element
+  ////////////////////////////////////
+  RealD mass = 0.1;
+  Real csw = 1.0;
+  WilsonCloverFermionR Dw(U, Grid, RBGrid, mass, csw, csw);
+  Dw.ImportGauge(U);
+  Dw.M(phi, Mphi);
+  ComplexD S = innerProduct(Mphi, Mphi); // Action : pdag MdagM p
+
+  // get the deriv of phidag MdagM phi with respect to "U"
+  LatticeGaugeField UdSdU(&Grid);
+  LatticeGaugeField tmp(&Grid);
+
+  ////////////////////////////////////////////
+  Dw.MDeriv(tmp, Mphi, phi, DaggerNo);
+  UdSdU = tmp;
+  Dw.MDeriv(tmp, phi, Mphi, DaggerYes);
+  UdSdU += tmp;
+  /////////////////////////////////////////////
+
+  ////////////////////////////////////
+  // Modify the gauge field a little
+  ////////////////////////////////////
+  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++)
+  {
+    // Traceless antihermitian momentum; gaussian in lie alg
+    SU3::GaussianFundamentalLieAlgebraMatrix(pRNG, mommu);
+    Hmom -= real(sum(trace(mommu * mommu)));
+    PokeIndex<LorentzIndex>(mom, mommu, mu);
+
+    parallel_for(int ss = 0; ss < mom._grid->oSites(); ss++)
+    {
+      Uprime[ss]._internal[mu] = ProjectOnGroup(Exponentiate(mom[ss]._internal[mu], dt, 12) * U[ss]._internal[mu]);
+    }
+  }
+
+  std::cout << GridLogMessage << "Initial mom hamiltonian is " << Hmom << std::endl;
+
+  // New action
+  Dw.ImportGauge(Uprime);
+  Dw.M(phi, MphiPrime);
+  ComplexD Sprime = innerProduct(MphiPrime, MphiPrime);
+
+  //////////////////////////////////////////////
+  // Use derivative to estimate dS
+  //////////////////////////////////////////////
+
+  LatticeComplex dS(&Grid);
+  dS = zero;
+  LatticeComplex dSmom(&Grid);
+  dSmom = zero;
+  LatticeComplex dSmom2(&Grid);
+  dSmom2 = zero;
+
+  for (int mu = 0; mu < Nd; 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" << 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 << " 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 << " dsdumu + dag  " << norm2(mommu) << std::endl;
+    std::cout << "" << std::endl;
+  }
+  /////////////////////////////////////////////////////
+
+  for (int mu = 0; mu < Nd; mu++)
+  {
+    forcemu = PeekIndex<LorentzIndex>(UdSdU, mu);
+    mommu = PeekIndex<LorentzIndex>(mom, mu);
+
+    // Update PF action density
+    dS = dS + trace(mommu * forcemu) * dt;
+
+    dSmom = dSmom - trace(mommu * forcemu) * dt;
+    dSmom2 = dSmom2 - trace(forcemu * forcemu) * (0.25 * dt * dt);
+
+    // Update mom action density
+    mommu = mommu + forcemu * (dt * 0.5);
+
+    Hmomprime -= real(sum(trace(mommu * mommu)));
+  }
+
+  ComplexD dSpred = sum(dS);
+  ComplexD dSm = sum(dSmom);
+  ComplexD dSm2 = sum(dSmom2);
+
+  std::cout << GridLogMessage << "Initial mom hamiltonian is " << Hmom << std::endl;
+  std::cout << GridLogMessage << "Final   mom hamiltonian is " << Hmomprime << std::endl;
+  std::cout << GridLogMessage << "Delta   mom hamiltonian is " << Hmomprime - Hmom << std::endl;
+
+  std::cout << GridLogMessage << " S      " << S << std::endl;
+  std::cout << GridLogMessage << " Sprime " << Sprime << std::endl;
+  std::cout << GridLogMessage << "dS (S' - S)          :" << Sprime - S << 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;
+
+  std::cout << GridLogMessage << "Total dS    " << Hmomprime - Hmom + Sprime - S << std::endl;
+
+  assert(fabs(real(Sprime - S - dSpred)) < 1.0);
+
+  std::cout << GridLogMessage << "Done" << std::endl;
+  Grid_finalize();
+}

tests/hadrons/Test_hadrons_2AS_spectrum.cc

@@ -0,0 +1,168 @@
+/*******************************************************************************
+ Grid physics library, www.github.com/paboyle/Grid
+ 
+ Source file: tests/hadrons/Test_hadrons_spectrum.cc
+ 
+ Copyright (C) 2015
+ 
+ Author: Antonin Portelli <antonin.portelli@me.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.
+ *******************************************************************************/
+
+#include <Grid/Hadrons/Application.hpp>
+
+using namespace Grid;
+using namespace Hadrons;
+
+
+    BEGIN_HADRONS_NAMESPACE
+    BEGIN_MODULE_NAMESPACE(MFermion)
+    MODULE_REGISTER_NS(GaugeProp2AS, TGaugeProp<WilsonTwoIndexAntiSymmetricImplR>, MFermion);
+    END_MODULE_NAMESPACE
+    BEGIN_MODULE_NAMESPACE(MSource)
+    MODULE_REGISTER_NS(Point2AS, TPoint<WilsonTwoIndexAntiSymmetricImplR>, MSource);
+    END_MODULE_NAMESPACE
+    BEGIN_MODULE_NAMESPACE(MContraction)
+    MODULE_REGISTER_NS(Meson2AS, ARG(TMeson<WilsonTwoIndexAntiSymmetricImplR, WilsonTwoIndexAntiSymmetricImplR>), MContraction);
+//    MODULE_REGISTER_NS(BaryonMultirep, ARG(TBaryon<FIMPL, FIMPL, FIMPL>), MContraction);
+    END_MODULE_NAMESPACE
+    BEGIN_MODULE_NAMESPACE(MSink)
+    MODULE_REGISTER_NS(ScalarPoint2AS, TPoint<WilsonTwoIndexAntiSymmetricImplR>, MSink);
+    END_MODULE_NAMESPACE
+    BEGIN_MODULE_NAMESPACE(MSolver)
+    MODULE_REGISTER_NS(RBPrecCG2AS, TRBPrecCG<WilsonTwoIndexAntiSymmetricImplR>, MSolver);
+    END_MODULE_NAMESPACE
+    BEGIN_MODULE_NAMESPACE(MAction)
+    MODULE_REGISTER_NS(WilsonClover2AS, TWilsonClover<WilsonTwoIndexAntiSymmetricImplR>, MAction);
+    END_MODULE_NAMESPACE
+    END_HADRONS_NAMESPACE
+
+
+int main(int argc, char *argv[])
+{
+    // initialization //////////////////////////////////////////////////////////
+    Grid_init(&argc, &argv);
+    HadronsLogError.Active(GridLogError.isActive());
+    HadronsLogWarning.Active(GridLogWarning.isActive());
+    HadronsLogMessage.Active(GridLogMessage.isActive());
+    HadronsLogIterative.Active(GridLogIterative.isActive());
+    HadronsLogDebug.Active(GridLogDebug.isActive());
+    LOG(Message) << "Grid initialized" << std::endl;
+    // run setup ///////////////////////////////////////////////////////////////
+    Application              application;
+    std::vector<std::string> flavour = {"l", "s"};
+    std::vector<double>      mass    = {-0.01, -0.04};
+    double                   csw     = 1.0;
+    // global parameters
+    Application::GlobalPar globalPar;
+    globalPar.trajCounter.start = 1500;
+    globalPar.trajCounter.end   = 1520;
+    globalPar.trajCounter.step  = 20;
+    globalPar.seed              = "1 2 3 4";
+    application.setPar(globalPar);
+    // gauge field
+    application.createModule<MGauge::Unit>("gauge");
+    MSource::Point2AS::Par ptPar;
+    ptPar.position = "0 0 0 0";
+    application.createModule<MSource::Point2AS>("pt", ptPar);
+    // sink
+    MSink::ScalarPoint2AS::Par sinkPar;
+    sinkPar.mom = "0 0 0";
+    application.createModule<MSink::ScalarPoint2AS>("sink", sinkPar);
+    
+    // set fermion boundary conditions to be periodic space, antiperiodic time.
+    std::string boundary = "1 1 1 -1";
+
+    for (unsigned int i = 0; i < flavour.size(); ++i)
+    {
+        // actions
+        MAction::WilsonClover2AS::Par actionPar;
+        actionPar.gauge = "gauge";
+        actionPar.mass  = mass[i];
+        actionPar.csw_r = csw;
+	    actionPar.csw_t = csw;
+        actionPar.clover_anisotropy.isAnisotropic= false;
+        actionPar.clover_anisotropy.t_direction  = Nd-1 ;
+        actionPar.clover_anisotropy.xi_0         = 1.0  ;
+        actionPar.clover_anisotropy.nu           = 1.0  ;
+	    actionPar.boundary = boundary;
+        application.createModule<MAction::WilsonClover2AS>("WilsonClover2AS_" + flavour[i], actionPar);
+        
+        // solvers
+        MSolver::RBPrecCG2AS::Par solverPar;
+        solverPar.action   = "WilsonClover2AS_" + flavour[i];
+        solverPar.residual = 1.0e-8;
+        application.createModule<MSolver::RBPrecCG2AS>("CG_" + flavour[i],
+                                                    solverPar);
+        
+        // propagators
+        MFermion::GaugeProp2AS::Par quarkPar;
+        quarkPar.solver = "CG_" + flavour[i];
+        quarkPar.source = "pt";
+        application.createModule<MFermion::GaugeProp2AS>("Qpt_" + flavour[i], quarkPar);
+        quarkPar.source = "z2";
+        application.createModule<MFermion::GaugeProp2AS>("QZ2_" + flavour[i], quarkPar);
+    }
+    for (unsigned int i = 0; i < flavour.size(); ++i)
+    for (unsigned int j = i; j < flavour.size(); ++j)
+    {
+        MContraction::Meson2AS::Par mesPar;
+        
+        mesPar.output  = "mesons2AS/pt_" + flavour[i] + flavour[j];
+        mesPar.q1      = "Qpt_" + flavour[i];
+        mesPar.q2      = "Qpt_" + flavour[j];
+        mesPar.gammas  = "all";
+        mesPar.sink    = "sink";
+        application.createModule<MContraction::Meson2AS>("meson_pt_"
+                                                      + flavour[i] + flavour[j],
+                                                      mesPar);
+     
+     //   mesPar.output  = "mesons2AS/Z2_" + flavour[i] + flavour[j];
+     //   mesPar.q1      = "QZ2_" + flavour[i];
+     //   mesPar.q2      = "QZ2_" + flavour[j];
+     //   mesPar.gammas  = "all";
+     //   mesPar.sink    = "sink";
+     //   application.createModule<MContraction::Meson2AS>("meson_Z2_"
+     //                                                 + flavour[i] + flavour[j],
+     //                                                 mesPar);
+    }
+    for (unsigned int i = 0; i < flavour.size(); ++i)
+    for (unsigned int j = i; j < flavour.size(); ++j)
+    for (unsigned int k = j; k < flavour.size(); ++k)
+    {
+        MContraction::Baryon::Par barPar;
+        
+        barPar.output = "baryons/pt_" + flavour[i] + flavour[j] + flavour[k];
+        barPar.q1     = "Qpt_" + flavour[i];
+        barPar.q2     = "Qpt_" + flavour[j];
+        barPar.q3     = "Qpt_" + flavour[k];
+        application.createModule<MContraction::Baryon>(
+            "baryon_pt_" + flavour[i] + flavour[j] + flavour[k], barPar);
+    }
+    
+    // execution
+    application.saveParameterFile("spectrum.xml");
+    application.run();
+    
+    // epilogue
+    LOG(Message) << "Grid is finalizing now" << std::endl;
+    Grid_finalize();
+    
+    return EXIT_SUCCESS;
+}

tests/hadrons/Test_hadrons_quark.cc

@@ -130,8 +130,8 @@ int main(int argc, char **argv)
     for (int c = 0; c < Nc; ++c)
     {
         ref = prop;
-        PropToFerm(ferm, prop, s, c);
+        PropToFerm<WilsonImplR>(ferm, prop, s, c);
-        FermToProp(prop, ferm, s, c);
+        FermToProp<WilsonImplR>(prop, ferm, s, c);
 
         std::cout << "Spin = " << s << ", Colour = " << c << std::endl;
         ref -= prop;

tests/hadrons/Test_hadrons_wilsonFund.cc

@@ -0,0 +1,157 @@
+/*******************************************************************************
+ Grid physics library, www.github.com/paboyle/Grid
+ 
+ Source file: tests/hadrons/Test_hadrons_spectrum.cc
+ 
+ Copyright (C) 2015
+ 
+ Author: Antonin Portelli <antonin.portelli@me.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.
+ *******************************************************************************/
+
+#include <Grid/Hadrons/Application.hpp>
+
+using namespace Grid;
+using namespace Hadrons;
+
+int main(int argc, char *argv[])
+{
+    // initialization //////////////////////////////////////////////////////////
+    Grid_init(&argc, &argv);
+    HadronsLogError.Active(GridLogError.isActive());
+    HadronsLogWarning.Active(GridLogWarning.isActive());
+    HadronsLogMessage.Active(GridLogMessage.isActive());
+    HadronsLogIterative.Active(GridLogIterative.isActive());
+    HadronsLogDebug.Active(GridLogDebug.isActive());
+    LOG(Message) << "Grid initialized" << std::endl;
+    
+    // run setup ///////////////////////////////////////////////////////////////
+    Application              application;
+    std::vector<std::string> flavour = {"l"};
+    std::vector<double>      mass    = {-0.1};
+    double                   csw     = 0.0;
+    
+    // global parameters
+    Application::GlobalPar globalPar;
+  
+    globalPar.trajCounter.start = 309;
+    globalPar.trajCounter.end   = 310;
+    globalPar.trajCounter.step  = 1;
+
+    globalPar.seed              = "1 2 3 4";
+
+    application.setPar(globalPar);
+    // gauge field
+    application.createModule<MIO::LoadNersc>("gauge");
+   
+    // sources
+    //MSource::Z2::Par z2Par;
+    //z2Par.tA = 0;
+    //z2Par.tB = 0;
+    //application.createModule<MSource::Z2>("z2", z2Par);
+    MSource::Point::Par ptPar;
+    ptPar.position = "0 0 0 0";
+    application.createModule<MSource::Point>("pt", ptPar);
+    // sink
+    MSink::Point::Par sinkPar;
+    sinkPar.mom = "0 0 0";
+    application.createModule<MSink::ScalarPoint>("sink", sinkPar);
+    
+    // set fermion boundary conditions to be periodic space, antiperiodic time.
+    std::string boundary = "1 1 1 -1";
+
+    for (unsigned int i = 0; i < flavour.size(); ++i)
+    {
+        // actions
+        MAction::WilsonClover::Par actionPar;
+        actionPar.gauge = "gauge";
+        actionPar.mass  = mass[i];
+        actionPar.boundary = boundary;
+        actionPar.csw_r = csw;
+        actionPar.csw_t = csw;
+
+        // !!!!! Check if Anisotropy works  !!!!!
+        actionPar.clover_anisotropy.isAnisotropic= false;
+        actionPar.clover_anisotropy.t_direction  = 3    ;   // Explicit for D=4
+        actionPar.clover_anisotropy.xi_0         = 1.0  ;
+        actionPar.clover_anisotropy.nu           = 1.0  ;
+
+        application.createModule<MAction::WilsonClover>("WilsonClover_" + flavour[i], actionPar);
+        
+        // solvers
+        MSolver::RBPrecCG::Par solverPar;
+        solverPar.action   = "WilsonClover_" + flavour[i];
+        solverPar.residual = 1.0e-8;
+        application.createModule<MSolver::RBPrecCG>("CG_" + flavour[i],
+                                                    solverPar);
+        
+        // propagators
+        MFermion::GaugeProp::Par quarkPar;
+        quarkPar.solver = "CG_" + flavour[i];
+        quarkPar.source = "pt";
+        application.createModule<MFermion::GaugeProp>("Qpt_" + flavour[i], quarkPar);
+ //       quarkPar.source = "z2";
+ //       application.createModule<MFermion::GaugeProp>("QZ2_" + flavour[i], quarkPar);
+    }
+    for (unsigned int i = 0; i < flavour.size(); ++i)
+    for (unsigned int j = i; j < flavour.size(); ++j)
+    {
+        MContraction::Meson::Par mesPar;
+        
+        mesPar.output  = "Fund_mesons/pt_" + flavour[i] + flavour[j];
+        mesPar.q1      = "Qpt_" + flavour[i];
+        mesPar.q2      = "Qpt_" + flavour[j];
+        mesPar.gammas  = "all";
+        mesPar.sink    = "sink";
+        application.createModule<MContraction::Meson>("meson_pt_"
+                                                      + flavour[i] + flavour[j],
+                                                      mesPar);
+   //     mesPar.output  = "mesons/Z2_" + flavour[i] + flavour[j];
+   //     mesPar.q1      = "QZ2_" + flavour[i];
+   //     mesPar.q2      = "QZ2_" + flavour[j];
+   //     mesPar.gammas  = "all";
+   //     mesPar.sink    = "sink";
+   //     application.createModule<MContraction::Meson>("meson_Z2_"
+   //                                                   + flavour[i] + flavour[j],
+   //                                                   mesPar);
+    }
+    for (unsigned int i = 0; i < flavour.size(); ++i)
+    for (unsigned int j = i; j < flavour.size(); ++j)
+    for (unsigned int k = j; k < flavour.size(); ++k)
+    {
+        MContraction::Baryon::Par barPar;
+        
+        barPar.output = "Fund_baryons/pt_" + flavour[i] + flavour[j] + flavour[k];
+        barPar.q1     = "Qpt_" + flavour[i];
+        barPar.q2     = "Qpt_" + flavour[j];
+        barPar.q3     = "Qpt_" + flavour[k];
+        application.createModule<MContraction::Baryon>(
+            "baryon_pt_" + flavour[i] + flavour[j] + flavour[k], barPar);
+    }
+    
+    // execution
+    application.saveParameterFile("WilsonClover_spectrum.xml");
+    application.run();
+    
+    // epilogue
+    LOG(Message) << "Grid is finalizing now" << std::endl;
+    Grid_finalize();
+    
+    return EXIT_SUCCESS;
+}

tests/hmc/Test_hmc_EOWilsonCloverFermionGauge.cc

@@ -0,0 +1,139 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid
+
+Source file: ./tests/Test_hmc_WilsonFermionGauge.cc
+
+Copyright (C) 2016
+
+Author: Guido Cossu <guido.cossu@ed.ac.uk>
+
+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 */
+#include <Grid/Grid.h>
+
+int main(int argc, char **argv) {
+  using namespace Grid;
+  using namespace Grid::QCD;
+
+  Grid_init(&argc, &argv);
+  int threads = GridThread::GetThreads();
+  // here make a routine to print all the relevant information on the run
+  std::cout << GridLogMessage << "Grid is setup to use " << threads << " threads" << std::endl;
+
+   // Typedefs to simplify notation
+  typedef GenericHMCRunner<MinimumNorm2> HMCWrapper;  // Uses the default minimum norm
+  typedef WilsonImplR FermionImplPolicy;
+  typedef WilsonCloverFermionR FermionAction;
+  typedef typename FermionAction::FermionField FermionField;
+
+
+  //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
+  HMCWrapper TheHMC;
+
+  // Grid from the command line
+  TheHMC.Resources.AddFourDimGrid("gauge");
+  // Possibile to create the module by hand 
+  // hardcoding parameters or using a Reader
+
+
+  // Checkpointer definition
+  CheckpointerParameters CPparams;  
+  CPparams.config_prefix = "ckpoint_lat";
+  CPparams.rng_prefix = "ckpoint_rng";
+  CPparams.saveInterval = 5;
+  CPparams.format = "IEEE64BIG";
+  
+  TheHMC.Resources.LoadNerscCheckpointer(CPparams);
+
+  RNGModuleParameters RNGpar;
+  RNGpar.serial_seeds = "1 2 3 4 5";
+  RNGpar.parallel_seeds = "6 7 8 9 10";
+  TheHMC.Resources.SetRNGSeeds(RNGpar);
+
+  // Construct observables
+  // here there is too much indirection 
+  typedef PlaquetteMod<HMCWrapper::ImplPolicy> PlaqObs;
+  TheHMC.Resources.AddObservable<PlaqObs>();
+  //////////////////////////////////////////////
+
+  /////////////////////////////////////////////////////////////
+  // Collect actions, here use more encapsulation
+  // need wrappers of the fermionic classes 
+  // that have a complex construction
+  // standard
+  RealD beta = 5.6 ;
+  WilsonGaugeActionR Waction(beta);
+  
+  // temporarily need a gauge field
+  auto GridPtr = TheHMC.Resources.GetCartesian();
+  auto GridRBPtr = TheHMC.Resources.GetRBCartesian();
+
+  LatticeGaugeField U(GridPtr);
+
+  Real mass = 0.01;
+  Real csw = 1.0;
+
+  FermionAction FermOp(U, *GridPtr, *GridRBPtr, mass, csw);
+
+  ConjugateGradient<FermionField> CG(1.0e-8, 2000);
+
+  TwoFlavourEvenOddPseudoFermionAction<FermionImplPolicy> Nf2(FermOp, CG, CG);
+
+    // Set smearing (true/false), default: false
+  Nf2.is_smeared = false;
+
+
+    // Collect actions
+  ActionLevel<HMCWrapper::Field> Level1(1);
+  Level1.push_back(&Nf2);
+
+  ActionLevel<HMCWrapper::Field> Level2(4);
+  Level2.push_back(&Waction);
+
+  TheHMC.TheAction.push_back(Level1);
+  TheHMC.TheAction.push_back(Level2);
+  /////////////////////////////////////////////////////////////
+
+  /*
+    double rho = 0.1;  // smearing parameter
+    int Nsmear = 2;    // number of smearing levels
+    Smear_Stout<HMCWrapper::ImplPolicy> Stout(rho);
+    SmearedConfiguration<HMCWrapper::ImplPolicy> SmearingPolicy(
+        UGrid, Nsmear, Stout);
+  */
+
+  // HMC parameters are serialisable 
+  TheHMC.Parameters.MD.MDsteps = 20;
+  TheHMC.Parameters.MD.trajL   = 1.0;
+
+  TheHMC.ReadCommandLine(argc, argv); // these can be parameters from file
+  TheHMC.Run();  // no smearing
+  // TheHMC.Run(SmearingPolicy); // for smearing
+
+  Grid_finalize();
+
+} // main
+
+
+
+
+
+
+

tests/hmc/Test_hmc_WC2ASFG_Production.cc

@@ -0,0 +1,213 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid
+
+Source file: ./tests/Test_hmc_WilsonFermionGauge.cc
+
+Copyright (C) 2017
+
+Author: Guido Cossu <guido.cossu@ed.ac.uk>
+
+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 */
+#include <Grid/Grid.h>
+
+
+namespace Grid{
+  struct FermionParameters: Serializable {
+    GRID_SERIALIZABLE_CLASS_MEMBERS(FermionParameters,
+            double, mass,
+            double, csw,
+				    double, StoppingCondition,
+				    int, MaxCGIterations,
+				    bool, ApplySmearing);
+  };
+
+  
+  struct WilsonCloverHMCParameters: Serializable {
+  GRID_SERIALIZABLE_CLASS_MEMBERS(WilsonCloverHMCParameters,
+				  double, gauge_beta,
+				  FermionParameters, WilsonClover)
+
+  template <class ReaderClass >
+  WilsonCloverHMCParameters(Reader<ReaderClass>& Reader){
+    read(Reader, "Action", *this);
+    }
+  };
+
+  struct SmearingParameters: Serializable {
+    GRID_SERIALIZABLE_CLASS_MEMBERS(SmearingParameters,
+				    double, rho,
+				    Integer, Nsmear)
+
+    template <class ReaderClass >
+    SmearingParameters(Reader<ReaderClass>& Reader){
+      read(Reader, "StoutSmearing", *this);
+    }
+
+  };
+  
+  
+}
+
+int main(int argc, char **argv)
+{
+  using namespace Grid;
+  using namespace Grid::QCD;
+
+  typedef Representations< FundamentalRepresentation, TwoIndexAntiSymmetricRepresentation > TheRepresentations;  
+
+  Grid_init(&argc, &argv);
+  int threads = GridThread::GetThreads();
+  // here make a routine to print all the relevant information on the run
+  std::cout << GridLogMessage << "Grid is setup to use " << threads << " threads" << std::endl;
+
+  // Typedefs to simplify notation
+  typedef GenericHMCRunnerHirep<TheRepresentations, MinimumNorm2> HMCWrapper; // Uses the default minimum norm
+  typedef WilsonTwoIndexAntiSymmetricImplR FermionImplPolicy; // gauge field implemetation for the pseudofermions
+  typedef WilsonCloverTwoIndexAntiSymmetricFermionR FermionAction; // type of lattice fermions (Wilson, DW, ...)
+  typedef typename FermionAction::FermionField FermionField;
+  typedef Grid::JSONReader Serialiser;
+
+  //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
+  HMCWrapper TheHMC;
+
+  // Grid from the command line
+  TheHMC.ReadCommandLine(argc, argv); 
+  if (TheHMC.ParameterFile.empty()){
+    std::cout << "Input file not specified."
+              << "Use --ParameterFile option in the command line.\nAborting" 
+              << std::endl;
+    exit(1);
+  }
+  Serialiser Reader(TheHMC.ParameterFile);
+  WilsonCloverHMCParameters MyParams(Reader);  
+
+  // Apply smearing to the fermionic action
+  bool ApplySmearing = MyParams.WilsonClover.ApplySmearing;
+
+  TheHMC.Resources.AddFourDimGrid("gauge");
+
+  // Checkpointer definition
+  CheckpointerParameters CPparams(Reader);
+  
+  /*
+  CPparams.config_prefix = "ckpoint_lat";
+  CPparams.rng_prefix = "ckpoint_rng";
+  CPparams.saveInterval = 5;
+  CPparams.format = "IEEE64BIG";
+  */
+  
+  TheHMC.Resources.LoadNerscCheckpointer(CPparams);
+
+  RNGModuleParameters RNGpar(Reader);
+  /*
+  RNGpar.serial_seeds = "1 2 3 4 5";
+  RNGpar.parallel_seeds = "6 7 8 9 10";
+  TheHMC.Resources.SetRNGSeeds(RNGpar);
+  */
+  TheHMC.Resources.SetRNGSeeds(RNGpar);
+
+  // Construct observables
+  typedef PlaquetteMod<HMCWrapper::ImplPolicy> PlaqObs;
+  TheHMC.Resources.AddObservable<PlaqObs>();
+
+  typedef PolyakovMod<HMCWrapper::ImplPolicy> PolyakovObs;
+  TheHMC.Resources.AddObservable<PolyakovObs>();
+
+  //typedef TopologicalChargeMod<HMCWrapper::ImplPolicy> QObs;
+  //TopologyObsParameters TopParams(Reader);
+  //TheHMC.Resources.AddObservable<QObs>(TopParams);
+  //////////////////////////////////////////////
+
+  /////////////////////////////////////////////////////////////
+  // Collect actions, here use more encapsulation
+  // need wrappers of the fermionic classes
+  // that have a complex construction
+  // standard
+  
+  //RealD beta = 5.6;
+  WilsonGaugeActionR Waction(MyParams.gauge_beta);
+
+  auto GridPtr   = TheHMC.Resources.GetCartesian();
+  auto GridRBPtr = TheHMC.Resources.GetRBCartesian();
+
+  // temporarily need a gauge field
+  TwoIndexAntiSymmetricRepresentation::LatticeField U(GridPtr);
+
+  //Real mass = 0.01;
+  //Real csw = 1.0;
+
+  Real mass = MyParams.WilsonClover.mass;
+  Real csw = MyParams.WilsonClover.csw;
+
+  std::cout << "mass and csw" << mass << " and " << csw << std::endl; 
+
+  FermionAction FermOp(U, *GridPtr, *GridRBPtr, mass, csw, csw);
+  ConjugateGradient<FermionField> CG(MyParams.WilsonClover.StoppingCondition, MyParams.WilsonClover.MaxCGIterations);
+  TwoFlavourPseudoFermionAction<FermionImplPolicy> Nf2(FermOp, CG, CG);
+
+  // Set smearing (true/false), default: false
+  Nf2.is_smeared = ApplySmearing;
+
+  // Collect actions
+  ActionLevel<HMCWrapper::Field, TheRepresentations> Level1(1);
+  Level1.push_back(&Nf2);
+
+  ActionLevel<HMCWrapper::Field, TheRepresentations> Level2(4);
+  Level2.push_back(&Waction);
+
+  TheHMC.TheAction.push_back(Level1);
+  TheHMC.TheAction.push_back(Level2);
+  /////////////////////////////////////////////////////////////
+
+
+  /*
+    double rho = 0.1;  // smearing parameter
+    int Nsmear = 2;    // number of smearing levels
+    Smear_Stout<HMCWrapper::ImplPolicy> Stout(rho);
+    SmearedConfiguration<HMCWrapper::ImplPolicy> SmearingPolicy(
+        UGrid, Nsmear, Stout);
+  */
+
+  // HMC parameters are serialisable
+
+  TheHMC.Parameters.initialize(Reader);
+  //TheHMC.Parameters.MD.MDsteps = 20;
+  //TheHMC.Parameters.MD.trajL = 1.0;
+
+  if (ApplySmearing){
+    SmearingParameters SmPar(Reader);
+    //double rho = 0.1;  // smearing parameter
+    //int Nsmear = 3;    // number of smearing levels
+    Smear_Stout<HMCWrapper::ImplPolicy> Stout(SmPar.rho);
+    SmearedConfiguration<HMCWrapper::ImplPolicy> SmearingPolicy(GridPtr, SmPar.Nsmear, Stout);
+    TheHMC.Run(SmearingPolicy); // for smearing
+  } else {
+    TheHMC.Run();  // no smearing
+  }
+
+  //TheHMC.ReadCommandLine(argc, argv); // these can be parameters from file
+  //TheHMC.Run();                       // no smearing
+  // TheHMC.Run(SmearingPolicy); // for smearing
+
+  Grid_finalize();
+
+} // main
+

tests/hmc/Test_hmc_WC2SFG_Production.cc

@@ -0,0 +1,212 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid
+
+Source file: ./tests/Test_hmc_WilsonFermionGauge.cc
+
+Copyright (C) 2017
+
+Author: Guido Cossu <guido.cossu@ed.ac.uk>
+
+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 */
+#include <Grid/Grid.h>
+
+
+namespace Grid{
+  struct FermionParameters: Serializable {
+    GRID_SERIALIZABLE_CLASS_MEMBERS(FermionParameters,
+            double, mass,
+            double, csw,
+				    double, StoppingCondition,
+				    int, MaxCGIterations,
+				    bool, ApplySmearing);
+  };
+
+  
+  struct WilsonCloverHMCParameters: Serializable {
+  GRID_SERIALIZABLE_CLASS_MEMBERS(WilsonCloverHMCParameters,
+				  double, gauge_beta,
+				  FermionParameters, WilsonClover)
+
+  template <class ReaderClass >
+  WilsonCloverHMCParameters(Reader<ReaderClass>& Reader){
+    read(Reader, "Action", *this);
+    }
+  };
+
+  struct SmearingParameters: Serializable {
+    GRID_SERIALIZABLE_CLASS_MEMBERS(SmearingParameters,
+				    double, rho,
+				    Integer, Nsmear)
+
+    template <class ReaderClass >
+    SmearingParameters(Reader<ReaderClass>& Reader){
+      read(Reader, "StoutSmearing", *this);
+    }
+
+  };
+  
+  
+}
+
+int main(int argc, char **argv)
+{
+  using namespace Grid;
+  using namespace Grid::QCD;
+
+  typedef Representations< FundamentalRepresentation, TwoIndexSymmetricRepresentation > TheRepresentations;  
+
+  Grid_init(&argc, &argv);
+  int threads = GridThread::GetThreads();
+  // here make a routine to print all the relevant information on the run
+  std::cout << GridLogMessage << "Grid is setup to use " << threads << " threads" << std::endl;
+
+  // Typedefs to simplify notation
+  typedef GenericHMCRunnerHirep<TheRepresentations, MinimumNorm2> HMCWrapper; // Uses the default minimum norm
+  typedef WilsonTwoIndexSymmetricImplR FermionImplPolicy; // gauge field implemetation for the pseudofermions
+  typedef WilsonCloverTwoIndexSymmetricFermionR FermionAction; // type of lattice fermions (Wilson, DW, ...)
+  typedef typename FermionAction::FermionField FermionField;
+  typedef Grid::JSONReader Serialiser;
+
+  //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
+  HMCWrapper TheHMC;
+
+  // Grid from the command line
+  TheHMC.ReadCommandLine(argc, argv); 
+  if (TheHMC.ParameterFile.empty()){
+    std::cout << "Input file not specified."
+              << "Use --ParameterFile option in the command line.\nAborting" 
+              << std::endl;
+    exit(1);
+  }
+  Serialiser Reader(TheHMC.ParameterFile);
+  WilsonCloverHMCParameters MyParams(Reader);  
+
+  // Apply smearing to the fermionic action
+  bool ApplySmearing = MyParams.WilsonClover.ApplySmearing;
+
+  TheHMC.Resources.AddFourDimGrid("gauge");
+
+  // Checkpointer definition
+  CheckpointerParameters CPparams(Reader);
+  
+  /*
+  CPparams.config_prefix = "ckpoint_lat";
+  CPparams.rng_prefix = "ckpoint_rng";
+  CPparams.saveInterval = 5;
+  CPparams.format = "IEEE64BIG";
+  */
+  
+  TheHMC.Resources.LoadNerscCheckpointer(CPparams);
+
+  RNGModuleParameters RNGpar(Reader);
+  /*
+  RNGpar.serial_seeds = "1 2 3 4 5";
+  RNGpar.parallel_seeds = "6 7 8 9 10";
+  TheHMC.Resources.SetRNGSeeds(RNGpar);
+  */
+  TheHMC.Resources.SetRNGSeeds(RNGpar);
+
+  // Construct observables
+  typedef PlaquetteMod<HMCWrapper::ImplPolicy> PlaqObs;
+  TheHMC.Resources.AddObservable<PlaqObs>();
+
+  typedef PolyakovMod<HMCWrapper::ImplPolicy> PolyakovObs;
+  TheHMC.Resources.AddObservable<PolyakovObs>();
+
+  //typedef TopologicalChargeMod<HMCWrapper::ImplPolicy> QObs;
+  //TopologyObsParameters TopParams(Reader);
+  //TheHMC.Resources.AddObservable<QObs>(TopParams);
+  //////////////////////////////////////////////
+
+  /////////////////////////////////////////////////////////////
+  // Collect actions, here use more encapsulation
+  // need wrappers of the fermionic classes
+  // that have a complex construction
+  // standard
+  
+  //RealD beta = 5.6;
+  WilsonGaugeActionR Waction(MyParams.gauge_beta);
+
+  auto GridPtr   = TheHMC.Resources.GetCartesian();
+  auto GridRBPtr = TheHMC.Resources.GetRBCartesian();
+
+  // temporarily need a gauge field
+  TwoIndexSymmetricRepresentation::LatticeField U(GridPtr);
+
+  //Real mass = 0.01;
+  //Real csw = 1.0;
+
+  Real mass = MyParams.WilsonClover.mass;
+  Real csw = MyParams.WilsonClover.csw;
+
+  std::cout << "mass and csw" << mass << " and " << csw << std::endl; 
+
+  FermionAction FermOp(U, *GridPtr, *GridRBPtr, mass, csw, csw);
+  ConjugateGradient<FermionField> CG(MyParams.WilsonClover.StoppingCondition, MyParams.WilsonClover.MaxCGIterations);
+  TwoFlavourPseudoFermionAction<FermionImplPolicy> Nf2(FermOp, CG, CG);
+
+  // Set smearing (true/false), default: false
+  Nf2.is_smeared = ApplySmearing;
+
+  // Collect actions
+  ActionLevel<HMCWrapper::Field, TheRepresentations> Level1(1);
+  Level1.push_back(&Nf2);
+
+  ActionLevel<HMCWrapper::Field, TheRepresentations> Level2(4);
+  Level2.push_back(&Waction);
+
+  TheHMC.TheAction.push_back(Level1);
+  TheHMC.TheAction.push_back(Level2);
+  /////////////////////////////////////////////////////////////
+
+
+  /*
+    double rho = 0.1;  // smearing parameter
+    int Nsmear = 2;    // number of smearing levels
+    Smear_Stout<HMCWrapper::ImplPolicy> Stout(rho);
+    SmearedConfiguration<HMCWrapper::ImplPolicy> SmearingPolicy(
+        UGrid, Nsmear, Stout);
+  */
+
+  // HMC parameters are serialisable
+
+  TheHMC.Parameters.initialize(Reader);
+  //TheHMC.Parameters.MD.MDsteps = 20;
+  //TheHMC.Parameters.MD.trajL = 1.0;
+
+  if (ApplySmearing){
+    SmearingParameters SmPar(Reader);
+    //double rho = 0.1;  // smearing parameter
+    //int Nsmear = 3;    // number of smearing levels
+    Smear_Stout<HMCWrapper::ImplPolicy> Stout(SmPar.rho);
+    SmearedConfiguration<HMCWrapper::ImplPolicy> SmearingPolicy(GridPtr, SmPar.Nsmear, Stout);
+    TheHMC.Run(SmearingPolicy); // for smearing
+  } else {
+    TheHMC.Run();  // no smearing
+  }
+
+  //TheHMC.ReadCommandLine(argc, argv); // these can be parameters from file
+  //TheHMC.Run();                       // no smearing
+  // TheHMC.Run(SmearingPolicy); // for smearing
+
+  Grid_finalize();
+
+} // main

tests/hmc/Test_hmc_WCFG_Production.cc

@@ -0,0 +1,210 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid
+
+Source file: ./tests/Test_hmc_WilsonFermionGauge.cc
+
+Copyright (C) 2017
+
+Author: Guido Cossu <guido.cossu@ed.ac.uk>
+
+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 */
+#include <Grid/Grid.h>
+
+
+namespace Grid{
+  struct FermionParameters: Serializable {
+    GRID_SERIALIZABLE_CLASS_MEMBERS(FermionParameters,
+            double, mass,
+            double, csw,
+				    double, StoppingCondition,
+				    int, MaxCGIterations,
+				    bool, ApplySmearing);
+  };
+
+  
+  struct WilsonCloverHMCParameters: Serializable {
+  GRID_SERIALIZABLE_CLASS_MEMBERS(WilsonCloverHMCParameters,
+				  double, gauge_beta,
+				  FermionParameters, WilsonClover)
+
+  template <class ReaderClass >
+  WilsonCloverHMCParameters(Reader<ReaderClass>& Reader){
+    read(Reader, "Action", *this);
+    }
+  };
+
+  struct SmearingParameters: Serializable {
+    GRID_SERIALIZABLE_CLASS_MEMBERS(SmearingParameters,
+				    double, rho,
+				    Integer, Nsmear)
+
+    template <class ReaderClass >
+    SmearingParameters(Reader<ReaderClass>& Reader){
+      read(Reader, "StoutSmearing", *this);
+    }
+
+  };
+  
+  
+}
+
+int main(int argc, char **argv)
+{
+  using namespace Grid;
+  using namespace Grid::QCD;
+
+  Grid_init(&argc, &argv);
+  int threads = GridThread::GetThreads();
+  // here make a routine to print all the relevant information on the run
+  std::cout << GridLogMessage << "Grid is setup to use " << threads << " threads" << std::endl;
+
+  // Typedefs to simplify notation
+  typedef GenericHMCRunner<MinimumNorm2> HMCWrapper; // Uses the default minimum norm
+  typedef WilsonImplR FermionImplPolicy;
+  typedef WilsonCloverFermionR FermionAction;
+  typedef typename FermionAction::FermionField FermionField;
+  typedef Grid::JSONReader Serialiser;
+
+  //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
+  HMCWrapper TheHMC;
+
+  // Grid from the command line
+  TheHMC.ReadCommandLine(argc, argv); 
+  if (TheHMC.ParameterFile.empty()){
+    std::cout << "Input file not specified."
+              << "Use --ParameterFile option in the command line.\nAborting" 
+              << std::endl;
+    exit(1);
+  }
+  Serialiser Reader(TheHMC.ParameterFile);
+  WilsonCloverHMCParameters MyParams(Reader);  
+
+  // Apply smearing to the fermionic action
+  bool ApplySmearing = MyParams.WilsonClover.ApplySmearing;
+
+  TheHMC.Resources.AddFourDimGrid("gauge");
+
+  // Checkpointer definition
+  CheckpointerParameters CPparams(Reader);
+  
+  /*
+  CPparams.config_prefix = "ckpoint_lat";
+  CPparams.rng_prefix = "ckpoint_rng";
+  CPparams.saveInterval = 5;
+  CPparams.format = "IEEE64BIG";
+  */
+  
+  TheHMC.Resources.LoadNerscCheckpointer(CPparams);
+
+  RNGModuleParameters RNGpar(Reader);
+  /*
+  RNGpar.serial_seeds = "1 2 3 4 5";
+  RNGpar.parallel_seeds = "6 7 8 9 10";
+  TheHMC.Resources.SetRNGSeeds(RNGpar);
+  */
+  TheHMC.Resources.SetRNGSeeds(RNGpar);
+
+  // Construct observables
+  typedef PlaquetteMod<HMCWrapper::ImplPolicy> PlaqObs;
+  TheHMC.Resources.AddObservable<PlaqObs>();
+
+  typedef PolyakovMod<HMCWrapper::ImplPolicy> PolyakovObs;
+  TheHMC.Resources.AddObservable<PolyakovObs>();
+
+  //typedef TopologicalChargeMod<HMCWrapper::ImplPolicy> QObs;
+  //TopologyObsParameters TopParams(Reader);
+  //TheHMC.Resources.AddObservable<QObs>(TopParams);
+  //////////////////////////////////////////////
+
+  /////////////////////////////////////////////////////////////
+  // Collect actions, here use more encapsulation
+  // need wrappers of the fermionic classes
+  // that have a complex construction
+  // standard
+  
+  //RealD beta = 5.6;
+  WilsonGaugeActionR Waction(MyParams.gauge_beta);
+
+  auto GridPtr   = TheHMC.Resources.GetCartesian();
+  auto GridRBPtr = TheHMC.Resources.GetRBCartesian();
+
+  // temporarily need a gauge field
+  LatticeGaugeField U(GridPtr);
+
+  //Real mass = 0.01;
+  //Real csw = 1.0;
+
+  Real mass = MyParams.WilsonClover.mass;
+  Real csw = MyParams.WilsonClover.csw;
+
+  std::cout << "mass and csw" << mass << " and " << csw << std::endl; 
+
+  FermionAction FermOp(U, *GridPtr, *GridRBPtr, mass, csw, csw);
+  ConjugateGradient<FermionField> CG(MyParams.WilsonClover.StoppingCondition, MyParams.WilsonClover.MaxCGIterations);
+  TwoFlavourPseudoFermionAction<FermionImplPolicy> Nf2(FermOp, CG, CG);
+
+  // Set smearing (true/false), default: false
+  Nf2.is_smeared = ApplySmearing;
+
+  // Collect actions
+  ActionLevel<HMCWrapper::Field> Level1(1);
+  Level1.push_back(&Nf2);
+
+  ActionLevel<HMCWrapper::Field> Level2(4);
+  Level2.push_back(&Waction);
+
+  TheHMC.TheAction.push_back(Level1);
+  TheHMC.TheAction.push_back(Level2);
+  /////////////////////////////////////////////////////////////
+
+
+  /*
+    double rho = 0.1;  // smearing parameter
+    int Nsmear = 2;    // number of smearing levels
+    Smear_Stout<HMCWrapper::ImplPolicy> Stout(rho);
+    SmearedConfiguration<HMCWrapper::ImplPolicy> SmearingPolicy(
+        UGrid, Nsmear, Stout);
+  */
+
+  // HMC parameters are serialisable
+
+  TheHMC.Parameters.initialize(Reader);
+  //TheHMC.Parameters.MD.MDsteps = 20;
+  //TheHMC.Parameters.MD.trajL = 1.0;
+
+  if (ApplySmearing){
+    SmearingParameters SmPar(Reader);
+    //double rho = 0.1;  // smearing parameter
+    //int Nsmear = 3;    // number of smearing levels
+    Smear_Stout<HMCWrapper::ImplPolicy> Stout(SmPar.rho);
+    SmearedConfiguration<HMCWrapper::ImplPolicy> SmearingPolicy(GridPtr, SmPar.Nsmear, Stout);
+    TheHMC.Run(SmearingPolicy); // for smearing
+  } else {
+    TheHMC.Run();  // no smearing
+  }
+
+  //TheHMC.ReadCommandLine(argc, argv); // these can be parameters from file
+  //TheHMC.Run();                       // no smearing
+  // TheHMC.Run(SmearingPolicy); // for smearing
+
+  Grid_finalize();
+
+} // main

tests/hmc/Test_hmc_WCMixedRepFG_Production.cc

@@ -0,0 +1,224 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid
+
+Source file: ./tests/Test_hmc_WilsonAdjointFermionGauge.cc
+
+Copyright (C) 2015
+
+Author: Peter Boyle <pabobyle@ph.ed.ac.uk>
+Author: Peter Boyle <paboyle@ph.ed.ac.uk>
+Author: neo <cossu@post.kek.jp>
+Author: paboyle <paboyle@ph.ed.ac.uk>
+
+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 */
+#include "Grid/Grid.h"
+
+
+namespace Grid{
+  struct FermionParameters: Serializable {
+    GRID_SERIALIZABLE_CLASS_MEMBERS(FermionParameters,
+            double, mass,
+            double, csw,
+				    double, StoppingCondition,
+				    int, MaxCGIterations,
+				    bool, ApplySmearing);
+  };
+
+  struct WilsonCloverHMCParameters: Serializable {
+  GRID_SERIALIZABLE_CLASS_MEMBERS(WilsonCloverHMCParameters,
+				  double, gauge_beta,
+          FermionParameters, WilsonCloverFund, 
+          FermionParameters, WilsonCloverAS)
+
+  template <class ReaderClass >
+  WilsonCloverHMCParameters(Reader<ReaderClass>& Reader){
+    read(Reader, "Action", *this);
+    }
+  };
+
+  struct SmearingParameters: Serializable {
+    GRID_SERIALIZABLE_CLASS_MEMBERS(SmearingParameters,
+				    double, rho,
+				    Integer, Nsmear)
+
+    template <class ReaderClass >
+    SmearingParameters(Reader<ReaderClass>& Reader){
+      read(Reader, "StoutSmearing", *this);
+    }
+
+  }; 
+}
+
+
+int main(int argc, char **argv) {
+  using namespace Grid;
+  using namespace Grid::QCD;
+
+  // Here change the allowed (higher) representations
+  typedef Representations< FundamentalRepresentation, TwoIndexAntiSymmetricRepresentation> TheRepresentations;
+
+  Grid_init(&argc, &argv);
+  int threads = GridThread::GetThreads();
+  // here make a routine to print all the relevant information on the run
+  std::cout << GridLogMessage << "Grid is setup to use " << threads << " threads" << std::endl;
+
+   // Typedefs to simplify notation
+  typedef GenericHMCRunnerHirep<TheRepresentations, MinimumNorm2> HMCWrapper;
+
+  typedef WilsonImplR FundImplPolicy;
+  typedef WilsonCloverFermionR FundFermionAction; 
+  typedef typename FundFermionAction::FermionField FundFermionField;
+
+  typedef WilsonTwoIndexAntiSymmetricImplR ASymmImplPolicy; 
+  typedef WilsonCloverTwoIndexAntiSymmetricFermionR ASymmFermionAction; 
+  typedef typename ASymmFermionAction::FermionField ASymmFermionField;
+
+  typedef Grid::JSONReader Serialiser;
+  //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
+  HMCWrapper TheHMC;
+  
+    // Grid from the command line
+    TheHMC.ReadCommandLine(argc, argv); 
+    if (TheHMC.ParameterFile.empty()){
+      std::cout << "Input file not specified."
+                << "Use --ParameterFile option in the command line.\nAborting" 
+                << std::endl;
+      exit(1);
+    }
+    Serialiser Reader(TheHMC.ParameterFile);
+    WilsonCloverHMCParameters MyParams(Reader);  
+  
+    // Apply smearing to the fermionic action
+    bool ApplySmearingFund = MyParams.WilsonCloverFund.ApplySmearing;
+    bool ApplySmearingAS = MyParams.WilsonCloverAS.ApplySmearing;
+    
+
+    TheHMC.Resources.AddFourDimGrid("gauge");
+  
+    // Checkpointer definition
+    CheckpointerParameters CPparams(Reader);
+    
+    /*
+    CPparams.config_prefix = "ckpoint_lat";
+    CPparams.rng_prefix = "ckpoint_rng";
+    CPparams.saveInterval = 5;
+    CPparams.format = "IEEE64BIG";
+    */
+    
+    TheHMC.Resources.LoadNerscCheckpointer(CPparams);
+  
+    RNGModuleParameters RNGpar(Reader);
+    /*
+    RNGpar.serial_seeds = "1 2 3 4 5";
+    RNGpar.parallel_seeds = "6 7 8 9 10";
+    TheHMC.Resources.SetRNGSeeds(RNGpar);
+    */
+    TheHMC.Resources.SetRNGSeeds(RNGpar);
+  
+    // Construct observables
+    typedef PlaquetteMod<HMCWrapper::ImplPolicy> PlaqObs;
+    TheHMC.Resources.AddObservable<PlaqObs>();
+  
+    typedef PolyakovMod<HMCWrapper::ImplPolicy> PolyakovObs;
+    TheHMC.Resources.AddObservable<PolyakovObs>();
+  
+    typedef TopologicalChargeMod<HMCWrapper::ImplPolicy> QObs;
+    TopologyObsParameters TopParams(Reader);
+    TheHMC.Resources.AddObservable<QObs>(TopParams);
+    //////////////////////////////////////////////
+  
+    /////////////////////////////////////////////////////////////
+    // Collect actions, here use more encapsulation
+    // need wrappers of the fermionic classes
+    // that have a complex construction
+    // standard
+    
+    //RealD beta = 5.6;
+    WilsonGaugeActionR Waction(MyParams.gauge_beta);
+    
+      auto GridPtr   = TheHMC.Resources.GetCartesian();
+      auto GridRBPtr = TheHMC.Resources.GetRBCartesian();
+    
+      // temporarily need a gauge field
+      FundamentalRepresentation::LatticeField UF(GridPtr);
+      TwoIndexAntiSymmetricRepresentation::LatticeField UAS(GridPtr);
+
+
+      Real Fundmass = MyParams.WilsonCloverFund.mass;
+      Real Fundcsw = MyParams.WilsonCloverFund.csw;
+      Real ASmass = MyParams.WilsonCloverAS.mass;
+      Real AScsw = MyParams.WilsonCloverAS.csw;
+
+      
+
+  std::cout << "Fund: mass and csw" << Fundmass << " and " << Fundcsw << std::endl; 
+  std::cout << "AS  : mass and csw" << ASmass << " and " << AScsw << std::endl; 
+  
+  
+  FundFermionAction FundFermOp(UF, *GridPtr, *GridRBPtr, Fundmass, Fundcsw, Fundcsw);
+  ConjugateGradient<FundFermionField> CG_Fund(MyParams.WilsonCloverFund.StoppingCondition, MyParams.WilsonCloverFund.MaxCGIterations);
+  TwoFlavourPseudoFermionAction<FundImplPolicy> Nf2_Fund(FundFermOp, CG_Fund, CG_Fund);
+
+  ASymmFermionAction ASFermOp(UAS, *GridPtr, *GridRBPtr, ASmass, AScsw, AScsw);
+  ConjugateGradient<ASymmFermionField> CG_AS(MyParams.WilsonCloverAS.StoppingCondition, MyParams.WilsonCloverAS.MaxCGIterations);
+  TwoFlavourPseudoFermionAction<ASymmImplPolicy> Nf2_AS(ASFermOp, CG_AS, CG_AS);
+
+  Nf2_Fund.is_smeared = ApplySmearingFund;
+  Nf2_AS.is_smeared   = ApplySmearingAS;
+  
+
+  // Collect actions
+  ActionLevel<HMCWrapper::Field, TheRepresentations > Level1(1);
+  Level1.push_back(&Nf2_Fund);
+  Level1.push_back(&Nf2_AS);
+
+
+  ActionLevel<HMCWrapper::Field, TheRepresentations > Level2(4);
+  Level2.push_back(&Waction);
+
+  TheHMC.TheAction.push_back(Level1);
+  TheHMC.TheAction.push_back(Level2);
+
+  TheHMC.Parameters.initialize(Reader);
+  //TheHMC.Parameters.MD.MDsteps = 20;
+  //TheHMC.Parameters.MD.trajL = 1.0;
+/*
+  if (ApplySmearingFund || ApplySmearingAS){
+    SmearingParameters SmPar(Reader);
+    //double rho = 0.1;  // smearing parameter
+    //int Nsmear = 3;    // number of smearing levels
+    Smear_Stout<HMCWrapper::ImplPolicy> Stout(SmPar.rho);
+    SmearedConfiguration<HMCWrapper::ImplPolicy> SmearingPolicy(GridPtr, SmPar.Nsmear, Stout);
+    TheHMC.Run(SmearingPolicy); // for smearing
+  } else {
+    TheHMC.Run();  // no smearing
+  }
+*/
+  TheHMC.Run(); 
+
+
+  //TheHMC.ReadCommandLine(argc, argv); // these can be parameters from file
+  //TheHMC.Run();                       // no smearing
+  // TheHMC.Run(SmearingPolicy); // for smearing
+
+  Grid_finalize();
+
+} // main

tests/hmc/Test_hmc_WCadjFG_Production.cc

@@ -0,0 +1,213 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid
+
+Source file: ./tests/Test_hmc_WilsonFermionGauge.cc
+
+Copyright (C) 2017
+
+Author: Guido Cossu <guido.cossu@ed.ac.uk>
+
+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 */
+#include <Grid/Grid.h>
+
+
+namespace Grid{
+  struct FermionParameters: Serializable {
+    GRID_SERIALIZABLE_CLASS_MEMBERS(FermionParameters,
+            double, mass,
+            double, csw,
+				    double, StoppingCondition,
+				    int, MaxCGIterations,
+				    bool, ApplySmearing);
+  };
+
+  
+  struct WilsonCloverHMCParameters: Serializable {
+  GRID_SERIALIZABLE_CLASS_MEMBERS(WilsonCloverHMCParameters,
+				  double, gauge_beta,
+				  FermionParameters, WilsonClover)
+
+  template <class ReaderClass >
+  WilsonCloverHMCParameters(Reader<ReaderClass>& Reader){
+    read(Reader, "Action", *this);
+    }
+  };
+
+  struct SmearingParameters: Serializable {
+    GRID_SERIALIZABLE_CLASS_MEMBERS(SmearingParameters,
+				    double, rho,
+				    Integer, Nsmear)
+
+    template <class ReaderClass >
+    SmearingParameters(Reader<ReaderClass>& Reader){
+      read(Reader, "StoutSmearing", *this);
+    }
+
+  };
+  
+  
+}
+
+int main(int argc, char **argv)
+{
+  using namespace Grid;
+  using namespace Grid::QCD;
+
+  typedef Representations< FundamentalRepresentation, AdjointRepresentation > TheRepresentations;  
+
+  Grid_init(&argc, &argv);
+  int threads = GridThread::GetThreads();
+  // here make a routine to print all the relevant information on the run
+  std::cout << GridLogMessage << "Grid is setup to use " << threads << " threads" << std::endl;
+
+  // Typedefs to simplify notation
+  typedef GenericHMCRunnerHirep<TheRepresentations, MinimumNorm2> HMCWrapper; // Uses the default minimum norm
+  typedef WilsonAdjImplR FermionImplPolicy; // gauge field implemetation for the pseudofermions
+  typedef WilsonCloverAdjFermionR FermionAction; // type of lattice fermions (Wilson, DW, ...)
+  typedef typename FermionAction::FermionField FermionField;
+  typedef Grid::JSONReader Serialiser;
+
+  //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
+  HMCWrapper TheHMC;
+
+  // Grid from the command line
+  TheHMC.ReadCommandLine(argc, argv); 
+  if (TheHMC.ParameterFile.empty()){
+    std::cout << "Input file not specified."
+              << "Use --ParameterFile option in the command line.\nAborting" 
+              << std::endl;
+    exit(1);
+  }
+  Serialiser Reader(TheHMC.ParameterFile);
+  WilsonCloverHMCParameters MyParams(Reader);  
+
+  // Apply smearing to the fermionic action
+  bool ApplySmearing = MyParams.WilsonClover.ApplySmearing;
+
+  TheHMC.Resources.AddFourDimGrid("gauge");
+
+  // Checkpointer definition
+  CheckpointerParameters CPparams(Reader);
+  
+  /*
+  CPparams.config_prefix = "ckpoint_lat";
+  CPparams.rng_prefix = "ckpoint_rng";
+  CPparams.saveInterval = 5;
+  CPparams.format = "IEEE64BIG";
+  */
+  
+  TheHMC.Resources.LoadNerscCheckpointer(CPparams);
+
+  RNGModuleParameters RNGpar(Reader);
+  /*
+  RNGpar.serial_seeds = "1 2 3 4 5";
+  RNGpar.parallel_seeds = "6 7 8 9 10";
+  TheHMC.Resources.SetRNGSeeds(RNGpar);
+  */
+  TheHMC.Resources.SetRNGSeeds(RNGpar);
+
+  // Construct observables
+  typedef PlaquetteMod<HMCWrapper::ImplPolicy> PlaqObs;
+  TheHMC.Resources.AddObservable<PlaqObs>();
+
+  typedef PolyakovMod<HMCWrapper::ImplPolicy> PolyakovObs;
+  TheHMC.Resources.AddObservable<PolyakovObs>();
+
+  typedef TopologicalChargeMod<HMCWrapper::ImplPolicy> QObs;
+  TopologyObsParameters TopParams(Reader);
+  TheHMC.Resources.AddObservable<QObs>(TopParams);
+  //////////////////////////////////////////////
+
+  /////////////////////////////////////////////////////////////
+  // Collect actions, here use more encapsulation
+  // need wrappers of the fermionic classes
+  // that have a complex construction
+  // standard
+  
+  //RealD beta = 5.6;
+  WilsonGaugeActionR Waction(MyParams.gauge_beta);
+
+  auto GridPtr   = TheHMC.Resources.GetCartesian();
+  auto GridRBPtr = TheHMC.Resources.GetRBCartesian();
+
+  // temporarily need a gauge field
+  AdjointRepresentation::LatticeField U(GridPtr);
+
+  //Real mass = 0.01;
+  //Real csw = 1.0;
+
+  Real mass = MyParams.WilsonClover.mass;
+  Real csw = MyParams.WilsonClover.csw;
+
+  std::cout << "mass and csw" << mass << " and " << csw << std::endl; 
+
+  FermionAction FermOp(U, *GridPtr, *GridRBPtr, mass, csw, csw);
+  ConjugateGradient<FermionField> CG(MyParams.WilsonClover.StoppingCondition, MyParams.WilsonClover.MaxCGIterations);
+  TwoFlavourPseudoFermionAction<FermionImplPolicy> Nf2(FermOp, CG, CG);
+
+  // Set smearing (true/false), default: false
+  Nf2.is_smeared = ApplySmearing;
+
+  // Collect actions
+  ActionLevel<HMCWrapper::Field, TheRepresentations> Level1(1);
+  Level1.push_back(&Nf2);
+
+  ActionLevel<HMCWrapper::Field, TheRepresentations> Level2(4);
+  Level2.push_back(&Waction);
+
+  TheHMC.TheAction.push_back(Level1);
+  TheHMC.TheAction.push_back(Level2);
+  /////////////////////////////////////////////////////////////
+
+
+  /*
+    double rho = 0.1;  // smearing parameter
+    int Nsmear = 2;    // number of smearing levels
+    Smear_Stout<HMCWrapper::ImplPolicy> Stout(rho);
+    SmearedConfiguration<HMCWrapper::ImplPolicy> SmearingPolicy(
+        UGrid, Nsmear, Stout);
+  */
+
+  // HMC parameters are serialisable
+
+  TheHMC.Parameters.initialize(Reader);
+  //TheHMC.Parameters.MD.MDsteps = 20;
+  //TheHMC.Parameters.MD.trajL = 1.0;
+
+  if (ApplySmearing){
+    SmearingParameters SmPar(Reader);
+    //double rho = 0.1;  // smearing parameter
+    //int Nsmear = 3;    // number of smearing levels
+    Smear_Stout<HMCWrapper::ImplPolicy> Stout(SmPar.rho);
+    SmearedConfiguration<HMCWrapper::ImplPolicy> SmearingPolicy(GridPtr, SmPar.Nsmear, Stout);
+    TheHMC.Run(SmearingPolicy); // for smearing
+  } else {
+    TheHMC.Run();  // no smearing
+  }
+
+  //TheHMC.ReadCommandLine(argc, argv); // these can be parameters from file
+  //TheHMC.Run();                       // no smearing
+  // TheHMC.Run(SmearingPolicy); // for smearing
+
+  Grid_finalize();
+
+} // main
+

tests/hmc/Test_hmc_WG_Production.cc

@@ -0,0 +1,114 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid
+
+Source file: ./tests/Test_hmc_WilsonFermionGauge.cc
+
+Copyright (C) 2015
+
+Author: Guido Cossu <guido.cossu@ed.ac.uk>
+
+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 */
+#include <Grid/Grid.h>
+
+namespace Grid{
+  struct ActionParameters: Serializable {
+    GRID_SERIALIZABLE_CLASS_MEMBERS(ActionParameters,
+				    double, beta)
+
+
+    template <class ReaderClass >
+    ActionParameters(Reader<ReaderClass>& Reader){
+      read(Reader, "Action", *this);
+    }
+
+  };
+
+}
+
+
+int main(int argc, char **argv) {
+  using namespace Grid;
+  using namespace Grid::QCD;
+
+  Grid_init(&argc, &argv);
+  GridLogLayout();
+
+  // Typedefs to simplify notation
+  typedef GenericHMCRunner<MinimumNorm2> HMCWrapper;  // Uses the default minimum norm
+  HMCWrapper TheHMC;
+  typedef Grid::JSONReader       Serialiser;
+
+  // Grid from the command line
+  TheHMC.Resources.AddFourDimGrid("gauge");
+  TheHMC.ReadCommandLine(argc, argv); // these can be parameters from file
+  // Reader, file should come from command line
+  if (TheHMC.ParameterFile.empty()){
+    std::cout << "Input file not specified."
+	      << "Use --ParameterFile option in the command line.\nAborting"
+	      << std::endl;
+    exit(1);
+  }
+  Serialiser Reader(TheHMC.ParameterFile);
+
+
+
+  // Checkpointer definition
+  CheckpointerParameters CPparams(Reader);
+  TheHMC.Resources.LoadNerscCheckpointer(CPparams);
+
+  RNGModuleParameters RNGpar(Reader);
+  TheHMC.Resources.SetRNGSeeds(RNGpar);
+
+  // Construct observables
+  // here there is too much indirection
+  typedef PlaquetteMod<HMCWrapper::ImplPolicy> PlaqObs;
+  typedef TopologicalChargeMod<HMCWrapper::ImplPolicy> QObs;
+  TheHMC.Resources.AddObservable<PlaqObs>();
+  TopologyObsParameters TopParams(Reader);
+  TheHMC.Resources.AddObservable<QObs>(TopParams);
+  //////////////////////////////////////////////
+
+  /////////////////////////////////////////////////////////////
+  // Collect actions, here use more encapsulation
+  // need wrappers of the fermionic classes
+  // that have a complex construction
+  // standard
+  ActionParameters WilsonPar(Reader);
+  //RealD beta = 6.4 ;
+  WilsonGaugeActionR Waction(WilsonPar.beta);
+
+  ActionLevel<HMCWrapper::Field> Level1(1);
+  Level1.push_back(&Waction);
+  //Level1.push_back(WGMod.getPtr());
+  TheHMC.TheAction.push_back(Level1);
+  /////////////////////////////////////////////////////////////
+
+  // HMC parameters are serialisable
+  TheHMC.Parameters.initialize(Reader);
+
+  //TheHMC.Parameters.MD.MDsteps = 17;
+  //TheHMC.Parameters.MD.trajL   = 1.0;
+
+  TheHMC.Run();  // no smearing
+
+  Grid_finalize();
+
+} // main

tests/hmc/Test_hmc_WilsonCloverFermionGauge.cc

@@ -0,0 +1,129 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid
+
+Source file: ./tests/Test_hmc_WilsonFermionGauge.cc
+
+Copyright (C) 2017
+
+Author: Guido Cossu <guido.cossu@ed.ac.uk>
+
+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 */
+#include <Grid/Grid.h>
+
+int main(int argc, char **argv)
+{
+  using namespace Grid;
+  using namespace Grid::QCD;
+
+  Grid_init(&argc, &argv);
+  int threads = GridThread::GetThreads();
+  // here make a routine to print all the relevant information on the run
+  std::cout << GridLogMessage << "Grid is setup to use " << threads << " threads" << std::endl;
+
+  // Typedefs to simplify notation
+  typedef GenericHMCRunner<MinimumNorm2> HMCWrapper; // Uses the default minimum norm
+  typedef WilsonImplR FermionImplPolicy;
+  typedef WilsonCloverFermionR FermionAction;
+  typedef typename FermionAction::FermionField FermionField;
+
+  //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
+  HMCWrapper TheHMC;
+
+  // Grid from the command line
+  TheHMC.Resources.AddFourDimGrid("gauge");
+
+  // Checkpointer definition
+  CheckpointerParameters CPparams;
+  CPparams.config_prefix = "ckpoint_lat";
+  CPparams.rng_prefix = "ckpoint_rng";
+  CPparams.saveInterval = 5;
+  CPparams.format = "IEEE64BIG";
+
+  TheHMC.Resources.LoadNerscCheckpointer(CPparams);
+
+  RNGModuleParameters RNGpar;
+  RNGpar.serial_seeds = "1 2 3 4 5";
+  RNGpar.parallel_seeds = "6 7 8 9 10";
+  TheHMC.Resources.SetRNGSeeds(RNGpar);
+
+  // Construct observables
+  typedef PlaquetteMod<HMCWrapper::ImplPolicy> PlaqObs;
+  TheHMC.Resources.AddObservable<PlaqObs>();
+
+  typedef PolyakovMod<HMCWrapper::ImplPolicy> PolyakovObs;
+  TheHMC.Resources.AddObservable<PolyakovObs>();
+  //////////////////////////////////////////////
+
+  /////////////////////////////////////////////////////////////
+  // Collect actions, here use more encapsulation
+  // need wrappers of the fermionic classes
+  // that have a complex construction
+  // standard
+  RealD beta = 5.6;
+  WilsonGaugeActionR Waction(beta);
+
+  auto GridPtr = TheHMC.Resources.GetCartesian();
+  auto GridRBPtr = TheHMC.Resources.GetRBCartesian();
+
+  // temporarily need a gauge field
+  LatticeGaugeField U(GridPtr);
+
+  Real mass = 0.01;
+  Real csw = 1.0;
+
+  FermionAction FermOp(U, *GridPtr, *GridRBPtr, mass, csw);
+  ConjugateGradient<FermionField> CG(1.0e-8, 5000);
+
+  TwoFlavourPseudoFermionAction<FermionImplPolicy> Nf2(FermOp, CG, CG);
+
+  // Set smearing (true/false), default: false
+  Nf2.is_smeared = false;
+
+  // Collect actions
+  ActionLevel<HMCWrapper::Field> Level1(1);
+  Level1.push_back(&Nf2);
+
+  ActionLevel<HMCWrapper::Field> Level2(4);
+  Level2.push_back(&Waction);
+
+  TheHMC.TheAction.push_back(Level1);
+  TheHMC.TheAction.push_back(Level2);
+  /////////////////////////////////////////////////////////////
+
+  /*
+    double rho = 0.1;  // smearing parameter
+    int Nsmear = 2;    // number of smearing levels
+    Smear_Stout<HMCWrapper::ImplPolicy> Stout(rho);
+    SmearedConfiguration<HMCWrapper::ImplPolicy> SmearingPolicy(
+        UGrid, Nsmear, Stout);
+  */
+
+  // HMC parameters are serialisable
+  TheHMC.Parameters.MD.MDsteps = 20;
+  TheHMC.Parameters.MD.trajL = 1.0;
+
+  TheHMC.ReadCommandLine(argc, argv); // these can be parameters from file
+  TheHMC.Run();                       // no smearing
+  // TheHMC.Run(SmearingPolicy); // for smearing
+
+  Grid_finalize();
+
+} // main

tests/lanczos/Test_WCMultiRep_lanczos.cc

@@ -0,0 +1,178 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid
+
+Source file: ./tests/Test_dwf_lanczos.cc
+
+Copyright (C) 2015
+
+Author: Peter Boyle <paboyle@ph.ed.ac.uk>
+
+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 */
+#include <Grid/Grid.h>
+
+using namespace std;
+using namespace Grid;
+using namespace Grid::QCD;
+
+
+//typedef WilsonCloverFermionR FermionOp;
+//typedef typename WilsonFermionR::FermionField FermionField;
+
+typedef WilsonImplR FundImplPolicy;
+typedef WilsonCloverFermionR FundFermionAction; 
+typedef typename FundFermionAction::FermionField FundFermionField;
+
+typedef WilsonTwoIndexAntiSymmetricImplR ASymmImplPolicy; 
+typedef WilsonCloverTwoIndexAntiSymmetricFermionR ASymmFermionAction; 
+typedef typename ASymmFermionAction::FermionField ASymmFermionField;
+
+
+RealD AllZero(RealD x) { return 0.; }
+
+int main(int argc, char** argv) {
+  Grid_init(&argc, &argv);
+
+  GridCartesian* UGrid = SpaceTimeGrid::makeFourDimGrid(
+      GridDefaultLatt(), GridDefaultSimd(Nd, vComplex::Nsimd()),
+      GridDefaultMpi());
+  GridRedBlackCartesian* UrbGrid =
+      SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid);
+  GridCartesian* FGrid = UGrid;
+  GridRedBlackCartesian* FrbGrid = UrbGrid;
+  printf("UGrid=%p UrbGrid=%p FGrid=%p FrbGrid=%p\n", UGrid, UrbGrid, FGrid,
+         FrbGrid);
+
+  std::vector<int> seeds4({1, 2, 3, 4});
+  std::vector<int> seeds5({5, 6, 7, 8});
+  GridParallelRNG RNG5(FGrid);
+  RNG5.SeedFixedIntegers(seeds5);
+  GridParallelRNG RNG4(UGrid);
+  RNG4.SeedFixedIntegers(seeds4);
+  GridParallelRNG RNG5rb(FrbGrid);
+  RNG5.SeedFixedIntegers(seeds5);
+
+  GridParallelRNG          pRNG(UGrid); 
+  GridSerialRNG            sRNG; 
+
+  FundamentalRepresentation::LatticeField Umu(UGrid);
+  
+  TwoIndexAntiSymmetricRepresentation HiRep(UGrid);
+  TwoIndexAntiSymmetricRepresentation::LatticeField UmuAS(UGrid);
+
+  
+  CheckpointerParameters CPparams;
+  
+  CPparams.config_prefix = "ckpoint_lat";
+  CPparams.rng_prefix = "ckpoint_rng";
+  CPparams.format = "IEEE64BIG";
+
+//NerscHmcCheckpointer<PeriodicGimplR> Checkpoint(std::string("ckpoint_lat"),
+ //                                                 std::string("ckpoint_rng"), 1);
+ 
+NerscHmcCheckpointer<PeriodicGimplR> Checkpoint(CPparams);
+
+  int CNFGSTART=1;
+  int CNFGEND=2;
+  int CNFGSTEP=1;
+
+  Real Fundmass = -0.1;
+  Real Fundcsw  =  1.0;
+  Real ASmass   = -0.1;
+  Real AScsw    =  1.0;
+    
+  std::cout << "Fund: mass and csw" << Fundmass << " and " << Fundcsw << std::endl; 
+  std::cout << "AS  : mass and csw" << ASmass << " and " << AScsw << std::endl; 
+
+  const int Nstop = 30;
+  const int Nk = 40;
+  const int Np = 40;
+  const int Nm = Nk + Np;
+  const int MaxIt = 10000;
+  RealD resid = 1.0e-8;
+
+    for (int cnfg=CNFGSTART;cnfg<=CNFGEND;cnfg+=CNFGSTEP){
+      Checkpoint.CheckpointRestore(cnfg,Umu, sRNG, pRNG);
+
+  //SU4::HotConfiguration(RNG4, Umu); // temporary, then read.
+  
+  HiRep.update_representation(Umu);
+  UmuAS = HiRep.U;
+
+  FundFermionAction FundFermOp(Umu,*FGrid,*FrbGrid, Fundmass, Fundcsw, Fundcsw);
+  MdagMLinearOperator<FundFermionAction,FundFermionField> HermOpFund(FundFermOp); /// <-----
+  
+  ASymmFermionAction ASFermOp(UmuAS,*FGrid,*FrbGrid, ASmass, AScsw, AScsw);
+  MdagMLinearOperator<ASymmFermionAction,ASymmFermionField> HermOpAS(ASFermOp); /// <-----
+  
+  std::vector<double> Coeffs{0, -1.};
+  Polynomial<FundFermionField> FundPolyX(Coeffs);
+  //Chebyshev<FundFermionField> FundCheb(0.0, 10., 12);
+  
+  FunctionHermOp<FundFermionField> FundPolyXOp(FundPolyX,HermOpFund);
+  PlainHermOp<FundFermionField>    FundOp     (HermOpFund);
+
+  ImplicitlyRestartedLanczos<FundFermionField> IRL_Fund(FundOp, FundPolyXOp, Nstop, Nk, Nm,
+                                               resid, MaxIt);
+  
+  Polynomial<ASymmFermionField> ASPolyX(Coeffs);
+  //Chebyshev<ASymmFermionField> ASCheb(0.0, 10., 12);
+
+  FunctionHermOp<ASymmFermionField> ASPolyXOp(ASPolyX,HermOpAS);
+  PlainHermOp<ASymmFermionField>    ASOp     (HermOpAS);
+
+  ImplicitlyRestartedLanczos<ASymmFermionField> IRL_AS(ASOp, ASPolyXOp, Nstop, Nk, Nm,
+                                               resid, MaxIt);
+                                               
+  std::vector<RealD> Fundeval(Nm);
+  std::vector<RealD> ASeval(Nm);
+
+  FundFermionField Fundsrc(FGrid);
+  ASymmFermionField   ASsrc(FGrid);
+  
+  gaussian(RNG5, Fundsrc);
+  gaussian(RNG5, ASsrc);
+
+  std::vector<FundFermionField> Fundevec(Nm, FGrid);
+  std::vector<ASymmFermionField>   ASevec(Nm, FGrid);
+  
+  for (int i = 0; i < 1; i++) {
+    std::cout << i << " / " << Nm << "Fund: grid pointer " << Fundevec[i]._grid
+              << std::endl;
+  };
+  for (int i = 0; i < 1; i++) {
+    std::cout << i << " / " << Nm << "AS: grid pointer " << ASevec[i]._grid
+              << std::endl;
+  };
+  
+  int FundNconv, ASNconv;
+  IRL_Fund.calc(Fundeval, Fundevec, Fundsrc, FundNconv);
+  IRL_AS.calc(ASeval, ASevec, ASsrc, ASNconv);
+
+      for (int i=0;i<FundNconv;i++){
+      std::cout << "Fund: eval[" << i << "] = " << Fundeval[i] << std::endl;
+    }  
+    for (int i=0;i<ASNconv;i++){
+      std::cout << "2Index: eval[" << i << "] = " << ASeval[i] << std::endl;
+    }  
+    }
+
+  Grid_finalize();
+}

tests/qdpxx/Makefile.am

@@ -1,4 +1,5 @@
 AM_CXXFLAGS += `chroma-config --cxxflags`
-AM_LDFLAGS  += `chroma-config --ldflags` `chroma-config --libs`
+AM_LDFLAGS  += `chroma-config --ldflags`
+LIBS += `chroma-config --libs`
 
 include Make.inc

tests/qdpxx/Makefile.am.qdpxx

@@ -1,6 +1,7 @@
 # additional include paths necessary to compile the C++ library
 
 AM_CXXFLAGS = -I$(top_srcdir)/include `chroma-config --cxxflags`
-AM_LDFLAGS = -L$(top_builddir)/lib `chroma-config --ldflags` `chroma-config --libs`
+AM_LDFLAGS = -L$(top_builddir)/lib `chroma-config --ldflags` 
+AM_LIBS = `chroma-config --libs`
 
 include Make.inc

tests/qdpxx/Test_qdpxx_wilson.cc

@@ -0,0 +1,519 @@
+/*************************************************************************************
+
+    Grid physics library, www.github.com/paboyle/Grid 
+
+    Source file: ./tests/qdpxx/Test_qdpxx_wilson.cc
+
+    Copyright (C) 2017
+
+    Author: Guido Cossu <guido.cossu@ed.ac.uk>
+
+    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 */
+
+#include <Grid/Grid.h>
+#include <chroma.h>
+#include <actions/ferm/invert/syssolver_linop_cg_array.h>
+#include <actions/ferm/invert/syssolver_linop_aggregate.h>
+
+// Mass
+double mq = 0.1;
+
+// Define Wilson Types
+typedef Grid::QCD::WilsonImplR::FermionField FermionField;
+typedef Grid::QCD::LatticeGaugeField GaugeField;
+
+enum ChromaAction
+{
+  Wilson,      // Wilson
+  WilsonClover // CloverFermions
+};
+
+namespace Chroma
+{
+
+class ChromaWrapper
+{
+public:
+  typedef multi1d<LatticeColorMatrix> U;
+  typedef LatticeFermion T4;
+
+  static void ImportGauge(GaugeField &gr,
+                          QDP::multi1d<QDP::LatticeColorMatrix> &ch)
+  {
+    Grid::QCD::LorentzColourMatrix LCM;
+    Grid::Complex cc;
+    QDP::ColorMatrix cm;
+    QDP::Complex c;
+
+    std::vector<int> x(4);
+    QDP::multi1d<int> cx(4);
+    std::vector<int> gd = gr._grid->GlobalDimensions();
+
+    for (x[0] = 0; x[0] < gd[0]; x[0]++)
+    {
+      for (x[1] = 0; x[1] < gd[1]; x[1]++)
+      {
+        for (x[2] = 0; x[2] < gd[2]; x[2]++)
+        {
+          for (x[3] = 0; x[3] < gd[3]; x[3]++)
+          {
+            cx[0] = x[0];
+            cx[1] = x[1];
+            cx[2] = x[2];
+            cx[3] = x[3];
+            Grid::peekSite(LCM, gr, x);
+
+            for (int mu = 0; mu < 4; mu++)
+            {
+              for (int i = 0; i < 3; i++)
+              {
+                for (int j = 0; j < 3; j++)
+                {
+                  cc = LCM(mu)()(i, j);
+                  c = QDP::cmplx(QDP::Real(real(cc)), QDP::Real(imag(cc)));
+                  QDP::pokeColor(cm, c, i, j);
+                }
+              }
+              QDP::pokeSite(ch[mu], cm, cx);
+            }
+          }
+        }
+      }
+    }
+  }
+
+  static void ExportGauge(GaugeField &gr,
+                          QDP::multi1d<QDP::LatticeColorMatrix> &ch)
+  {
+    Grid::QCD::LorentzColourMatrix LCM;
+    Grid::Complex cc;
+    QDP::ColorMatrix cm;
+    QDP::Complex c;
+
+    std::vector<int> x(4);
+    QDP::multi1d<int> cx(4);
+    std::vector<int> gd = gr._grid->GlobalDimensions();
+
+    for (x[0] = 0; x[0] < gd[0]; x[0]++)
+    {
+      for (x[1] = 0; x[1] < gd[1]; x[1]++)
+      {
+        for (x[2] = 0; x[2] < gd[2]; x[2]++)
+        {
+          for (x[3] = 0; x[3] < gd[3]; x[3]++)
+          {
+            cx[0] = x[0];
+            cx[1] = x[1];
+            cx[2] = x[2];
+            cx[3] = x[3];
+
+            for (int mu = 0; mu < 4; mu++)
+            {
+              for (int i = 0; i < 3; i++)
+              {
+                for (int j = 0; j < 3; j++)
+                {
+                  cm = QDP::peekSite(ch[mu], cx);
+                  c = QDP::peekColor(cm, i, j);
+                  cc = Grid::Complex(toDouble(real(c)), toDouble(imag(c)));
+                  LCM(mu)
+                  ()(i, j) = cc;
+                }
+              }
+            }
+            Grid::pokeSite(LCM, gr, x);
+          }
+        }
+      }
+    }
+  }
+
+  // Specific for Wilson Fermions
+  static void ImportFermion(Grid::QCD::LatticeFermion &gr,
+                            QDP::LatticeFermion &ch)
+  {
+    Grid::QCD::SpinColourVector F;
+    Grid::Complex c;
+
+    QDP::Fermion cF;
+    QDP::SpinVector cS;
+    QDP::Complex cc;
+
+    std::vector<int> x(4); // explicit 4d fermions in Grid
+    QDP::multi1d<int> cx(4);
+    std::vector<int> gd = gr._grid->GlobalDimensions();
+
+    for (x[0] = 0; x[0] < gd[0]; x[0]++)
+    {
+      for (x[1] = 0; x[1] < gd[1]; x[1]++)
+      {
+        for (x[2] = 0; x[2] < gd[2]; x[2]++)
+        {
+          for (x[3] = 0; x[3] < gd[3]; x[3]++)
+          {
+            cx[0] = x[0];
+            cx[1] = x[1];
+            cx[2] = x[2];
+            cx[3] = x[3];
+
+            Grid::peekSite(F, gr, x);
+
+            for (int j = 0; j < 3; j++)
+            {
+              for (int sp = 0; sp < 4; sp++)
+              {
+
+                c = F()(sp)(j);
+
+                cc = QDP::cmplx(QDP::Real(real(c)), QDP::Real(imag(c)));
+
+                QDP::pokeSpin(cS, cc, sp);
+              }
+              QDP::pokeColor(cF, cS, j);
+            }
+            QDP::pokeSite(ch, cF, cx);
+          }
+        }
+      }
+    }
+  }
+
+  // Specific for 4d Wilson fermions
+  static void ExportFermion(Grid::QCD::LatticeFermion &gr,
+                            QDP::LatticeFermion &ch)
+  {
+    Grid::QCD::SpinColourVector F;
+    Grid::Complex c;
+
+    QDP::Fermion cF;
+    QDP::SpinVector cS;
+    QDP::Complex cc;
+
+    std::vector<int> x(4); // 4d fermions
+    QDP::multi1d<int> cx(4);
+    std::vector<int> gd = gr._grid->GlobalDimensions();
+
+    for (x[0] = 0; x[0] < gd[0]; x[0]++)
+    {
+      for (x[1] = 0; x[1] < gd[1]; x[1]++)
+      {
+        for (x[2] = 0; x[2] < gd[2]; x[2]++)
+        {
+          for (x[3] = 0; x[3] < gd[3]; x[3]++)
+          {
+            cx[0] = x[0];
+            cx[1] = x[1];
+            cx[2] = x[2];
+            cx[3] = x[3];
+
+            cF = QDP::peekSite(ch, cx);
+            for (int sp = 0; sp < 4; sp++)
+            {
+              for (int j = 0; j < 3; j++)
+              {
+                cS = QDP::peekColor(cF, j);
+                cc = QDP::peekSpin(cS, sp);
+                c = Grid::Complex(QDP::toDouble(QDP::real(cc)),
+                                  QDP::toDouble(QDP::imag(cc)));
+                F()
+                (sp)(j) = c;
+              }
+            }
+            Grid::pokeSite(F, gr, x);
+          }
+        }
+      }
+    }
+  }
+
+  static Handle<Chroma::UnprecLinearOperator<T4, U, U>> GetLinOp(U &u, ChromaAction params)
+  {
+    QDP::Real _mq(mq);
+    QDP::multi1d<int> bcs(QDP::Nd);
+
+    // Boundary conditions
+    bcs[0] = bcs[1] = bcs[2] = bcs[3] = 1;
+
+    if (params == Wilson)
+    {
+
+      Chroma::WilsonFermActParams p;
+      p.Mass = _mq;
+      AnisoParam_t _apar;
+      _apar.anisoP = true;
+      _apar.t_dir = 3; // in 4d
+      _apar.xi_0 = 2.0;
+      _apar.nu = 1.0;
+      p.anisoParam = _apar;
+
+      Chroma::Handle<Chroma::FermBC<T4, U, U>> fbc(new Chroma::SimpleFermBC<T4, U, U>(bcs));
+      Chroma::Handle<Chroma::CreateFermState<T4, U, U>> cfs(new Chroma::CreateSimpleFermState<T4, U, U>(fbc));
+      Chroma::UnprecWilsonFermAct S_f(cfs, p);
+      Chroma::Handle<Chroma::FermState<T4, U, U>> ffs(S_f.createState(u));
+      return S_f.linOp(ffs);
+    }
+
+    if (params == WilsonClover)
+    {
+      Chroma::CloverFermActParams p;
+      p.Mass = _mq;
+      p.clovCoeffR = QDP::Real(1.0);
+      p.clovCoeffT = QDP::Real(2.0);
+      p.u0 = QDP::Real(1.0);
+      AnisoParam_t _apar;
+      _apar.anisoP = true;
+      _apar.t_dir = 3; // in 4d
+      _apar.xi_0 = 2.0;
+      _apar.nu = 1.0;
+      p.anisoParam = _apar;
+
+      Chroma::Handle<Chroma::FermBC<T4, U, U>> fbc(new Chroma::SimpleFermBC<T4, U, U>(bcs));
+      Chroma::Handle<Chroma::CreateFermState<T4, U, U>> cfs(new Chroma::CreateSimpleFermState<T4, U, U>(fbc));
+      Chroma::UnprecCloverFermAct S_f(cfs, p);
+      Chroma::Handle<Chroma::FermState<T4, U, U>> ffs(S_f.createState(u));
+      return S_f.linOp(ffs);
+    }
+  }
+};
+} // namespace Chroma
+
+void calc_chroma(ChromaAction action, GaugeField &lat, FermionField &src, FermionField &res, int dag)
+{
+  QDP::multi1d<QDP::LatticeColorMatrix> u(4);
+  Chroma::ChromaWrapper::ImportGauge(lat, u);
+
+  QDP::LatticeFermion check;
+  QDP::LatticeFermion result;
+  QDP::LatticeFermion psi;
+
+  Chroma::ChromaWrapper::ImportFermion(src, psi);
+
+  for (int mu = 0; mu < 4; mu++)
+  {
+    std::cout << "Imported Gauge norm [" << mu << "] " << QDP::norm2(u[mu]) << std::endl;
+  }
+  std::cout << "Imported Fermion norm " << QDP::norm2(psi) << std::endl;
+
+  typedef QDP::LatticeFermion T;
+  typedef QDP::multi1d<QDP::LatticeColorMatrix> U;
+
+  auto linop = Chroma::ChromaWrapper::GetLinOp(u, action);
+
+  printf("Calling Chroma Linop\n");
+  fflush(stdout);
+
+  if (dag)
+    (*linop)(check, psi, Chroma::MINUS);
+  else
+    (*linop)(check, psi, Chroma::PLUS);
+
+  printf("Called Chroma Linop\n");
+  fflush(stdout);
+
+  // std::cout << "Calling Chroma Linop " << std::endl;
+  // linop->evenEvenLinOp(tmp, psi, isign);
+  // check[rb[0]] = tmp;
+  // linop->oddOddLinOp(tmp, psi, isign);
+  // check[rb[1]] = tmp;
+  // linop->evenOddLinOp(tmp, psi, isign);
+  // check[rb[0]] += tmp;
+  // linop->oddEvenLinOp(tmp, psi, isign);
+  // check[rb[1]] += tmp;
+
+  Chroma::ChromaWrapper::ExportFermion(res, check);
+}
+
+void make_gauge(GaugeField &Umu, FermionField &src)
+{
+  using namespace Grid;
+  using namespace Grid::QCD;
+
+  std::vector<int> seeds4({1, 2, 3, 4});
+
+  Grid::GridCartesian *UGrid = (Grid::GridCartesian *)Umu._grid;
+  Grid::GridParallelRNG RNG4(UGrid);
+  RNG4.SeedFixedIntegers(seeds4);
+  Grid::QCD::SU3::HotConfiguration(RNG4, Umu);
+
+  // Fermion field
+  Grid::gaussian(RNG4, src);
+  /*
+  Grid::QCD::SpinColourVector F;
+  Grid::Complex c;
+
+  
+
+  std::vector<int> x(4); // 4d fermions
+  std::vector<int> gd = src._grid->GlobalDimensions();
+
+  for (x[0] = 0; x[0] < gd[0]; x[0]++)
+  {
+    for (x[1] = 0; x[1] < gd[1]; x[1]++)
+    {
+      for (x[2] = 0; x[2] < gd[2]; x[2]++)
+      {
+        for (x[3] = 0; x[3] < gd[3]; 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)|| (sp==3)) && (j==2))
+              {
+                c = Grid::Complex(1.0, 0.0);
+                F()(sp)(j) = c;
+              }
+            }
+          }
+          Grid::pokeSite(F, src, x);
+          
+        }
+      }
+    }
+  }
+  */
+}
+
+void calc_grid(ChromaAction action, Grid::QCD::LatticeGaugeField &Umu, Grid::QCD::LatticeFermion &src, Grid::QCD::LatticeFermion &res, int dag)
+{
+  using namespace Grid;
+  using namespace Grid::QCD;
+
+  Grid::GridCartesian *UGrid = (Grid::GridCartesian *)Umu._grid;
+  Grid::GridRedBlackCartesian *UrbGrid = Grid::QCD::SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid);
+
+  Grid::RealD _mass = mq;
+
+  if (action == Wilson)
+  {
+    WilsonAnisotropyCoefficients anis;
+    anis.isAnisotropic = true;
+    anis.t_direction = 3;
+    anis.xi_0 = 2.0;
+    anis.nu = 1.0;
+    WilsonImplParams iParam;
+    Grid::QCD::WilsonFermionR Wf(Umu, *UGrid, *UrbGrid, _mass, iParam, anis);
+
+    std::cout << Grid::GridLogMessage << " Calling Grid Wilson Fermion multiply " << std::endl;
+
+    if (dag)
+      Wf.Mdag(src, res);
+    else
+      Wf.M(src, res);
+    return;
+  }
+
+  if (action == WilsonClover)
+  {
+    Grid::RealD _csw_r = 1.0;
+    Grid::RealD _csw_t = 2.0;
+    WilsonAnisotropyCoefficients anis;
+    anis.isAnisotropic = true;
+    anis.t_direction = 3;
+    anis.xi_0 = 2.0;
+    anis.nu = 1.0;
+    WilsonImplParams CloverImplParam;
+    Grid::QCD::WilsonCloverFermionR Wf(Umu, *UGrid, *UrbGrid, _mass, _csw_r, _csw_t, anis, CloverImplParam);
+    Wf.ImportGauge(Umu);
+
+    std::cout << Grid::GridLogMessage << " Calling Grid Wilson Clover Fermion multiply " << std::endl;
+
+    if (dag)
+      Wf.Mdag(src, res);
+    else
+      Wf.M(src, res);
+    return;
+  }
+
+  assert(0);
+}
+
+int main(int argc, char **argv)
+{
+
+  /********************************************************
+   * Setup QDP
+   *********************************************************/
+  Chroma::initialize(&argc, &argv);
+  Chroma::WilsonTypeFermActs4DEnv::registerAll();
+
+  /********************************************************
+   * Setup Grid
+   *********************************************************/
+  Grid::Grid_init(&argc, &argv);
+  Grid::GridCartesian *UGrid = Grid::QCD::SpaceTimeGrid::makeFourDimGrid(Grid::GridDefaultLatt(),
+                                                                         Grid::GridDefaultSimd(Grid::QCD::Nd, Grid::vComplex::Nsimd()),
+                                                                         Grid::GridDefaultMpi());
+
+  std::vector<int> gd = UGrid->GlobalDimensions();
+  QDP::multi1d<int> nrow(QDP::Nd);
+  for (int mu = 0; mu < 4; mu++)
+    nrow[mu] = gd[mu];
+
+  QDP::Layout::setLattSize(nrow);
+  QDP::Layout::create();
+
+  GaugeField Ug(UGrid);
+  FermionField src(UGrid);
+  FermionField res_chroma(UGrid);
+  FermionField res_grid(UGrid);
+  FermionField only_wilson(UGrid);
+  FermionField difference(UGrid);
+
+  std::vector<ChromaAction> ActionList({Wilson, WilsonClover});
+  std::vector<std::string> ActionName({"Wilson", "WilsonClover"});
+
+  {
+
+    for (int i = 0; i < ActionList.size(); i++)
+    {
+      std::cout << "*****************************" << std::endl;
+      std::cout << "Action " << ActionName[i] << std::endl;
+      std::cout << "*****************************" << std::endl;
+      make_gauge(Ug, src); // fills the gauge field and the fermion field with random numbers
+
+      for (int dag = 0; dag < 2; dag++)
+      {
+
+        {
+
+          std::cout << "Dag =  " << dag << std::endl;
+
+          calc_chroma(ActionList[i], Ug, src, res_chroma, dag);
+
+          // Remove the normalisation of Chroma Gauge links ????????
+          std::cout << "Norm of Chroma " << ActionName[i] << " multiply " << Grid::norm2(res_chroma) << std::endl;
+          calc_grid(ActionList[i], Ug, src, res_grid, dag);
+
+          std::cout << "Norm of gauge " << Grid::norm2(Ug) << std::endl;
+
+          std::cout << "Norm of Grid " << ActionName[i] << " multiply " << Grid::norm2(res_grid) << std::endl;
+
+          difference = res_chroma - res_grid;
+          std::cout << "Norm of difference " << Grid::norm2(difference) << std::endl;
+        }
+      }
+
+      std::cout << "Finished test " << std::endl;
+
+      Chroma::finalize();
+    }
+  }
+}