Merge branch 'feature/hadrons' into feature/qed-fvol

# Conflicts:
#	extras/Hadrons/Modules.hpp
#	extras/Hadrons/Modules/MGauge/StochEm.cc
#	extras/Hadrons/modules.inc

extras/Hadrons/Environment.cc

@@ -270,7 +270,7 @@ int Environment::getObjectModule(const std::string name) const
 
 unsigned int Environment::getObjectLs(const unsigned int address) const
 {
-    if (hasObject(address))
+    if (hasCreatedObject(address))
     {
         return object_[address].Ls;
     }

extras/Hadrons/Exceptions.hpp

@@ -37,7 +37,6 @@ See the full license in the file "LICENSE" in the top level distribution directo
 #define SRC_LOC std::string(__FUNCTION__) + " at " + std::string(__FILE__) + ":"\
                 + std::to_string(__LINE__)
 #define HADRON_ERROR(exc, msg)\
-LOG(Error) << msg << std::endl;\
 throw(Exceptions::exc(msg, SRC_LOC));
 
 #define DECL_EXC(name, base) \

extras/Hadrons/GeneticScheduler.hpp

@@ -130,7 +130,7 @@ void GeneticScheduler<V, T>::nextGeneration(void)
     {
         initPopulation();
     }
-    LOG(Debug) << "Starting population:\n" << *this << std::endl;
+    //LOG(Debug) << "Starting population:\n" << *this << std::endl;
     
     // random mutations
     //PARALLEL_FOR_LOOP
@@ -138,7 +138,7 @@ void GeneticScheduler<V, T>::nextGeneration(void)
     {
         doMutation();
     }
-    LOG(Debug) << "After mutations:\n" << *this << std::endl;
+    //LOG(Debug) << "After mutations:\n" << *this << std::endl;
     
     // mating
     //PARALLEL_FOR_LOOP
@@ -146,14 +146,14 @@ void GeneticScheduler<V, T>::nextGeneration(void)
     {
         doCrossover();
     }
-    LOG(Debug) << "After mating:\n" << *this << std::endl;
+    //LOG(Debug) << "After mating:\n" << *this << std::endl;
     
     // grim reaper
     auto it = population_.begin();
     
     std::advance(it, par_.popSize);
     population_.erase(it, population_.end());
-    LOG(Debug) << "After grim reaper:\n" << *this << std::endl;
+    //LOG(Debug) << "After grim reaper:\n" << *this << std::endl;
 }
 
 // evolution steps /////////////////////////////////////////////////////////////

extras/Hadrons/Global.cc

@@ -37,20 +37,38 @@ HadronsLogger Hadrons::HadronsLogWarning(1,"Warning");
 HadronsLogger Hadrons::HadronsLogMessage(1,"Message");
 HadronsLogger Hadrons::HadronsLogIterative(1,"Iterative");
 HadronsLogger Hadrons::HadronsLogDebug(1,"Debug");
+HadronsLogger Hadrons::HadronsLogIRL(1,"IRL");
 
 void Hadrons::initLogger(void)
 {
-    auto w = std::string("Hadrons").length();
+    auto w  = std::string("Hadrons").length();
+    int  cw = 8;
+
+
     GridLogError.setTopWidth(w);
     GridLogWarning.setTopWidth(w);
     GridLogMessage.setTopWidth(w);
     GridLogIterative.setTopWidth(w);
     GridLogDebug.setTopWidth(w);
+    GridLogIRL.setTopWidth(w);
+    GridLogError.setChanWidth(cw);
+    GridLogWarning.setChanWidth(cw);
+    GridLogMessage.setChanWidth(cw);
+    GridLogIterative.setChanWidth(cw);
+    GridLogDebug.setChanWidth(cw);
+    GridLogIRL.setChanWidth(cw);
     HadronsLogError.Active(GridLogError.isActive());
     HadronsLogWarning.Active(GridLogWarning.isActive());
     HadronsLogMessage.Active(GridLogMessage.isActive());
     HadronsLogIterative.Active(GridLogIterative.isActive());
     HadronsLogDebug.Active(GridLogDebug.isActive());
+    HadronsLogIRL.Active(GridLogIRL.isActive());
+    HadronsLogError.setChanWidth(cw);
+    HadronsLogWarning.setChanWidth(cw);
+    HadronsLogMessage.setChanWidth(cw);
+    HadronsLogIterative.setChanWidth(cw);
+    HadronsLogDebug.setChanWidth(cw);
+    HadronsLogIRL.setChanWidth(cw);
 }
 
 // type utilities //////////////////////////////////////////////////////////////

extras/Hadrons/Global.hpp

@@ -58,6 +58,9 @@ using Grid::operator<<;
 #ifndef FIMPL
 #define FIMPL WilsonImplR
 #endif
+#ifndef ZFIMPL
+#define ZFIMPL ZWilsonImplR
+#endif
 #ifndef SIMPL
 #define SIMPL ScalarImplCR
 #endif
@@ -111,6 +114,7 @@ extern HadronsLogger HadronsLogWarning;
 extern HadronsLogger HadronsLogMessage;
 extern HadronsLogger HadronsLogIterative;
 extern HadronsLogger HadronsLogDebug;
+extern HadronsLogger HadronsLogIRL;
 
 void initLogger(void);
 
@@ -180,6 +184,18 @@ typedef XmlWriter ResultWriter;
 #define RESULT_FILE_NAME(name) \
 name + "." + std::to_string(vm().getTrajectory()) + "." + resultFileExt
 
+// default Schur convention
+
+#ifndef HADRONS_DEFAULT_SCHUR 
+#define HADRONS_DEFAULT_SCHUR DiagMooee
+#endif
+#define _HADRONS_SCHUR_OP_(conv) Schur##conv##Operator
+#define HADRONS_SCHUR_OP(conv) _HADRONS_SCHUR_OP_(conv)
+#define HADRONS_DEFAULT_SCHUR_OP HADRONS_SCHUR_OP(HADRONS_DEFAULT_SCHUR)
+#define _HADRONS_SCHUR_SOLVE_(conv) SchurRedBlack##conv##Solve
+#define HADRONS_SCHUR_SOLVE(conv) _HADRONS_SCHUR_SOLVE_(conv)
+#define HADRONS_DEFAULT_SCHUR_SOLVE HADRONS_SCHUR_SOLVE(HADRONS_DEFAULT_SCHUR)
+
 END_HADRONS_NAMESPACE
 
 #include <Grid/Hadrons/Exceptions.hpp>

extras/Hadrons/HadronsXmlRun.cc

@@ -56,14 +56,26 @@ int main(int argc, char *argv[])
     Grid_init(&argc, &argv);
     
     // execution
-    Application application(parameterFileName);
+    try
-    
-    application.parseParameterFile(parameterFileName);
-    if (!scheduleFileName.empty())
     {
-        application.loadSchedule(scheduleFileName);
+        Application application(parameterFileName);
+        
+        application.parseParameterFile(parameterFileName);
+        if (!scheduleFileName.empty())
+        {
+            application.loadSchedule(scheduleFileName);
+        }
+        application.run();
+    }
+    catch (const std::exception& e)
+    {
+        LOG(Error) << "FATAL ERROR -- Exception " << typeName(&typeid(e)) << std::endl;
+        LOG(Error) << e.what() << std::endl;
+        LOG(Error) << "Aborting program" << std::endl;
+        Grid_finalize();
+
+        return EXIT_FAILURE;
     }
-    application.run();
     
     // epilogue
     LOG(Message) << "Grid is finalizing now" << std::endl;

extras/Hadrons/HadronsXmlSchedule.cc

@@ -1,65 +0,0 @@
-/*************************************************************************************
-
-Grid physics library, www.github.com/paboyle/Grid 
-
-Source file: extras/Hadrons/HadronsXmlSchedule.cc
-
-Copyright (C) 2015-2018
-
-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 */
-
-#include <Grid/Hadrons/Application.hpp>
-
-using namespace Grid;
-using namespace QCD;
-using namespace Hadrons;
-
-int main(int argc, char *argv[])
-{
-    // parse command line
-    std::string parameterFileName, scheduleFileName;
-    
-    if (argc < 3)
-    {
-        std::cerr << "usage: " << argv[0] << " <parameter file> <schedule output> [Grid options]";
-        std::cerr << std::endl;
-        std::exit(EXIT_FAILURE);
-    }
-    parameterFileName = argv[1];
-    scheduleFileName  = argv[2];
-    
-    // initialization
-    Grid_init(&argc, &argv);
-    
-    // execution
-    Application application;
-    
-    application.parseParameterFile(parameterFileName);
-    application.schedule();
-    application.printSchedule();
-    application.saveSchedule(scheduleFileName);
-    
-    // epilogue
-    LOG(Message) << "Grid is finalizing now" << std::endl;
-    Grid_finalize();
-    
-    return EXIT_SUCCESS;
-}

extras/Hadrons/LanczosUtils.hpp

@@ -0,0 +1,115 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: extras/Hadrons/LanczosUtils.hpp
+
+Copyright (C) 2015-2018
+
+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_LanczosUtils_hpp_
+#define Hadrons_LanczosUtils_hpp_
+
+#include <Grid/Hadrons/Global.hpp>
+#include <Grid/algorithms/iterative/LocalCoherenceLanczos.h>
+
+BEGIN_HADRONS_NAMESPACE
+
+// Lanczos type
+#ifndef HADRONS_DEFAULT_LANCZOS_NBASIS
+#define HADRONS_DEFAULT_LANCZOS_NBASIS 60
+#endif
+
+template <typename T>
+struct EigenPack
+{
+    typedef T VectorType;
+    std::vector<RealD> eval;
+    std::vector<T>     evec;
+    
+    EigenPack(void) = default;
+
+    EigenPack(const size_t size, GridBase *grid)
+    {
+        resize(size, grid);
+    }
+
+    void resize(const size_t size, GridBase *grid)
+    {
+        eval.resize(size);
+        evec.resize(size, grid);
+    }
+
+    void read(const std::string fileStem)
+    {
+        std::string     evecFilename = fileStem + "_evec.bin";
+        std::string     evalFilename = fileStem + "_eval.xml";
+        emptyUserRecord record;
+        ScidacReader    binReader;
+        XmlReader       xmlReader(evalFilename);
+
+        LOG(Message) << "Reading " << evec.size() << " eigenvectors from '" 
+                     << evecFilename << "'" << std::endl;
+        binReader.open(evecFilename);
+        for(int k = 0; k < evec.size(); ++k) 
+        {
+            binReader.readScidacFieldRecord(evec[k], record);
+        }
+        binReader.close();
+        LOG(Message) << "Reading " << eval.size() << " eigenvalues from '" 
+                     << evalFilename << "'" << std::endl;
+        Grid::read(xmlReader, "evals", eval);
+    }
+
+    void write(const std::string fileStem)
+    {
+        std::string     evecFilename = fileStem + "_evec.bin";
+        std::string     evalFilename = fileStem + "_eval.xml";
+        emptyUserRecord record;
+        ScidacWriter    binWriter;
+        XmlWriter       xmlWriter(evalFilename);
+
+        LOG(Message) << "Writing " << evec.size() << " eigenvectors to '" 
+                     << evecFilename << "'" << std::endl;
+        binWriter.open(fileStem + "_evec.bin");
+        for(int k = 0; k < evec.size(); ++k) 
+        {
+            binWriter.writeScidacFieldRecord(evec[k], record);
+        }
+        binWriter.close();
+        LOG(Message) << "Writing " << eval.size() << " eigenvalues to '" 
+                     << evalFilename << "'" << std::endl;
+        Grid::write(xmlWriter, "evals", eval);
+    }
+};
+
+template <typename FImpl>
+using FineEigenPack = EigenPack<typename FImpl::FermionField>;
+
+template <typename FImpl, int nBasis>
+using CoarseEigenPack = EigenPack<
+    typename LocalCoherenceLanczos<typename FImpl::SiteSpinor, 
+                                   typename FImpl::SiteComplex, 
+                                   nBasis>::CoarseField>;
+
+END_HADRONS_NAMESPACE
+
+#endif // Hadrons_LanczosUtils_hpp_

extras/Hadrons/Makefile.am

@@ -1,5 +1,5 @@
 lib_LIBRARIES = libHadrons.a
-bin_PROGRAMS  = HadronsXmlRun HadronsXmlSchedule
+bin_PROGRAMS  = HadronsXmlRun
 
 include modules.inc
 
@@ -21,6 +21,7 @@ nobase_libHadrons_a_HEADERS = \
 	GeneticScheduler.hpp      \
 	Global.hpp                \
 	Graph.hpp                 \
+	LanczosUtils.hpp          \
 	Module.hpp                \
 	Modules.hpp               \
 	ModuleFactory.hpp         \
@@ -28,6 +29,3 @@ nobase_libHadrons_a_HEADERS = \
 
 HadronsXmlRun_SOURCES = HadronsXmlRun.cc
 HadronsXmlRun_LDADD   = libHadrons.a -lGrid
-
-HadronsXmlSchedule_SOURCES = HadronsXmlSchedule.cc
-HadronsXmlSchedule_LDADD   = libHadrons.a -lGrid

