Merge pull request #136 from pretidav/feature/clover

Feature/clover

extras/Hadrons/Modules.hpp

@@ -1,25 +1,27 @@
-#include <Grid/Hadrons/Modules/MAction/DWF.hpp>
-#include <Grid/Hadrons/Modules/MAction/Wilson.hpp>
-#include <Grid/Hadrons/Modules/MContraction/Baryon.hpp>
-#include <Grid/Hadrons/Modules/MContraction/DiscLoop.hpp>
-#include <Grid/Hadrons/Modules/MContraction/Gamma3pt.hpp>
-#include <Grid/Hadrons/Modules/MContraction/Meson.hpp>
-#include <Grid/Hadrons/Modules/MContraction/WeakHamiltonian.hpp>
-#include <Grid/Hadrons/Modules/MContraction/WeakHamiltonianEye.hpp>
-#include <Grid/Hadrons/Modules/MContraction/WeakHamiltonianNonEye.hpp>
-#include <Grid/Hadrons/Modules/MContraction/WeakNeutral4ptDisc.hpp>
-#include <Grid/Hadrons/Modules/MFermion/GaugeProp.hpp>
-#include <Grid/Hadrons/Modules/MGauge/Load.hpp>
-#include <Grid/Hadrons/Modules/MGauge/Random.hpp>
-#include <Grid/Hadrons/Modules/MGauge/StochEm.hpp>
-#include <Grid/Hadrons/Modules/MGauge/Unit.hpp>
 #include <Grid/Hadrons/Modules/MLoop/NoiseLoop.hpp>
 #include <Grid/Hadrons/Modules/MScalar/ChargedProp.hpp>
-#include <Grid/Hadrons/Modules/MScalar/FreeProp.hpp>
 #include <Grid/Hadrons/Modules/MScalar/Scalar.hpp>
+#include <Grid/Hadrons/Modules/MScalar/FreeProp.hpp>
+#include <Grid/Hadrons/Modules/MSource/Wall.hpp>
+#include <Grid/Hadrons/Modules/MSource/SeqGamma.hpp>
+#include <Grid/Hadrons/Modules/MSource/Point.hpp>
+#include <Grid/Hadrons/Modules/MSource/Z2.hpp>
+#include <Grid/Hadrons/Modules/MFermion/GaugeProp.hpp>
+#include <Grid/Hadrons/Modules/MContraction/Meson.hpp>
+#include <Grid/Hadrons/Modules/MContraction/WeakHamiltonianNonEye.hpp>
+#include <Grid/Hadrons/Modules/MContraction/WeakHamiltonianEye.hpp>
+#include <Grid/Hadrons/Modules/MContraction/DiscLoop.hpp>
+#include <Grid/Hadrons/Modules/MContraction/Baryon.hpp>
+#include <Grid/Hadrons/Modules/MContraction/Gamma3pt.hpp>
+#include <Grid/Hadrons/Modules/MContraction/WeakNeutral4ptDisc.hpp>
+#include <Grid/Hadrons/Modules/MContraction/WeakHamiltonian.hpp>
 #include <Grid/Hadrons/Modules/MSink/Point.hpp>
 #include <Grid/Hadrons/Modules/MSolver/RBPrecCG.hpp>
-#include <Grid/Hadrons/Modules/MSource/Point.hpp>
+#include <Grid/Hadrons/Modules/MGauge/StochEm.hpp>
-#include <Grid/Hadrons/Modules/MSource/SeqGamma.hpp>
+#include <Grid/Hadrons/Modules/MGauge/FundtoHirep.hpp>
-#include <Grid/Hadrons/Modules/MSource/Wall.hpp>
+#include <Grid/Hadrons/Modules/MGauge/Unit.hpp>
-#include <Grid/Hadrons/Modules/MSource/Z2.hpp>
+#include <Grid/Hadrons/Modules/MGauge/Load.hpp>
+#include <Grid/Hadrons/Modules/MGauge/Random.hpp>
+#include <Grid/Hadrons/Modules/MAction/WilsonClover.hpp>
+#include <Grid/Hadrons/Modules/MAction/DWF.hpp>
+#include <Grid/Hadrons/Modules/MAction/Wilson.hpp>

