Minor changes

extras/Hadrons/Modules.hpp

@@ -1,25 +1,26 @@
-#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/MLoop/NoiseLoop.hpp>
+#include <Grid/Hadrons/Modules/MFermion/GaugeProp.hpp>
 #include <Grid/Hadrons/Modules/MContraction/WeakHamiltonian.hpp>
+#include <Grid/Hadrons/Modules/MContraction/Meson.hpp>
+#include <Grid/Hadrons/Modules/MContraction/DiscLoop.hpp>
 #include <Grid/Hadrons/Modules/MContraction/WeakHamiltonianEye.hpp>
+#include <Grid/Hadrons/Modules/MContraction/Baryon.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/MContraction/Gamma3pt.hpp>
+#include <Grid/Hadrons/Modules/MSource/Z2.hpp>
+#include <Grid/Hadrons/Modules/MSource/SeqGamma.hpp>
+#include <Grid/Hadrons/Modules/MSource/Point.hpp>
+#include <Grid/Hadrons/Modules/MSource/Wall.hpp>
+#include <Grid/Hadrons/Modules/MSource/Laplacian.hpp>
+#include <Grid/Hadrons/Modules/MSolver/RBPrecCG.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/MAction/DWF.hpp>
+#include <Grid/Hadrons/Modules/MAction/Wilson.hpp>
+#include <Grid/Hadrons/Modules/MGauge/StochEm.hpp>
+#include <Grid/Hadrons/Modules/MGauge/Unit.hpp>
+#include <Grid/Hadrons/Modules/MGauge/Random.hpp>
+#include <Grid/Hadrons/Modules/MGauge/Load.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/MSource/SeqGamma.hpp>
-#include <Grid/Hadrons/Modules/MSource/Wall.hpp>
-#include <Grid/Hadrons/Modules/MSource/Z2.hpp>

extras/Hadrons/Modules/MSource/Laplacian.hpp