extras/Hadrons/Modules.hpp

@@ -7,7 +7,9 @@ Source file: extras/Hadrons/Modules.hpp
 Copyright (C) 2015-2018
 
 Author: Antonin Portelli <antonin.portelli@me.com>
+Author: Guido Cossu <guido.cossu@ed.ac.uk>
 Author: Lanny91 <andrew.lawson@gmail.com>
+Author: pretidav <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
@@ -43,12 +45,13 @@ See the full license in the file "LICENSE" in the top level distribution directo
 #include <Grid/Hadrons/Modules/MSource/SeqConserved.hpp>
 #include <Grid/Hadrons/Modules/MSink/Smear.hpp>
 #include <Grid/Hadrons/Modules/MSink/Point.hpp>
+#include <Grid/Hadrons/Modules/MSolver/LocalCoherenceLanczos.hpp>
 #include <Grid/Hadrons/Modules/MSolver/RBPrecCG.hpp>
 #include <Grid/Hadrons/Modules/MGauge/Unit.hpp>
 #include <Grid/Hadrons/Modules/MGauge/Random.hpp>
+#include <Grid/Hadrons/Modules/MGauge/FundtoHirep.hpp>
 #include <Grid/Hadrons/Modules/MGauge/StochEm.hpp>
 #include <Grid/Hadrons/Modules/MGauge/UnitEm.hpp>
-#include <Grid/Hadrons/Modules/MGauge/FundtoHirep.hpp>
 #include <Grid/Hadrons/Modules/MUtilities/TestSeqGamma.hpp>
 #include <Grid/Hadrons/Modules/MUtilities/TestSeqConserved.hpp>
 #include <Grid/Hadrons/Modules/MLoop/NoiseLoop.hpp>
@@ -60,9 +63,14 @@ See the full license in the file "LICENSE" in the top level distribution directo
 #include <Grid/Hadrons/Modules/MAction/DWF.hpp>
 #include <Grid/Hadrons/Modules/MAction/Wilson.hpp>
 #include <Grid/Hadrons/Modules/MAction/WilsonClover.hpp>
+#include <Grid/Hadrons/Modules/MAction/ZMobiusDWF.hpp>
 #include <Grid/Hadrons/Modules/MScalarSUN/Div.hpp>
 #include <Grid/Hadrons/Modules/MScalarSUN/TrMag.hpp>
+#include <Grid/Hadrons/Modules/MScalarSUN/EMT.hpp>
 #include <Grid/Hadrons/Modules/MScalarSUN/TwoPoint.hpp>
 #include <Grid/Hadrons/Modules/MScalarSUN/TrPhi.hpp>
+#include <Grid/Hadrons/Modules/MScalarSUN/Utils.hpp>
+#include <Grid/Hadrons/Modules/MScalarSUN/TransProj.hpp>
+#include <Grid/Hadrons/Modules/MScalarSUN/TrKinetic.hpp>
 #include <Grid/Hadrons/Modules/MIO/LoadNersc.hpp>
 #include <Grid/Hadrons/Modules/MIO/LoadBinary.hpp>

extras/Hadrons/Modules/MAction/WilsonClover.hpp

@@ -2,12 +2,13 @@
 
 Grid physics library, www.github.com/paboyle/Grid 
 
-Source file: extras/Hadrons/Modules/MAction/Wilson.hpp
+Source file: extras/Hadrons/Modules/MAction/WilsonClover.hpp
 
-Copyright (C) 2015
+Copyright (C) 2015-2018
-Copyright (C) 2016
 
 Author: Antonin Portelli <antonin.portelli@me.com>
+Author: Guido Cossu <guido.cossu@ed.ac.uk>
+Author: pretidav <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

extras/Hadrons/Modules/MAction/ZMobiusDWF.hpp