extras/Hadrons/Modules/MAction/WilsonClover.hpp

@@ -0,0 +1,142 @@
+/*************************************************************************************
+
+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);
+}
+
+// execution ///////////////////////////////////////////////////////////////////
+template <typename FImpl>
+void TWilsonClover<FImpl>::execute()
+{
+    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      = *env().template getObject<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);
+    FMat *fMatPt = new WilsonCloverFermion<FImpl>(U, grid, gridRb, par().mass,
+						  par().csw_r,
+						  par().csw_t,
+					      par().clover_anisotropy,
+						  implParams);
+    env().setObject(getName(), fMatPt);
+}
+
+END_MODULE_NAMESPACE
+
+END_HADRONS_NAMESPACE
+
+#endif // Hadrons_WilsonClover_hpp_

extras/Hadrons/Modules/MContraction/Baryon.hpp

@@ -112,7 +112,9 @@ void TBaryon<FImpl1, FImpl2, FImpl3>::execute(void)
                  << " quarks '" << par().q1 << "', '" << par().q2 << "', and '"
                  << par().q3 << "'" << std::endl;
     
-    CorrWriter             writer(par().output);
+    std::string output_name = par().output + "." + std::to_string(env().getTrajectory());
+    
+    CorrWriter             writer(output_name);
     PropagatorField1      &q1 = *env().template getObject<PropagatorField1>(par().q1);
     PropagatorField2      &q2 = *env().template getObject<PropagatorField2>(par().q2);
     PropagatorField3      &q3 = *env().template getObject<PropagatorField3>(par().q2);

extras/Hadrons/Modules/MContraction/Meson.hpp

@@ -165,8 +165,10 @@ void TMeson<FImpl1, FImpl2>::execute(void)
     LOG(Message) << "Computing meson contractions '" << getName() << "' using"
                  << " quarks '" << par().q1 << "' and '" << par().q2 << "'"
                  << std::endl;
+
+    std::string output_name = par().output + "." + std::to_string(env().getTrajectory());
     
-    CorrWriter             writer(par().output);
+    CorrWriter             writer(output_name);
     std::vector<TComplex>  buf;
     std::vector<Result>    result;
     Gamma                  g5(Gamma::Algebra::Gamma5);

extras/Hadrons/Modules/MFermion/GaugeProp.hpp

@@ -43,7 +43,6 @@ private:
 };
 
 MODULE_REGISTER_NS(GaugeProp, TGaugeProp<FIMPL>, MFermion);
-
 /******************************************************************************
  *                      TGaugeProp implementation                             *
  ******************************************************************************/