@@ -0,0 +1,153 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: extras/Hadrons/Modules/MSource/Laplacian.hpp
+
+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 */
+
+#ifndef Hadrons_MSource_Laplacian_hpp_
+#define Hadrons_MSource_Laplacian_hpp_
+
+#include <Grid/Hadrons/Global.hpp>
+#include <Grid/Hadrons/Module.hpp>
+#include <Grid/Hadrons/ModuleFactory.hpp>
+
+BEGIN_HADRONS_NAMESPACE
+
+/*
+ 
+ Laplacian smearing source
+ -----------------------------
+ 
+ * options:
+ - source: name of source object to be smeared (string)
+ - N: number of steps (integer)
+ - alpha: smearing parameter (real)
+ 
+ */
+
+/******************************************************************************
+ *                          Laplace smearing operator                         *
+ ******************************************************************************/
+BEGIN_MODULE_NAMESPACE(MSource)
+
+class LaplacianPar : Serializable
+{
+  public:
+    GRID_SERIALIZABLE_CLASS_MEMBERS(LaplacianPar,
+                                    std::string, source,
+                                    std::string, gauge,
+                                    unsigned int, N,
+                                    double, alpha);
+};
+
+template <typename FImpl>
+class TLaplacian : public Module<LaplacianPar>
+{
+  public:
+    FERM_TYPE_ALIASES(FImpl, );
+
+  public:
+    // constructor
+    TLaplacian(const std::string name);
+    // destructor
+    virtual ~TLaplacian(void) = default;
+    // dependency relation
+    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(LaplaceSmearing, TLaplacian<FIMPL>, MSource);
+
+/******************************************************************************
+ *                       TLaplacian template implementation                   *
+ ******************************************************************************/
+// constructor /////////////////////////////////////////////////////////////////
+template <typename FImpl>
+TLaplacian<FImpl>::TLaplacian(const std::string name)
+    : Module<LaplacianPar>(name)
+{
+}
+
+// dependencies/products ///////////////////////////////////////////////////////
+template <typename FImpl>
+std::vector<std::string> TLaplacian<FImpl>::getInput(void)
+{
+    std::vector<std::string> in = {par().source, par().gauge};
+
+    return in;
+}
+
+template <typename FImpl>
+std::vector<std::string> TLaplacian<FImpl>::getOutput(void)
+{
+    std::vector<std::string> out = {getName()};
+
+    return out;
+}
+
+// setup ///////////////////////////////////////////////////////////////////////
+template <typename FImpl>
+void TLaplacian<FImpl>::setup(void)
+{
+    env().template registerLattice<PropagatorField>(getName());
+}
+
+// execution ///////////////////////////////////////////////////////////////////
+template <typename FImpl>
+void TLaplacian<FImpl>::execute(void)
+{
+
+    FermionField source(env().getGrid()), tmp(env().getGrid());
+    PropagatorField &SmrSrc = *env().template createLattice<PropagatorField>(getName());
+    PropagatorField &fullSrc = *env().template getObject<PropagatorField>(par().source);
+    auto &U      = *env().template getObject<LatticeGaugeField>(par().gauge);
+    Laplacian<FImpl> LaplaceOperator(env().getGrid());
+    LaplaceOperator.ImportGauge(U);
+    double prefactor = par().alpha / (double)(par().N);
+
+    for (unsigned int s = 0; s < Ns; ++s)
+    {
+        for (unsigned int c = 0; c < Nc; ++c)
+        {
+            PropToFerm(source, fullSrc, s, c);
+            for (int smr = 0; smr < par().N; ++smr)
+            {
+                LaplaceOperator.M(source, tmp);
+                source += prefactor * tmp;
+            }
+            FermToProp(SmrSrc, source, s, c);
+        }
+    }
+}
+
+END_MODULE_NAMESPACE
+
+END_HADRONS_NAMESPACE
+
+#endif // Hadrons_MSource_Z2_hpp_

extras/Hadrons/modules.inc

@@ -1,38 +1,39 @@
 modules_cc =\
-  Modules/MContraction/WeakHamiltonianEye.cc \
   Modules/MContraction/WeakHamiltonianNonEye.cc \
   Modules/MContraction/WeakNeutral4ptDisc.cc \
-  Modules/MGauge/Load.cc \
+  Modules/MContraction/WeakHamiltonianEye.cc \
+  Modules/MScalar/FreeProp.cc \
+  Modules/MScalar/ChargedProp.cc \
+  Modules/MGauge/Unit.cc \
   Modules/MGauge/Random.cc \
   Modules/MGauge/StochEm.cc \
-  Modules/MGauge/Unit.cc \
-  Modules/MScalar/ChargedProp.cc \
-  Modules/MScalar/FreeProp.cc
+  Modules/MGauge/Load.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/MLoop/NoiseLoop.hpp \
+  Modules/MFermion/GaugeProp.hpp \
   Modules/MContraction/WeakHamiltonian.hpp \
+  Modules/MContraction/Meson.hpp \
+  Modules/MContraction/DiscLoop.hpp \
   Modules/MContraction/WeakHamiltonianEye.hpp \
+  Modules/MContraction/Baryon.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/MContraction/Gamma3pt.hpp \
+  Modules/MSource/Z2.hpp \
+  Modules/MSource/SeqGamma.hpp \
+  Modules/MSource/Point.hpp \
+  Modules/MSource/Wall.hpp \
+  Modules/MSource/Laplacian.hpp \
+  Modules/MSolver/RBPrecCG.hpp \
   Modules/MScalar/ChargedProp.hpp \
   Modules/MScalar/FreeProp.hpp \
   Modules/MScalar/Scalar.hpp \
-  Modules/MSink/Point.hpp \
-  Modules/MSolver/RBPrecCG.hpp \
-  Modules/MSource/Point.hpp \
-  Modules/MSource/SeqGamma.hpp \
-  Modules/MSource/Wall.hpp \
-  Modules/MSource/Z2.hpp
+  Modules/MAction/DWF.hpp \
+  Modules/MAction/Wilson.hpp \
+  Modules/MGauge/StochEm.hpp \
+  Modules/MGauge/Unit.hpp \
+  Modules/MGauge/Random.hpp \
+  Modules/MGauge/Load.hpp \
+  Modules/MSink/Point.hpp
 

lib/qcd/action/Action.h

@@ -47,4 +47,8 @@ Author: paboyle <paboyle@ph.ed.ac.uk>
 ////////////////////////////////////////
 #include <Grid/qcd/action/pseudofermion/PseudoFermion.h>
 
+////////////////////////////////////////////////////////////////////////
+// Laplacian on fermion fields
+////////////////////////////////////////////////////////////////////////
+#include <Grid/qcd/utils/CovariantLaplacian.h>
 #endif

lib/qcd/action/ActionCore.h

@@ -53,7 +53,7 @@ directory
 // Utility functions
 ////////////////////////////////////////////
 #include <Grid/qcd/utils/Metric.h>
-#include <Grid/qcd/utils/CovariantLaplacian.h>
+#include <Grid/qcd/utils/CovariantAdjointLaplacian.h>
 
 
 

lib/qcd/action/fermion/DomainWallEOFAFermion.cc

@@ -60,11 +60,11 @@ namespace QCD {
         Approx::zolotarev_free(zdata);
     }
 
-    /***************************************************************
-    /* Additional EOFA operators only called outside the inverter.
-    /* Since speed is not essential, simple axpby-style
-    /* implementations should be fine.
-    /***************************************************************/
+    /*
+     Additional EOFA operators only called outside the inverter.
+     Since speed is not essential, simple axpby-style
+     implementations should be fine.
+    */
     template<class Impl>
     void DomainWallEOFAFermion<Impl>::Omega(const FermionField& psi, FermionField& Din, int sign, int dag)
     {
@@ -115,9 +115,9 @@ namespace QCD {
         return(norm2(chi));
     }
 
-    /********************************************************************
-    /* Performance critical fermion operators called inside the inverter
-    /********************************************************************/
+    
+    // Performance critical fermion operators called inside the inverter
+    
 
     template<class Impl>
     void DomainWallEOFAFermion<Impl>::M5D(const FermionField& psi, FermionField& chi)

lib/qcd/action/fermion/MobiusEOFAFermion.cc

@@ -77,11 +77,11 @@ namespace QCD {
       }
     }
 
-    /***************************************************************
-    /* Additional EOFA operators only called outside the inverter.
-    /* Since speed is not essential, simple axpby-style
-    /* implementations should be fine.
-    /***************************************************************/
+    /*
+     Additional EOFA operators only called outside the inverter.
+     Since speed is not essential, simple axpby-style
+     implementations should be fine.
+    */
     template<class Impl>
     void MobiusEOFAFermion<Impl>::Omega(const FermionField& psi, FermionField& Din, int sign, int dag)
     {
@@ -193,9 +193,9 @@ namespace QCD {
       return(norm2(chi));
     }
 
-    /********************************************************************
-    /* Performance critical fermion operators called inside the inverter
-    /********************************************************************/
+    
+    // Performance critical fermion operators called inside the inverter
+    
 
     template<class Impl>
     void MobiusEOFAFermion<Impl>::M5D(const FermionField& psi, FermionField& chi)

lib/qcd/action/pseudofermion/EvenOddSchurDifferentiable.h

@@ -38,7 +38,7 @@ namespace Grid{
     //         (Moe Moo)    (Moe Mee^-1    1 )   (0   Moo-Moe Mee^-1 Meo)  (0   1         )
     //
     // Determinant is det of middle factor
-    // This assumes Mee is indept of U.
+    // NOTICE: This assumes Mee is indept of U in computing the derivative
     //
     template<class Impl>
     class SchurDifferentiableOperator :  public SchurDiagMooeeOperator<FermionOperator<Impl>,typename Impl::FermionField> 
@@ -77,7 +77,7 @@ namespace Grid{
           //  X^dag Der_oe MeeInv Meo Y
           // Use Mooee as nontrivial but gauge field indept
           this->_Mat.Meooe   (V,tmp1);      // odd->even -- implicit -0.5 factor to be applied
-	  this->_Mat.MooeeInv(tmp1,tmp2);   // even->even 
+	        this->_Mat.MooeeInv(tmp1,tmp2);   // even->even 
           this->_Mat.MoeDeriv(ForceO,U,tmp2,DaggerNo);
           //  Accumulate X^dag M_oe MeeInv Der_eo Y
           this->_Mat.MeooeDag   (U,tmp1);    // even->odd -- implicit -0.5 factor to be applied

lib/qcd/utils/CovariantAdjointLaplacian.h

@@ -0,0 +1,209 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid
+
+Source file: ./lib/qcd/action/scalar/CovariantAdjointLaplacian.h
+
+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 */
+
+#ifndef COVARIANT_ADJOINT_LAPLACIAN_H
+#define COVARIANT_ADJOINT_LAPLACIAN_H
+
+namespace Grid
+{
+namespace QCD
+{
+
+struct LaplacianParams : Serializable
+{
+  GRID_SERIALIZABLE_CLASS_MEMBERS(LaplacianParams,
+                                  RealD, lo,
+                                  RealD, hi,
+                                  int, MaxIter,
+                                  RealD, tolerance,
+                                  int, degree,
+                                  int, precision);
+
+  // constructor
+  LaplacianParams(RealD lo = 0.0,
+                  RealD hi = 1.0,
+                  int maxit = 1000,
+                  RealD tol = 1.0e-8,
+                  int degree = 10,
+                  int precision = 64)
+      : lo(lo),
+        hi(hi),
+        MaxIter(maxit),
+        tolerance(tol),
+        degree(degree),
+        precision(precision){};
+};
+
+////////////////////////////////////////////////////////////
+// Laplacian operator L on adjoint fields
+//
+// phi: adjoint field
+// L: D_mu^dag D_mu
+//
+// L phi(x) = Sum_mu [ U_mu(x)phi(x+mu)U_mu(x)^dag +
+//                     U_mu(x-mu)^dag phi(x-mu)U_mu(x-mu)
+//                     -2phi(x)]
+//
+// Operator designed to be encapsulated by
+// an HermitianLinearOperator<.. , ..>
+////////////////////////////////////////////////////////////
+
+template <class Impl>
+class LaplacianAdjointField : public Metric<typename Impl::Field>
+{
+  OperatorFunction<typename Impl::Field> &Solver;
+  LaplacianParams param;
+  MultiShiftFunction PowerHalf;
+  MultiShiftFunction PowerInvHalf;
+
+public:
+  INHERIT_GIMPL_TYPES(Impl);
+
+  LaplacianAdjointField(GridBase *grid, OperatorFunction<GaugeField> &S, LaplacianParams &p, const RealD k = 1.0)
+      : U(Nd, grid), Solver(S), param(p), kappa(k)
+  {
+    AlgRemez remez(param.lo, param.hi, param.precision);
+    std::cout << GridLogMessage << "Generating degree " << param.degree << " for x^(1/2)" << std::endl;
+    remez.generateApprox(param.degree, 1, 2);
+    PowerHalf.Init(remez, param.tolerance, false);
+    PowerInvHalf.Init(remez, param.tolerance, true);
+  };
+
+  void Mdir(const GaugeField &, GaugeField &, int, int) { assert(0); }
+  void Mdiag(const GaugeField &, GaugeField &) { assert(0); }
+
+  void ImportGauge(const GaugeField &_U)
+  {
+    for (int mu = 0; mu < Nd; mu++)
+    {
+      U[mu] = PeekIndex<LorentzIndex>(_U, mu);
+    }
+  }
+
+  void M(const GaugeField &in, GaugeField &out)
+  {
+    // in is an antihermitian matrix
+    // test
+    //GaugeField herm = in + adj(in);
+    //std::cout << "AHermiticity: " << norm2(herm) << std::endl;
+
+    GaugeLinkField tmp(in._grid);
+    GaugeLinkField tmp2(in._grid);
+    GaugeLinkField sum(in._grid);
+
+    for (int nu = 0; nu < Nd; nu++)
+    {
+      sum = zero;
+      GaugeLinkField in_nu = PeekIndex<LorentzIndex>(in, nu);
+      GaugeLinkField out_nu(out._grid);
+      for (int mu = 0; mu < Nd; mu++)
+      {
+        tmp = U[mu] * Cshift(in_nu, mu, +1) * adj(U[mu]);
+        tmp2 = adj(U[mu]) * in_nu * U[mu];
+        sum += tmp + Cshift(tmp2, mu, -1) - 2.0 * in_nu;
+      }
+      out_nu = (1.0 - kappa) * in_nu - kappa / (double(4 * Nd)) * sum;
+      PokeIndex<LorentzIndex>(out, out_nu, nu);
+    }
+  }
+
+  void MDeriv(const GaugeField &in, GaugeField &der)
+  {
+    // in is anti-hermitian
+    RealD factor = -kappa / (double(4 * Nd));
+
+    for (int mu = 0; mu < Nd; mu++)
+    {
+      GaugeLinkField der_mu(der._grid);
+      der_mu = zero;
+      for (int nu = 0; nu < Nd; nu++)
+      {
+        GaugeLinkField in_nu = PeekIndex<LorentzIndex>(in, nu);
+        der_mu += U[mu] * Cshift(in_nu, mu, 1) * adj(U[mu]) * in_nu;
+      }
+      // the minus sign comes by using the in_nu instead of the
+      // adjoint in the last multiplication
+      PokeIndex<LorentzIndex>(der, -2.0 * factor * der_mu, mu);
+    }
+  }
+
+  // separating this temporarily
+  void MDeriv(const GaugeField &left, const GaugeField &right,
+              GaugeField &der)
+  {
+    // in is anti-hermitian
+    RealD factor = -kappa / (double(4 * Nd));
+
+    for (int mu = 0; mu < Nd; mu++)
+    {
+      GaugeLinkField der_mu(der._grid);
+      der_mu = zero;
+      for (int nu = 0; nu < Nd; nu++)
+      {
+        GaugeLinkField left_nu = PeekIndex<LorentzIndex>(left, nu);
+        GaugeLinkField right_nu = PeekIndex<LorentzIndex>(right, nu);
+        der_mu += U[mu] * Cshift(left_nu, mu, 1) * adj(U[mu]) * right_nu;
+        der_mu += U[mu] * Cshift(right_nu, mu, 1) * adj(U[mu]) * left_nu;
+      }
+      PokeIndex<LorentzIndex>(der, -factor * der_mu, mu);
+    }
+  }
+
+  void Minv(const GaugeField &in, GaugeField &inverted)
+  {
+    HermitianLinearOperator<LaplacianAdjointField<Impl>, GaugeField> HermOp(*this);
+    Solver(HermOp, in, inverted);
+  }
+
+  void MSquareRoot(GaugeField &P)
+  {
+    GaugeField Gp(P._grid);
+    HermitianLinearOperator<LaplacianAdjointField<Impl>, GaugeField> HermOp(*this);
+    ConjugateGradientMultiShift<GaugeField> msCG(param.MaxIter, PowerHalf);
+    msCG(HermOp, P, Gp);
+    P = Gp;
+  }
+
+  void MInvSquareRoot(GaugeField &P)
+  {
+    GaugeField Gp(P._grid);
+    HermitianLinearOperator<LaplacianAdjointField<Impl>, GaugeField> HermOp(*this);
+    ConjugateGradientMultiShift<GaugeField> msCG(param.MaxIter, PowerInvHalf);
+    msCG(HermOp, P, Gp);
+    P = Gp;
+  }
+
+private:
+  RealD kappa;
+  std::vector<GaugeLinkField> U;
+};
+
+} // QCD
+} // Grid
+#endif

lib/qcd/utils/CovariantLaplacian.h

@@ -4,7 +4,7 @@ Grid physics library, www.github.com/paboyle/Grid
 
 Source file: ./lib/qcd/action/scalar/CovariantLaplacian.h
 
-Copyright (C) 2016
+Copyright (C) 2017
 
 Author: Guido Cossu <guido.cossu@ed.ac.uk>
 
@@ -30,168 +30,57 @@ directory
 #ifndef COVARIANT_LAPLACIAN_H
 #define COVARIANT_LAPLACIAN_H
 
-namespace Grid {
-namespace QCD {
-
-struct LaplacianParams : Serializable {
-  GRID_SERIALIZABLE_CLASS_MEMBERS(LaplacianParams, 
-                                  RealD, lo, 
-                                  RealD, hi, 
-                                  int,   MaxIter, 
-                                  RealD, tolerance, 
-                                  int,   degree, 
-                                  int,   precision);
-  
-  // constructor 
-  LaplacianParams(RealD lo      = 0.0, 
-                  RealD hi      = 1.0, 
-                  int maxit     = 1000,
-                  RealD tol     = 1.0e-8, 
-                  int degree    = 10,
-                  int precision = 64)
-    : lo(lo),
-      hi(hi),
-      MaxIter(maxit),
-      tolerance(tol),
-      degree(degree),
-      precision(precision){};
-};
-
-
+namespace Grid
+{
+namespace QCD
+{
 
 ////////////////////////////////////////////////////////////
-// Laplacian operator L on adjoint fields
+// Laplacian operator L on fermion fields
 //
-// phi: adjoint field
-// L: D_mu^dag D_mu
+// phi: fermion field
 //
-// L phi(x) = Sum_mu [ U_mu(x)phi(x+mu)U_mu(x)^dag + 
-//                     U_mu(x-mu)^dag phi(x-mu)U_mu(x-mu)
-//                     -2phi(x)]
+// L phi(x) = Sum_mu [ U_mu(x) phi(x+mu) + U_mu(x-mu) phi(x-mu) - 2phi(x)]
 //
 // Operator designed to be encapsulated by
 // an HermitianLinearOperator<.. , ..>
 ////////////////////////////////////////////////////////////
 
+// has to inherit from a fermion implementation
 template <class Impl>
-class LaplacianAdjointField: public Metric<typename Impl::Field> {
-  OperatorFunction<typename Impl::Field> &Solver;
-  LaplacianParams param;
-  MultiShiftFunction PowerHalf;    
-  MultiShiftFunction PowerInvHalf;    
+class Laplacian
+{
+public:
+  INHERIT_IMPL_TYPES(Impl);
 
- public:
-  INHERIT_GIMPL_TYPES(Impl);
+  // add a bool to smear only in the spatial directions
+  Laplacian(GridBase *grid, bool spatial = false)
+      : U(Nd, grid), spatial_laplacian(spatial){};
 
-  LaplacianAdjointField(GridBase* grid, OperatorFunction<GaugeField>& S, LaplacianParams& p, const RealD k = 1.0)
-      : U(Nd, grid), Solver(S), param(p), kappa(k){
-        AlgRemez remez(param.lo,param.hi,param.precision);
-        std::cout<<GridLogMessage << "Generating degree "<<param.degree<<" for x^(1/2)"<<std::endl;
-        remez.generateApprox(param.degree,1,2);
-        PowerHalf.Init(remez,param.tolerance,false);
-        PowerInvHalf.Init(remez,param.tolerance,true);
-        
+  void Mdir(const FermionField &, FermionField &, int, int) { assert(0); }
+  void Mdiag(const FermionField &, FermionField &) { assert(0); }
 
-      };
-
-  void Mdir(const GaugeField&, GaugeField&, int, int){ assert(0);}
-  void Mdiag(const GaugeField&, GaugeField&){ assert(0);}
-
-  void ImportGauge(const GaugeField& _U) {
-    for (int mu = 0; mu < Nd; mu++) {
+  void ImportGauge(const GaugeField &_U)
+  {
+    for (int mu = 0; mu < Nd; mu++)
       U[mu] = PeekIndex<LorentzIndex>(_U, mu);
-    }
   }
 
-  void M(const GaugeField& in, GaugeField& out) {
-    // in is an antihermitian matrix
-    // test
-    //GaugeField herm = in + adj(in);
-    //std::cout << "AHermiticity: " << norm2(herm) << std::endl;
+  void M(const FermionField &in, FermionField &out)
+  {
+    int dims = spatial_laplacian ? (Nd - 1) : Nd;
 
-    GaugeLinkField tmp(in._grid);
-    GaugeLinkField tmp2(in._grid);
-    GaugeLinkField sum(in._grid);
-
-    for (int nu = 0; nu < Nd; nu++) {
-      sum = zero;
-      GaugeLinkField in_nu = PeekIndex<LorentzIndex>(in, nu);
-      GaugeLinkField out_nu(out._grid);
-      for (int mu = 0; mu < Nd; mu++) {
-        tmp = U[mu] * Cshift(in_nu, mu, +1) * adj(U[mu]);
-        tmp2 = adj(U[mu]) * in_nu * U[mu];
-        sum += tmp + Cshift(tmp2, mu, -1) - 2.0 * in_nu;
-      }
-      out_nu = (1.0 - kappa) * in_nu - kappa / (double(4 * Nd)) * sum;
-      PokeIndex<LorentzIndex>(out, out_nu, nu);
-    }
+    out = -2.0 * dims * in;
+    // eventually speed up with the stencil operator, if necessary
+    for (int mu = 0; mu < dims; mu++)
+      out += Impl::CovShiftForward(U[mu], mu, in) + Impl::CovShiftBackward(U[mu], mu, in);
   }
 
-  void MDeriv(const GaugeField& in, GaugeField& der) {
-    // in is anti-hermitian
-    RealD factor = -kappa / (double(4 * Nd));
-    
-    for (int mu = 0; mu < Nd; mu++){
-      GaugeLinkField der_mu(der._grid);
-      der_mu = zero;
-      for (int nu = 0; nu < Nd; nu++){
-        GaugeLinkField in_nu = PeekIndex<LorentzIndex>(in, nu);
-        der_mu += U[mu] * Cshift(in_nu, mu, 1) * adj(U[mu]) * in_nu;
-      }
-      // the minus sign comes by using the in_nu instead of the
-      // adjoint in the last multiplication
-      PokeIndex<LorentzIndex>(der,  -2.0 * factor * der_mu, mu);
-    } 
-  }
-
-  // separating this temporarily
-  void MDeriv(const GaugeField& left, const GaugeField& right,
-              GaugeField& der) {
-    // in is anti-hermitian
-    RealD factor = -kappa / (double(4 * Nd));
-
-    for (int mu = 0; mu < Nd; mu++) {
-      GaugeLinkField der_mu(der._grid);
-      der_mu = zero;
-      for (int nu = 0; nu < Nd; nu++) {
-        GaugeLinkField left_nu = PeekIndex<LorentzIndex>(left, nu);
-        GaugeLinkField right_nu = PeekIndex<LorentzIndex>(right, nu);
-        der_mu += U[mu] * Cshift(left_nu, mu, 1) * adj(U[mu]) * right_nu;
-        der_mu += U[mu] * Cshift(right_nu, mu, 1) * adj(U[mu]) * left_nu;
-      }
-      PokeIndex<LorentzIndex>(der, -factor * der_mu, mu);
-    }
-  }
-
-  void Minv(const GaugeField& in, GaugeField& inverted){
-    HermitianLinearOperator<LaplacianAdjointField<Impl>,GaugeField> HermOp(*this);
-    Solver(HermOp, in, inverted);
-  }
-
-  void MSquareRoot(GaugeField& P){
-    GaugeField Gp(P._grid);
-    HermitianLinearOperator<LaplacianAdjointField<Impl>,GaugeField> HermOp(*this);
-    ConjugateGradientMultiShift<GaugeField> msCG(param.MaxIter,PowerHalf);
-    msCG(HermOp,P,Gp);
-    P = Gp; 
-  }
-
-  void MInvSquareRoot(GaugeField& P){
-    GaugeField Gp(P._grid);
-    HermitianLinearOperator<LaplacianAdjointField<Impl>,GaugeField> HermOp(*this);
-    ConjugateGradientMultiShift<GaugeField> msCG(param.MaxIter,PowerInvHalf);
-    msCG(HermOp,P,Gp);
-    P = Gp; 
-  }
-
-
-
- private:
-  RealD kappa;
+private:
+  bool spatial_laplacian;
   std::vector<GaugeLinkField> U;
-};
-
-}
-}
+}; // Laplacian
 
+} // QCD
+} // Grid
 #endif

tests/core/Test_laplacian.cc

@@ -0,0 +1,105 @@
+    /*************************************************************************************
+
+    Grid physics library, www.github.com/paboyle/Grid
+
+    Source file: ./tests/Test_laplacian.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>
+
+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;
+
+  GridParallelRNG          pRNG(&Grid);
+  pRNG.SeedFixedIntegers(std::vector<int>({45,12,81,9}));
+
+
+  std::vector<int> point({0,0,0,0});
+
+  LatticeFermion src   (&Grid); //random(pRNG,src);
+  SpinColourVectorD Sp;
+  for (unsigned int s = 0; s < Ns; ++s)
+      for (unsigned int c = 0; c < Nc; ++c)
+        Sp()(s)(c) = 1;
+
+  src = zero;
+  pokeSite(Sp,src,point);
+
+  LatticeFermion result(&Grid); result=zero;
+  LatticeFermion tmp(&Grid); tmp=zero;
+
+  // Gauge configuration
+  LatticeGaugeField Umu(&Grid); SU3::HotConfiguration(pRNG,Umu);
+
+  std::cout<<GridLogMessage<<"=============================================================="<<std::endl;
+  std::cout<<GridLogMessage<<"= Testing the laplacian operator on a point source            "<<std::endl;
+  std::cout<<GridLogMessage<<"=============================================================="<<std::endl;
+
+  Laplacian<WilsonImplR> LaplaceOperator(src._grid);
+  LaplaceOperator.ImportGauge(Umu);
+  LaplaceOperator.M(src, result);
+
+  std::cout << "Source vector" << std::endl;
+  std::cout << src << std::endl;
+
+  std::cout << "Result vector" << std::endl;
+  std::cout << result << std::endl;
+
+  std::cout<<GridLogMessage<<"=============================================================="<<std::endl;
+  std::cout<<GridLogMessage<<"= Testing the laplacian smearing operator on a point source   "<<std::endl;
+  std::cout<<GridLogMessage<<"=============================================================="<<std::endl;
+
+  LatticeFermion smeared  (&Grid); smeared = src;
+  for (int smr = 0; smr < 10; ++smr)
+  {
+      LaplaceOperator.M(smeared, tmp);
+      smeared += 0.1*tmp;
+  }
+
+  std::cout << "Smeared vector" << std::endl;
+  std::cout << smeared << std::endl;
+
+  // Norm of vector
+  LatticeComplex smr_norm(&Grid);
+  smr_norm = localNorm2(smeared);
+  std::cout << "Smeared vector norm" << std::endl;
+  std::cout << smr_norm << std::endl;
+
+
+
+  Grid_finalize();
+}

tests/hadrons/Test_hadrons_meson_3pt_laplacian.cc

@@ -0,0 +1,195 @@
+/*******************************************************************************
+ Grid physics library, www.github.com/paboyle/Grid
+ 
+ Source file: tests/hadrons/Test_hadrons_meson_3pt.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", "s", "c1", "c2", "c3"};
+    std::vector<double>      mass    = {.01, .04, .2  , .25 , .3  };
+    unsigned int             nt      = GridDefaultLatt()[Tp];
+    
+    // global parameters
+    Application::GlobalPar globalPar;
+    globalPar.trajCounter.start    = 1500;
+    globalPar.trajCounter.end      = 1520;
+    globalPar.trajCounter.step     = 20;
+    globalPar.seed                 = "1 2 3 4";
+    globalPar.genetic.maxGen       = 1000;
+    globalPar.genetic.maxCstGen    = 200;
+    globalPar.genetic.popSize      = 20;
+    globalPar.genetic.mutationRate = .1;
+    application.setPar(globalPar);
+
+    
+    // gauge field
+    application.createModule<MGauge::Unit>("gauge");
+    
+    // set fermion boundary conditions to be periodic space, antiperiodic time.
+    std::string boundary = "1 1 1 -1";
+
+    // sink
+    MSink::Point::Par sinkPar;
+    sinkPar.mom = "0 0 0";
+    application.createModule<MSink::ScalarPoint>("sink", sinkPar);
+    for (unsigned int i = 0; i < flavour.size(); ++i)
+    {
+        // actions
+        MAction::DWF::Par actionPar;
+        actionPar.gauge = "gauge";
+        actionPar.Ls    = 12;
+        actionPar.M5    = 1.8;
+        actionPar.mass  = mass[i];
+        actionPar.boundary = boundary;
+        application.createModule<MAction::DWF>("DWF_" + flavour[i], actionPar);
+        
+        // solvers
+        MSolver::RBPrecCG::Par solverPar;
+        solverPar.action   = "DWF_" + flavour[i];
+        solverPar.residual = 1.0e-8;
+        application.createModule<MSolver::RBPrecCG>("CG_" + flavour[i],
+                                                    solverPar);
+    }
+    for (unsigned int t = 0; t < nt; t += 1)
+    {
+        std::string                           srcName;
+        std::string                           lapName;
+        std::vector<std::string>              qName;
+        std::vector<std::vector<std::string>> seqName;
+        
+        // Z2 source
+        MSource::Z2::Par z2Par;
+        z2Par.tA = t;
+        z2Par.tB = t;
+        srcName  = "z2_" + std::to_string(t);
+        application.createModule<MSource::Z2>(srcName, z2Par);
+
+        // Example of smearing of the source 
+        MSource::LaplaceSmearing::Par LapPar;
+        LapPar.N = 10;
+        LapPar.alpha = 0.1;
+        LapPar.source = srcName;
+        LapPar.gauge = "gauge";
+        lapName = "z2smr_" + std::to_string(t);
+        application.createModule<MSource::LaplaceSmearing>(lapName, LapPar);
+
+        for (unsigned int i = 0; i < flavour.size(); ++i)
+        {
+            // sequential sources
+            MSource::SeqGamma::Par seqPar;
+            qName.push_back("QZ2_" + flavour[i] + "_" + std::to_string(t));
+            seqPar.q   = qName[i];
+            seqPar.tA  = (t + nt/4) % nt;
+            seqPar.tB  = (t + nt/4) % nt;
+            seqPar.mom = "1. 0. 0. 0.";
+            seqName.push_back(std::vector<std::string>(Nd));
+            for (unsigned int mu = 0; mu < Nd; ++mu)
+            {
+                seqPar.gamma   = 0x1 << mu;
+                seqName[i][mu] = "G" + std::to_string(seqPar.gamma)
+                                 + "_" + std::to_string(seqPar.tA) + "-"
+                                 + qName[i];
+                application.createModule<MSource::SeqGamma>(seqName[i][mu], seqPar);
+            }
+            
+            // propagators
+            MFermion::GaugeProp::Par quarkPar;
+            quarkPar.solver = "CG_" + flavour[i];
+            quarkPar.source = srcName;
+            application.createModule<MFermion::GaugeProp>(qName[i], quarkPar);
+            for (unsigned int mu = 0; mu < Nd; ++mu)
+            {
+                quarkPar.source = seqName[i][mu];
+                seqName[i][mu]  = "Q_" + flavour[i] + "-" + seqName[i][mu];
+                application.createModule<MFermion::GaugeProp>(seqName[i][mu], quarkPar);
+            }
+        }
+        
+        // contractions
+        MContraction::Meson::Par mesPar;
+        for (unsigned int i = 0; i < flavour.size(); ++i)
+        for (unsigned int j = i; j < flavour.size(); ++j)
+        {
+            mesPar.output = "mesons/Z2_" + flavour[i] + flavour[j];
+            mesPar.q1     = qName[i];
+            mesPar.q2     = qName[j];
+            mesPar.gammas = "all";
+            mesPar.sink   = "sink";
+            application.createModule<MContraction::Meson>("meson_Z2_"
+                                                          + std::to_string(t)
+                                                          + "_"
+                                                          + flavour[i]
+                                                          + flavour[j],
+                                                          mesPar);
+        }
+        for (unsigned int i = 0; i < flavour.size(); ++i)
+        for (unsigned int j = 0; j < flavour.size(); ++j)
+        for (unsigned int mu = 0; mu < Nd; ++mu)
+        {
+            MContraction::Meson::Par mesPar;
+            
+            mesPar.output = "3pt/Z2_" + flavour[i] + flavour[j] + "_"
+                            + std::to_string(mu);
+            mesPar.q1     = qName[i];
+            mesPar.q2     = seqName[j][mu];
+            mesPar.gammas = "all";
+            mesPar.sink   = "sink";
+            application.createModule<MContraction::Meson>("3pt_Z2_"
+                                                          + std::to_string(t)
+                                                          + "_"
+                                                          + flavour[i]
+                                                          + flavour[j]
+                                                          + "_"
+                                                          + std::to_string(mu),
+                                                          mesPar);
+        }
+    }
+    
+    // execution
+    application.saveParameterFile("meson3pt.xml");
+    application.run();
+    
+    // epilogue
+    LOG(Message) << "Grid is finalizing now" << std::endl;
+    Grid_finalize();
+    
+    return EXIT_SUCCESS;
+}