@@ -0,0 +1,143 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: extras/Hadrons/Modules/MAction/ZMobiusDWF.hpp
+
+Copyright (C) 2015-2018
+
+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_ZMobiusDWF_hpp_
+#define Hadrons_MAction_ZMobiusDWF_hpp_
+
+#include <Grid/Hadrons/Global.hpp>
+#include <Grid/Hadrons/Module.hpp>
+#include <Grid/Hadrons/ModuleFactory.hpp>
+
+BEGIN_HADRONS_NAMESPACE
+
+/******************************************************************************
+ *                         ZMobiusDWF                                 *
+ ******************************************************************************/
+BEGIN_MODULE_NAMESPACE(MAction)
+
+class ZMobiusDWFPar: Serializable
+{
+public:
+    GRID_SERIALIZABLE_CLASS_MEMBERS(ZMobiusDWFPar,
+                                    std::string                      , gauge,
+                                    unsigned int                     , Ls,
+                                    double                           , mass,
+                                    double                           , M5,
+                                    double                           , b,
+                                    double                           , c,
+                                    std::vector<std::complex<double>>, omega,
+                                    std::string                      , boundary);
+};
+
+template <typename FImpl>
+class TZMobiusDWF: public Module<ZMobiusDWFPar>
+{
+public:
+    FGS_TYPE_ALIASES(FImpl,);
+public:
+    // constructor
+    TZMobiusDWF(const std::string name);
+    // destructor
+    virtual ~TZMobiusDWF(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(ZMobiusDWF, TZMobiusDWF<ZFIMPL>, MAction);
+
+/******************************************************************************
+ *                 TZMobiusDWF implementation                             *
+ ******************************************************************************/
+// constructor /////////////////////////////////////////////////////////////////
+template <typename FImpl>
+TZMobiusDWF<FImpl>::TZMobiusDWF(const std::string name)
+: Module<ZMobiusDWFPar>(name)
+{}
+
+// dependencies/products ///////////////////////////////////////////////////////
+template <typename FImpl>
+std::vector<std::string> TZMobiusDWF<FImpl>::getInput(void)
+{
+    std::vector<std::string> in = {par().gauge};
+    
+    return in;
+}
+
+template <typename FImpl>
+std::vector<std::string> TZMobiusDWF<FImpl>::getOutput(void)
+{
+    std::vector<std::string> out = {getName()};
+    
+    return out;
+}
+
+// setup ///////////////////////////////////////////////////////////////////////
+template <typename FImpl>
+void TZMobiusDWF<FImpl>::setup(void)
+{
+    LOG(Message) << "Setting up z-Mobius domain wall fermion matrix with m= "
+                 << par().mass << ", M5= " << par().M5 << ", Ls= " << par().Ls 
+                 << ", b= " << par().b << ", c= " << par().c
+                 << " using gauge field '" << par().gauge << "'"
+                 << std::endl;
+    LOG(Message) << "Omegas: " << std::endl;
+    for (unsigned int i = 0; i < par().omega.size(); ++i)
+    {
+        LOG(Message) << "  omega[" << i << "]= " << par().omega[i] << std::endl;
+    }
+    LOG(Message) << "Fermion boundary conditions: " << par().boundary
+                 << std::endl;
+
+    env().createGrid(par().Ls);
+    auto &U    = envGet(LatticeGaugeField, par().gauge);
+    auto &g4   = *env().getGrid();
+    auto &grb4 = *env().getRbGrid();
+    auto &g5   = *env().getGrid(par().Ls);
+    auto &grb5 = *env().getRbGrid(par().Ls);
+    auto omega = par().omega;
+    std::vector<Complex> boundary = strToVec<Complex>(par().boundary);
+    typename ZMobiusFermion<FImpl>::ImplParams implParams(boundary);
+    envCreateDerived(FMat, ZMobiusFermion<FImpl>, getName(), par().Ls, U, g5,
+                     grb5, g4, grb4, par().mass, par().M5, omega,
+                     par().b, par().c, implParams);
+}
+
+// execution ///////////////////////////////////////////////////////////////////
+template <typename FImpl>
+void TZMobiusDWF<FImpl>::execute(void)
+{}
+
+END_MODULE_NAMESPACE
+
+END_HADRONS_NAMESPACE
+
+#endif // Hadrons_MAction_ZMobiusDWF_hpp_

extras/Hadrons/Modules/MFermion/GaugeProp.hpp

@@ -7,7 +7,9 @@ Source file: extras/Hadrons/Modules/MFermion/GaugeProp.hpp
 Copyright (C) 2015-2018
 
 Author: Antonin Portelli <antonin.portelli@me.com>
+Author: Guido Cossu <guido.cossu@ed.ac.uk>
 Author: Lanny91 <andrew.lawson@gmail.com>
+Author: pretidav <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

extras/Hadrons/Modules/MGauge/FundtoHirep.hpp

@@ -4,11 +4,10 @@ Grid physics library, www.github.com/paboyle/Grid
 
 Source file: extras/Hadrons/Modules/MGauge/FundtoHirep.hpp
 
-Copyright (C) 2015
+Copyright (C) 2015-2018
-Copyright (C) 2016
 
-Author: David Preti <david.preti@to.infn.it>
+Author: Antonin Portelli <antonin.portelli@me.com>
-	Guido Cossu <guido.cossu@ed.ac.uk>
+Author: pretidav <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

extras/Hadrons/Modules/MGauge/StochEm.hpp

@@ -7,6 +7,7 @@ Source file: extras/Hadrons/Modules/MGauge/StochEm.hpp
 Copyright (C) 2015-2018
 
 Author: Antonin Portelli <antonin.portelli@me.com>
+Author: Vera Guelpers <vmg1n14@soton.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

extras/Hadrons/Modules/MIO/LoadNersc.cc

@@ -71,6 +71,4 @@ void TLoadNersc::execute(void)
 
     auto &U = envGet(LatticeGaugeField, getName());
     NerscIO::readConfiguration(U, header, fileName);
-    LOG(Message) << "NERSC header:" << std::endl;
-    dump_meta_data(header, LOG(Message));
 }

extras/Hadrons/Modules/MScalarSUN/Div.hpp

@@ -31,18 +31,18 @@ See the full license in the file "LICENSE" in the top level distribution directo
 #include <Grid/Hadrons/Global.hpp>
 #include <Grid/Hadrons/Module.hpp>
 #include <Grid/Hadrons/ModuleFactory.hpp>
+#include <Grid/Hadrons/Modules/MScalarSUN/Utils.hpp>
 
 BEGIN_HADRONS_NAMESPACE
 
 /******************************************************************************
- *                         Div                                 *
+ *                       Divergence of a vector field                         *
  ******************************************************************************/
 BEGIN_MODULE_NAMESPACE(MScalarSUN)
 
 class DivPar: Serializable
 {
 public:
-    GRID_SERIALIZABLE_ENUM(DiffType, undef, forward, 1, backward, 2, central, 3);
     GRID_SERIALIZABLE_CLASS_MEMBERS(DivPar,
                                     std::vector<std::string>, op,
                                     DiffType,                 type,
@@ -59,8 +59,8 @@ public:
     {
     public:
         GRID_SERIALIZABLE_CLASS_MEMBERS(Result,
-                                        DivPar::DiffType, type,
+                                        DiffType, type,
-                                        Complex,          value);
+                                        Complex,  value);
     };
 public:
     // constructor
@@ -83,7 +83,7 @@ MODULE_REGISTER_NS(DivSU5, TDiv<ScalarNxNAdjImplR<5>>, MScalarSUN);
 MODULE_REGISTER_NS(DivSU6, TDiv<ScalarNxNAdjImplR<6>>, MScalarSUN);
 
 /******************************************************************************
- *                 TDiv implementation                             *
+ *                           TDiv implementation                              *
  ******************************************************************************/
 // constructor /////////////////////////////////////////////////////////////////
 template <typename SImpl>
@@ -135,18 +135,7 @@ void TDiv<SImpl>::execute(void)
     for (unsigned int mu = 0; mu < nd; ++mu)
     {
         auto &op = envGet(ComplexField, par().op[mu]);
-        switch(par().type)
+        dmuAcc(div, op, mu, par().type);
-        {
-            case DivPar::DiffType::backward:
-                div += op - Cshift(op, mu, -1);
-                break;
-            case DivPar::DiffType::forward:
-                div += Cshift(op, mu, 1) - op;
-                break;
-            case DivPar::DiffType::central:
-                div += 0.5*(Cshift(op, mu, 1) - Cshift(op, mu, -1));
-                break;
-        }
     }
     if (!par().output.empty())
     {

extras/Hadrons/Modules/MScalarSUN/EMT.hpp

@@ -0,0 +1,181 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: extras/Hadrons/Modules/MScalarSUN/EMT.hpp
+
+Copyright (C) 2015-2018
+
+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_MScalarSUN_EMT_hpp_
+#define Hadrons_MScalarSUN_EMT_hpp_
+
+#include <Grid/Hadrons/Global.hpp>
+#include <Grid/Hadrons/Module.hpp>
+#include <Grid/Hadrons/ModuleFactory.hpp>
+#include <Grid/Hadrons/Modules/MScalarSUN/Utils.hpp>
+
+BEGIN_HADRONS_NAMESPACE
+
+/******************************************************************************
+ *                         Energy-momentum tensor                             *
+ ******************************************************************************/
+BEGIN_MODULE_NAMESPACE(MScalarSUN)
+
+class EMTPar: Serializable
+{
+public:
+    GRID_SERIALIZABLE_CLASS_MEMBERS(EMTPar,
+                                    std::string, kinetic,
+                                    std::string, phiPow,
+                                    std::string, improvement,
+                                    double     , m2,
+                                    double     , lambda,
+                                    double     , g,
+                                    double     , xi,
+                                    std::string, output);
+};
+
+template <typename SImpl>
+class TEMT: public Module<EMTPar>
+{
+public:
+    typedef typename SImpl::Field        Field;
+    typedef typename SImpl::ComplexField ComplexField;
+public:
+    // constructor
+    TEMT(const std::string name);
+    // destructor
+    virtual ~TEMT(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(EMTSU2, TEMT<ScalarNxNAdjImplR<2>>, MScalarSUN);
+MODULE_REGISTER_NS(EMTSU3, TEMT<ScalarNxNAdjImplR<3>>, MScalarSUN);
+MODULE_REGISTER_NS(EMTSU4, TEMT<ScalarNxNAdjImplR<4>>, MScalarSUN);
+MODULE_REGISTER_NS(EMTSU5, TEMT<ScalarNxNAdjImplR<5>>, MScalarSUN);
+MODULE_REGISTER_NS(EMTSU6, TEMT<ScalarNxNAdjImplR<6>>, MScalarSUN);
+
+/******************************************************************************
+ *                           TEMT implementation                              *
+ ******************************************************************************/
+// constructor /////////////////////////////////////////////////////////////////
+template <typename SImpl>
+TEMT<SImpl>::TEMT(const std::string name)
+: Module<EMTPar>(name)
+{}
+
+// dependencies/products ///////////////////////////////////////////////////////
+template <typename SImpl>
+std::vector<std::string> TEMT<SImpl>::getInput(void)
+{
+    std::vector<std::string> in;
+    
+    for (unsigned int mu = 0; mu < env().getNd(); ++mu)
+    for (unsigned int nu = mu; nu < env().getNd(); ++nu)
+    {
+        in.push_back(varName(par().kinetic, mu, nu));
+        in.push_back(varName(par().improvement, mu, nu));
+    }
+    in.push_back(varName(par().phiPow, 2));
+    in.push_back(varName(par().phiPow, 4));
+
+    return in;
+}
+
+template <typename SImpl>
+std::vector<std::string> TEMT<SImpl>::getOutput(void)
+{
+    std::vector<std::string> out;
+    
+    for (unsigned int mu = 0; mu < env().getNd(); ++mu)
+    for (unsigned int nu = mu; nu < env().getNd(); ++nu)
+    {
+        out.push_back(varName(getName(), mu, nu));
+    }
+
+    return out;
+}
+
+// setup ///////////////////////////////////////////////////////////////////////
+template <typename SImpl>
+void TEMT<SImpl>::setup(void)
+{
+    for (unsigned int mu = 0; mu < env().getNd(); ++mu)
+    for (unsigned int nu = mu; nu < env().getNd(); ++nu)
+    {
+        envCreateLat(ComplexField, varName(getName(), mu, nu));
+    }
+    envTmpLat(ComplexField, "sumkin");
+}
+
+// execution ///////////////////////////////////////////////////////////////////
+template <typename SImpl>
+void TEMT<SImpl>::execute(void)
+{
+    LOG(Message) << "Computing energy-momentum tensor" << std::endl;
+    LOG(Message) << "  kinetic terms: '" << par().kinetic << "'" << std::endl;
+    LOG(Message) << "      tr(phi^n): '" << par().phiPow << "'" << std::endl;
+    LOG(Message) << "    improvement: '" << par().improvement << "'" << std::endl;
+    LOG(Message) << "            m^2= " << par().m2 << std::endl;
+    LOG(Message) << "         lambda= " << par().lambda << std::endl;
+    LOG(Message) << "              g= " << par().g << std::endl;
+    LOG(Message) << "             xi= " << par().xi << std::endl;
+
+    const unsigned int N       = SImpl::Group::Dimension;
+    auto               &trphi2 = envGet(ComplexField, varName(par().phiPow, 2));
+    auto               &trphi4 = envGet(ComplexField, varName(par().phiPow, 4));
+
+    envGetTmp(ComplexField, sumkin);
+    sumkin = zero;
+    for (unsigned int mu = 0; mu < env().getNd(); ++mu)
+    {
+        auto &trkin = envGet(ComplexField, varName(par().kinetic, mu, mu));
+
+        sumkin += trkin;
+    }
+    for (unsigned int mu = 0; mu < env().getNd(); ++mu)
+    for (unsigned int nu = mu; nu < env().getNd(); ++nu)
+    {
+        auto &out   = envGet(ComplexField, varName(getName(), mu, nu));
+        auto &trkin = envGet(ComplexField, varName(par().kinetic, mu, nu));
+        auto &imp   = envGet(ComplexField, varName(par().improvement, mu, nu));
+
+        out = 2.*trkin + par().xi*imp;
+        if (mu == nu)
+        {
+            out -= sumkin + par().m2*trphi2 + par().lambda*trphi4;
+        }
+        out *= N/par().g;
+    }
+}
+
+END_MODULE_NAMESPACE
+
+END_HADRONS_NAMESPACE
+
+#endif // Hadrons_MScalarSUN_EMT_hpp_

extras/Hadrons/Modules/MScalarSUN/TrKinetic.hpp

@@ -0,0 +1,170 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: extras/Hadrons/Modules/MScalarSUN/TrKinetic.hpp
+
+Copyright (C) 2015-2018
+
+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_MScalarSUN_TrKinetic_hpp_
+#define Hadrons_MScalarSUN_TrKinetic_hpp_
+
+#include <Grid/Hadrons/Global.hpp>
+#include <Grid/Hadrons/Module.hpp>
+#include <Grid/Hadrons/ModuleFactory.hpp>
+#include <Grid/Hadrons/Modules/MScalarSUN/Utils.hpp>
+
+BEGIN_HADRONS_NAMESPACE
+
+/******************************************************************************
+ *                         Trace of kinetic term                              *
+ ******************************************************************************/
+BEGIN_MODULE_NAMESPACE(MScalarSUN)
+
+class TrKineticPar: Serializable
+{
+public:
+    GRID_SERIALIZABLE_CLASS_MEMBERS(TrKineticPar,
+                                    std::string,  field,
+                                    DiffType,     type,
+                                    std::string,  output);
+};
+
+template <typename SImpl>
+class TTrKinetic: public Module<TrKineticPar>
+{
+public:
+    typedef typename SImpl::Field        Field;
+    typedef typename SImpl::ComplexField ComplexField;
+    class Result: Serializable
+    {
+    public:
+        GRID_SERIALIZABLE_CLASS_MEMBERS(Result,
+                                        std::string, op,
+                                        Complex    , value);
+    };
+public:
+    // constructor
+    TTrKinetic(const std::string name);
+    // destructor
+    virtual ~TTrKinetic(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(TrKineticSU2, TTrKinetic<ScalarNxNAdjImplR<2>>, MScalarSUN);
+MODULE_REGISTER_NS(TrKineticSU3, TTrKinetic<ScalarNxNAdjImplR<3>>, MScalarSUN);
+MODULE_REGISTER_NS(TrKineticSU4, TTrKinetic<ScalarNxNAdjImplR<4>>, MScalarSUN);
+MODULE_REGISTER_NS(TrKineticSU5, TTrKinetic<ScalarNxNAdjImplR<5>>, MScalarSUN);
+MODULE_REGISTER_NS(TrKineticSU6, TTrKinetic<ScalarNxNAdjImplR<6>>, MScalarSUN);
+
+/******************************************************************************
+ *                      TTrKinetic implementation                             *
+ ******************************************************************************/
+// constructor /////////////////////////////////////////////////////////////////
+template <typename SImpl>
+TTrKinetic<SImpl>::TTrKinetic(const std::string name)
+: Module<TrKineticPar>(name)
+{}
+
+// dependencies/products ///////////////////////////////////////////////////////
+template <typename SImpl>
+std::vector<std::string> TTrKinetic<SImpl>::getInput(void)
+{
+    std::vector<std::string> in = {par().field};
+    
+    return in;
+}
+
+template <typename SImpl>
+std::vector<std::string> TTrKinetic<SImpl>::getOutput(void)
+{
+    std::vector<std::string> out ;
+
+    for (unsigned int mu = 0; mu < env().getNd(); ++mu)
+    for (unsigned int nu = mu; nu < env().getNd(); ++nu)
+    {
+        out.push_back(varName(getName(), mu, nu));
+    }
+    
+    return out;
+}
+
+// setup ///////////////////////////////////////////////////////////////////////
+template <typename SImpl>
+void TTrKinetic<SImpl>::setup(void)
+{
+    for (unsigned int mu = 0; mu < env().getNd(); ++mu)
+    for (unsigned int nu = mu; nu < env().getNd(); ++nu)
+    {
+        envCreateLat(ComplexField, varName(getName(), mu, nu));
+    }
+    envTmp(std::vector<Field>, "der", 1, env().getNd(), env().getGrid());
+}
+
+// execution ///////////////////////////////////////////////////////////////////
+template <typename SImpl>
+void TTrKinetic<SImpl>::execute(void)
+{
+    LOG(Message) << "Computing tr(d_mu phi*d_nu phi) using " << par().type
+                 << " derivative" << std::endl; 
+
+    std::vector<Result> result;
+    auto                &phi = envGet(Field, par().field);
+
+    envGetTmp(std::vector<Field>, der);
+    for (unsigned int mu = 0; mu < env().getNd(); ++mu)
+    {
+        dmu(der[mu], phi, mu, par().type);
+    }
+    for (unsigned int mu = 0; mu < env().getNd(); ++mu)
+    for (unsigned int nu = mu; nu < env().getNd(); ++nu)
+    {
+        auto &out = envGet(ComplexField, varName(getName(), mu, nu));
+
+        out = -trace(der[mu]*der[nu]);
+        if (!par().output.empty())
+        {
+            Result r;
+
+            r.op    = "tr(d_" + std::to_string(mu) + "phi*d_" 
+                      + std::to_string(nu) + "phi)";
+            r.value = TensorRemove(sum(out));
+            result.push_back(r);
+        }
+    }
+    if (result.size() > 0)
+    {
+        saveResult(par().output, "trkinetic", result);
+    }
+}
+
+END_MODULE_NAMESPACE
+
+END_HADRONS_NAMESPACE
+
+#endif // Hadrons_MScalarSUN_TrKinetic_hpp_

extras/Hadrons/Modules/MScalarSUN/TrMag.hpp

@@ -31,11 +31,12 @@ See the full license in the file "LICENSE" in the top level distribution directo
 #include <Grid/Hadrons/Global.hpp>
 #include <Grid/Hadrons/Module.hpp>
 #include <Grid/Hadrons/ModuleFactory.hpp>
+#include <Grid/Hadrons/Modules/MScalarSUN/Utils.hpp>
 
 BEGIN_HADRONS_NAMESPACE
 
 /******************************************************************************
- *                       Module to compute tr(mag^n)                          *
+ *                     Trace of powers of the magnetisation                   *
  ******************************************************************************/
 BEGIN_MODULE_NAMESPACE(MScalarSUN)
 
@@ -117,7 +118,7 @@ template <typename SImpl>
 void TTrMag<SImpl>::execute(void)
 {
     LOG(Message) << "Computing tr(mag^n) for n even up to " << par().maxPow
-                 << "..." << std::endl;
+                 << std::endl;
 
     std::vector<Result> result;
     auto                &phi = envGet(Field, par().field);

extras/Hadrons/Modules/MScalarSUN/TrPhi.hpp

@@ -31,11 +31,12 @@ See the full license in the file "LICENSE" in the top level distribution directo
 #include <Grid/Hadrons/Global.hpp>
 #include <Grid/Hadrons/Module.hpp>
 #include <Grid/Hadrons/ModuleFactory.hpp>
+#include <Grid/Hadrons/Modules/MScalarSUN/Utils.hpp>
 
 BEGIN_HADRONS_NAMESPACE
 
 /******************************************************************************
- *                         Module to compute tr(phi^n)                        *
+ *                      Trace of powers of a scalar field                     *
  ******************************************************************************/
 BEGIN_MODULE_NAMESPACE(MScalarSUN)
 
@@ -73,9 +74,6 @@ public:
     virtual void setup(void);
     // execution
     virtual void execute(void);
-private:
-    // output name generator
-    std::string outName(const unsigned int n);
 };
 
 MODULE_REGISTER_NS(TrPhiSU2, TTrPhi<ScalarNxNAdjImplR<2>>, MScalarSUN);
@@ -109,7 +107,7 @@ std::vector<std::string> TTrPhi<SImpl>::getOutput(void)
 
     for (unsigned int n = 2; n <= par().maxPow; n += 2)
     {
-        out.push_back(outName(n));
+        out.push_back(varName(getName(), n));
     }
     
     return out;
@@ -127,7 +125,7 @@ void TTrPhi<SImpl>::setup(void)
     envTmpLat(Field, "buf");
     for (unsigned int n = 2; n <= par().maxPow; n += 2)
     {
-        envCreateLat(ComplexField, outName(n));
+        envCreateLat(ComplexField, varName(getName(), n));
     }
 }
 
@@ -136,7 +134,7 @@ template <typename SImpl>
 void TTrPhi<SImpl>::execute(void)
 {
     LOG(Message) << "Computing tr(phi^n) for n even up to " << par().maxPow
-                 << "..." << std::endl; 
+                 << std::endl; 
 
     std::vector<Result> result;
     auto                &phi = envGet(Field, par().field);
@@ -147,7 +145,7 @@ void TTrPhi<SImpl>::execute(void)
     phi2 = -phi*phi; 
     for (unsigned int n = 2; n <= par().maxPow; n += 2)
     {
-        auto &phin = envGet(ComplexField, outName(n));
+        auto &phin = envGet(ComplexField, varName(getName(), n));
 
         buf  = buf*phi2;
         phin = trace(buf);
@@ -166,13 +164,6 @@ void TTrPhi<SImpl>::execute(void)
     }
 }
 
-// output name generator ///////////////////////////////////////////////////////
-template <typename SImpl>
-std::string TTrPhi<SImpl>::outName(const unsigned int n)
-{
-    return getName() + "_" + std::to_string(n);
-}
-
 END_MODULE_NAMESPACE
 
 END_HADRONS_NAMESPACE

extras/Hadrons/Modules/MScalarSUN/TransProj.hpp

@@ -0,0 +1,185 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: extras/Hadrons/Modules/MScalarSUN/TransProj.hpp
+
+Copyright (C) 2015-2018
+
+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_MScalarSUN_TransProj_hpp_
+#define Hadrons_MScalarSUN_TransProj_hpp_
+
+#include <Grid/Hadrons/Global.hpp>
+#include <Grid/Hadrons/Module.hpp>
+#include <Grid/Hadrons/ModuleFactory.hpp>
+#include <Grid/Hadrons/Modules/MScalarSUN/Utils.hpp>
+
+BEGIN_HADRONS_NAMESPACE
+
+/******************************************************************************
+ *                         Transverse projection                              *
+ ******************************************************************************/
+BEGIN_MODULE_NAMESPACE(MScalarSUN)
+
+class TransProjPar: Serializable
+{
+public:
+    GRID_SERIALIZABLE_CLASS_MEMBERS(TransProjPar,
+                                    std::string,  op,
+                                    DiffType,     type,
+                                    std::string,  output);
+};
+
+template <typename SImpl>
+class TTransProj: public Module<TransProjPar>
+{
+public:
+    typedef typename SImpl::Field        Field;
+    typedef typename SImpl::ComplexField ComplexField;
+    class Result: Serializable
+    {
+    public:
+        GRID_SERIALIZABLE_CLASS_MEMBERS(Result,
+                                        std::string, op,
+                                        Complex    , value);
+    };
+public:
+    // constructor
+    TTransProj(const std::string name);
+    // destructor
+    virtual ~TTransProj(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(TransProjSU2, TTransProj<ScalarNxNAdjImplR<2>>, MScalarSUN);
+MODULE_REGISTER_NS(TransProjSU3, TTransProj<ScalarNxNAdjImplR<3>>, MScalarSUN);
+MODULE_REGISTER_NS(TransProjSU4, TTransProj<ScalarNxNAdjImplR<4>>, MScalarSUN);
+MODULE_REGISTER_NS(TransProjSU5, TTransProj<ScalarNxNAdjImplR<5>>, MScalarSUN);
+MODULE_REGISTER_NS(TransProjSU6, TTransProj<ScalarNxNAdjImplR<6>>, MScalarSUN);
+
+/******************************************************************************
+ *                        TTransProj implementation                           *
+ ******************************************************************************/
+// constructor /////////////////////////////////////////////////////////////////
+template <typename SImpl>
+TTransProj<SImpl>::TTransProj(const std::string name)
+: Module<TransProjPar>(name)
+{}
+
+// dependencies/products ///////////////////////////////////////////////////////
+template <typename SImpl>
+std::vector<std::string> TTransProj<SImpl>::getInput(void)
+{
+    std::vector<std::string> in = {par().op};
+    
+    return in;
+}
+
+template <typename SImpl>
+std::vector<std::string> TTransProj<SImpl>::getOutput(void)
+{
+    std::vector<std::string> out;
+    
+    for (unsigned int mu = 0; mu < env().getNd(); ++mu)
+    for (unsigned int nu = mu; nu < env().getNd(); ++nu)
+    {
+        out.push_back(varName(getName(), mu, nu));
+    }
+    
+    return out;
+}
+
+// setup ///////////////////////////////////////////////////////////////////////
+template <typename SImpl>
+void TTransProj<SImpl>::setup(void)
+{
+    for (unsigned int mu = 0; mu < env().getNd(); ++mu)
+    for (unsigned int nu = mu; nu < env().getNd(); ++nu)
+    {
+        envCreateLat(ComplexField, varName(getName(), mu, nu));
+    }
+    envTmpLat(ComplexField, "buf1");
+    envTmpLat(ComplexField, "buf2");
+    envTmpLat(ComplexField, "lap");
+}
+
+// execution ///////////////////////////////////////////////////////////////////
+template <typename SImpl>
+void TTransProj<SImpl>::execute(void)
+{
+    LOG(Message) << "Computing (delta_mu,nu d^2 - d_mu*d_nu)*op using " 
+                 << par().type << " derivatives and op= '" << par().op 
+                 << "'" << std::endl; 
+
+    std::vector<Result> result;
+    auto                &op = envGet(ComplexField, par().op);
+
+    envGetTmp(ComplexField, buf1);
+    envGetTmp(ComplexField, buf2);
+    envGetTmp(ComplexField, lap);
+    lap = zero;
+    for (unsigned int mu = 0; mu < env().getNd(); ++mu)
+    {
+        dmu(buf1, op, mu, par().type);
+        dmu(buf2, buf1, mu, par().type);
+        lap += buf2;
+    }
+    for (unsigned int mu = 0; mu < env().getNd(); ++mu)
+    for (unsigned int nu = mu; nu < env().getNd(); ++nu)
+    {
+        auto &out = envGet(ComplexField, varName(getName(), mu, nu));
+        dmu(buf1, op, mu, par().type);
+        dmu(buf2, buf1, nu, par().type);
+        out = -buf2;
+        if (mu == nu)
+        {
+            out += lap;
+        }
+        if (!par().output.empty())
+        {
+            Result r;
+
+            r.op    = "(delta_" + std::to_string(mu) + "," + std::to_string(nu)
+                      + " d^2 - d_" + std::to_string(mu) + "*d_" 
+                      + std::to_string(nu) + ")*op";
+            r.value = TensorRemove(sum(out));
+            result.push_back(r);
+        }
+    }
+    if (result.size() > 0)
+    {
+        saveResult(par().output, "transproj", result);
+    }
+}
+
+END_MODULE_NAMESPACE
+
+END_HADRONS_NAMESPACE
+
+#endif // Hadrons_MScalarSUN_TransProj_hpp_

extras/Hadrons/Modules/MScalarSUN/TwoPoint.hpp

@@ -31,6 +31,7 @@ See the full license in the file "LICENSE" in the top level distribution directo
 #include <Grid/Hadrons/Global.hpp>
 #include <Grid/Hadrons/Module.hpp>
 #include <Grid/Hadrons/ModuleFactory.hpp>
+#include <Grid/Hadrons/Modules/MScalarSUN/Utils.hpp>
 
 BEGIN_HADRONS_NAMESPACE
 
@@ -87,7 +88,7 @@ MODULE_REGISTER_NS(TwoPointSU5, TTwoPoint<ScalarNxNAdjImplR<5>>, MScalarSUN);
 MODULE_REGISTER_NS(TwoPointSU6, TTwoPoint<ScalarNxNAdjImplR<6>>, MScalarSUN);
 
 /******************************************************************************
- *                 TTwoPoint implementation                             *
+ *                       TTwoPoint implementation                             *
  ******************************************************************************/
 // constructor /////////////////////////////////////////////////////////////////
 template <typename SImpl>

extras/Hadrons/Modules/MScalarSUN/Utils.hpp

@@ -0,0 +1,107 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: extras/Hadrons/Modules/MScalarSUN/Utils.hpp
+
+Copyright (C) 2015-2018
+
+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_MScalarSUN_Utils_hpp_
+#define Hadrons_MScalarSUN_Utils_hpp_
+
+#include <Grid/Hadrons/Global.hpp>
+#include <Grid/Hadrons/Module.hpp>
+
+BEGIN_HADRONS_NAMESPACE
+
+BEGIN_MODULE_NAMESPACE(MScalarSUN)
+
+GRID_SERIALIZABLE_ENUM(DiffType, undef, forward, 1, backward, 2, central, 3);
+
+template <typename Field>
+inline void dmu(Field &out, const Field &in, const unsigned int mu, const DiffType type)
+{
+    auto & env = Environment::getInstance();
+
+    if (mu >= env.getNd())
+    {
+        HADRON_ERROR(Range, "Derivative direction out of range");
+    }
+    switch(type)
+    {
+        case DiffType::backward:
+            out = in - Cshift(in, mu, -1);
+            break;
+        case DiffType::forward:
+            out = Cshift(in, mu, 1) - in;
+            break;
+        case DiffType::central:
+            out = 0.5*(Cshift(in, mu, 1) - Cshift(in, mu, -1));
+            break;
+        default:
+            HADRON_ERROR(Argument, "Derivative type invalid");
+            break;
+    }
+}
+
+template <typename Field>
+inline void dmuAcc(Field &out, const Field &in, const unsigned int mu, const DiffType type)
+{
+    auto & env = Environment::getInstance();
+
+    if (mu >= env.getNd())
+    {
+        HADRON_ERROR(Range, "Derivative direction out of range");
+    }
+    switch(type)
+    {
+        case DiffType::backward:
+            out += in - Cshift(in, mu, -1);
+            break;
+        case DiffType::forward:
+            out += Cshift(in, mu, 1) - in;
+            break;
+        case DiffType::central:
+            out += 0.5*(Cshift(in, mu, 1) - Cshift(in, mu, -1));
+            break;
+        default:
+            HADRON_ERROR(Argument, "Derivative type invalid");
+            break;
+    }
+}
+
+inline std::string varName(const std::string name, const unsigned int mu)
+{
+    return name + "_" + std::to_string(mu);
+}
+
+inline std::string varName(const std::string name, const unsigned int mu, 
+                           const unsigned int nu)
+{
+    return name + "_" + std::to_string(mu) + "_" + std::to_string(nu);
+}
+
+END_MODULE_NAMESPACE
+
+END_HADRONS_NAMESPACE
+
+#endif // Hadrons_MScalarSUN_Utils_hpp_

extras/Hadrons/Modules/MSolver/LocalCoherenceLanczos.hpp

@@ -0,0 +1,249 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: extras/Hadrons/Modules/MSolver/LocalCoherenceLanczos.hpp
+
+Copyright (C) 2015-2018
+
+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_MSolver_LocalCoherenceLanczos_hpp_
+#define Hadrons_MSolver_LocalCoherenceLanczos_hpp_
+
+#include <Grid/Hadrons/Global.hpp>
+#include <Grid/Hadrons/Module.hpp>
+#include <Grid/Hadrons/ModuleFactory.hpp>
+#include <Grid/Hadrons/LanczosUtils.hpp>
+
+BEGIN_HADRONS_NAMESPACE
+
+/******************************************************************************
+ *                    Local coherence Lanczos eigensolver                     *
+ *****************************************************************************/
+BEGIN_MODULE_NAMESPACE(MSolver)
+
+class LocalCoherenceLanczosPar: Serializable
+{
+public:
+    GRID_SERIALIZABLE_CLASS_MEMBERS(LocalCoherenceLanczosPar,
+                                    std::string,   action,
+                                    int,           doFine,
+                                    int,           doCoarse,
+                                    LanczosParams, fineParams,
+                                    LanczosParams, coarseParams,
+                                    ChebyParams,   smoother,
+                                    RealD,         coarseRelaxTol,
+                                    std::string,   blockSize,
+                                    std::string,   output);
+};
+
+template <typename FImpl, int nBasis>
+class TLocalCoherenceLanczos: public Module<LocalCoherenceLanczosPar>
+{
+public:
+    FERM_TYPE_ALIASES(FImpl,);
+    typedef LocalCoherenceLanczos<typename FImpl::SiteSpinor, 
+                                  typename FImpl::SiteComplex, 
+                                  nBasis>                LCL;
+    typedef FineEigenPack<FImpl>                         FinePack;
+    typedef CoarseEigenPack<FImpl, nBasis>               CoarsePack; 
+    typedef HADRONS_DEFAULT_SCHUR_OP<FMat, FermionField> SchurFMat;
+public:
+    // constructor
+    TLocalCoherenceLanczos(const std::string name);
+    // destructor
+    virtual ~TLocalCoherenceLanczos(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);
+private:
+    void makeCoarseGrid(void);
+private:
+    std::vector<int>                       coarseDim_;
+    int                                    Ls_, cLs_{1};
+    std::unique_ptr<GridCartesian>         coarseGrid4_{nullptr};
+    std::unique_ptr<GridCartesian>         coarseGrid_{nullptr};
+    std::unique_ptr<GridRedBlackCartesian> coarseGrid4Rb_{nullptr};
+    std::unique_ptr<GridRedBlackCartesian> coarseGridRb_{nullptr};
+    std::string                            fineName_, coarseName_;
+};
+
+MODULE_REGISTER_NS(LocalCoherenceLanczos, 
+    ARG(TLocalCoherenceLanczos<FIMPL, HADRONS_DEFAULT_LANCZOS_NBASIS>), 
+    MSolver);
+MODULE_REGISTER_NS(ZLocalCoherenceLanczos, 
+    ARG(TLocalCoherenceLanczos<ZFIMPL, HADRONS_DEFAULT_LANCZOS_NBASIS>), 
+    MSolver);
+
+/******************************************************************************
+ *                 TLocalCoherenceLanczos implementation                      *
+ ******************************************************************************/
+// constructor /////////////////////////////////////////////////////////////////
+template <typename FImpl, int nBasis>
+TLocalCoherenceLanczos<FImpl, nBasis>::TLocalCoherenceLanczos(const std::string name)
+: Module<LocalCoherenceLanczosPar>(name)
+{
+    fineName_   = getName() + "_fine";
+    coarseName_ = getName() + "_coarse";
+}
+
+// dependencies/products ///////////////////////////////////////////////////////
+template <typename FImpl, int nBasis>
+std::vector<std::string> TLocalCoherenceLanczos<FImpl, nBasis>::getInput(void)
+{
+    std::vector<std::string> in = {par().action};
+    
+    return in;
+}
+
+template <typename FImpl, int nBasis>
+std::vector<std::string> TLocalCoherenceLanczos<FImpl, nBasis>::getOutput(void)
+{
+    std::vector<std::string> out = {fineName_, coarseName_};
+    
+    return out;
+}
+
+// setup ///////////////////////////////////////////////////////////////////////
+template <typename FImpl, int nBasis>
+void TLocalCoherenceLanczos<FImpl, nBasis>::makeCoarseGrid(void)
+{
+    int              nd        = env().getNd();
+    std::vector<int> blockSize = strToVec<int>(par().blockSize);
+    auto             fineDim   = env().getDim();
+
+    Ls_ = env().getObjectLs(par().action);
+    env().createGrid(Ls_);
+    coarseDim_.resize(nd);
+    for (int d = 0; d < coarseDim_.size(); d++)
+    {
+        coarseDim_[d] = fineDim[d]/blockSize[d];
+        if (coarseDim_[d]*blockSize[d] != fineDim[d])
+        {
+            HADRON_ERROR(Size, "Fine dimension " + std::to_string(d) 
+                         + " (" + std::to_string(fineDim[d]) 
+                         + ") not divisible by coarse dimension ("
+                         + std::to_string(coarseDim_[d]) + ")"); 
+        }
+    }
+    if (blockSize.size() > nd)
+    {
+        cLs_ = Ls_/blockSize[nd];
+        if (cLs_*blockSize[nd] != Ls_)
+        {
+            HADRON_ERROR(Size, "Fine Ls (" + std::to_string(Ls_) 
+                         + ") not divisible by coarse Ls ("
+                         + std::to_string(cLs_) + ")");
+        }
+    }
+    if (Ls_ > 1)
+    {
+        coarseGrid4_.reset(SpaceTimeGrid::makeFourDimGrid(
+            coarseDim_, GridDefaultSimd(nd, vComplex::Nsimd()),
+            GridDefaultMpi()));
+        coarseGrid4Rb_.reset(SpaceTimeGrid::makeFourDimRedBlackGrid(coarseGrid4_.get()));
+        coarseGrid_.reset(SpaceTimeGrid::makeFiveDimGrid(cLs_, coarseGrid4_.get()));
+        coarseGridRb_.reset(SpaceTimeGrid::makeFiveDimRedBlackGrid(cLs_, coarseGrid4_.get()));
+    }
+    else
+    {
+        coarseGrid_.reset(SpaceTimeGrid::makeFourDimGrid(
+            coarseDim_, GridDefaultSimd(nd, vComplex::Nsimd()),
+            GridDefaultMpi()));
+        coarseGridRb_.reset(SpaceTimeGrid::makeFourDimRedBlackGrid(coarseGrid_.get()));
+    }
+}
+
+template <typename FImpl, int nBasis>
+void TLocalCoherenceLanczos<FImpl, nBasis>::setup(void)
+{
+    LOG(Message) << "Setting up local coherence Lanczos eigensolver for"
+                 << " action '" << par().action << "' (" << nBasis
+                 << " eigenvectors)..." << std::endl;
+    
+    if (!coarseGrid_)
+    {
+        makeCoarseGrid();
+    }
+    LOG(Message) << "Coarse grid: " << coarseGrid_->GlobalDimensions() << std::endl;
+    envCreate(FinePack, fineName_, Ls_, par().fineParams.Nm, env().getRbGrid(Ls_));
+    envCreate(CoarsePack, coarseName_, Ls_, par().coarseParams.Nm, coarseGridRb_.get());
+    auto &fine   = envGet(FinePack, fineName_);
+    auto &coarse = envGet(CoarsePack, coarseName_);
+    envTmp(SchurFMat, "mat", Ls_, envGet(FMat, par().action));
+    envGetTmp(SchurFMat, mat);
+    envTmp(LCL, "solver", Ls_, env().getRbGrid(Ls_), coarseGridRb_.get(), mat, 
+           Odd, fine.evec, coarse.evec, fine.eval, coarse.eval);
+}
+
+// execution ///////////////////////////////////////////////////////////////////
+template <typename FImpl, int nBasis>
+void TLocalCoherenceLanczos<FImpl, nBasis>::execute(void)
+{
+    auto &finePar   = par().fineParams;
+    auto &coarsePar = par().coarseParams;
+    auto &fine      = envGet(FinePack, fineName_);
+    auto &coarse    = envGet(CoarsePack, coarseName_);
+
+    envGetTmp(LCL, solver);
+    if (par().doFine)
+    {
+        LOG(Message) << "Performing fine grid IRL -- Nstop= " 
+                     << finePar.Nstop << ", Nk= " << finePar.Nk << ", Nm= " 
+                     << finePar.Nm << std::endl;
+        solver.calcFine(finePar.Cheby, finePar.Nstop, finePar.Nk, finePar.Nm,
+                        finePar.resid,finePar.MaxIt, finePar.betastp, 
+                        finePar.MinRes);
+        solver.testFine(finePar.resid*100.0);
+        LOG(Message) << "Orthogonalising" << std::endl;
+        solver.Orthogonalise();
+        if (!par().output.empty())
+        {
+            fine.write(par().output + "_fine");
+        }
+    }
+    if (par().doCoarse)
+    {
+        LOG(Message) << "Performing coarse grid IRL -- Nstop= " 
+                     << coarsePar.Nstop << ", Nk= " << coarsePar.Nk << ", Nm= " 
+                     << coarsePar.Nm << std::endl;
+        solver.calcCoarse(coarsePar.Cheby, par().smoother, par().coarseRelaxTol,
+			              coarsePar.Nstop, coarsePar.Nk, coarsePar.Nm, 
+                          coarsePar.resid, coarsePar.MaxIt, coarsePar.betastp, 
+                          coarsePar.MinRes);
+        solver.testCoarse(coarsePar.resid*100.0, par().smoother, 
+                          par().coarseRelaxTol);
+        if (!par().output.empty())
+        {
+            coarse.write(par().output + "_coarse");
+        }
+    }
+}
+
+END_MODULE_NAMESPACE
+
+END_HADRONS_NAMESPACE
+
+#endif // Hadrons_MSolver_LocalCoherenceLanczos_hpp_

extras/Hadrons/Modules/MSolver/RBPrecCG.hpp

@@ -43,9 +43,10 @@ BEGIN_MODULE_NAMESPACE(MSolver)
 class RBPrecCGPar: Serializable
 {
 public:
-    GRID_SERIALIZABLE_CLASS_MEMBERS(RBPrecCGPar,
+    GRID_SERIALIZABLE_CLASS_MEMBERS(RBPrecCGPar ,
-                                    std::string, action,
+                                    std::string    , action,
-                                    double     , residual);
+                                    unsigned int   , maxIteration,
+                                    double         , residual);
 };
 
 template <typename FImpl>
@@ -69,7 +70,8 @@ protected:
     virtual void execute(void);
 };
 
-MODULE_REGISTER_NS(RBPrecCG, TRBPrecCG<FIMPL>, MSolver);
+MODULE_REGISTER_NS(RBPrecCG,  TRBPrecCG<FIMPL>, MSolver);
+MODULE_REGISTER_NS(ZRBPrecCG, TRBPrecCG<ZFIMPL>, MSolver);
 
 /******************************************************************************
  *                      TRBPrecCG template implementation                     *
@@ -117,14 +119,16 @@ void TRBPrecCG<FImpl>::setup(void)
     auto &mat   = envGet(FMat, par().action);
     auto solver = [&mat, this](FermionField &sol, const FermionField &source)
     {
-        ConjugateGradient<FermionField>           cg(par().residual, 10000);
+        ConjugateGradient<FermionField>           cg(par().residual, 
-        SchurRedBlackDiagMooeeSolve<FermionField> schurSolver(cg);
+                                                     par().maxIteration);
+        HADRONS_DEFAULT_SCHUR_SOLVE<FermionField> schurSolver(cg);
         
         schurSolver(mat, source, sol);
     };
     envCreate(SolverFn, getName(), Ls, solver);
 }
 
+
 // execution ///////////////////////////////////////////////////////////////////
 template <typename FImpl>
 void TRBPrecCG<FImpl>::execute(void)

extras/Hadrons/Modules/MSource/SeqConserved.hpp

@@ -8,6 +8,7 @@ Copyright (C) 2015-2018
 
 Author: Antonin Portelli <antonin.portelli@me.com>
 Author: Lanny91 <andrew.lawson@gmail.com>
+Author: Vera Guelpers <vmg1n14@soton.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
@@ -38,9 +39,11 @@ BEGIN_HADRONS_NAMESPACE
 
 /*
  
- Sequential source
+ Sequential source with insertion of conserved current. 
+ Additionally optional insertion of a photon field A_\mu(x).
  -----------------------------
- * src_x = q_x * theta(x_3 - tA) * theta(tB - x_3) * J_mu * exp(i x.mom)
+ * src_x = sum_{mu=mu_min}^{mu_max} 
+     q_x * theta(x_3 - tA) * theta(tB - x_3) * J_mu * exp(i x.mom) (* A_\mu(x))
  
  * options:
  - q: input propagator (string)
@@ -48,8 +51,10 @@ BEGIN_HADRONS_NAMESPACE
  - tA: begin timeslice (integer)
  - tB: end timesilce (integer)
  - curr_type: type of conserved current to insert (Current)
- - mu: Lorentz index of current to insert (integer)
+ - mu_min: begin Lorentz Index (integer)
+ - mu_max: end Lorentz Index (integer)
  - mom: momentum insertion, space-separated float sequence (e.g ".1 .2 1. 0.")
+ - photon: optional photon field (string)
  
  */
 
@@ -67,8 +72,10 @@ public:
                                     unsigned int, tA,
                                     unsigned int, tB,
                                     Current,      curr_type,
-                                    unsigned int, mu,
+                                    unsigned int, mu_min,
-                                    std::string,  mom);
+                                    unsigned int, mu_max,
+                                    std::string,  mom,
+                                    std::string,  photon);
 };
 
 template <typename FImpl>
@@ -76,6 +83,8 @@ class TSeqConserved: public Module<SeqConservedPar>
 {
 public:
     FERM_TYPE_ALIASES(FImpl,);
+public:
+    typedef PhotonR::GaugeField     EmField;
 public:
     // constructor
     TSeqConserved(const std::string name);
@@ -89,10 +98,14 @@ protected:
     virtual void setup(void);
     // execution
     virtual void execute(void);
+private:
+    bool        SeqhasPhase_{false}; 
+    std::string SeqmomphName_;
 };
 
 MODULE_REGISTER_NS(SeqConserved, TSeqConserved<FIMPL>, MSource);
 
+
 /******************************************************************************
  *                      TSeqConserved implementation                          *
  ******************************************************************************/
@@ -100,6 +113,7 @@ MODULE_REGISTER_NS(SeqConserved, TSeqConserved<FIMPL>, MSource);
 template <typename FImpl>
 TSeqConserved<FImpl>::TSeqConserved(const std::string name)
 : Module<SeqConservedPar>(name)
+, SeqmomphName_ (name + "_Seqmomph")
 {}
 
 // dependencies/products ///////////////////////////////////////////////////////
@@ -107,7 +121,8 @@ template <typename FImpl>
 std::vector<std::string> TSeqConserved<FImpl>::getInput(void)
 {
     std::vector<std::string> in = {par().q, par().action};
-    
+    if (!par().photon.empty()) in.push_back(par().photon);
+        
     return in;
 }
 
@@ -116,7 +131,7 @@ std::vector<std::string> TSeqConserved<FImpl>::getOutput(void)
 {
     std::vector<std::string> out = {getName()};
     
-    return out;
+   return out;
 }
 
 // setup ///////////////////////////////////////////////////////////////////////
@@ -125,6 +140,10 @@ void TSeqConserved<FImpl>::setup(void)
 {
     auto Ls_ = env().getObjectLs(par().action);
     envCreateLat(PropagatorField, getName(), Ls_);
+    envTmpLat(PropagatorField, "src_tmp");
+    envCacheLat(LatticeComplex, SeqmomphName_);
+    envTmpLat(LatticeComplex, "coor");
+    envTmpLat(LatticeComplex, "latt_compl");
 }
 
 // execution ///////////////////////////////////////////////////////////////////
@@ -134,27 +153,79 @@ void TSeqConserved<FImpl>::execute(void)
     if (par().tA == par().tB)
     {
         LOG(Message) << "Generating sequential source with conserved "
-                     << par().curr_type << " current insertion (mu = " 
+                     << par().curr_type << " current at " 
-                     << par().mu << ") at " << "t = " << par().tA << std::endl;
+		     << "t = " << par().tA << " summed over the indices " 
+		     << par().mu_min << " <= mu <= " << par().mu_max 
+		     << std::endl;
     }
     else
     {
         LOG(Message) << "Generating sequential source with conserved "
-                     << par().curr_type << " current insertion (mu = " 
+                     << par().curr_type << " current for " 
-                     << par().mu << ") for " << par().tA << " <= t <= " 
+                     << par().tA << " <= t <= " 
-                     << par().tB << std::endl;
+                     << par().tB << " summed over the indices " 
+		     << par().mu_min << " <= mu <= " << par().mu_max
+	             << std::endl;
     }
     auto &src = envGet(PropagatorField, getName());
+    envGetTmp(PropagatorField, src_tmp);
+    src_tmp = src;
     auto &q   = envGet(PropagatorField, par().q);
     auto &mat = envGet(FMat, par().action);
+    envGetTmp(LatticeComplex, latt_compl);
 
-    std::vector<Real> mom = strToVec<Real>(par().mom);
+    src = zero;
-    mat.SeqConservedCurrent(q, src, par().curr_type, par().mu, 
+
-                            mom, par().tA, par().tB);
+    //exp(ipx)
+    auto &mom_phase = envGet(LatticeComplex, SeqmomphName_);
+    if (!SeqhasPhase_)
+    {    
+        std::vector<Real> mom = strToVec<Real>(par().mom);
+        mom_phase = zero;
+        Complex           i(0.0,1.0);
+        envGetTmp(LatticeComplex, coor);
+        for(unsigned int mu = 0; mu < env().getNd(); mu++)
+        {
+            LatticeCoordinate(coor, mu);
+            mom_phase = mom_phase + (mom[mu]/env().getGrid()->_fdimensions[mu])*coor;
+        }
+        mom_phase = exp((Real)(2*M_PI)*i*mom_phase);
+        SeqhasPhase_ = true;
+    }
+    LOG(Message) << "Inserting momentum " << strToVec<Real>(par().mom) << std::endl;
+
+
+
+    if (!par().photon.empty())    	
+    {
+	 LOG(Message) << "Inserting the stochastic photon field " << par().photon << std::endl;
+    }
+
+    for(unsigned int mu=par().mu_min;mu<=par().mu_max;mu++)
+    {
+        if (!par().photon.empty())    	
+        {
+	    //Get the stochastic photon field, if required
+            auto &stoch_photon = envGet(EmField,  par().photon);
+    	    latt_compl =  PeekIndex<LorentzIndex>(stoch_photon, mu) * mom_phase;
+        }
+        else
+        {
+            latt_compl = mom_phase;
+        } 
+
+    	mat.SeqConservedCurrent(q, src_tmp, par().curr_type, mu, 
+                             par().tA, par().tB, latt_compl);
+	src += src_tmp;
+
+    }	
+
+ 
 }
 
+
 END_MODULE_NAMESPACE
 
 END_HADRONS_NAMESPACE
 
-#endif // Hadrons_SeqConserved_hpp_
+#endif // Hadrons_MSource_SeqConserved_hpp_

extras/Hadrons/VirtualMachine.cc

@@ -111,6 +111,7 @@ void VirtualMachine::pushModule(VirtualMachine::ModPt &pt)
             {
                 // output does not exists, add it
                 env().addObject(out, address);
+                module_[address].output.push_back(env().getObjectAddress(out));
             }
             else
             {
@@ -296,12 +297,65 @@ void VirtualMachine::makeModuleGraph(void)
     {
         for (auto &in: module_[m].input)
         {
-            graph.addEdge(env().getObjectModule(in), m);
+            int min = env().getObjectModule(in);
+
+            if (min < 0)
+            {
+                HADRON_ERROR(Definition, "object with address " 
+                             + std::to_string(in) 
+                             + " is not produced by any module");
+            }
+            else
+            {
+                graph.addEdge(min, m);
+            }
         }
     }
     graph_ = graph;
 }
 
+// dump GraphViz graph /////////////////////////////////////////////////////////
+void VirtualMachine::dumpModuleGraph(std::ostream &out)
+{
+    makeModuleGraph();
+    out << "digraph hadrons {" << std::endl;
+    out << "node [shape=record, fontname=\"Courier\", fontsize=\"11\"];" << std::endl;
+    out << "graph [fontname = \"Courier\", fontsize=\"11\"];" << std::endl;
+    out << "edge [fontname = \"Courier\", fontsize=\"11\"];"<< std::endl;
+    for (unsigned int m = 0; m < module_.size(); ++m)
+    {
+
+    }
+    for (unsigned int m = 0; m < module_.size(); ++m)
+    {
+        for (auto &in: module_[m].input)
+        {
+            int min = env().getObjectModule(in);
+
+            out << min << " -> " << m << " [ label = \""
+                << env().getObjectName(in) << "\" ];" << std::endl;
+        }
+    }
+    for (unsigned int m = 0; m < module_.size(); ++m)
+    {
+        out <<  m << " [ label = \"{<f0> " << getModule(m)->getRegisteredName()
+            << " |<f1> " << getModuleName(m) << "}\" ];" << std::endl;
+    }
+    out << "}\n" << std::endl;
+}
+
+void VirtualMachine::dumpModuleGraph(void)
+{
+    dumpModuleGraph(std::cout);
+}
+
+void VirtualMachine::dumpModuleGraph(const std::string filename)
+{
+    std::ofstream f(filename);
+
+    dumpModuleGraph(f);
+}
+
 // memory profile //////////////////////////////////////////////////////////////
 const VirtualMachine::MemoryProfile & VirtualMachine::getMemoryProfile(void)
 {
@@ -424,11 +478,17 @@ void VirtualMachine::memoryProfile(const unsigned int address)
         cleanEnvironment();
         for (auto &in: m->getInput())
         {
-            memoryProfile(env().getObjectModule(in));
+            if (!env().hasCreatedObject(in))
+            {
+                memoryProfile(env().getObjectModule(in));
+            }
         }
         for (auto &ref: m->getReference())
         {
-            memoryProfile(env().getObjectModule(ref));
+            if (!env().hasCreatedObject(ref))
+            {
+                memoryProfile(env().getObjectModule(ref));
+            }
         }
         m->setup();
         updateProfile(address);
@@ -532,7 +592,7 @@ VirtualMachine::Program VirtualMachine::schedule(const GeneticPar &par)
     gen = 0;
     do
     {
-        LOG(Debug) << "Generation " << gen << ":" << std::endl;
+        //LOG(Debug) << "Generation " << gen << ":" << std::endl;
         scheduler.nextGeneration();
         if (gen != 0)
         {

extras/Hadrons/VirtualMachine.hpp

@@ -84,7 +84,7 @@ private:
         const std::type_info      *type{nullptr};
         std::string               name;
         ModPt                     data{nullptr};
-        std::vector<unsigned int> input;
+        std::vector<unsigned int> input, output;
         size_t                    maxAllocated;
     };
 public:
@@ -120,6 +120,10 @@ public:
     void                printContent(void) const;
     // module graph (could be a const reference if topoSort was const)
     Graph<unsigned int> getModuleGraph(void);
+    // dump GraphViz graph
+    void                dumpModuleGraph(std::ostream &out);
+    void                dumpModuleGraph(void);
+    void                dumpModuleGraph(const std::string filename);
     // memory profile
     const MemoryProfile &getMemoryProfile(void);
     // garbage collector

extras/Hadrons/modules.inc

@@ -31,12 +31,13 @@ modules_hpp =\
   Modules/MSource/SeqConserved.hpp \
   Modules/MSink/Smear.hpp \
   Modules/MSink/Point.hpp \
+  Modules/MSolver/LocalCoherenceLanczos.hpp \
   Modules/MSolver/RBPrecCG.hpp \
   Modules/MGauge/Random.hpp \
+  Modules/MGauge/FundtoHirep.hpp \
   Modules/MGauge/StochEm.hpp \
   Modules/MGauge/Unit.hpp \
   Modules/MGauge/UnitEm.hpp \
-  Modules/MGauge/FundtoHirep.hpp \
   Modules/MUtilities/TestSeqGamma.hpp \
   Modules/MUtilities/TestSeqConserved.hpp \
   Modules/MLoop/NoiseLoop.hpp \
@@ -48,10 +49,15 @@ modules_hpp =\
   Modules/MAction/DWF.hpp \
   Modules/MAction/Wilson.hpp \
   Modules/MAction/WilsonClover.hpp \
+  Modules/MAction/ZMobiusDWF.hpp \
   Modules/MScalarSUN/Div.hpp \
   Modules/MScalarSUN/TrMag.hpp \
+  Modules/MScalarSUN/EMT.hpp \
   Modules/MScalarSUN/TwoPoint.hpp \
   Modules/MScalarSUN/TrPhi.hpp \
+  Modules/MScalarSUN/Utils.hpp \
+  Modules/MScalarSUN/TransProj.hpp \
+  Modules/MScalarSUN/TrKinetic.hpp \
   Modules/MIO/LoadNersc.hpp \
   Modules/MIO/LoadBinary.hpp
 

lib/algorithms/iterative/Deflation.h

@@ -59,7 +59,7 @@ public:
     assert(evec.size()==eval.size());
     auto N = evec.size();
     for (int i=0;i<N;i++) {
-      Field& tmp = evec[i];
+      const Field& tmp = evec[i];
       axpy(guess,TensorRemove(innerProduct(tmp,src)) / eval[i],tmp,guess);
     }
   }
@@ -86,9 +86,9 @@ public:
     int N = (int)evec_coarse.size();
     CoarseField src_coarse(evec_coarse[0]._grid);
     CoarseField guess_coarse(evec_coarse[0]._grid);    guess_coarse = zero;
-    blockProject(src,src_coarse,subspace);    
+    blockProject(src_coarse,src,subspace);    
     for (int i=0;i<N;i++) {
-      CoarseField & tmp = evec_coarse[i];
+      const CoarseField & tmp = evec_coarse[i];
       axpy(guess_coarse,TensorRemove(innerProduct(tmp,src_coarse)) / eval_coarse[i],tmp,guess_coarse);
     }
     blockPromote(guess_coarse,guess,subspace);

lib/algorithms/iterative/LocalCoherenceLanczos.h

@@ -332,7 +332,7 @@ public:
     // create a smoother and see if we can get a cheap convergence test and smooth inside the IRL
     //////////////////////////////////////////////////////////////////////////////////////////////////
     Chebyshev<FineField>                          ChebySmooth(cheby_smooth);
-    ProjectedFunctionHermOp<Fobj,CComplex,nbasis> ChebyOp (ChebySmooth,_FineOp,_subspace);
+    ProjectedFunctionHermOp<Fobj,CComplex,nbasis> ChebyOp (ChebySmooth,_FineOp,subspace);
     ImplicitlyRestartedLanczosSmoothedTester<Fobj,CComplex,nbasis> ChebySmoothTester(ChebyOp,ChebySmooth,_FineOp,subspace,relax);
 
     for(int k=0;k<evec_coarse.size();k++){
@@ -373,14 +373,14 @@ public:
 		  RealD MaxIt, RealD betastp, int MinRes)
   {
     Chebyshev<FineField>                          Cheby(cheby_op);
-    ProjectedHermOp<Fobj,CComplex,nbasis>         Op(_FineOp,_subspace);
+    ProjectedHermOp<Fobj,CComplex,nbasis>         Op(_FineOp,subspace);
-    ProjectedFunctionHermOp<Fobj,CComplex,nbasis> ChebyOp (Cheby,_FineOp,_subspace);
+    ProjectedFunctionHermOp<Fobj,CComplex,nbasis> ChebyOp (Cheby,_FineOp,subspace);
     //////////////////////////////////////////////////////////////////////////////////////////////////
     // create a smoother and see if we can get a cheap convergence test and smooth inside the IRL
     //////////////////////////////////////////////////////////////////////////////////////////////////
 
     Chebyshev<FineField>                                           ChebySmooth(cheby_smooth);
-    ImplicitlyRestartedLanczosSmoothedTester<Fobj,CComplex,nbasis> ChebySmoothTester(ChebyOp,ChebySmooth,_FineOp,_subspace,relax);
+    ImplicitlyRestartedLanczosSmoothedTester<Fobj,CComplex,nbasis> ChebySmoothTester(ChebyOp,ChebySmooth,_FineOp,subspace,relax);
 
     evals_coarse.resize(Nm);
     evec_coarse.resize(Nm,_CoarseGrid);

lib/log/Log.h

@@ -86,7 +86,7 @@ protected:
   Colours &Painter;
   int active;
   int timing_mode;
-  int topWidth{-1};
+  int topWidth{-1}, chanWidth{-1};
   static int timestamp;
   std::string name, topName;
   std::string COLOUR;
@@ -126,6 +126,7 @@ public:
     }
   }
   void setTopWidth(const int w) {topWidth = w;}
+  void setChanWidth(const int w) {chanWidth = w;}
 
   friend std::ostream& operator<< (std::ostream& stream, Logger& log){
 
@@ -136,7 +137,12 @@ public:
         stream << std::setw(log.topWidth);
       }
       stream << log.topName << log.background()<< " : ";
-      stream << log.colour() <<  std::left << log.name << log.background() << " : ";
+      stream << log.colour() <<  std::left;
+      if (log.chanWidth > 0)
+      {
+        stream << std::setw(log.chanWidth);
+      }
+      stream << log.name << log.background() << " : ";
       if ( log.timestamp ) {
 	log.StopWatch->Stop();
 	GridTime now = log.StopWatch->Elapsed();

lib/qcd/action/fermion/FermionOperator.h

@@ -125,9 +125,9 @@ namespace Grid {
                                        PropagatorField &q_out,
                                        Current curr_type,
                                        unsigned int mu,
-                                       std::vector<Real> mom,
                                        unsigned int tmin, 
-                                       unsigned int tmax)=0;
+                                       unsigned int tmax,
+				       ComplexField &lattice_cmplx)=0;
     };
 
   }

lib/qcd/action/fermion/ImprovedStaggeredFermion.cc

@@ -407,17 +407,19 @@ void ImprovedStaggeredFermion<Impl>::ContractConservedCurrent(PropagatorField &q
 }
 
 template <class Impl>
-void ImprovedStaggeredFermion<Impl>::SeqConservedCurrent(PropagatorField &q_in,
+void ImprovedStaggeredFermion<Impl>::SeqConservedCurrent(PropagatorField &q_in, 
-                                                         PropagatorField &q_out,
+                                              PropagatorField &q_out,
-                                                         Current curr_type,
+                                              Current curr_type,
-                                                         unsigned int mu, 
+                                              unsigned int mu,
-                                                         std::vector<Real> mom,
+                                              unsigned int tmin, 
-                                                         unsigned int tmin,
+                                              unsigned int tmax,
-                                                         unsigned int tmax)
+					      ComplexField &lattice_cmplx)
 {
     assert(0);
+
 }
 
+
 FermOpStaggeredTemplateInstantiate(ImprovedStaggeredFermion);
 
   //AdjointFermOpTemplateInstantiate(ImprovedStaggeredFermion);

lib/qcd/action/fermion/ImprovedStaggeredFermion.h

@@ -166,13 +166,13 @@ class ImprovedStaggeredFermion : public StaggeredKernels<Impl>, public ImprovedS
                                 PropagatorField &q_out,
                                 Current curr_type,
                                 unsigned int mu);