@@ -103,7 +102,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;
@@ -112,12 +111,12 @@ void TGaugeProp<FImpl>::execute(void)
         {
             if (Ls_ == 1)
             {
-                PropToFerm(source, fullSrc, s, c);
+               PropToFerm<FImpl>(source, fullSrc, s, c);
             }
             else
             {
                 source = zero;
-                PropToFerm(tmp, fullSrc, s, c);
+                PropToFerm<FImpl>(tmp, fullSrc, s, c);
                 InsertSlice(tmp, source, 0, 0);
                 InsertSlice(tmp, source, Ls_-1, 0);
                 axpby_ssp_pplus(source, 0., source, 1., source, 0, 0);
@@ -133,12 +132,12 @@ 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)
         {
@@ -148,7 +147,7 @@ void TGaugeProp<FImpl>::execute(void)
             axpby_ssp_pminus(sol, 0., sol, 1., sol, 0, 0);
             axpby_ssp_pplus(sol, 1., sol, 1., sol, 0, Ls_-1);
             ExtractSlice(tmp, sol, 0, 0);
-            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/MGauge/Load.cc

@@ -66,7 +66,8 @@ void TLoad::setup(void)
 void TLoad::execute(void)
 {
     FieldMetaData  header;
-    std::string fileName = par().file + "."
+
+    std::string fileName = par().file + "ckpoint_lat."
                            + std::to_string(env().getTrajectory());
     
     LOG(Message) << "Loading NERSC configuration from file '" << fileName

extras/Hadrons/modules.inc

@@ -1,38 +1,41 @@
 modules_cc =\
+  Modules/MScalar/ChargedProp.cc \
+  Modules/MScalar/FreeProp.cc \
   Modules/MContraction/WeakHamiltonianEye.cc \
   Modules/MContraction/WeakHamiltonianNonEye.cc \
   Modules/MContraction/WeakNeutral4ptDisc.cc \
-  Modules/MGauge/Load.cc \
-  Modules/MGauge/Random.cc \
   Modules/MGauge/StochEm.cc \
   Modules/MGauge/Unit.cc \
-  Modules/MScalar/ChargedProp.cc \
+  Modules/MGauge/Load.cc \
-  Modules/MScalar/FreeProp.cc
+  Modules/MGauge/FundtoHirep.cc \
+  Modules/MGauge/Random.cc
 
 modules_hpp =\
-  Modules/MAction/DWF.hpp \
-  Modules/MAction/Wilson.hpp \
-  Modules/MContraction/Baryon.hpp \
-  Modules/MContraction/DiscLoop.hpp \
-  Modules/MContraction/Gamma3pt.hpp \
-  Modules/MContraction/Meson.hpp \
-  Modules/MContraction/WeakHamiltonian.hpp \
-  Modules/MContraction/WeakHamiltonianEye.hpp \
-  Modules/MContraction/WeakHamiltonianNonEye.hpp \
-  Modules/MContraction/WeakNeutral4ptDisc.hpp \
-  Modules/MFermion/GaugeProp.hpp \
-  Modules/MGauge/Load.hpp \
-  Modules/MGauge/Random.hpp \
-  Modules/MGauge/StochEm.hpp \
-  Modules/MGauge/Unit.hpp \
   Modules/MLoop/NoiseLoop.hpp \
   Modules/MScalar/ChargedProp.hpp \
-  Modules/MScalar/FreeProp.hpp \
   Modules/MScalar/Scalar.hpp \
+  Modules/MScalar/FreeProp.hpp \
+  Modules/MSource/Wall.hpp \
+  Modules/MSource/SeqGamma.hpp \
+  Modules/MSource/Point.hpp \
+  Modules/MSource/Z2.hpp \
+  Modules/MFermion/GaugeProp.hpp \
+  Modules/MContraction/Meson.hpp \
+  Modules/MContraction/WeakHamiltonianNonEye.hpp \
+  Modules/MContraction/WeakHamiltonianEye.hpp \
+  Modules/MContraction/DiscLoop.hpp \
+  Modules/MContraction/Baryon.hpp \
+  Modules/MContraction/Gamma3pt.hpp \
+  Modules/MContraction/WeakNeutral4ptDisc.hpp \
+  Modules/MContraction/WeakHamiltonian.hpp \
   Modules/MSink/Point.hpp \
   Modules/MSolver/RBPrecCG.hpp \
-  Modules/MSource/Point.hpp \
+  Modules/MGauge/StochEm.hpp \
-  Modules/MSource/SeqGamma.hpp \
+  Modules/MGauge/FundtoHirep.hpp \
-  Modules/MSource/Wall.hpp \
+  Modules/MGauge/Unit.hpp \
-  Modules/MSource/Z2.hpp
+  Modules/MGauge/Load.hpp \
+  Modules/MGauge/Random.hpp \
+  Modules/MAction/WilsonClover.hpp \
+  Modules/MAction/DWF.hpp \
+  Modules/MAction/Wilson.hpp
 

lib/qcd/QCD.h

@@ -421,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)
+      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);
             }
@@ -437,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);
             }

lib/qcd/action/fermion/Fermion.h

@@ -106,6 +106,10 @@ 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;
@@ -116,6 +120,19 @@ 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

@@ -1004,6 +1004,10 @@ typedef WilsonImpl<vComplex,  TwoIndexSymmetricRepresentation, CoeffReal > Wilso
 typedef WilsonImpl<vComplexF, TwoIndexSymmetricRepresentation, CoeffReal > WilsonTwoIndexSymmetricImplF;  // Float
 typedef WilsonImpl<vComplexD, TwoIndexSymmetricRepresentation, CoeffReal > WilsonTwoIndexSymmetricImplD;  // Double
  
+typedef WilsonImpl<vComplex,  TwoIndexAntiSymmetricRepresentation, CoeffReal > WilsonTwoIndexAntiSymmetricImplR;   // Real.. whichever prec
+typedef WilsonImpl<vComplexF, TwoIndexAntiSymmetricRepresentation, CoeffReal > WilsonTwoIndexAntiSymmetricImplF;  // Float
+typedef WilsonImpl<vComplexD, TwoIndexAntiSymmetricRepresentation, CoeffReal > WilsonTwoIndexAntiSymmetricImplD;  // Double
+
 typedef DomainWallVec5dImpl<vComplex ,Nc, CoeffReal> DomainWallVec5dImplR; // Real.. whichever prec
 typedef DomainWallVec5dImpl<vComplexF,Nc, CoeffReal> DomainWallVec5dImplF; // Float
 typedef DomainWallVec5dImpl<vComplexD,Nc, CoeffReal> DomainWallVec5dImplD; // Double

lib/qcd/action/fermion/WilsonCloverFermion.cc

@@ -235,9 +235,9 @@ void WilsonCloverFermion<Impl>::MeeDeriv(GaugeField &mat, const FermionField &U,
   assert(0); // not implemented yet
 }
 
-FermOpTemplateInstantiate(WilsonCloverFermion); // now only for the fundamental representation
+FermOpTemplateInstantiate(WilsonCloverFermion);
-//AdjointFermOpTemplateInstantiate(WilsonCloverFermion);
+AdjointFermOpTemplateInstantiate(WilsonCloverFermion);
-//TwoIndexFermOpTemplateInstantiate(WilsonCloverFermion);
+TwoIndexFermOpTemplateInstantiate(WilsonCloverFermion);
 //GparityFermOpTemplateInstantiate(WilsonCloverFermion);
 }
 }

lib/qcd/action/fermion/WilsonFermion.h

@@ -44,7 +44,8 @@ class WilsonFermionStatic {
   static const int npoint = 8;
 };
 
-struct WilsonAnisotropyCoefficients{
+ struct WilsonAnisotropyCoefficients: Serializable
+ {
   GRID_SERIALIZABLE_CLASS_MEMBERS(WilsonAnisotropyCoefficients,
   bool, isAnisotropic,
   int, t_direction,

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/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
   //////////////////////////////////////////////////

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_wilsonFund.cc

@@ -0,0 +1,160 @@
+/*******************************************************************************
+ 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 = 1;
+    globalPar.trajCounter.end   = 2;
+    globalPar.trajCounter.step  = 1;
+
+    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<MGauge::Load>("gauge");
+   //application.createModule<MGauge::Unit>("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_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_WilsonCloverFermionGauge.cc

@@ -67,6 +67,9 @@ int main(int argc, char **argv)
   // Construct observables
   typedef PlaquetteMod<HMCWrapper::ImplPolicy> PlaqObs;
   TheHMC.Resources.AddObservable<PlaqObs>();
+
+  typedef PolyakovMod<HMCWrapper::ImplPolicy> PolyakovObs;
+  TheHMC.Resources.AddObservable<PolyakovObs>();
   //////////////////////////////////////////////
 
   /////////////////////////////////////////////////////////////

tests/lanczos/Test_WCMultiRep_lanczos.cc

@@ -0,0 +1,170 @@
+/*************************************************************************************
+
+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);
+  ImplicitlyRestartedLanczos<FundFermionField> IRL_Fund(HermOpFund, FundPolyX, Nstop, Nk, Nm,
+                                               resid, MaxIt);
+  
+  Polynomial<ASymmFermionField> ASPolyX(Coeffs);
+  Chebyshev<ASymmFermionField> ASCheb(0.0, 10., 12);
+  ImplicitlyRestartedLanczos<ASymmFermionField> IRL_AS(HermOpAS, ASPolyX, 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();
+}