-  void SeqConservedCurrent(PropagatorField &q_in,
+  void SeqConservedCurrent(PropagatorField &q_in, 
                            PropagatorField &q_out,
-                           Current curr_type,
+                           Current curr_type, 
-                           unsigned int mu, 
+                           unsigned int mu,
-                           std::vector<Real> mom,
+                           unsigned int tmin, 
-                           unsigned int tmin,
+                           unsigned int tmax,
-                           unsigned int tmax);
+			   ComplexField &lattice_cmplx);
 };
 
 typedef ImprovedStaggeredFermion<StaggeredImplF> ImprovedStaggeredFermionF;

lib/qcd/action/fermion/ImprovedStaggeredFermion5D.cc

@@ -419,15 +419,16 @@ void ImprovedStaggeredFermion5D<Impl>::ContractConservedCurrent(PropagatorField
 }
 
 template <class Impl>
-void ImprovedStaggeredFermion5D<Impl>::SeqConservedCurrent(PropagatorField &q_in,
+void ImprovedStaggeredFermion5D<Impl>::SeqConservedCurrent(PropagatorField &q_in, 
-                                                          PropagatorField &q_out,
+                                              PropagatorField &q_out,
-                                                          Current curr_type,
+                                              Current curr_type,
-                                                          unsigned int mu, 
+                                              unsigned int mu,
-                                                          std::vector<Real> mom,
+                                              unsigned int tmin, 
-                                                          unsigned int tmin,
+                                              unsigned int tmax,
-                                                          unsigned int tmax)
+					      ComplexField &lattice_cmplx)
 {
     assert(0);
+
 }
 
 FermOpStaggeredTemplateInstantiate(ImprovedStaggeredFermion5D);

lib/qcd/action/fermion/ImprovedStaggeredFermion5D.h

@@ -178,13 +178,13 @@ namespace QCD {
                                   PropagatorField &q_out,
                                   Current curr_type,
                                   unsigned int mu);
-    void SeqConservedCurrent(PropagatorField &q_in,
+    void SeqConservedCurrent(PropagatorField &q_in, 
                              PropagatorField &q_out,
-                             Current curr_type,
+                             Current curr_type, 
-                             unsigned int mu, 
+                             unsigned int mu,
-                             std::vector<Real> mom,
+                             unsigned int tmin, 
-                             unsigned int tmin,
+                             unsigned int tmax,
-                             unsigned int tmax);
+                 	     ComplexField &lattice_cmplx);
   };
 
 }}

lib/qcd/action/fermion/WilsonFermion.cc

@@ -407,40 +407,30 @@ void WilsonFermion<Impl>::ContractConservedCurrent(PropagatorField &q_in_1,
     }
 }
 
+
 template <class Impl>
 void WilsonFermion<Impl>::SeqConservedCurrent(PropagatorField &q_in, 
                                               PropagatorField &q_out,
                                               Current curr_type,
                                               unsigned int mu,
-                                              std::vector<Real> mom,
                                               unsigned int tmin, 
-                                              unsigned int tmax)
+                                              unsigned int tmax,
+					      ComplexField &lattice_cmplx)
 {
     conformable(_grid, q_in._grid);
     conformable(_grid, q_out._grid);
-    Lattice<iSinglet<Simd>> ph(_grid), coor(_grid);
-    ComplexD i(0.0,1.0);
     PropagatorField tmpFwd(_grid), tmpBwd(_grid), tmp(_grid);
     unsigned int tshift = (mu == Tp) ? 1 : 0;
     unsigned int LLt    = GridDefaultLatt()[Tp];
 
-    // Momentum projection
-    ph = zero;
-    for(unsigned int mu = 0; mu < Nd - 1; mu++)
-    {
-        LatticeCoordinate(coor, mu);
-        ph = ph + mom[mu]*coor*((1./(_grid->_fdimensions[mu])));
-    }
-    ph = exp((RealD)(2*M_PI)*i*ph);
-
     q_out = zero;
     LatticeInteger coords(_grid);
     LatticeCoordinate(coords, Tp);
 
     // Need q(x + mu) and q(x - mu).
     tmp = Cshift(q_in, mu, 1);
-    tmpFwd = tmp*ph;
+    tmpFwd = tmp*lattice_cmplx;
-    tmp = ph*q_in;
+    tmp = lattice_cmplx*q_in;
     tmpBwd = Cshift(tmp, mu, -1);
 
     parallel_for (unsigned int sU = 0; sU < Umu._grid->oSites(); ++sU)
@@ -475,6 +465,8 @@ void WilsonFermion<Impl>::SeqConservedCurrent(PropagatorField &q_in,
                                                 Umu, sU, mu, t_mask);
         }
     }
+
+
 }
 
 FermOpTemplateInstantiate(WilsonFermion);

lib/qcd/action/fermion/WilsonFermion.h

@@ -174,13 +174,13 @@ class WilsonFermion : public WilsonKernels<Impl>, public WilsonFermionStatic {
                                 PropagatorField &q_out,
                                 Current curr_type,
                                 unsigned int mu);
-  void SeqConservedCurrent(PropagatorField &q_in,
+  void SeqConservedCurrent(PropagatorField &q_in, 
-                           PropagatorField &q_out,
+                             PropagatorField &q_out,
-                           Current curr_type,
+                             Current curr_type, 
-                           unsigned int mu, 
+                             unsigned int mu,
-                           std::vector<Real> mom,
+                             unsigned int tmin, 
-                           unsigned int tmin,
+                             unsigned int tmax,
-                           unsigned int tmax);
+			     ComplexField &lattice_cmplx);
 };
 
 typedef WilsonFermion<WilsonImplF> WilsonFermionF;

lib/qcd/action/fermion/WilsonFermion5D.cc

@@ -779,92 +779,89 @@ void WilsonFermion5D<Impl>::ContractConservedCurrent(PropagatorField &q_in_1,
 }
 
 
+
 template <class Impl>
 void WilsonFermion5D<Impl>::SeqConservedCurrent(PropagatorField &q_in, 
                                                 PropagatorField &q_out,
                                                 Current curr_type, 
                                                 unsigned int mu,
-                                                std::vector<Real> mom,
                                                 unsigned int tmin, 
-                                                unsigned int tmax)
+                                                unsigned int tmax,
+						ComplexField &lattice_cmplx)
 {
     conformable(q_in._grid, FermionGrid());
     conformable(q_in._grid, q_out._grid);
-    Lattice<iSinglet<Simd>> ph(FermionGrid()), coor(FermionGrid());
+    PropagatorField tmp(GaugeGrid()),tmp2(GaugeGrid());
-    PropagatorField tmpFwd(FermionGrid()), tmpBwd(FermionGrid()),
-                    tmp(FermionGrid());
-    ComplexD i(0.0, 1.0);
     unsigned int tshift = (mu == Tp) ? 1 : 0;
     unsigned int LLs = q_in._grid->_rdimensions[0];
     unsigned int LLt    = GridDefaultLatt()[Tp];
 
-    // Momentum projection.
-    ph = zero;
-    for(unsigned int nu = 0; nu < Nd - 1; nu++)
-    {
-        // Shift coordinate lattice index by 1 to account for 5th dimension.
-        LatticeCoordinate(coor, nu + 1);
-        ph = ph + mom[nu]*coor*((1./(_FourDimGrid->_fdimensions[nu])));
-    }
-    ph = exp((RealD)(2*M_PI)*i*ph);
-
     q_out = zero;
     LatticeInteger coords(_FourDimGrid);
     LatticeCoordinate(coords, Tp);
 
-    // Need q(x + mu, s) and q(x - mu, s). 5D lattice so shift 4D coordinate mu
-    // by one.
-    tmp = Cshift(q_in, mu + 1, 1);
-    tmpFwd = tmp*ph;
-    tmp = ph*q_in;
-    tmpBwd = Cshift(tmp, mu + 1, -1);
 
-    parallel_for (unsigned int sU = 0; sU < Umu._grid->oSites(); ++sU)
+    for (unsigned int s = 0; s < LLs; ++s)
     {
-        // Compute the sequential conserved current insertion only if our simd
+        bool axial_sign = ((curr_type == Current::Axial) && (s < (LLs / 2)));
-        // object contains a timeslice we need.
+	bool tadpole_sign = (curr_type == Current::Tadpole);
-        vInteger t_mask   = ((coords._odata[sU] >= tmin) &&
+	bool switch_sgn = tadpole_sign || axial_sign;
-                             (coords._odata[sU] <= tmax));
-        Integer timeSlices = Reduce(t_mask);
 
-        if (timeSlices > 0)
+
-        {
+        //forward direction: Need q(x + mu, s)*A(x)
-            unsigned int sF = sU * LLs;
+        ExtractSlice(tmp2, q_in, s, 0);  //q(x,s) 
-            for (unsigned int s = 0; s < LLs; ++s)
+        tmp = Cshift(tmp2, mu, 1);	 //q(x+mu,s)
+        tmp2 = tmp*lattice_cmplx;	 //q(x+mu,s)*A(x)	
+
+    	parallel_for (unsigned int sU = 0; sU < Umu._grid->oSites(); ++sU)
+    	{
+            // Compute the sequential conserved current insertion only if our simd
+            // object contains a timeslice we need.
+            vInteger t_mask   = ((coords._odata[sU] >= tmin) &&
+                	         (coords._odata[sU] <= tmax));
+            Integer timeSlices = Reduce(t_mask);
+
+            if (timeSlices > 0)
             {
-                bool axial_sign = ((curr_type == Current::Axial) && (s < (LLs / 2)));
+		unsigned int sF = sU * LLs + s;
-                Kernels::SeqConservedCurrentSiteFwd(tmpFwd._odata[sF], 
+                Kernels::SeqConservedCurrentSiteFwd(tmp2._odata[sU], 
-                                                    q_out._odata[sF], Umu, sU,
+                                              q_out._odata[sF], Umu, sU,
-                                                    mu, t_mask, axial_sign);
+                                              mu, t_mask, switch_sgn);
-                ++sF;
             }
+
         }
 
-        // Repeat for backward direction.
+        //backward direction: Need q(x - mu, s)*A(x-mu)
-        t_mask     = ((coords._odata[sU] >= (tmin + tshift)) && 
+        ExtractSlice(tmp2, q_in, s, 0);  //q(x,s)
-                      (coords._odata[sU] <= (tmax + tshift)));
+        tmp = lattice_cmplx*tmp2;	 //q(x,s)*A(x)
+        tmp2 = Cshift(tmp, mu, -1);	 //q(x-mu,s)*A(x-mu,s)
 
-	//if tmax = LLt-1 (last timeslice) include timeslice 0 if the time is shifted (mu=3)	
+    	parallel_for (unsigned int sU = 0; sU < Umu._grid->oSites(); ++sU)
-	unsigned int t0 = 0;
+    	{
-	if((tmax==LLt-1) && (tshift==1)) t_mask = (t_mask || (coords._odata[sU] == t0 ));
+            vInteger  t_mask     = ((coords._odata[sU] >= (tmin + tshift)) && 
+                   	  	    (coords._odata[sU] <= (tmax + tshift)));
 
-        timeSlices = Reduce(t_mask);
+	    //if tmax = LLt-1 (last timeslice) include timeslice 0 if the time is shifted (mu=3)	
+	    unsigned int t0 = 0;
+	    if((tmax==LLt-1) && (tshift==1)) t_mask = (t_mask || (coords._odata[sU] == t0 ));
 
-        if (timeSlices > 0)
+            Integer timeSlices = Reduce(t_mask);
-        {
+
-            unsigned int sF = sU * LLs;
+            if (timeSlices > 0)
-            for (unsigned int s = 0; s < LLs; ++s)
             {
-                bool axial_sign = ((curr_type == Current::Axial) && (s < (LLs / 2)));
+		unsigned int sF = sU * LLs + s; 
-                Kernels::SeqConservedCurrentSiteBwd(tmpBwd._odata[sF], 
+        	Kernels::SeqConservedCurrentSiteBwd(tmp2._odata[sU], 
-                                                    q_out._odata[sF], Umu, sU,
+                                             q_out._odata[sF], Umu, sU,
-                                                    mu, t_mask, axial_sign);
+                                             mu, t_mask, axial_sign);
-                ++sF;
             }
-        }
+	}
     }
 }
 
+
+
+
+
 FermOpTemplateInstantiate(WilsonFermion5D);
 GparityFermOpTemplateInstantiate(WilsonFermion5D);
   

lib/qcd/action/fermion/WilsonFermion5D.h

@@ -222,13 +222,13 @@ namespace QCD {
                                   PropagatorField &q_out,
                                   Current curr_type, 
                                   unsigned int mu);
-    void SeqConservedCurrent(PropagatorField &q_in,
+    void SeqConservedCurrent(PropagatorField &q_in, 
                              PropagatorField &q_out,
-                             Current curr_type,
+                             Current curr_type, 
                              unsigned int mu,
-                             std::vector<Real> mom,
+                             unsigned int tmin, 
-                             unsigned int tmin,
+                             unsigned int tmax,
-                             unsigned int tmax);
+			     ComplexField &lattice_cmplx);
   };
 
 }}

tests/hadrons/Test_QED.cc

@@ -0,0 +1,260 @@
+/*******************************************************************************
+ Grid physics library, www.github.com/paboyle/Grid
+ 
+ Source file: tests/hadrons/Test_hadrons_spectrum.cc
+ 
+ Copyright (C) 2015-2018
+ 
+ Author: Antonin Portelli <antonin.portelli@me.com>
+ Author: Vera Guelpers    <v.m.guelpers@soton.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.
+ *******************************************************************************/
+
+#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 = {"h"}; //{"l", "s", "c1", "c2", "c3"};
+    std::vector<double>      mass    = {.2}; //{.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";
+    application.setPar(globalPar);
+    // gauge field
+    application.createModule<MGauge::Unit>("gauge");
+    // pt source
+    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";
+
+    //stochastic photon field
+    MGauge::StochEm::Par photonPar;
+    photonPar.gauge = PhotonR::Gauge::feynman;
+    photonPar.zmScheme = PhotonR::ZmScheme::qedL;
+    application.createModule<MGauge::StochEm>("ph_field", photonPar);
+
+
+
+    for (unsigned int i = 0; i < flavour.size(); ++i)
+    {
+        // actions
+        MAction::DWF::Par actionPar;
+        actionPar.gauge = "gauge";
+        actionPar.Ls    = 8;
+        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);
+        
+        // propagators
+        MFermion::GaugeProp::Par quarkPar;
+        quarkPar.solver = "CG_" + flavour[i];
+        quarkPar.source = "pt";
+        application.createModule<MFermion::GaugeProp>("Qpt_" + flavour[i],
+							 quarkPar);
+
+
+	//seq sources with tadpole insertion
+        MSource::SeqConserved::Par seqPar_T;
+        seqPar_T.q         = "Qpt_" + flavour[i] + "_5d";
+        seqPar_T.action    = "DWF_" + flavour[i];
+        seqPar_T.tA        = 0;
+        seqPar_T.tB        = nt-1;
+        seqPar_T.curr_type = Current::Tadpole;
+	seqPar_T.mu_min	   = 0;
+	seqPar_T.mu_max	   = 3;
+        seqPar_T.mom       = "0. 0. 0. 0.";
+        application.createModule<MSource::SeqConserved>("Qpt_" + flavour[i]
+							 + "_seq_T", seqPar_T);
+        // seq propagator with tadpole insertion
+        MFermion::GaugeProp::Par quarkPar_seq_T;
+        quarkPar_seq_T.solver = "CG_" + flavour[i];
+        quarkPar_seq_T.source = "Qpt_" + flavour[i] + "_seq_T";
+        application.createModule<MFermion::GaugeProp>("Qpt_" + flavour[i]
+							+ "_seq_T" + flavour[i],
+							quarkPar_seq_T);
+
+
+
+	//seq sources with conserved vector and photon insertion
+        MSource::SeqConserved::Par seqPar_V;
+        seqPar_V.q         = "Qpt_" + flavour[i] + "_5d";
+        seqPar_V.action    = "DWF_" + flavour[i];
+        seqPar_V.tA        = 0;
+        seqPar_V.tB        = nt-1;
+        seqPar_V.curr_type = Current::Vector;
+	seqPar_V.mu_min	   = 0;
+	seqPar_V.mu_max	   = 3;
+        seqPar_V.mom       = "0. 0. 0. 0.";
+	seqPar_V.photon	   = "ph_field";
+        application.createModule<MSource::SeqConserved>("Qpt_" + flavour[i] 
+						    + "_seq_V_ph", seqPar_V);
+        // seq propagator with conserved vector and photon insertion
+        MFermion::GaugeProp::Par quarkPar_seq_V;
+        quarkPar_seq_V.solver = "CG_" + flavour[i];
+        quarkPar_seq_V.source = "Qpt_" + flavour[i] + "_seq_V_ph";
+        application.createModule<MFermion::GaugeProp>("Qpt_" + flavour[i] 
+						+ "_seq_V_ph_" + flavour[i], 
+							quarkPar_seq_V);
+
+
+
+	//double seq sources with conserved vector and photon insertion
+	//(for self energy)
+        MSource::SeqConserved::Par seqPar_VV;
+        seqPar_VV.q         = "Qpt_" + flavour[i] + "_seq_V_ph_" 
+				+ flavour[i] + "_5d";
+        seqPar_VV.action    = "DWF_" + flavour[i];
+        seqPar_VV.tA        = 0;
+        seqPar_VV.tB        = nt-1;
+        seqPar_VV.curr_type = Current::Vector;
+	seqPar_VV.mu_min    = 0;
+	seqPar_VV.mu_max    = 3;
+        seqPar_VV.mom       = "0. 0. 0. 0.";
+	seqPar_VV.photon    = "ph_field";
+        application.createModule<MSource::SeqConserved>("Qpt_" + flavour[i] 
+						+ "_seq_V_ph" + flavour[i] 
+						+ "_seq_V_ph", seqPar_VV);
+        //double seq propagator with conserved vector and photon insertion
+        MFermion::GaugeProp::Par quarkPar_seq_VV;
+        quarkPar_seq_VV.solver = "CG_" + flavour[i];
+        quarkPar_seq_VV.source = "Qpt_" + flavour[i] + "_seq_V_ph" 
+						+ flavour[i] + "_seq_V_ph";
+        application.createModule<MFermion::GaugeProp>("Qpt_" + flavour[i] 
+						+ "_seq_V_ph_" + flavour[i] 
+						+ "_seq_V_ph_" + flavour[i], 
+							quarkPar_seq_VV);
+
+
+
+    }
+    for (unsigned int i = 0; i < flavour.size(); ++i)
+    for (unsigned int j = i; j < flavour.size(); ++j)
+    {
+        //2pt function contraction
+        MContraction::Meson::Par mesPar;
+        mesPar.output  = "QED/pt_" + flavour[i] + flavour[j];
+        mesPar.q1      = "Qpt_" + flavour[i];
+        mesPar.q2      = "Qpt_" + flavour[j];
+        mesPar.gammas  = "<Gamma5 Gamma5>";
+        mesPar.sink    = "sink";
+        application.createModule<MContraction::Meson>("meson_pt_"
+                                                      + flavour[i] + flavour[j],
+                                                      mesPar);
+
+
+
+        //tadpole contraction
+        MContraction::Meson::Par mesPar_seq_T;
+        mesPar_seq_T.output  = "QED/tadpole_pt_" + flavour[i] + "_T_" 
+					+ flavour[i] + "__" + flavour[j];
+        mesPar_seq_T.q1      = "Qpt_" + flavour[i] + "_seq_T" + flavour[i];
+        mesPar_seq_T.q2      = "Qpt_" + flavour[j];
+        mesPar_seq_T.gammas  = "<Gamma5 Gamma5>";
+        mesPar_seq_T.sink    = "sink";
+        application.createModule<MContraction::Meson>("meson_tadpole_pt_" +
+                                                      flavour[i] + "_seq_T" 
+						      + flavour[i] + flavour[j],
+                                                      mesPar_seq_T);
+
+
+
+        //photon exchange contraction
+        MContraction::Meson::Par mesPar_seq_E;
+        mesPar_seq_E.output  = "QED/exchange_pt_" + flavour[i] + "_V_ph_" 
+				+ flavour[i] + "__" + flavour[j] + "_V_ph_"
+				+ flavour[j];
+        mesPar_seq_E.q1      = "Qpt_" + flavour[i] + "_seq_V_ph_" + flavour[i];
+        mesPar_seq_E.q2      = "Qpt_" + flavour[j] + "_seq_V_ph_" + flavour[j];
+        mesPar_seq_E.gammas  = "<Gamma5 Gamma5>";
+        mesPar_seq_E.sink    = "sink";
+        application.createModule<MContraction::Meson>("meson_exchange_pt_" 
+					+ flavour[i] + "_seq_V_ph_" + flavour[i] 
+					+ flavour[j] + "_seq_V_ph_" + flavour[j],
+                                                      mesPar_seq_E);
+
+
+
+        //self energy contraction
+        MContraction::Meson::Par mesPar_seq_S;
+        mesPar_seq_S.output  = "QED/selfenergy_pt_" + flavour[i] + "_V_ph_" 
+				+ flavour[i] + "_V_ph_" + flavour[i] + "__" 
+				+  flavour[j];
+        mesPar_seq_S.q1      = "Qpt_" + flavour[i] + "_seq_V_ph_" + flavour[i] 
+				+ "_seq_V_ph_" + flavour[i];
+        mesPar_seq_S.q2      = "Qpt_" + flavour[j];
+        mesPar_seq_S.gammas  = "<Gamma5 Gamma5>";
+        mesPar_seq_S.sink    = "sink";
+        application.createModule<MContraction::Meson>("meson_selfenergy_pt_" 
+						    + flavour[i] + "_seq_V_ph_" 
+						    + flavour[i] + "_seq_V_ph_" 
+						    + flavour[i] + flavour[j],
+                                                       mesPar_seq_S);
+
+    }
+
+
+
+    
+    // execution
+    application.saveParameterFile("QED.xml");
+    application.run();
+    
+    // epilogue
+    LOG(Message) << "Grid is finalizing now" << std::endl;
+    Grid_finalize();
+    
+    return EXIT_SUCCESS;
+}

tests/hadrons/Test_hadrons.hpp

@@ -263,7 +263,8 @@ inline void makeConservedSequentialSource(Application &application,
         seqPar.tA        = tS;
         seqPar.tB        = tS;
         seqPar.curr_type = curr;
-        seqPar.mu        = mu;
+        seqPar.mu_min    = mu;
+        seqPar.mu_min    = mu;
         seqPar.mom       = mom;
         application.createModule<MSource::SeqConserved>(srcName, seqPar);
     }